Glacial refugia and postglacial expansion of the alpine–prealpine plant species Polygala chamaebuxus

Abstract The shrubby milkwort (Polygala chamaebuxus L.) is widely distributed in the Alps, but occurs also in the lower mountain ranges of Central Europe such as the Franconian Jura or the Bohemian uplands. Populations in these regions may either originate from glacial survival or from postglacial recolonization. In this study, we analyzed 30 populations of P. chamaebuxus from the whole distribution range using AFLP (Amplified Fragment Length Polymorphism) analysis to identify glacial refugia and to illuminate the origin of P. chamaebuxus in the lower mountain ranges of Central Europe. Genetic variation and the number of rare fragments within populations were highest in populations from the central part of the distribution range, especially in the Southern Alps (from the Tessin Alps and the Prealps of Lugano to the Triglav Massiv) and in the middle part of the northern Alps. These regions may have served, in accordance with previous studies, as long‐term refugia for the glacial survival of the species. The geographic pattern of genetic variation, as revealed by analysis of molecular variance, Bayesian cluster analysis and a PopGraph genetic network was, however, only weak. Instead of postglacial recolonization from only few long‐term refugia, which would have resulted in deeper genetic splits within the data set, broad waves of postglacial expansion from several short‐term isolated populations in the center to the actual periphery of the distribution range seem to be the scenario explaining the observed pattern of genetic variation most likely. The populations from the lower mountain ranges in Central Europe were more closely related to the populations from the southwestern and northern than from the nearby eastern Alps. Although glacial survival in the Bohemian uplands cannot fully be excluded, P. chamaebuxus seems to have immigrated postglacially from the southwestern or central‐northern parts of the Alps into these regions during the expansion of the pine forests in the early Holocene.

However, the ecological requirements of plant species have a strong impact on their glacial and postglacial history and different hypotheses about the migration and survival of plant species during Quaternary can therefore be proposed for species with different ecological preferences (Holderegger & Thiel-Egenter, 2009;Vargas, 2003).
Many temperate species, originally occurring in Central Europe, became extinct during the Quaternary ice ages and retreated to southern refugia and survived glacial maxima. In contrast, high-alpine species even persisted in central refugia on ice-free mountain tops, so called "nunataks." Less cold resistant alpine species survived either in refugia at the periphery of the Alps or may have migrated to lowland areas.
Knowledge about the vegetation of these lowlands between the Scandinavian and Alpine ice sheet in Central Europe during glaciation is yet scarce. There are stratigraphic records of pollen and macrofossils for Salix herbacea, Betula nana, Dryas octopetala, or Koeningia islandica, whereas dwarf shrubs counting among Ericaceae played an unexpectedly subordinate role (Burga, Klötzli, & Grabherr, 2004;Lang, 1994).
The shrubby milkwort (Polygala chamaebuxus) is an endemic European species with a remarkably broad ecological niche and a wide distribution range including the Alps but also Central European mountain ranges like the Franconian Jura or the Bohemian uplands. There, it occurs mainly in pine forests and on rocky mountain slopes. In the study presented here, we tried to illuminate the origin of the species in these lower mountain regions. More specifically, our aim was (i) to identify glacial refugia of P. chamaebuxus and (ii) to analyze whether the populations of the species in the low mountain ranges can be attributed rather to glacial survival or to postglacial immigration.  (Calvo, Hantson, & Paiva, 2014;Lorite, Peňas, Benito, Caňadas, & Valle, 2010). In addition, the subgenus includes one species which is restricted to northwestern Africa: P. munbyana Boiss. & Reut.

