Introduction

Quaternary climate oscillations have caused entire plant zones to shift (Huntley and Birks 1983; Comes and Kadereit 1998; Tzedakis et al. 2013). The geographical ranges of many species have expanded and contracted in accordance with the climatic changes (Bennett et al. 1991; Hewitt 2004). The glaciated areas covered the lowlands of northern and central Europe and the high mountain ranges in the south, mainly the Alps and the Pyrenees, as well as part of the Carpathians (Ehlers and Gibbard 2004; Birks 2008). Contrary to what is expected, conditions during the last glaciation in Europe might not have given the high-mountain, alpine-adapted species, which are cold tolerant but strictly habitat-limited, the opportunity to enlarge their ranges (Stewart et al. 2010). Cooling might have increased the area that was available to some species and allowed them to expand. This is probably the case with some tundra plants, which had a larger range in the periglacial tundra than in the present period (e.g. Skrede et al. 2006). The relatively mild climate in the lowlands south of the Alps during the glacial periods allowed tree species to grow and persist but these conditions were likely not suitable for alpine species (Ravazzi 2002; Schönswetter et al. 2005; Stewart et al. 2010). Extensive areas to the north and east of the Alps also were likely not habitable for many high elevation plants because of the steppe vegetation that prevailed there under the cold and dry climate during the glacial period (Grichuk 1992; Schönswetter et al. 2005). Along with the expansion of the glaciers in the cold Quaternary phases, the habitats of the European high-mountain plants were becoming increasingly limited. Consequently, another type of glacial refugia must have existed to support their survival during these periods.

In the period preceding the wide use of molecular markers, our knowledge regarding the refugia that could host alpine species was largely speculative. Schönswetter et al. (2005) defined the potential refugial zones in the Alps based on the paleoenvironmental and geological data and subsequently tested the existence of these zones by comparison with the phylogeographic patterns of 12 mountain species. This combined approach allowed major refugial areas to be identified (Fig. 1). However, some doubts still remain about the existence of refugia in the middle section of the alpine massif, which was on ice-free mountaintops that likely protruded above the glaciers. Only organisms adapted to severe high-elevation conditions could survive there. This so-called ‘nunatak hypothesis’ has gained some support with molecular evidence suggesting the glacial persistence of some high-mountain plant and animal species in nunatak refugia (Schönswetter and Schneeweiss 2019). In particular, some long-lived plant species, frequently diploid, that had adapted to cold conditions and with a high ability for the clonal growth might have survived glaciation periods in the peripheral refugia or on nunataks (Favarger 1972; Nève and Verlaque 2010). In terms of genetic patterns, isolated nunatak populations should presumably be characterized by clear differentiation and distinctive diversity that may manifest with private alleles/haplotypes that are restricted to refugial populations (Schönswetter and Schneeweiss 2019). However, these genetic signatures of nunatak survival may be difficult to detect because of the genetic swamping that is likely to occur during (re)colonization process (Gabrielsen et al. 1997; Schneeweiss and Schönswetter 2011; Schönswetter and Schneeweiss 2019).

Fig. 1
figure 1

Sampled populations of Salix serpillifolia (black dots) and refugial areas in the Alps, as defined by Schönswetter et al. (2005): I–IV: southern-Alpine (brown areas), V–VII: northern-Alpine (blue areas), VIII–IX: central alpine (purple areas); red asterisks: sites with macrofossils of Salix retusa agg

Full size image

The fossil record of the Alpine flora is scarce due to the lack of conditions suitable for the accumulation of macroremains at the highest elevations (Lang 1994). The available evidence does not provide definite support for continuous in situ existence of high-elevation species throughout the time with the greatest ice-sheet area. Paleodistribution models can fill this gap by supplying an independent test of the predictions of genetic models. Modelling of past climates and species distributions may provide complementary insight into their complex histories (Gavin et al. 2014). Species distribution models (SDMs) assume that the climatic niche of a species can be based on present-day observations and that it remains the same over the time, which enables the climate-based reconstruction of the potential distributional range of the species in the past using paleoclimatic data (Wiens and Graham 2005). However, since SDMs relay on climatic variables only, the results gained via SDMs should be tested against the other sources including genetic data and fossil records.

Salix serpillifolia Scop. is a very resilient species and was described by Ellenberg (1988) as “the only woody plant to feel at home in the region of eternal snow”. This slow-growing, dioecious, diploid, procumbent dwarf willow is endemic to the Alps (Dostál 1989; Rechinger 1992; Skvortsov 1999; Kosiński et al. 2019). It inhabits severe alpine vegetation zone and even ledges in the subnival belt, up to 3200 m a.s.l., where it is found mainly on a base-rich substrate in the habitats of grassland communities, screes and rock crevices. Despite its small size, it is supposed that it can live for over 100 years (Ellenberg 1988). It forms prostrate mats on exposed mountain slopes and ridges and produces creeping stems with the capacity of rooting, and thereby, clonal growth (personal observations). Fragmentation of plant material and re-rooting further down the slope increases the probability of its survival in a harsh and unstable environment.

