Genetic variability of Linnaea borealis – remnant of Eastern-Southeastern Alpine populations indicates a strong bottleneck and in situ glacial survival

Introduction

The current distribution areas of plants were shaped by various environmental changes during and after the glaciation period. Rapid fluctuations in environmental conditions shifted the boundaries of colonized areas (Hewitt 1996) and caused the extreme fragmentation and disjunct distribution of many plant species (Jermakowicz et al. 2017). Plant populations that occur at the edge of their distributional range or, in the case of disjunct distribution, which is caused by habitat rarity, are prone to environmental and biotic stressors. Plants that are scattered and isolated are more susceptible to variability in climate, habitat and biogenetic factors (Kołos et al. 2015). This applies to isolated remnants of e.g. Arctic-Alpine species, resulting from shifts in their distribution range during glacial and inter-glacial periods. However, it also applies to boreal-alpine species that survived climate change during the Holocene and now only occur in a few places in Central Europe (Svenning et al. 2008, Sundberg 2014, Kołos et al. 2015, Jermakowicz et al. 2017). In recent decades, the conservation value of relict populations has been emphasized as they are adapted to stressful ecological conditions, which could be important for future distribution range shifts owing to climate change (Vogler and Reisch 2013).

In a rapidly changing environment, many plant species have become rare due to habitat loss and fragmentation and are restricted to very small, isolated populations with an increased risk of local or regional extinction (Jermakowicz et al. 2017). Apart from immediate extinction, the most severe consequence of habitat fragmentation is the loss of genetic diversity. This occurs due to random genetic drift, increased inbreeding, and reduced gene flow between populations. This leads to a loss of offspring and of population fitness and at the same time reduces adaptability to environmental changes (Willi et al. 2006).

Long distances between isolated populations further reduce pollen transfer and limit pollination (Wilcock and Neiland 2002). Many authors (Scobie and Wilcock 2009, Gaudeul et al. 2019, Kikowska et al. 2022) have reported a decline or even complete failure of sexual reproduction in small populations of rare plant species. All these reproduction problems are even greater in species that rely on cross-pollination (Wróblewska 2013). A typical example is the twinflower, Linnaea borealis L. (Caprifoliaceae), a circumboreal plant species that is already in severe decline throughout Europe (except in its boreal range in Scandinavia). However, its geographical distribution ranges across the Northern Hemisphere, occurring from Scotland and northern Europe through Russia to Siberia, northern Asia to Kamchatka and Japan, northern China and Mongolia, and from Alaska and Canada to Greenland (Alm 2006). This clonal, self-incompatible plant is a shade-tolerant, evergreen dwarf shrub. In the natural Swedish populations, L. borealis reproduces vegetatively, mainly by stolons, and forms large clones up to 10 m long (Eriksson 1988, Niva 2003). In Scotland, low seed production is attributed to reproductive isolation caused by the lack of availability of compatible mates within populations and limited pollen dispersal between them (Ross and Roi 1990, Wilcock and Neiland 2002, Scobie and Wilcock 2009, Wiberg et al. 2016).