| Species description
Based on karyological and palynological studies (Merxmüller & Heubl, 1983), it was suggested that P. munbyana (2n = 14) belongs to the diploid level, P. webbiana, P. balansae, and P. vayredae are tetraploids with 2n = 28, whereas hyperhexaploidy (2n = 44) was found in P. chamaebuxus. Karyotype analysis revealed that P. chamaebuxus developed most probably by autopolyploidy from P. vayredae or the African P. webbiana or by allopolyploidy of these species. The evolution of the group concerned seems to have taken place in the southwestern Mediterranean and to have continued on the Iberian way as far as the Alps and Central Europe (Merxmüller & Heubl, 1983).
In contrast to the Iberian taxa which are narrow endemics, P. chamaebuxus L. has the largest and northernmost distribution range of all members. It occurs in the Alps, the northern Apennine, the northern parts of the Dinaric Mountains, and in parts of the prealpine moraine landscape as well as some in low mountain ranges like such as Jurassic mountains, the Bavarian Forest, the Fichtelgebirge, and the Bohemian uplands (Sebald, Seybold, Philippi, & Wörz, 1998). A white flowered form of P. chamaebuxus occurs, most probably, over the whole distribution range, whereas a red flowered form (var grandiflora Gaudin; var rhodoptera Ball) can only be found in the cantons of Graubünden and Tessin and down the Apennine (Meusel, Jäger, Rauschert, & Weinert, 1978).
Polygala chamaebuxus is a 5-to 30-cm-high dwarf shrub. Full flowering occurs in spring and early summer. The species is, like the closely related species P. vayredae (Castro, Loureiro, Ferrero, Silveira, & Navarro, 2013;Castro, Silveira, & Navarro, 2008), insectpollinated, allogamous, and self-incompatible (Hegi, 1986;Jauch, 1917). Polygala chamaebuxus exhibits a broad ecological range. It grows in open forests, mainly pine woods, among rocks and mountain slopes. According to phytosociological classification, this taxon is together with Erica carnea a characteristic element of the order Erico-Pinetalia. In the Alps, it reaches up to 2,650 m above sea level in Graubünden (Braun-Blanquet & Rübel, 1932), and at Monte Baldo, it can be found from 80 m above sea level up to 2,100 m altitude (Prosser, Bertolli, & Festi, 2009). It grows predominantly on calcareous soil types but also some populations on more acidic soils have been reported. Polygala chamaebuxus is a medium shade plant and the light supply seems to be one of the most important factors, which is strongly influenced by the surrounding vegetation (Gauckler, 1938). Therefore, it occurs predominantly in sparse pine woods, dry oak forests, as well as on calcareous low-nutrient meadows (Sebald et al., 1998).

| Study design and sampling of plant material
For the study presented here, plant material was sampled from 30 populations (Table 1, Figure 1) covering continuously almost the entire range of P. chamaebuxus. When possible, within populations, ten samples were taken with a minimum distance of ten meters following a transect to avoid double sampling of the same individual.

| AFLP analysis
For AFLPs, the DNA was extracted from the dried sampling material following the CTAB protocol from Rogers and Bendich (1994) adapted by Reisch and Kellermeier (2007). After photometrical measurement of the concentration, solutions were diluted with water to 7.8 ng/ μl and were subsequently used for AFLPs, which were conducted in T A B L E 1 Geographic location of the studied Polygala chamaebuxus populations with number, population code, name of the location as well as geographic longitude (Long.), latitude (Lat.) and altitude.
Populations were numbered across the distribution range from west to east and north to south accordance with the protocol of Beckmann Coulter as described before (Bylebyl, Poschlod, & Reisch, 2008;Reisch, 2008).
After an initial screening of 30 primer combinations, three of them were chosen for the subsequent selective PCR reaction using labeled and fivefold (D4) with 1× TE 0.1 buffer for AFLP, while the D3 products remained undiluted. Subsequently, 5 μl of each of the diluted PCR products of a given sample was pooled and added to a mixture of 2 μl sodium acetate (3 mol/L, pH 5.2), 2 μl Na 2 EDTA (100 mmol/L, pH 8), and 1 μl glycogen (Roche). DNA was precipitated in a 1.5-ml tube by adding 60 μl of 96% ethanol (−20°C) and 20-min centrifugation at 14,000 × g at 4°C.
The supernatant was poured off, and the pellet was washed by adding 200 μl 76% ethanol (−20°C) and centrifugation at the latter conditions.
The pelleted DNA was vacuum dried in a vacuum concentrator.
Subsequently, the pellet was dissolved in a mixture of 24.8 μl Sample Loading Solution (SLS, Beckman Coulter) and 0.2 μl CEQ Size Standard 400 (Beckman Coulter) and subsequently selective PCR products were separated by capillary gel electrophoresis on an automated sequencer (CEQ 8000, Beckmann Coulter).
Results were examined using the CEQ 8000 software (Beckman Coulter) and analyzed using the software Bionumerics 6.6 (Applied Maths, Kortrijk, Belgium). In order to assess the reproducibility of the scored fragments, about 10% (29 samples) of all analyzed samples were repeated and the genotyping error rate (Bonin et al., 2004) was estimated, which was 4.8%.