In the present study, we investigated the spatial genetic structure of the Alpine endemic woody plant, S. serpillifolia. So far, most studies dedicated to alpine or subalpine calcicolous species has been focused on the herbaceous plants (e.g. Stehlik et al. 2002a, b; Stehlik 2002; Naciri and Gaudeul 2007; Schönswetter et al. 2009; Schneeweiss and Schönswetter 2010; Schönswetter and Schneeweiss 2019). Alpine zones are important areas of global biodiversity because they host a high diversity of cold-adapted specialists, including many endemic species. Nevertheless, global warming is expected to cause an upward shifting of plant species in the high mountain systems and, consequently, the habitats available for the alpine and nival plants may be significantly limited leading to extinction of some of them (Rangwala and Miller 2012; Steinbauer et al. 2018; Asse et al. 2018). Thus, population genetic data gained for alpine species may be useful not only in determining historical events but also in predicting the ability of taxa to respond to current climate change.

In this study, genetic analyses were used with the aim of answering the following questions: (1) What is the pattern of geographic distribution of the genetic diversity of S. serpillifolia? (2) Is there any congruence between this pattern and those reported for other European obligate high-mountain plant species? (3) Is there evidence for nunatak survival, which is expected to be plausible for this species given its current habitat requirements?

Materials and methods

Sampling, DNA extraction and genotyping

Seventeen populations of S. serpillifolia were sampled across the Alps (Table 1; Fig. 1). Leaves were collected from plants that were separated by at least 20 m, as willows from the S. retusa group reproduce vegetatively by the rooting of the stems (Rechinger 1957; Hörandl et al. 2002; Koblížek 2006). Samples of these populations were deposited at the Herbarium of the Institute of Dendrology, Polish Academy of Sciences (KOR). Genomic DNA was extracted from dried leaf tissues following a CTAB protocol (Doyle and Doyle 1987).

Table 1 Sampled populations of Salix serpillifolia; N—number of sampled individuals (nSSR/cpSSR)
Full size table

Twelve populations with a sufficient number of sampled individuals were used for the screening of polymorphisms at the nuclear microsatellite loci (nSSRs) (Table 1). Initially, a total of 30 nSSRs markers originally designed for Salix sp. (Barker et al. 2003; Stamati et al. 2003) and Populus sp. (Tuskan et al. 2004; Carvalho et al. 2010) were tested. Then, a set of ten markers was selected (Online Resource 1), as they gave repeatable, high-quality amplification products; these were gSIMCT024, gSIMCT035, gSIMCT041, gSIMCT052 (Stamati et al. 2003), ORPM-028, ORPM-030, ORPM-202 (Tuskan et al. 2004), ASP322 (Carvalho et al. 2010), GCPM 1255, and GCPM 1812 (Zeng et al. 2016). The PCR was performed in two multiplexes in a total volume of 10 μl using Qiagen Multiplex PCR Kit (Qiagen, Hilden, Germany) according to the manufacturer’s instructions. Amplification was conducted according to the thermal profile presented in Tuskan et al. (2004) and Stamati et al. (2003) (Online Resource 1).

All sampled populations were screened for their chloroplast microsatellite polymorphisms (cpSSRs) using three primer pairs, CSU01, CSU05 and CSU07, developed for S. reinii (Lian et al. 2003). The samples were amplified in PCR reactions in a total volume of 10 µl using Qiagen Multiplex PCR Kits (Qiagen, Hilden, Germany) under conditions recommended by the manufacturer.

PCR products of both nSSRs and cpSSRs were separated on an Applied Biosystems 3130 Genetic Analyzer (Thermo Fisher Scientific, Waltham, USA) with an internal size standard, GeneScan™ 500 LIZ®. Genotypes were scored manually using GeneMapper 4.0 software (Thermo Fisher Scientific, Waltham, USA). Raw2Gen software (I.J. Chybicki, Kazimierz Wielki University, Poland, unpubl. data) was used to bin the observed nSSR and cpSSR allele sizes into representative discrete alleles according to the method of Idury and Cardon (1997).

Data analysis

Nuclear SSR genetic diversity and differentiation

