Monk's rhubarb, Rumex alpinus L. (R. alpinus), is a perennial plant native to the mountains of Central and Southern Europe. Currently, the distribution of R. alpinus has been partly affected by its utilization as a vegetable and a medicinal herb. In the mountains of the Czech Republic, it is considered an invasive plant, probably introduced into the Krkonoše Mountains by colonists from the Alps.

This study's main aim was to verify whether R. alpinus was introduced into the Krkonoše Mountains by alpine colonists or whether it was anthropogenically introduced from the Carpathians. Furthermore, the genetic structure of native and introduced populations of R. alpinus was determined.

For the evaluation of genetic structure, 417 samples of R. alpinus were collected from the Alps, Carpathians, Balkan, Pyrenees, and Czech Mountains. In total, 12 simple sequence repeat (SSR) markers were applied.

The results of AMOVA showed a high 60% variation within populations, 27% variation among groups, and 13% among the population within groups. The overall unbiased gene diversity was high (ĥ = 0.55). The higher level of genetic differentiation among populations (F_ST = 0.35; p < 0.01) indicated restricted gene flow between populations. Compared to native populations, limited genetic variability was observed in the nonnative populations. It was concluded that local adaptation, low gene exchange, and genetic drift affected the genetic diversity of nonnative R. alpinus.

The results support a genetic link between Alpine and Czech genotypes of R. alpinus, while the Carpathians genotypes corresponded to the Balkan genotype.

Description of localities

Between 2017 and 2020, plant samples of R. alpinus were collected from different locations in Europe (Table 1, Figure 1). According to Professor Klimeš research´s focused on ancestors of the original alpine colonists and individual families, R. alpinus was collected in the exact places (Krkonoše) where the colonists lived and R. alpinus probably occurred. Equally, individual plants of R. alpinus were collected in Tyrol and Styria in Austria. As control samples, R. alpinus samples were collected from other mountain localities in the Alps (Lombardy, and Graubünden) and the Pyrenees. Finally, to determine the true origin of R. alpinus were collected samples representing plant populations in the East and West Carpathians, and then in the mountainous regions of the Balkans. Populations from other mountainous regions of the Czech Republic (Jizera and Eagle Mountains) were collected for comparison and possible exclusion of other origins of R. alpinus.  

In total, 417 individual leaf samples were collected, representing 31 populations with 9–20 individual plants per population. The 329 plant samples were collected in the European Mountains, Alps, West and East Carpathians, Pyrenees, and Balkan Peninsula, which represent areas where R. alpinus is a native species.

In the Czech Mountains, represented by Eagle Mountains, Jizera Mountains, and Krkonoše (Giant) Mountains where R. alpinus is considered a nonnative/invasive plant (Kopecký, 1973; Lokvenc, 1978; Kubát, 1990; Šťastná et al., 2010; Pyšek et al., 2012), overall88 plant samples were collected (Table 1 and Figure 1).

Population sampling and DNA extraction

At each locality, leave samples were collected in plants separated approximately 200 m apart to avoid collecting the same plant because R. alpinus also reproduces vegetatively and to ensure a representative sample of the population. Samples were dried on silica gel and stored in a collection at the Faculty of Environmental Sciences, Czech University of Life Sciences Prague.

Total genomic DNA was extracted from silica gel-dried leaves of R. alpinus in two repetitions using the protocol of Doyle & Doyle (1978) with a minor modification that involved the addition of 10 mg of polyvinylpyrrolidone (PVPP) (Carl Roth, Germany) and 5 µl of 10 mg/µl RNase A (Thermo Scientific, Czech Republic) during the initial phase prior to incubation. The quality and yield of isolated DNA were assessed on a 0.8% agarose gel in 1× TAE buffer. The concentration and quality of DNA were measured using a spectrophotometer UVS-99/UVIS Drop (Avans Biotech, Taiwan). The extracted DNA samples were diluted to the concentration of 20 ng/μl for subsequent analysis and stored at - 20°C.

Microsatellite (SSR) analysis