| Statistical analysis
Using the resulting binary matrix, genetic variation within populations was determined applying the program PopGene 1.32 (Yeh, Yang, Boyles, Ye, & Mao, 1997) as percentage of polymorphic bands F I G U R E 1 Genetic variation within the studied populations, measured as AMOVA-derived SSWP/n − 1 values (SSWP) and rarity index (DW). Circle diameter and color indicate the degree of genetic variation. The dotted line marks the area with high levels of genetic variation and rarity within populations in the center of the distribution range PB and Nei's gene diversity H = 1 − Σ(p i )². Additionally, we calculated rarity as frequency down weighted markers (DW) for each population  with AFLPdat in R (Ehrich, 2006). Therefore, we randomly chose eight individuals per population in five iterations.
A Bayesian cluster analysis using 10,000 Markov chain Monte Carlo (MCMC) simulations was computed with 20 iterations per K = 1-31 and a burning period of 10,000 with the software Structure 2.3.3 (Pritchard, Stephens, & Donelly, 2000). The most probable number of classes was calculated (Evanno, Regnaut, & Goudet, 2005), and the mean probability of the individuals of each population to be assigned to the respective classes was calculated over all 20 repeats for the most probable number of classes.
Furthermore, a nonhierarchical AMOVA was carried out with GenAlEx 6.41 (Peakall & Smouse, 2006) based on pairwise Euclidian distances to assess the variation within and among populations. This also yielded pairwise PhiPT values as well as the SSWP value (sum of squares within population) for each population. Dividing the latter value through the number of individuals reduced by one, provided the sample size-independent measure of variation SSWP/(n − 1).
A Mantel test was performed to analyze whether the genetic distances and the geographic distances between populations were correlated (Mantel, 1967).
Finally, we used PopGraph (Dyer & Nason, 2004)   However, populations from the northeastern part of the distribution range were mainly assigned to one group, while the populations from the southwest and the southeast were more frequently classified in two other groups.
In a nonhierarchical analysis of molecular variance (AMOVA), only 16.5% of the total genetic variation was found among all populations while 83.5% were detected within populations (Table 3). The overall Φ PT was therefore 0.17. Variation between the groups detected in the Bayesian cluster analysis was significant but with only 3% very low. Similarly, molecular variance between the northeastern group on the one hand and the southeastern and southwestern group on the other hand was only 4% and, therefore, also very low. A Mantel test showed a significant correlation of the genetic variation between populations obtained from the AMOVA (Φ PT ) and the respective geographic distance between populations (r = .570, p < .001).
In the PopGraph genetic, network populations were highly interconnected ( Figure 3). However, the populations from the northern group detected in the Bayesian cluster analysis were more closely related to the populations from the southwestern than to the populations from the southeastern group. One of the most variable populations also containing a higher number of rare fragments (population LV) was completely separated from the network. T A B L E 3 Results of the conducted analyses of molecular variance (AMOVA). We calculated variation between all populations (1), between the three groups derived from the Bayesian cluster analysis (2) between the northern group and the western (3) and eastern group (4)

| Genetic variation of Polygala chamaebuxus in the context of life history traits
It has already been demonstrated that life history traits have a strong impact on genetic variation within and between populations. In particular, life span, frequency, and mating system are of outstanding importance for genetic variation (Nybom, 2004;Reisch & Bernhardt-Römermann, 2014). The genetic variation within populations of P. chamaebuxus observed in our study (H = 0.21) was comparable to the variation recently reported for other long-lived, common, and outcrossing plant species (H = 0.20) using AFLPs (Reisch & Bernhardt-Römermann, 2014). The results of our study match, from this point of view, the findings of the preceding reviews.
In contrast to our expectations, we observed, however, only a low level of genetic variation between populations of P. chamaebuxus.
Previously, for long-lived, common, and outcrossing plant species, a mean Φ PT of 0.20-0.34 was reported (Reisch & Bernhardt-Römermann, 2014). As genetic variation depends on life history traits, the comparison of single species with differing traits is always delicate. Nevertheless, many alpine species exhibited even higher levels of genetic differentiation (Schönswetter et al., 2004;Vogler & Reisch, 2013). With a Φ PT of only 0.17 between all populations across the whole distribution range, P. chamaebuxus exhibited only a weak geographic pattern of genetic variation. This suggests a comparatively short period of isolation during the glaciations and rather broad waves of postglacial recolonization as discussed more detailed below.