Linkage disequilibrium among the nSSR loci for each population was estimated using FSTAT 2.9.4 with 1500 permutations (Goudet 1995, 2003). Deviation from Hardy–Weinberg equilibrium was tested using an exact test implemented in GenePop v. 4.0 (Raymond and Rousset 1995). The multilocus parameters of genetic diversity for populations, including average number of alleles (A), private alleles per population (PA), average number of effective alleles (AE), allelic richness (AR), expected heterozygosity (HE), observed heterozygosity (HO) and inbreeding coefficient (FIS) were estimated using GenAlEx v.6.2 (Peakall and Smouse 2012). The null alleles frequency (Null) and multilocus inbreeding coefficient including ‘null alleles’ correction (FISNull) were estimated following an individual inbreeding model (IIM) implemented in INEst 2.2 (Chybicki and Burczyk 2009). Using the Deviance Information Criterion (DIC), we compared the full model (‘nfb’; null alleles, and genotyping failures, FIS > 0) with a random mating model (‘nb’; null alleles and genotyping failures, under an assumption of FIS = 0) to assess the factors affecting the homozygosity level in populations. The analysis was run with 500,000 MCMC cycles, with every 200th cycle sampled and a burn-in of 50,000.

The overall genetic differentiation among populations, Wright’s fixation index (FST) (Wright 1931), was assessed using GenAlEx and the significance of FST was calculated with 10,000 permutations. Additionally, the ENA correction method excluding null alleles (FSTENA) was applied to re-calculate multilocus FST after 9999 replicates in FreeNA, as null alleles may overestimate the value of FST (Chapuis and Estoup 2007). Finally, pairwise population comparisons based on FST were performed in Arlequin v. 3.1 and tested with a permutation test.

Evidence of a recent bottleneck was evaluated using the M-ratio method (Garza and Williamson 2001) implemented in INEst 2.2. The mean ratio was calculated as M = k/r, where ‘k’ is the total number of alleles and ‘r’ is the overall range in allele size for each population. The M-ratio value is expected to be lower in populations that have experienced a bottleneck than in populations in mutation-drift equilibrium (Garza and Williamson 2001). The analysis was run under the two-phase mutation model (TPM) with the proportion of multi-step mutations (pg) set as 0.22 and the mean size of multi-step mutations (δg) as 3.1 according to recommendations of Peery et al. (2012). The reduction in the M-ratio was verified by comparing the mean observed M-ratio (M) to the simulated M-ratio generated under mutation-drift equilibrium (MReq). To evaluate the significance of a deficiency in the M-ratio, the Wilcoxon signed-rank test based on coalescent simulations consisting of 106 iterations was applied.

Chloroplast SSR genetic diversity

Haplotypes were defined as unique combinations of size variants across the cpSSRs. HAPLOTYPE ANALYSIS 1.05 (Eliades and Eliades 2009) was used to calculate within-population genetic diversity, number of different haplotypes (A), number of private haplotypes (PH), effective number of haplotypes (NE), haplotypic richness (HR), Nei’s index of genetic diversity (HE), mean genetic distance between individuals (\( D_{\text{sh}}^{2} \)) (Goldstein et al. 1995) and genetic differentiation (FST) among populations (Finkeldey and Murillo 1999).

Geographic patterns of differentiation

A model-based Bayesian clustering approach was used to investigate the geographic pattern of genetic differentiation at nSSRs among populations using STRUCTURE 2.3.3 (Pritchard et al. 2000). The analysis was performed based on a non-spatial admixture model with correlated allele frequencies without prior information of population location. Ten independent repetitions for each number of groups (K) ranging from 1 to 13 were performed with a burn-in of 105 steps, followed by 2 × 105 MCMC iterations. The alignment of the results across ten replicates of analyses was assessed using CLUMPAK (Kopelman et al. 2015) and the best number of clusters was determined according to Evanno’s ΔK method (Evanno et al. 2005). To infer the population genetic structure at linked cpSSRs, a spatial genetic mixture model implemented in BAPS was used (Corander et al. 2008). The optimal number of K clusters was used with the clustering of groups option in which populations with known coordinates are used as clustering units. After the testing stage, an analysis was run for ten replicates of K ranging between 1 and 18.

To identify potential genetic barriers among populations Monmonier’s maximum-difference algorithm implemented in BARRIER 2.2 was used (Manni et al. 2004). The analysis was conducted using Reynold’s distance matrix and the geographical coordinates of sample sites. The significance of barriers was evaluated by bootstrapping with 1000 replications of 1000 distance matrices generated by the SEQBOOT and GENDIST programmes of PHYLIP package 3.696 (Felsenstein 1989).