In the Alps, L. borealis is considered a relict species with some localities in France, Switzerland, Italy (Fornaciari 1964, Hegi 1966, Wilhalm 2010), Austria and Slovenia (the southeasternmost locality) (Wraber 1963). According to the literature, there should be three almost neighbouring sites in the Eastern-Southeastern Alps (Ernet in Franz 2011): in the Austrian part of Carinthia (1) and Styria (2), and in the Slovenian part of Styria (3). However, the species is still present at two sites (Slovenian and Austrian parts of Styria, but the third site (Carinthia (AT)) has not been confirmed since 1992 (Ernet and Franz 2011). All these sites are of great importance from a biogeographical point of view and raise a number of research questions related to the survival of L. borealis in this disjunct southeasternmost distribution. In this respect, the L. borealis can be considered a model species for the study of glacial relicts in the Alps. In addition, various conservation measures have already been implemented to guard against any further loss of its habitat and genetic diversity (Wróblewska 2013, Wiberg et al. 2016).

Our aim was therefore to investigate the genetic variability within and between two isolated (but spatially close) populations of L. borealis from the Eastern-Southeastern Alps (Julian Alps (SI) and Gurktaler Alps (AT)) and to compare this with the genetic variability of three populations from the northern range area (Sweden, Scotland and Russia). We used nine previously described microsatellite loci (A’Hara et al. 2011) to assess genetic variability. In this way, we attempted (1) to determine the extent to which the two isolated populations are genetically distinct and (2) to answer the question of whether the L. borealis of the Eastern-Southeastern Alps is a postglacial colonizer or a glacial relict. Our findings will provide essential information for the successful conservation strategies for the survival of this plant species in situ over a longer period of time, especially from the aspect of climate change.

Material and methods

Study area and field methods

In our study of the L. borealis population, we focused on two disjunct populations in the Eastern-Southeastern Alps, while the other three populations sampled across the distribution range served as control populations, for comparison (only). The Eastern-Southeastern Alps populations were in the Julian Alps, Soteska near Nomenj (495 m a.s.l), Slovenia (SI) and the second location was Dieslingsee in the Gurktaler Alps (1850 m a.s.l.), Styria, Austria (AT). Control populations were: (1) in a mixed coniferous forest (dominated by Scots pine) just outside the town of Bollnäs in central Sweden at an elevation of about 70 m a.s.l. (SE); (2) in the Cairngorms National Park in Scotland, United Kingdom, Great Britain (GB), and (3) in the Altay Mountains along the valley of the Aktru River, Russia (RU) (Fig. 1).

Fig. 1. Locations of study areas () of populations of Linnaea borealis. A – Soteska by Nomenj, Slovenia, B – Altay Mountains along the valley of the Aktru River, Russia, C – Dieslingsee, Austria, D – the city of Bollnäs, central Sweden, E – Great Britain.

The population in the Julian Alps is the southeasternmost locality of L. borealis in Europe (Wraber 1963). The twinflower population (Fig. 2) grows in stands of the association Rhodothamno-Rhododendretum hirsuti (Aichinger 1933) Br.‐Bl. & Sissingh in Br.‐Bl. et al. 1939 (Wraber 1963, Šilc and Čarni 2012). The location in the Gurktaler Alps is in the Styrian part of the mountain range and was first mentioned by Ernet and Franz (2011). At this location, stands with L. borealis belong to the association Rhododendretum ferruginei Rübel 1911. The two Southeast Alpine populations grow on slopes exposed to the north, where cold air emerges from under the rocks, especially in summer (Martinčič 1977).

Fig. 2. Habitat in Dieslingsee in the Gurktaler Alps, Styria, Austria (AT) (left) and Linnaea borealis in flower (right). Photo: N. Pipenbaher.

The plants were randomly sampled at four study sites (SI, AT, SE, RU) (Fig.1). To avoid selecting the same genet, we sampled shoots growing at least 3-5 m apart, as suggested by Wróblewska (2013). The small number of samples is correlated with study site size of the Eastern-Southeastern Alps population which it is around 0.2 ha. The shoots were stored in a cool place (4 °C) in the laboratory until DNA isolation. DNA samples from all populations were collected with permission from representative institutions. In this study, we used pre-isolated DNA samples from GB. A total of seventy-three samples of L. borealis were included in the study.

DNA isolation and microsatellite analysis

Total DNA was extracted from fresh young leaves using the CTAB protocol (Doyle and Doyle 1987). To approximately 2-3 square centimetres of fresh leaf tissue, one ml of preheated (68 °C) CTAB extraction buffer [2% (w/v) CTAB, cetyltriammonium bromide (Sigma, Germany), 1.4 NaCl, 20 mM EDTA, 100mM Tris-HCl (pH 8.0), 0.2% 2-mercaptoethanol] was added and well homogenized with a mortar and pestle and transferred to a 2 mL tube. Samples were incubated for 1.5 h at 68 °C in a water bath. After incubation, 600 μL of chloroform: isoamyl alcohol in a proportion of 24:1 were added, and the samples were thoroughly mixed. The mixtures were centrifuged at 11000 rev./min for 10 min. After centrifugation, the supernatant was transferred to a fresh 1.5 mL tube and the DNA was precipitated by the addition of 0.1 vol. of 3 M sodium acetate and 1 vol. of ice-cold isopropanol and kept at -20 °C for 30 min. Samples were again centrifuged at 11000 rev./min for 10 min. The pellet was washed in 70% ethanol for 20 min, air dried and rehydrated in 100 μL of TE buffer (10 mM Tris-HCl, 1 mM EDTA, pH 8.0). Two separate extractions were performed for each plant sample. The DNA concentration was estimated using a DNA fluorometer DQ 300 (Hoefer, Inc., Holliston, Massachusetts). For microsatellite analysis nine microsatellite loci developed by A’Hara et al. (2011) were used (On-line Suppl. Tab. 1). The number of primers sufficient for reliable variety identification depends on the nature and discriminating power of each primer. Cytogenetic study revealed that the number of chromosomes in cells of L. borealis is 2n = 32 (Packer 1964).

For the polymerase chain reaction (PCR) 5 μL multiplex PCR kit (Qiagen, Germany) 20 ng DNA, 0.5 μL of each primer (Sigma, Germany) and RNase-free water were used. PCR conditions consisted of a 15 min hot start at 94 °C; 40 cycles of denaturation at 94 °C for 30 s, annealing at 60.4 °C for 120 s and an extension step at 72 °C for 90 s. The reactions were terminated by incubation at 60 °C for 30 min. PCR was performed using a Whatman Biometra T-Gradient thermocycler (Goettingen, Germany). Capillary electrophoresis of PCR products was performed using a Beckman Coulter CEQ8000 sequencer, according to the manufacturer’s instructions. Fragment size analysis was performed using the integrated software (Fragment Analyzer). A fluorescently labelled size marker (Beckman Coulter DNA Size Standard Kit 400 bp) was used as a molecular weight reference.

Data analysis

For each microsatellite locus, the following genetic diversity parameters were assessed: the number of alleles per locus (n), allele frequencies, observed (H_o) and expected heterozygosity (H_e), and polymorphic information content (PIC) were calculated using the computer program Cervus 3.0.7 (Marshall et al. 2015). PIC is a measure of the quality of the informativeness of molecular markers. PIC values of more than 0.5 are considered highly informative for codominant markers such as SSRs (Serrote et al. 2020). In addition, we also evaluated the number of alleles present in only one genotype (private alleles). The observed (H_o) and expected heterozygosity (H_e) were also calculated for each population.

The genetic relationships among the 73 genotypes were evaluated using an unweighted Neighbor-Joining analysis with DARwin 6.0.18 software (Perrier and Jacquemoud-Collet 2005). Dissimilarities (30,000 bootstraps) were calculated using the Dice coefficient (Dice 1945), and the dendrogram was constructed.

The genetic structure of the collection was analysed using the Bayesian method performed with STRUCTURE V2.3.4 software (Pritchard et al. 2000). For each simulation, 20 independent runs were performed for each K from 1-20. A burn-in period of 100.000 followed by 750.000 MCMC (Markov Chain Monte Carlo) repetitions was performed. The Structure Harvester V0.6.94 application (Earl and vonHoldt 2012), implemented using the ΔK method described by Evanno et al. (2005) determined the most relevant K value. The samples with a threshold value of 90% were assigned to a specific group and the remaining samples were classified as admixed. The Q-matrix obtained by the STRUCTURE software, complemented with the ‘tidyverse’ R package (Wickham et al. 2019), was used to generate the bar plots. In addition, the following parameters were calculated with the latter software: average distances (expected heterozygosity) between individuals of the same group and the divergence of allele frequencies between populations (net nucleotide distance).

To assess the proportional differences in genetic variability with and between populations (Slovenia, Austria, Sweden, Great Britain, and Russia), a PERMANOVA analysis was performed using the packages VEGAN (Oksanen et al. 2006) and LATTICE (Deepayan 2008, Sarkar 2008) in R statistical environment (R Development Core Team 2021). We operated with the Euclidian distance parameter and 999 permutations. To plot the results, an NMDS approached was performed by applying ‘vegan’ R package (Venables and Ripley 2002, Oksanen et al. 2006). The betadisper function and the corresponding ANOVA was used to test for differences or similarities in data dispersion with and between the study sites.

Results

Molecular analysis: Microsatellite analysis revealed 86 alleles at 9 microsatellite loci. The number of alleles detected per locus (Tab. 1) ranged from 6 (locus D110 and D118) to 16 (locus B119), with an average of 10 alleles per locus.

Tab. 1. Allele sizes (bp) and allele frequencies (in parentheses) of the 73 Linnaea borealis genotypes at nine microsatellite loci (A5, A102, D110a, D7, D110, D118, A112, B119, C105).

Parameters of genetic variability were calculated for all loci (clones excluded) (Tab. 2). Observed heterozygosity (H_o) ranged from 0.426 (locus D110) to 0.957 (locus C105), with an average of 0.794. Expected heterozygosity (H_e) ranged from 0.604 (locus D110) to 0.906 (locus A112), with an average of 0.768. Differences between observed and expected heterozygosity were examined for all loci studied. The largest difference was observed at locus A5 (0.187) and the smallest at locus D7 (0.007). The average of observed (0.794) and expected (0.768) heterozygosity was quite similar. At 5 of 9 loci (A5, D110a, A112, B119, C105) the observed heterozygosity (H_o) was higher than the expected (H_e), but at four loci (A102, D7, D110 and D118) H_o was lower than H_e. The most informative locus for this group of genotypes was A112, with a polymorphic information content (PIC) of 0.898, and the least informative locus was D110a, with a PIC of 0.553.

Tab. 2. Parameters of genetic variability calculated for different microsatellite loci of the 47 (clones excluded) Linnaea borealis genotypes: n – number of alleles, H_o – observed heterozygosity, H_e – expected heterozygosity, PIC – polymorphic information content.

The parameters of genetic variability were also calculated for each of the 5 populations studied (Tab. 3). The results show the difference between the observed (H_o) and the expected heterozygosity (H_e) in all populations. The lowest value of H_e (0.560) was found in population RU and the highest (0.717) in population SE. The higher expected heterozygosity in the SE population can be explained by the higher number of individuals included (14 different genotypes).

Tab. 3. Parameters of genetic variability calculated for five populations of the 47 Linnaea borealis genotypes: number of genotypes, n – number of alleles, H_o – average observed heterozygosity, H_e – average expected heterozygosity. AT – Austria, SI – Slovenia, SE – Sweden, GB – Great Britain, RU – Russia.

A total of 12 private alleles were found for most loci: A5 and D110 (3 private alleles), A10 (2 private alleles), D110, D118, B119 and C105 (1 private allele). No private alleles were found for loci D7 and A112. Sample AT 1 stood out from the Austrian population with three private alleles. In addition, two private alleles were found in samples SE 12 and GB 10 from Sweden and UK, while five samples had only one private allele (SE 3, SE 6, SE10, GB 2 and GB 4). The population from Sweden not only contained the highest number of genotypes with private alleles, but also the highest number of private alleles per population.

Genetic relationships: The unweighted Neighbor-Joining dendrogram calculated from the microsatellite dataset of L. borealis divided the material into three main clusters (Fig. 3). The first main cluster comprised five subclusters. Subcluster 1 included all samples collected in Slovenia. The second subcluster consisted only of samples from Great Britain. The third and fourth subclusters consisted of samples from Great Britain, Russia and Sweden, and the fifth subcluster contained only samples from Sweden. The second cluster was divided into two subclusters, with the first subcluster containing all samples collected in Austria and the second subcluster containing samples from Sweden. The third cluster contained only samples collected in Russia.

Fig. 3. Neighbor-Joining dendrogram describing genetic relationships among 73 Linnaea borealis genotypes from five locations: AT – Austria, SI – Slovenia, SE – Sweden, GB – Great Britain, RU – Russia.

Although we avoided sampling the same specimens, only eight samples (out of 15) represented distinct genotypes in Slovenia and only six (out of 15) in Austria. In the Sweden population, 14 samples represented distinct genotypes, in Russia 6 out of 15, and in Great Britain 13 out of 13. Looking at the length of the horizontal branches in the dendrogram, it can be seen that, as expected, genetic variability was the highest in the boreal population and diversity was very low in the two Eastern-Southeastern Alpine populations. Surprisingly, the isolated populations from Austria and Slovenia, which are very close to each other, turned out to be completely different and showed a very high degree of variation between the populations.

Genetic structure: the results of the Bayesian analysis revealed that the most relevant number of groups (ΔK) for the analysed data was K = 3 (On-line Suppl. Fig. 1, On-line Suppl. Tab. 2), which divided the material into three groups (Fig. 4).

Fig. 4. Bar plot of Bayesian analysis results (K=3): group 1 (green), 13 genotypes from Austria; group 2 (red) 15 genotypes from Slovenia, and group 3 (blue) representing genotypes from different locations (GB, RU, and SE). Five genotypes were considered admixed.

Using a membership threshold value of > 90%, only 5 (7%) samples were considered admixed, while 68 (93%) samples belonged to identified groups (On-line Suppl. Tab. 3). When the distribution of groups in the material studied was observed, Group 1 (green) consisted of 13 genotypes, Group 2 of 15 genotypes and Group 3 of the highest number of genotypes (40). Group 1 contained material from Austria, Group 2 comprised genotypes from Slovenia, while Group 3 consisted of diverse material from different origins (Great Britain, Russia and Sweden). In the two admixed samples (AT 1 and AT 15) from Austria, the majority of the genome was assigned to Group 1, while the remaining samples (SE 1, SE 3 and SE 10) corresponded to Group 3. Allele frequency divergence between the three groups computed with STRUCTURE software revealed a relatively high proximity between Group 1 and Group 2 (0.0488), while Group 2 and Group 3 appeared more diverged (0.2134), although less than Group 1 and Group 3 (0.2547) (On-line Suppl. Tab. 4). In addition, allelic variation revealed the highest mean distances between individuals for Group 3, while the mean distances for individuals in Groups 1 and 2 were lower (On-line Suppl. Tab. 4).

The PERMANOVA analyses showed significant differences in genetic variability (P < α, α = 0.05) among the observed populations. However, the Euclidean distance-based NMDS plot (Fig. 5) revealed that the RU, SE and GB populations were genetically closer to each other than the AT and SI populations. This was also confirmed by the significant (P < α, α = 0.05) beta-dispersion test. Accordingly, genetic variability was significantly lower in the AT and SI populations. The results showed that the two observed Eastern-Southeastern Alpine populations also differ genetically (Fig. 5).

Fig. 5. Genetic variability with and between population at study sites. AT – Austria, SI – Slovenia, SE – Sweden, GB – Great Britain, RU – Russia.

Discussion

Analysis based on nine microsatellite loci, revealed very low genetic variability within the two isolated Eastern-Southeastern Alpine populations of L. borealis (SI and AT). This was to be expected as only a low number of specimens (individuals) were recorded at each location (sampling site) and therefore only a limited number of distinguishable genets were available. A low number of individuals could also be because of the small size of the remnant. In Slovenia the site of the population is around 0.2 ha, in Austria it is even smaller. However, a very high genetic variability was found between the two populations (SI and AT), even though they are, compared to other sampled populations, relatively closely positioned. This indicates a very low or no dispersal ability; nevertheless, the absence of suitable habitats should not be underestimated. This leads to the conclusion that the two spatially and reproductively isolated Eastern-Southeastern Alpine populations most likely represent glacial relicts, such a pattern having been observed in previous studies of glacial relict plants (e.g. Hensen et al. 2010, Jermakowicz et al. 2017). These researchers explained the high genetic variability between populations by consecutive bottlenecks and long-term reproductive isolation. In addition, Vogler and Reisch (2013) claimed that genetic variability may be high in the case of different populations of glacial relicts and low in the case of postglacial colonization. Low genetic diversity within the population might be also explained as a shift in mating system from outcrossing to selfing (Surina et al. 2023).

In a comprehensive study by Wróblewska (2013), a phylogeographic structure derived from plastid DNA and a moderate genome-wide diversity estimated from AFLP markers were found for L. borealis in Eurasia. Six haplotypes were identified in an area from Scotland to Altai and from Norway to Italy. However, Wróblewska (2013) concluded that although half of the populations studied were highly isolated, they still showed similar levels of genetic diversity across the geographic range. Furthermore, she stated that she had found no support for the hypothesis that a bottleneck and/or inbreeding. factors affecting genetic diversity, were associated with habitat fragmentation. The case of L. borealis clearly shows that it was not habitat fragmentation, but rapidly changed climatic conditions after the last glaciation that caused the narrow isolation of the species leading to glacial relict populations, like the one in Slovenia and Austria. Regarding the geographical distribution of Alpine populations, she only sampled in the Western Alps (Aosta, Italy; Les Allues, France) and the Central Alps (Zernez, Switzerland). However, according to Meusel and Jäger (2017), its distribution is not that rare. East of the line Innsbruck–Wipp Valley–Brenner Pass–Eisack Valley–Etsch Valley there are only scattered small populations, especially in the Hohe Tauern mountain range. From there, L. borealis has not been found in the rest of the Eastern Alps, with the exception of three localities on the south-eastern borders of the Alps (Ernet and Franz 2011). This could be the reason why four out of six haplotypes occur at the above-mentioned locations in the Western and Central Alps.

In any case, the warming after the end of the last glaciation led to a loss of suitable habitats, resulting in a more or less strong fragmentation of the Alpine population, which eventually led – at least in the Eastern Alps – to only a few isolated patches remaining in specific habitat conditions related to elevation, northern exposure or particular situations in terms of substrate conditions. Ernet and Franz (2011) mentioned the boulder-strewn slopes where cold air typically emerges from under the scree in summer. This is exactly the characteristic of the Austrian location in the Gurktaler Alps, but is even more pronounced at the Slovenian site, where at the extremely low elevation of only 495 m a.s.l., from holes in the rubble ice-cold air blows during the growing season. There is no direct sunlight on the steep northern slopes, which also increases soil and air humidity. According to Wraber (1963), therefore, no zonal beech forest has developed at this location, but a community dominated by Rhododendron hirsutum L., Pinus sylvestris L. and Larix decidua Miller.

To conclude, it is well known that the most severe negative impact of fragmentation is caused by genetic drift, increased inbreeding, and reduced gene flow among populations (Scobie and Wilcock 2009, Wiberg et al. 2016). This led to a reduced number of offspring and a reduction in population fitness (Scobie and Wilcock 2009). In the case of L. borealis which is a self-incompatible plant, the effect of isolation on reproductive success is even greater. For successful fertilization, self-incompatible plants like L. borealis require cross-pollination with a compatible mate, which does not belong to the same clone. It has been reported that in extreme cases where all remaining plants have the same alleles for genetic control of pollen-pistil interactions seed set can be reduced to zero (Demauro 1993). However, in small populations of many clonal, self-incompatible plants such as Rubus saxatilis L., Aster furcatus E.S.Burgess, Calystegia collina (Greene) Brummitt, Arnica montana L., Ranunculus reptans L. and Maianthemum bifolium (L.) F.W.Schmidt (Scobie and Wilcock 2009), the limited availability of partner plants is increasingly being identified as the cause of reproductive failure, the plants perhaps adapting to turning to selfing.

Assuming that the two isolated populations of L. borealis in the Eastern-Southeastern Alps are glacial relicts with considerably reduced genetic variability, further research on their viability and reproductive mating system is needed to ensure better conservation success.

Author contribution statement

All authors contributed to the study conception and design. All authors read and approved the final manuscript.

Data availability

The datasets generated during and/or analysed during the current study are available from the corresponding author on reasonable request.