| Glacial refugia and postglacial recolonization
Following our data, especially the high level of rarity, suggests long-term survival of P. chamaebuxus in the Southern Alps between Switzerland and Italy. This area has already been identified as refugium for other calcicolous, subalpine to lower alpine plant species in previous studies . Another putative refugium of P. chamaebuxus has probably been located in the middle part of the northern Alps, where we also observed a higher number of rare fragments. The occurrence of P. chamaebuxus along the northern margin of the Alps at least during the last interglacial (Eemian) has been proved by fossil evidence (Murr, 1926;Wettstein, 1892) and previous studies have already postulated glacial refugia at the northern edge of the Alps (Schönswetter, Stehlik, Holderegger, & Tribsch, 2005;Stehlik, 2003), which supports the assumption that P. chamaebuxus could have survived glaciations also in this region.
However, our results indicate rather a genetic continuum than deep genetic splits between populations of P. chamaebuxus, which may be a sign of a comparatively short period of isolation during the LGM. It is known that the strong glaciations of the Würm glaciation were limited to few periods of extreme cold climate with culmination during the LGM (Veit, 2002 range, which allows the species to grow under various climatic conditions and is even considered as cold germinator (Jäger, 2011). Polygala chamaebuxus may, for this reason, have been affected not that strongly by the glaciations like other highly specialized species. It is possible that the refugia described above were locations where the species

| Glacial survival in the lower mountain ranges or not?
The populations of P. chamaebuxus in the lower mountains of Central This assumption is supported by previous studies reporting glacial survival of forest-related plant species in cryptic refugia located in the lower Central European mountain ranges (Kramp et al., 2009;Michl et al., 2010;Slovák et al., 2012;Tyler, 2002), although some studies also revealed ambiguous results (Dvořáková et al., 2010;Fér et al., 2007). Kramp et al. (2009)  Similarly, it is assumed that the boreo-montane tall forb Cicerbita alpina survived glaciations in sheltered pockets with a humid climate in some parts of Central Europe (Michl et al., 2010) and that Cyclamen purpurascens may also have survived glaciations in prealpine northern refugia (Slovák et al., 2012). For the woodland grass Melica nutans, several independent "strongly restricted and isolated" refugia in Central Europe have been detected (Tyler, 2002 In the PopGraph genetic network, the populations from the lower mountain regions were more closely related to the populations from the western part than to the populations from the eastern part of the distribution range. This suggests that P. chamaebuxus may have immigrated postglacially from the southwestern or central-northern part of the Alps to the lower mountains of Central Europe. This migration process of P. chamaebuxus to the lower mountain regions may be associated with the expansion of pine forests after the last LGM. It is assumed that Pinus sylvestris survived glaciations on the Iberian and the Balkan Peninsula (Sinclair, Morman, & Ennos, 1999;Soranzo, Alia, Provan, & Powell, 2000;Wójkiewicz & Wachiowak, 2016). However, cryptic northern refugia have also been postulated for Scots pine (Kinloch, Westfall, & Forrest, 1986;Stewart & Lister, 2001), similar to the herbaceous forest species mentioned above. Whereas the Iberian populations are considered as relicts, Central Europe and Scandinavia were recolonized postglacially from the Balkan (Wójkiewicz & Wachiowak, 2016). From there, pine forests spread in the early postglacial phases and covered large parts of the alpine forelands and Central Europe (Lang, 1994). Polygala chamaebuxus is considered as a species typically for these early pine forests (Hardtke & Ihl, 2000) and still occurs today in this type of habitat (Gauckler, 1938). The widely distributed postglacial pine forests seem to have provided well conditions for a broad and continuous co-migration of P. chamaebuxus together with Scots pine toward the north. Migration could already have been started in the Late Glacial from 15,000 BP to 10,000 BP as pine and birch were already present in the Alps and the alpine forelands until about 8,000 BP when the continuous distribution of pine forests ended (Lang, 1994;Veit, 2002).
Similarly, the species seems to have migrated from the center of the distribution range to the eastern and western Alps. In this context, it is a remarkable finding of our study that the population from the Velebit in Croatia was more closely related to the population from the Apennine and westernward populations than to the populations from the nearby southeastern Alps. This observation was also made for Saxifraga paniculata in a previous study (Reisch, 2008) and seems to be linked to the desiccation of the Adriatic during glaciation, which seems to have alleviated migration processes.

ACKNOWLEDGMENTS
Special thanks go to Sabine Fischer for the design of the maps, Daniela Listl for assistance with PopGraph, Petra Schitko for help in the laboratory, Franziska Kaulfuß and Josef Simmel for fruitful discussions, and Peter Poschlod for his generous support.