The isolation-by-distance (IBD) analysis was performed in GenAlEx by testing the correlation between the matrix of pairwise Slatkin’s linearized FST values (Slatkin 1995) with ENA correction and the pairwise geographical distance matrix using the Mantel test, and the significance of IDB was assessed with 1000 permutations.

Species distribution modelling

To evaluate the effect of past and present climatic conditions on the geographic range of S. serpillifolia, species distribution modelling that was based on the maximum entropy approach implemented in MaxEnt 3.4.1 was applied (Phillips et al. 2004; Elith et al. 2011). The climatic information available in the WorldClim database (Hijmans et al. 2005) was used to construct the models of the potential distribution of this species at the present time (1960–1990; PRE) and at the last glacial maximum (c. 22,000 year BP; LGM). To construct the models for the current species range, 19 bioclimatic variables were used at a 30 arc-sec resolution (Hijmans et al. 2005) (Online Resource 2). The LGM layers were based on the Community Climate System Model (CCSM3) at a 2.5 arc-min resolution (Collins et al. 2006). Analyses were performed as a bootstrap with 10 replicates. The ‘random seed’ option was applied, and 20% of the data were used to be set aside as the test points. The maximum number of iterations was set to 10,000 and the convergence threshold to 0.00001. To evaluate the model accuracy, the receiver operating characteristic (ROC) curve and value of area under the curve (AUC) were used. An AUC value close to 0.5 indicates that the model is weak, and a value close to 1.0 suggests that the model has excellent performance. The list of species occurrence records that was used to model the range of S. serpillifolia was based on the data presented in Kosiński et al. (2017) and the ZOBODAT database (http://www.zobodat.at) (Online Resource 3).

Results

Genetic diversity

Nuclear microsatellites

Detailed information on the polymorphisms of the nSSR loci that were used in this study is presented in Online Resource 1. In total, 150 alleles were identified among 351 individuals at ten nSSR loci. SSR loci showed a wide level of variability ranging from 5 (gSIMCT041) to 29 (gSIMCT052) alleles with an average of 15 alleles per locus. Null alleles were detected in all loci, with an average frequency approaching 0.082 (Online Resource 1).

Genetic variation at the nSSR loci of the analysed populations of S. serpillifolia was relatively low, with the expected heterozygosity reaching an average HE = 0.595 and average number of alleles A = 6.758 (Table 2). The highest HE was found in population 12 (HE = 0.723) and the lowest in population 8 (HE = 0.481). Allelic richness (AR) was comparable in most of the populations, except for in populations 12 and 9, in which the highest values of AR were noted, reaching 8.711 and 8.191, respectively (Table 2). Private alleles were noted in all populations, but the highest number of private alleles (9 and 10) was observed in populations 9 and 12, respectively. Null alleles were observed in all populations, and their frequency ranged from 0.041 to 0.127 with an average of 0.086. A significant excess of homozygotes (p < 0.001) was detected in all populations and generally high values of inbreeding coefficient (FIS) were noted, ranging from 0.114 to 0.358 (0.250 on average). The FIS value that was estimated with a ‘null alleles’ correction (FISNull) was lower in most of the populations (0.069 on average), suggesting that the observed deficiency of heterozygotes may result not only from inbreeding but also from the presence of null alleles, especially in populations 1, 4, 8 and 14, where only null alleles had a substantial impact on FIS estimation (Table 2).

Table 2 Summary of genetic variability estimated for populations of Salix serpillifolia across ten nSSR loci; population abbreviations as in Table 1
Full size table

Estimations of the global level of genetic differentiation with and without the ENA correction were almost the same and attained values of FSTENA = 0.074 (CI 0.054–0.107) and FST = 0.075 (CI 0.056–0.103), respectively. The population pairwise FST ranged from 0.0086 (between populations 2 and 7) to 0.232 (between populations 1 and 11), and 63 of 65 comparisons were significant (p < 0.05); non-significant differentiation was observed between populations 2 and 7, and between populations 4 and 6 (data not shown).

The mean observed M-ratio ranged from 0.168 to 0.830. It was significantly (p < 0.05) lower than the simulated equilibrium distribution of the M-ratio (MReq) in five populations (Table 3). Specifically, evidence for a recent bottleneck was detected in populations from the Graian Alps (populations 3 and 5), Bernese Alps (populations 7 and 8) and the Italian part of the Central Eastern Alps (population 11).

Table 3 Results of the M-ratio test under the two-phase model (TPM) estimated for Salix serpillifolia populations across ten polymorphic nSSR loci; population abbreviations as in Table 1
Full size table

Chloroplast microsatellites

Chloroplast microsatellites displayed a high level of variation: the three cpSSR loci used in our study yielded four (CSU01), six (CSU05) and nine (CSU07) size variants differing by a single nucleotide. Together, the cpSSR size variants were combined into 53 haplotypes. Generally, the haplotypes occurred at low frequencies; only a single haplotype reached a frequency of 10.4% and three attained a frequency of 5–10%; the remaining haplotypes occurred at frequencies below 5%. A summary of the haplotypic diversity of populations is provided in Table 4. The highest total number of haplotypes per population was 10, and it was detected in populations 7 and 11 (Central and Eastern Alps, respectively), while population 17 was characterized by the presence of only two haplotypes, which was the lowest number noted in the whole dataset. Private haplotypes were noted in most of the populations (Table 4). In population 7, six out of the ten haplotypes detected were the private ones, which was the highest number noted. Generally, more private haplotypes were found in populations from the Western Alps than in populations from the Eastern Alps. In most populations, the effective number of haplotypes (NE) was not very different from the total number of haplotypes (A), i.e., 6.824 vs. 5.278 on average. High genetic diversity (HE) was noted and ranged from 0.436 (population 17) to 1.0 in population 15, in which five of the analysed individuals had their own distinct haplotypes. The mean intrapopulation genetic distance ranged from \( D_{\text{sh}}^{2} \) = 0.582 in population 17 up to 8.206 in population 10, with the average value amounting to 2.949. The overall genetic differentiation among populations was moderate and reached FST = 0.198. The highest values of pairwise FST were found between populations 6 and 17 (0.109) and populations 14 and 17 (0.101). The lowest differentiation was noted between populations 1 and 12 (0.006) (data not shown).

Table 4 Summary of genetic variability estimated for populations of Salix serpillifolia across cpSSR loci; population abbreviations as in Table 1
Full size table

Spatial genetic structure

A weak isolation by distance was revealed by the Mantel test based on the genetic (nSSRs) and geographic distance matrices (R2 = 0.13, p = 0.02) (Online Resource 4). We did not find a similar correlation for the cpSSR dataset.

Monmonier’s algorithm identified five genetic barriers among the analysed populations (Fig. 2a), but only the first three barriers received meaningful bootstrap support. These three strongest barriers were detected for the isolated population 1 from the Maritime Alps and the most north-easterly located population 11 from the Sarntal Alps. The distinct characteristics of populations 1 and 11 were repeatedly noted in the other analyses (Figs. 2, 3). The third barrier divided the Alps into western and eastern domains. The weakest barrier separated the most westerly situated population 3 from the remaining western Alpine populations.

Fig. 2
figure 2

Population genetic structure of Salix serpillifolia, as based on the nSSR loci; a proportion of membership of individuals from a given population in each of the three clusters as indicated by the Bayesian clustering analysis (STRUCTURE); red lines: genetic boundaries, as detected by BARRIER and obtained with Monmonier’s algorithm, which was based on Reynolds’ genetic distance (numbers represent the bootstrap values indicating support for the respective boundaries), b barplot proportion of membership of each individual in the three assumed clusters (K = 3), c chart of ΔK values. Refugial areas as explained in Fig. 1

Full size image
Fig. 3
figure 3

Population-based genetic structure of Salix serpillifolia, as determined by the cpSSR loci; a pie charts representing the average per cluster assignment values for all individuals from a given population based on the individual-based clustering algorithm in BAPS (K = 11) and the refugial areas in the Alps, as defined by Schönswetter et al. (2005a); b bar plot with the assignment analysis in BAPS (K = 11). Refugial areas as explained in Fig. 1

Full size image

The Bayesian clustering of the nSSRs (Fig. 2) revealed that the best supported number of clusters was K = 3 (Fig. 2c). Some geographic pattern of differentiation corresponding to the main regions of the Alps was detected (Fig. 2a, b). The most easterly situated populations 11 and 2 were mostly assigned to cluster I and they were the most distinct from the rest of the populations. In particular, population 11 from the Central Eastern Alps (the Sarntal Alps) showed only limited overlap with the gene pools of other inferred clusters. Most of the populations from the Western Alps (populations 2–8) and the single population from the Southern Limestone Alps (population 14) formed a separate and rather homogenous group belonging to cluster II with the high individual proportion of the memberships ranging from 66 to 98%. Only populations 3 and 14 from this group presented high admixture with other gene pools. Interestingly, the most northerly located population 9 (the Glarus Alps) showed a genetic mixture of the western and eastern alpine gene pools (cluster I and III, respectively) (Fig. 2a).

The BAPS clustering for the cpSSRs did not reveal any clear geographically-dependent pattern (Fig. 3). The most likely number of clusters inferred was K = 11, and most populations consisted of individuals showing genetic affinities with different clusters. Considering the major cluster only, some geographic specificity appeared and indicated the distinctiveness of the populations from the Western and Eastern Alps. However, both the western and eastern distribution domains were inhabited with populations showing a large intermixing of a few clusters that occurred in those areas (Fig. 3a, b). Interestingly, in some populations, there were individuals that represented exclusive BAPS clusters, which were not detected elsewhere. Especially, the monotypic population 6 and the polytypic populations 7 and 8, all three from the inner part of the massif, were characterized by unique genetic structure. In populations 7 and 8 some individuals formed own clusters, while the remaining ones were grouped into clusters found in the Western and Eastern Alps (populations 7 and 8, respectively).

Species distribution modelling

Predicted current distribution of S. serpillifolia estimated by the MaxEnt software primarily covers the Alps, with more favourable climatic conditions in the western part of the Alps (Fig. 4a). The model-based distribution of S. serpillifolia during the LGM is smaller than the contemporary one and shifted to the south-west: it covers the Apuan Alps in the Northern Apennines and the south-western parts of the Alps, including the Maritime, Ligurian, Cottian and Provence Alps (Fig. 4b). With lower suitability, this projected range includes the mountain areas of the central and southern parts of the Swiss Alps, and, to a lesser extent, areas of the Central Eastern Alps and the Southern Limestone Alps. The most important climatic factor that shapes the putative range of S. serpillifolia is the mean temperature of the warmest quarter (58.6%) (Online Resource 2).

Fig. 4
figure 4

MaxEnt-modelled geographic distributions under a current, and b Last Glacial Maximum climate conditions

Full size image

Discussion

Diversity and spatial genetic structure of Salix serpillifolia

Our study shows that S. serpillifolia harbours moderate genetic and allelic diversity at nuclear microsatellites (HE = 0.644; A = 10.933; at the species level, not presented in the results section), yet this diversity was comparable to previously observed estimations in the arctic-alpine willow species, S. lanata, S. lapponum and S. herbacea (HE = 0.706, 0.703, 0.527; A = 8.4, 12.4, 5.8, respectively) (Stamati et al. 2007). Considerable diversity was reported at chloroplast microsatellites, but they had lower allelic variability (HE = 0.844; A = 6.824.). However, the number of size variants that were obtained from the three cpSSR loci in S. serpillifolia was nearly two times higher than that reported earlier for S. reinii (19 and 10, respectively) (Lian et al. 2003).

According to M-ratio, five populations (3, 5, 7, 8 and 11) experienced a bottleneck effect. Generally, detection of demographic decline with M-ratio indicates that the bottleneck lasted several generations and population has recovered (Williamson-Natesan 2005). In the case of S. serpillifolia, these results suggest that population decline might have been related to the severe conditions prevailing during the LGM. As purely stochastic process, bottleneck leads to long-term reduction of effective population sizes. Consequently, it results in loss of allelic and gene diversity and may distort the inter-population differentiation pattern.

A significant among-population component of variation (FST = 7.4%) and weak yet significant IBD was noted at nSSRs (Online Resource 4), suggesting that the differentiation was rather drift-induced. As expected from the more variable cpSSRs that recruit from the non-recombining plastid genome, the among-population differentiation was high (FST = 19.8%).

The genetic barriers among populations detected using nSSRs correspond roughly to the general division of the genetic diversity of Alpine biota into four groups encompassing the South-Western, Western, Central and Eastern Alps (Schmitt 2017). Less pronounced yet still noticeable was the division into two parts that agrees with the genetic division assuming the western and eastern domain (Thiel-Egenter et al. 2011; Mosca et al. 2014). Among three STRUCTURE clusters, weak geographically dependent structure was inferred. Specifically, the homogenization of the diversity among populations at the Western Alps was apparent in which a single cluster predominated (cluster II); only populations 1 and 3 from this group displayed considerable admixture from other gene pools (cluster I). Slightly more differentiation and admixing were visible in populations from the eastern part of the Alps. Cluster III prevailed in the most easterly located populations 11 and 12, possibly indicating a south-eastern or even eastern colonization route (Fig. 2). However, denser sampling in the Eastern Alps would be highly recommended to make clear inferences about eastern Alpine refugia and their contribution in recolonization.

Possible peripheral LGM refugia for Salix serpillifolia

Many high-alpine species show clear edaphic limitations that may determine the possibility of survival and the genetic structure (Alvarez et al. 2009). Siliceous substrate dominates mainly in the inner Alps, which were covered with ice sheets during the glaciation period (Kelly et al. 2004). Consequently, wider possibilities for refugial persistence are predicted for calcicolous plants such as S. serpillifolia, for which the widespread territories in the southern and northern outskirts of the Alps were likely available (Tribsch et al. 2002; Tribsch and Schönswetter 2003; Tribsch 2004).

The south-western periphery of the Alps remained largely free of ice during the Pleistocene and was characterized by suitable calcareous habitats (Ehlers and Gibbard 2004; Naciri and Gaudeul 2007). This part of the Alps is frequently described as a refugial area for numerous herbaceous plant (e.g. Gaudeul et al. 2000; Stehlik et al. 2002b; Schönswetter et al. 2002; 2003a; 2004; Ehrich et al. 2007; Naciri and Gaudeul 2007; Paun et al. 2008; Rogivue et al. 2018). With considerable support, this part of the Alps was indicated by MaxEnt as the main area of S. serpillifolia occurrence during the LGM (Fig. 4b). Interestingly, MaxEnt indicated also more southerly located the Apuan Alps in the Northern Apennines as a potential area of the species survival of the glacial period, but there is no additional evidence in support of this possibility, therefore, the niche modelling can be here an artefact. Moreover, the applied model was based on bioclimatic variables without considering the other important factors affecting the distribution of species, such as substrate type or orographic complexity.

The genetic analysis also revealed a significant structure among the populations inhabiting the south-western part of the Alps, especially in the cpSSRs. The private haplotypes that were detected in our dataset were primarily concentrated in this area (Table 4). Specifically, population 1 located in the Maritime Alps was characterized with a distinct genetic makeup that was confirmed by both marker systems (nSSRs and cpSSRs), and this may be expected due to its long-term isolation. Low to moderate indicators of genetic diversity at the nSSRs and cpSSRs were noted in this population. A founder effect or long-term in situ persistence under isolation could result in the genetic impoverishment of population 1. However, which of those two processes was responsible for the reduction in genetic diversity cannot be fully explained with our genetic data. With regards to the other populations from this region, the location of populations 2 and 4 strictly overlap with the already outlined peripheral south-western refugium (Schönswetter et al. 2005), and their haplotypic composition further supports their genetic distinctiveness and thus likely refugial character.

Nuclear SSRs indicates that the LGM refugia for S. serpillifolia might have been located also in the southern part of the Eastern Alps. The genetic diversity was above average in most of the populations from this region that is one of the indicators of refugial status. Population 12 (the Dolomites), in particular, was characterized by the highest average number of alleles, private alleles and expected heterozygosity. This population, along with population 14, belongs to the area of southern-alpine peripheral refugium between Lake Como and the Dolomites, where many other species probably survived during the LGM (e.g. Schönswetter et al. 2002, 2003a, 2004; Tribsch et al. 2002; Naciri and Gaudeul 2007). Interestingly, population 9 from the Glarus Alps refers to the eastern-alpine populations. It is located in the Western Alps (near the border with the Eastern Alps), but more northerly and closer to putative northern alpine peripheral refugia. This population was second in terms of genetic diversity. Nevertheless, regarding S. serpillifolia, there was not a distinct north–south split between populations 12 and 9. Conversely, their genetic distance and structure were rather close, suggesting considerable postglacial gene flow in this longitudinal direction.

Although the easternmost parts of the Alps remained largely unglaciated during the LGM (Ivy-Ochs et al. 2008), this area was not indicated by MaxEnt as a potential refugium for S. serpillifolia. Nevertheless, it is acknowledged to be an important glacial refugium for many high-mountain plants (e.g. Gugerli et al. 2009; Mosca et al. 2012; Ronikier et al. 2008; Stehlik et al. 2001; Schönswetter et al. 2002, a, b, 2004; Tribsch et al. 2002; Mráz et al. 2007). Considerable genetic distinctiveness of population 11 (the Sarntal Alps), which is the most easterly located population in our study, may suggest that this part of the Alps might have harboured populations of this species during the LGM. On the other hand, the detected bottleneck effect may also be responsible for the obtained pattern of differentiation. Unfortunately, we did not sample any population of S. serpillifolia in the eastern peripheries of the Alpine range, which does not allow us to draw a final conclusion on the role of this area as a potential glacial refugium for this species.

Nunatak refugia for S. serpillifolia: yes or no?

The central part of the Swiss Alps was indicated by MaxEnt as climatically suitable for S. serpillifolia during the LGM. This area corresponds to the locations of the nunatak refugia, which increases the overall likelihood of our hypothesis concerning S. serpillifolia survival on nunataks (Stehlik 2000; Schönswetter et al. 2005) (Fig. 4b).

The past species occurrence may be inferred directly from the in situ fossils. To the best of our knowledge, there are only three fossil records of S. retusa agg. dated back to the Late Glacial (14,600–12,000 years BP) (Fig. 1) (Gaillard and Moulin 1989; Tinner et al. 1996; Hadorn et al. 2002; Tinner and Kaltenrieder 2005; Drescher-Schneider et al. 2007; Leesch et al. 2012). The example of one of them, the Gouillé Rion, which is situated at high elevation (2343 m a.s.l.) in the Swiss Central Alps, shows that S. retusa-type (including S. serpillifolia) willows might have appeared relatively quickly after the retreating of the ice sheet. This may be explained both by effective long-distance dispersal as well as by the sources of diaspores located nearby, i.e. nunatak survival. Nevertheless, these paleobotanical data originate from the late glacial period, and thus cannot directly confirm the nunatak survival of S. serpillifolia. Interestingly, the existence of nunatak refugia was recently suggested for several tree species (Carcaillet and Blarquez 2017, 2019; Carcaillet et al. 2018). The authors gave many examples to support the idea that the inner Alpine massifs (even above 2000 m a.s.l.) sheltered the glacial refugia for mountain and stone pines, larch and birch. Although this idea is inspiring and should be considered in our discussion, it still remains controversial and continues to be the subject of interesting debate (Finsinger et al. 2018).

Our results give an indication of possible nunatak survival of S. serpillifolia in populations from the inner part of the Alps, specifically from the central Swiss Alps, due to their location and the very distinct genetic characteristics demonstrated by cpSSRs (Fig. 3). Specifically, in our opinion, the best candidates for nunatak survival are monotypic populations 6 and 10 and the polytypic populations 7 and 8; the two latter show the combination of their exclusive clusters and an admixture from clusters occurring in the Western and Eastern Alps. Moreover, populations 6 and 10 are situated in the vicinity of potential nunatak refugial areas that were defined by Schönswetter et al. (2005), and population 7 is almost entirely composed of the private cpSSRs haplotypes. In this respect, population 5 in the western Alps, which is located at an elevation of 2100 m, also shows distinct genetic characteristics and could, therefore, be included in the category of the nunatak refugia. According to the nunatak hypothesis, small and isolated populations might have built up their genetic distinctiveness by the accumulation of new mutations, bottlenecks, inbreeding, and random fluctuations in haplotype frequencies (genetic drift effect) that would eventually lead to loss or fixation of some haplotypes (Ellstrand and Elam 1993; Frankham 1996; Hewitt 2004; Mee and Moore 2014). The latter explanation would be exemplified by populations 5 and 7, which both experienced bottlenecks that are likely responsible for their distinctiveness. Most of the populations that we suggest to be the nunatak survivors experienced demographic decline which is highly possible during harsh conditions of the LGM. Additionally, nunatak refugia need to be assumed to be rather spatially confined that also makes the occurrence of bottleneck more likely. Consequently, populations 5, 6, 7, 8 and 10 would outline the location of the nunatak refugia in the Alps that would extend those regions that were previously proposed by Schönswetter et al. (2005).

Chloroplast markers which are maternally inherited and efficiently dispersed by seeds, allow seed migration to be traced. However, the efficient dispersion of seeds may easily influence the original gene pool of the isolated populations by swamping their genetic distinctiveness. Willows are anemochoric, which is generally acknowledged to be an efficient mode of gene dispersal. Hence, three of the five populations designated here to be nunatak survivors (5, 7, and 8), present admixture from surrounding stands during the postglacial period.

Salix serpillifolia can reproduce clonally. During sampling, we intended to avoid ramets but we found two identical multilocus genotypes in four populations, which were separated by at least 20 m; this finding gives a very rough hint on the clonal spread of this species. It has been suggested, that this life history trait might have been crucial for forming refugia in northern latitudes (Bhagwat and Willis 2008; Dering et al. 2016; Westergaard et al. 2018) and would surely be so in the case of nunatak refugia. However, the increased prevalence of clonal spread may lead to the partial loss of genetic variability (Dering et al. 2016). Clonality may also result in the fixation of a few genotypes, due to its prevalence in size-confined refugial populations, giving them a distinct genetic composition.

Conclusions

Nunatak survival is still controversial topic in phylogeographic studies. In this work, we combined genetic structure analysis and species distribution modelling to verify the hypothesis on nunatak refugia for S. serpillifolia. We have found a certain genetic support for in situ glacial survival on Alpine nunataks in five populations from the inner part of the Swiss Alps. Nevertheless, this genetic signal is rather weak, since the genetic composition of the local refugial micropopulations was likely blurred by gene flow of periglacial genotypes during the Holocene.