Based on the test, 12 polymorphic primer pairs were selected from the 15 primer pairs according to Šurinová et al. (2018) and used (Table S1). DNA amplification was performed in 5 μl reactions consisting of 2.5 μl QIAGEN Multiplex PCR Master Mix; 0.125 μl of each M13-labelled forward, reverse, and fluorolabelled (NED™, PET®, 6-FAM™, VIC® – Table S1) M13 primers (10 μM each in initial volume); 20 ng of DNA dissolved in 0.5 μl TE buffer; and 1.625 μl dH2O. The PCR protocol was performed according to Schuelke (2000). PCR was performed using an Applied Biosystem Thermal Cycler (Applied Biosystems, USA) as follows: an initial denaturation step at 95°C for 15 min followed by 25 cycles of denaturation (95°C for 20 s), annealing (59°C for 30 s) and extension (72°C for 20 s), followed by 10 cycles of denaturation (95°C for 30 s), annealing (53°C for 45 s) and extension (72°C for 45 s) and a final extension at 72°C for 15 min. During the first 25 cycles, specific PCR products are produced, and in the following 10 cycles, the fluorescent M13 tag is ligated to the M13 forward primer. The quality of PCR products has been verified on 2% agarose gels. Three multiplexes were built (Table S1).

Fragment analyses were performed using capillary electrophoresis in an ABI PRISM 3500 Genetic Analyser automated sequencer (Applied Biosystems, USA). Electropherograms were analysed and scored using GeneMarker ver. 1.8 (SoftGenetics, USA).

Statistical analysis

Analysis of the molecular variance test (AMOVA) with 1 000 permutations was calculated in ARLEQUIN software ver. 3.5.2 (Excoffier & Lischer, 2010). The degree of genetic differentiation among populations was also evaluated using ARLEQUIN software using the distance matrix based on the fixation index (F_ST) generated by the program. Further, a distance matrix based on geographical distances was calculated for R. alpinus populations within R program (R development core team, 2021) version 4.0.3 using the routines in geosphere (Hijmans, 2021) library. These were subsequently logarithmically transformed and correlated with F_ST distance matrix using the Mantel test and 9999 permutations.

Nei’s genetic distance was employed to obtain a UPGMA dendrogram after 1,000 bootstrap samplings in TFPGA software (Miller, 1997).

The diversity indices for each population included the percentage of polymorphic loci, the average diversity of the loci using Nei's unbiased gene diversity ĥ (Nei, 1973), and the Shannon information index (Lewontin, 1972; Shannon & Weaver, 1949) were calculated using the POPGENE, version 1.32 (Yeh et al., 1999).

To assess the Hardy-Weinberg equilibrium, we used the ARLEQUIN software ver. 3.5.2 (Excoffier & Lischer, 2010). Was conducted an exact test using a Markov chain with a forecasted chain length of 1,000,000 and 100,000 dememorization steps (Guo & Thompson, 1992). Deviation from HWE was assessed at a significance level of p < 0.05. The results were interpreted according to established guidelines (Levene, 1949).

Another approach to studying population structure analysis is based on Bayesian statistics STRUCTURE, version 2.3.4 (Pritchard et al., 2000) was used to determine the genetic architecture of the R. alpinus populations. Ten independent runs of 1–20 groups (K = 1–20) were performed using the locprior model with admixture and correlated allele frequency (Falush et al., 2003; Hubisz et al., 2009) with the recommended 200,000 Markov chain iterations after a burn-in period of 100 000 iterations. The optimal value of K was estimated based on ln(K) and on the ΔK calculation, which considers the rate of change in the ln P(D) values among successive K runs to account for patterns of dispersal that are not homogeneous among populations (Evanno et al., 2005). The number (K) of clusters into which the sample data (X) were fitted with posterior probability Pr (X|K) was estimated using the same model with 1,000,000 Markov chain iterations after a burn-in period of 100,000 iterations (Evanno et al., 2005).

An exact test for population differentiation was calculated using the Tools for Population Genetic Analyses (TFPGA; version 1.3; Miller, 1997) with 100 000 recommended permutation steps.

To identify potential bottleneck events in the populations under investigation, we employed BOTTLENECK 1.2.02 software (Cornuet & Luikart, 1996; Piry et al., 1999) and heterozygosity excess resulting from population reduction was examined. We utilized three models of mutational equilibrium: the infinite allele model (IAM), the stepwise mutation model (SMM), and the two-phase mutation model (TPM), with the latter being the most appropriate for microsatellites. For the TPM, we employed the default settings, which assumed that 70% of mutations occur in a single step, with a variance of 30 among multiple steps. The significance of these models was assessed using a one-tailed Wilcoxon rank test, which is suitable for dataset analysis with less than 20 microsatellite loci (Piry et al., 1999). A population was deemed to have experienced a bottleneck event only if all three models produced significant results (p-value ≤ 0.05).