First evidence of post‐glacial contraction of Alpine endemics: Insights from Berardia subacaulis in the European Alps

Late Quaternary glaciations left an enduring imprint on the distribution of species and their genetic structure. The responses of plants endemic to the Alps can be summarized in three major demographic hypotheses: (i) post‐glacial expansion hypothesis; (ii) post‐glacial contraction hypothesis; and (iii) long‐term stability hypothesis. Here we test these hypotheses and reconstruct the time and extent of demographic responses of an endemic plant to the Late Quaternary climate dynamics.


| INTRODUC TI ON
Late Quaternary glaciations, in particular the latest glacial cycle (120-0 Kya; Seguinot et al., 2018), left an enduring imprint on biodiversity globally (Bennett & Provan, 2008;Hewitt, 2004;Svenning et al., 2015;Weigelt et al., 2016). The rapid and severe climatic fluctuations during the Late Quaternary induced marked environmental changes, often resulting in dramatic changes in species population sizes, genetic structure and distributional ranges (Bennett & Provan, 2008;Hewitt, 2004;Stewart et al., 2010). The effects of these climatic fluctuations varied according to latitude, topography and intrinsic features of each species, including their ecological requirements and dispersal capabilities (Hewitt, 2004;Stewart et al., 2010).
Therefore, different demographic responses and population genetic patterns should be expected in different species depending on the environmental conditions they experienced and their ecological preferences (Bennett & Provan, 2008;Nieto Feliner, 2014;Stewart et al., 2010).
Despite the different patterns expected in response to glaciations, in the European Alps (hereafter Alps) the majority of studies suggests that endemic species survived in refugia during the glaciations, regardless of their ecological requirements. In fact, the Alps were nearly entirely covered by ice (Ehlers et al., 2011) with numerous nunataks (unglaciated peaks) within the ice sheet and larger unglaciated areas at the periphery of the ice sheet (Schönswetter et al., 2005;Van Husen, 1997). Consequently, during glacial maxima, endemic plants contracted their distributional range, surviving in areas both within (i.e. nunataks) and outside (i.e. peripheral or lowland refugia) the ice sheet (Pawłowski, 1970;Schönswetter & Schneeweiss, 2019;Schönswetter et al., 2005;Tribsch, 2004;Tribsch & Schönswetter, 2003). When the glaciers retreated and new areas became available, these endemic species may have expanded their distributional range, but only moderately ('post-glacial expansion hypothesis'). The limited extent of post-glacial expansion of such endemics has been explained in terms of reduced colonization ability resulting from dispersal limitation (Dullinger et al., 2012;Essl et al., 2011;Svenning & Skov, 2004, 2007a, 2007b, lower reproductive output than widespread species (Lavergne et al., 2004), and loss of genetic diversity during glacial contraction (Ehrendorfer, 1965;Niklfeld, 1972). In the Alps, this hypothesis has been widely supported in several studies on alpine endemic species based on molecular analyses and species distribution models (Bettin et al., 2007;Schönswetter et al., 2002;Stehlik, 2003;Stehlik et al., 2002;Tribsch, 2004).
However, few studies suggest that, in regions of the Alps where the effects of Late Quaternary glaciations were less dramatic (e.g. Maritime and Ligurian Alps), endemics survived in situ via short altitudinal shifts Patsiou et al., 2014). This long-term stability ('long-term stability hypothesis') is thought to result from continued moisture availability and complex local topography (Bátori et al., 2017;Scherrer & Körner, 2011;Tzedakis et al., 2002) that allowed in situ survival via short-distance dispersal during climatic shifts, minimizing the extinction of populations and genotypes (Médail & Diadema, 2009;Ohlemüller et al., 2008;Suchan et al., 2019). Moreover, this long-term stability allowed the accumulation of species and especially of old endemic lineages in certain areas, including the southern European mountains (Médail & Diadema, 2009;Ronikier et al., 2012;Tzedakis et al., 2002).
Alternatively, some cold-adapted endemics attained larger distributions during glacial maxima by surviving and expanding in icefree areas (Birks, 2008). When glaciers retreated and temperatures increased, these cold-adapted endemics contracted their distributional ranges ('post-glacial contraction hypothesis'), thus mountain chains might serve as current refugial areas. To the best of our knowledge, the hypothesis of post-glacial contraction has never been demonstrated for endemics in the Alps, while it has been supported for endemic species in areas less affected by glaciations, such as mountains in the central Iberian Peninsula and the Japanese Alps (Ikeda et al., 2008;Ikeda & Setoguchi, 2007;Martín-Bravo et al., 2009;Peredo et al., 2009). In the central Iberian Peninsula and the Japanese Alps, given the limited ice extent during the Last Glacial Maximum (LGM), some cold-adapted endemics underwent range expansions during glacial maxima, and contracted their range into highaltitude refugia as mountain forests recovered during post-glacial (Ikeda et al., 2008;Martín-Bravo et al., 2009). The post-glacial contraction of such endemics has been mainly explained by an increase in competition of larger plants and/or their intolerance to warmer temperatures (Birks, 2008).
An ideal study area to test the three competing hypotheses, which, to the best of our knowledge, have never been explicitly tested before in Alpine endemics, are the South Western European Alps (hereafter SW Alps). Contrary to the rest of the Alps, in the SW Alps the effects of glaciations were less severe: in fact, the LGM (~24 Kya BP in the SW Alps, Seguinot et al., 2018) temperature depression was moderated by the Mediterranean Sea influence (Seguinot et al., 2018). The SW Alps are characterized by high local climatic heterogeneities, as a consequence of the close proximity of the Mediterranean and Alpine climates and high topographic heterogeneity (Casazza et al., 2005(Casazza et al., , 2008Fauquette et al., 2018). This climatic heterogeneity results in the variety of phylogeographical patterns detected in this region, where both the post-glacial expansion and the long-term stability (mainly by altitudinal shift) hypotheses have been proposed Médail & Diadema, 2009). Moreover, the limited extent of the glacial sheet in this area may have allowed cold-adapted species to expand during glaciation, as expected by the post-glacial contraction hypothesis, even if this K E Y W O R D S Approximate Bayesian Computation, endemism centre, genotyping by sequencing, past climate change, post-glacial contraction, species distribution models pattern has never been detected. Consequently, the SW Alps are an ideal study area for investigating all three hypotheses outlined above. Here, we chose Berardia subacaulis Vill., an Alpine endemic widely distributed in the SW Alps, to test the three alternative hypotheses. On the basis of its current distributional range, Susanna and Garcia-Jacas (2009) hypothesized that B. subacaulis survived the glacial period in refugia in the southern part of the SW Alps, seemingly supporting the post-glacial expansion hypothesis. Alternatively,  suggested that B. subacaulis likely survived in situ, in line with the long-term stability hypothesis, but these authors used paleoclimatic data at a coarse scale, hence they might have underestimated climate-change impacts compared to finer-scale data (Franklin et al., 2013). Moreover, the authors of the mentioned studies on B. subacaulis did not use molecular data or demographic inference to explicitly test alternative hypotheses of post-glacial colonization.
Here we test the three competing hypotheses described above Bayesian Computation framework. We first reconstruct spatial dynamics through the Late Quaternary climatic cycles using a highresolution paleoclimatic data set spanning the last 28 Ky (from Upper Pleniglacial to Holocene - Maiorano et al., 2013;Theodoridis et al., 2017), and quantify genetic variation in B. subacaulis using a reduced representation library approach. We finally integrate the results of the two independent approaches to reconstruct the recent demographic history of B. subacaulis. The results of this study should allow us to better understand the processes that shaped the distribution of Alpine endemic species and their responses to postglacial warming.

| Study species
Berardia subacaulis Vill. belongs to a monospecific genus endemic to the SW Alps (Figure 2), it is an element of the Tertiary paleo flora (Ozenda, 2009), which went almost completely extinct during Late Quaternary glaciations (Ozenda, 2009 and stony slopes at high altitudes (between 1700 and 2700 m), currently restricted to the SW Alps ( Figure 2). Capitula are mainly solitary, pollinated by a wide array of insects, whose visits remain scarce . The flowers are self-compatible but protandrous, favouring cross-fertilization first and allowing for self-fertilization when cross-fertilization does not occur . Pappus is characterized by short hairs, which combined with the acaulescence of the plant, seems to be less efficient. On the basis of the poorly efficient plume, we categorized the dispersal mode of B. subacaulis as trichometeorochory (dispersal distance of ~15 mt for the 99% of the seeds of a plant, classification of Vittoz & Engler, 2007). Information on the generation times of B. subacaulis is scarce. It is known that seed dormancy is absent and that from the germination it might occur from 40 to 80 days F I G U R E 1 Three main alternative hypotheses for the response of endemic species to climatic oscillations during the Late Quaternary: (i) 'post-glacial expansion hypothesis', characterized by expansion of population size and distributional range after the Last Glacial Maximum (LGM) (ii) 'post-glacial contraction hypothesis', characterized by reduction in population sizes and distributional range after the LGM, and (iii) 'long-term stability hypothesis', characterized by no or weak changes in population sizes and distributional range  Table S4), covering the entire distributional range of the species. Extractions were performed using a modified CTAB method (Doyle & Doyle, 1990) and concentrated to a minimum of 30 ng/μl for samples. DNA quality was visualized on 0.8% agarose gels and quantity was assessed using a QuBit 2.0 fluorometer (2.0, Life Technologies).
The extracted DNA was sent to the Cornell Institute for Genomic Diversity for genotyping-by-sequencing (GBS; Elshire et al., 2011). Individual DNA samples were digested with the restriction enzyme EcoT22I. Single-end sequencing was carried out using one lane of an Illumina HiSeq 2000, producing raw reads that were 101 bp long.

| SNPs discovery and genotyping
To demultiplex the raw data and recover the individual samples in the Illumina library, the sequences were processed with the program process_radtags in Stacks 1.35 (Catchen et al., 2011(Catchen et al., , 2013 and the single-end reads were then merged into a single file per individual. We used Trimmomatic v0.33 (Bolger et al., 2014) to trim low quality or uncalled bases from raw reads, remove potential adapter contamination and discard reads that were less than 85 bp (i.e. the length of the trimmed reads) long after low quality base trimming. Because of the lack of a reference genome for B. subacaulis, a de novo assembly of the filtered reads was performed with de_novo_map.pl in Stacks 1.35 using the following parameters: a minimum number of identical raw reads required to create a stack (m) 10; the number of mismatches allowed between loci when processing a single individual (M) 4; and the number of mismatches allowed between loci when building the catalog (n) 8.
At this stage, we eliminated four individuals, each one of them belonging to a different population, due to low sequence coverage.
Finally, we used the POPULATIONS program in Stacks 1.35 to create data matrices of sequences and single nucleotide polymorphisms (SNPs) in fasta and vcf format, respectively, for all subsequent molecular analyses.

| Population genetic structure
The assignment of individuals to population clusters was conducted using the program STRUCTURE version 2.3.4 (Pritchard et al., 2000). Since STRUCTURE assumes that the investigated loci are unlinked, we used the package Gonospy v0.1 (https://github. com/spyro stheo dorid is/gonospy) to randomly select only one SNP per locus and repeated the analysis five times to account for biases in the SNP selection process, as routinely done (Theodoridis et al., 2017). We used a model of admixture to determine the number of population clusters (K) with a burn-in of 50,000 and 100,000 iterations. We performed five replicate runs for each value of the number of population clusters, K, between 1 and 10. The average and standard deviation (SD) of the natural log probability across replicates were used to calculate deltaK (Evanno et al., 2005) using Structure Harvester (Earl & vonHoldt, 2012). For the best K value, the results of the five randomly selected SNP repetitions and five replicated STRUCTURE runs per repetition were combined using the 'greedy' algorithm within CLUMPP v1.1.2n (Jakobsson & Rosenberg, 2007) and plotted using DISTRUCT (Rosenberg, 2004).
To estimate partitioning of genetic variance among the groups determined with STRUCTURE, among populations within groups, within populations and within individuals we used an analysis of molecular variance (AMOVA; Excoffier et al., 1992) in Arlequin 3.5 (Excoffier & Lischer, 2010), where the total variance is partitioned into components analogous to F-statistics. AMOVA tests were conducted using 10,000 permutations between groups on the three matrices identifying two groups as the best partition (see STRUCTURE results). Because principal component analysis (PCA) is free from many of the population genetic assumptions underlying STRUCTURE (Gao et al., 2007;Jombart et al., 2009), we further performed a PCA to visualize the major axes of variation applied to the five random SNP matrices using the dudi.pca function implemented in 'ade4' package (Dray & Dufour, 2007) in R (R Development Core Team, 2008).
Finally, we used TreeMix to estimate the maximum likelihood tree for the populations of B. subacaulis. TreeMix uses as input F I G U R E 2 Study area in the South-Western Alps [Colour figure can be viewed at wileyonlinelibrary.com] a set of population allele frequencies and assumes biallelic, unlinked sites (Pickrell & Pritchard, 2012). Therefore, the analysis was repeated on each of the random SNP matrices (with only one random SNP per locus, excluding SNPs with more than two nucleotide states across populations) and for each SNP repetition we ran 1,000 bootstrap replicates. A consensus of the trees across all SNP repetitions and bootstrap replicates was computed in R using 'phangorn' (Schliep, 2011) and 'ape' (Paradis et al., 2004) packages.

| Genetic diversity statistics
The simplest measure of genetic diversity at a locus is the number of alleles (allelic richness) and a simple measure of genetic distinctiveness is the number of unique alleles in a population (private allelic richness) (Kalinowski, 2004). Following a rarefaction approach to account for sample size differences across lineages (Kalinowski, 2004), allelic richness and private allelic richness per locus were estimated in the genetic groups identified by STRUCTURE. Following the approach of Theodoridis et al. (2017), we used the full sequence of each locus as recovered in the Stacks analysis (fasta file). Loci present in at least 15 diploid individuals were kept, resulting in a minimum of 30 sequences per locus and per genetic group. Both statistics were calculated for sequence sample sizes (g) from 2 to 30 (1 to 15 diploid individuals). This analysis was conducted using custom scripts written in python (see Theodoridis et al., 2017).

| Species occurrence data and environmental layers
Occurrence data were retrieved from the database "SILENE" of the Conservatoire Botanique National Méditerranéen de Porquerolles (CBNMED) and "FLORE" of the Conservatoire Botanique National Alpin de Gap (CBNA). We gathered a total of 1942 occurrences for the species from CBNMED, from CBNA and from our field surveys.
To increase the accuracy of the distribution models, we used only data collected with a GPS instrument during the field surveys; a final data set of 1184 presence records was used in the analyses.
To train SDMs and estimate current climatically suitable habitats, we first extracted the current monthly mean temperature (tmean) and precipitation (prec) from the WorldClim data set (version 1.4) at 30-s (c. 1 km) spatial resolution (Hijmans et al., 2005). Using the monthly climatic values and the function biovars in the package 'dismo' (Hijmans et al., 2012), we derived a set of 17 climatic variables (Table S1). We then performed a pairwise Pearson correlation among bioclimatic predictors aimed at reducing multicollinearity and minimizing model overfitting, retaining only predictors that were not highly correlated (r ≤ |0.80|; see recommendation of Elith et al., 2006) and that are physiologically important for plant species.
Three bioclimatic variables were retained for further analyses: mean temperature of driest quarter (BIO9); mean temperature of warmest quarter (BIO10); and precipitation seasonality (BIO15). Because B.
subacaulis grows only on specific calcareous substrates, for all the time slices we added a layer reporting the presence/absence of suitable substrate, based on the global lithological map dataset GLiM (Hartmann & Moosdorf, 2012). The GLiM represents the rock types  Table S3).
To reconstruct the potential past distributional ranges of B. subacaulis, we used a paleoclimate dataset of monthly mean temperature and precipitation data from Maiorano et al. (2013) described in Singarayer and Valdes (2010) and downscaled it at 1 km spatial resolution using the delta method (see Patsiou et al., 2014 for the downscaling method of the paleoclimate dataset). These paleoclimate layers were used to derive the three bioclimatic variables described above. The final paleoclimate dataset covered a period of 28 Ky at 1 ka (from 1 to 21 Kya) and 4 ka intervals (from 21 to 28 Kya).
Because most modelling techniques implemented in BIOMOD2 require that species distribution data are represented as presence and absence, we generated ten different sets of pseudo-absences.
For each SDM technique, the number of pseudo-absences (Table   S2) was selected according to the recommendation of Barbet-Massin et al. (2012). The predictive performance of the models was evaluated for each pseudo-absence run by repeating a splitsample cross-validation 10 times, using a random subset (70%) of the initial dataset each time to calibrate the models, while the remaining 30% was used to evaluate the models. Three different evaluation measures were calculated: area under the receiver operating characteristic curve (AUC), true skill statistic (TSS) and Kappa statistic (Thuiller et al., 2009). To further transform continuous probability values from model outputs to binary presenceabsence we used the threshold that maximized TSS. For the final ensemble projections, we employed 'majority' combinations (i.e. at least three out of five techniques predicting suitable areas) to produce binomial maps of presence and absence through time. To exclude potential suitable areas that were covered by ice, we further filtered these presence-absence maps using reconstructed ice margins spanning a period from 28 to 15 Kya (Ehlers et al., 2011;Rolland et al., 2020).

| Estimation of demographic history
To test for alternative demographic responses of B. subacaulis to the Late Quaternary glaciations, we integrated the genomic data and the results of SDMs in a statistical phylogeography framework using the ABC procedure (Cornuet et al., 2014;Patiño et al., 2015;Ray et al., 2010;. We evaluated three  (Table 2). Current population sizes were based on census population sizes (personal observations in the field) and converted to effective population sizes with reference to population sizes of other perennial plants (Gossmann et al., 2010). The time interval for post-glacial demographic changes (expansion and contraction hypotheses) was inferred by visually identifying changes in the number of suitable cells predicted by the SDMs and according to geological evidence (Rolland et al., 2020) and vegetation reconstruction (Brisset et al., 2015). In fact, in the slopes (~2300 m a.s.l.) where B. subacaulis occurs, it has been showed that the main deglaciation took place at 15.3-14.2 Kya (Rolland et al., 2020) and that the meso-thermophilous flora in- We simulated the same amount of diploid loci that was used to estimate the genetic diversity statistics in B. subacaulis, maintaining the r = log N 0 N t t same size of the empirical data set (see the 'Genetic diversity statistics' section). For each simulated dataset, we calculated the same two summary statistics used in the genomic diversity analysis of the empirical data set: allelic richness and private allelic richness for the two genetic groups (as in François et al., 2008), for a fixed sample size of 15 individuals (30 sequences) as point estimate. We chose allelic richness and private allelic richness because they explicitly account for different sample sizes among genetic groups. Moreover, commonly used summary statistics of genetic diversity calculated from data generated by Reduced Representation Libraries (such as RADseq) are known to deviate considerably from real values, thus considered inappropriate for demographic inferences (Arnold et al., 2013). Nevertheless, we additionally calculated expected heterozygosity for each genetic group using the full sequence of each simulated locus following Nei (1973). To estimate the posterior probabilities of the five models (i.e. the fit between observed and simulated data under each model), we first calculated the Euclidean distance between the observed and simulated summary statistics (standardized across all simulations for each model), following Beaumont et al. (2002). We then adopted the standard rejection procedure proposed by Pritchard et al. (1999), as follows. For each of the five models, we retained the top 1% (i.e. 6000) simulations with the smallest Euclidean distances to the observed summary statistics. The retained simulations for each model were combined and the final 30,000 simulations were ordered by ascending Euclidean distances recomputed on summary statistics standardized with common mean and standard deviation (Patiño et al., 2015;Ray et al., 2010). The posterior probability of each model was then computed as the proportion of simulations executed under the respective model that were included in the set of the top 1% (i.e. 600) simulations closest to the observed data (Estoup et al., 2004;Miller et al., 2005). To estimate the parameters of the model with the higher posterior probability, we used the neural network approach implemented in R package 'abc' (Csilléry et al., 2012). Postsampling adjustment and parameter estimation were carried out using the 6000 top simulations. We report the mode of the posterior distribution and the 95% highest posterior density (HPD) interval for each parameter.

| Population genetic structure
For three out of the five replicates, STRUCTURE analyses identified K = 2 (Figure 3a, Figure S3) as the most probable number of genetic   Figure S4). Nevertheless, AMOVA shows that most of the genetic variance is due to differentiation among populations within groups and within individuals (Table 1, Table S6), indicating a relatively weak genetic differentiation between the two main groups.

| Genetic diversity statistics
We obtained estimates of allelic richness and private allelic richness   Table S5).

| Species distribution modelling
Projections of the potential distribution of B. subacaulis into the past suggested that its putative range strongly decreased as temperature increased ( Figure S1). Additionally, the predicted phases of persistence and contraction were consistent among the five techniques ( Figure S1). Importantly, all model techniques identified a peak in climatically suitable habitats during the LGM and a reduction in these habitats in the period between 24 and 11 Kya

| Estimation of demographic history
The results of ABC-based analyses supported the post-glacial contraction of B. subacaulis with a recent divergence between the genetic groups ( Figure

| DISCUSS ION
In this study, we investigated the response of the endemic B. subacaulis to the Late Quaternary climatic oscillation. Combining several independent lines of evidence in an ABC modeling framework, we were able to discern among three main competing hypotheses (i.e. post glacial contraction, post glacial expansion and long-term stability) and to provide, for the first time, empirical evidence of postglacial demographic contraction and a recent split between the two genetic groups for an endemic plant in the European Alps during the Late Quaternary. However, considering the high species diversity and environmental heterogeneity in the SW Alps, and that the present study is based only on a single species, further research is needed to draw general conclusions about the regional Quaternary history of endemic species in the Alps.
The statistical phylogeographic approach used in our study provided strong support for a recent divergence between the two genetic groups (Figure 4b, Table 2 Figure S1) showed a reduction in potentially suitable areas for B. subacaulis during the post-glacial climate warming, suggesting that B. subacaulis started to contract in this period, while it likely attained a larger distribution during the glaciations . In this scenario, the post-glacial isolation due to range contraction likely led to the recent divergence between the two main genetic groups.
In B. subacaulis we observed a weak genetic subdivision of populations, supported by the low number of genetic groups (the  Figure 3, Figures S2 and   S4). Usually, in a scenario of post-glacial range expansion, a strong phylogeographic structure is expected as result of high differentiation among refugial populations due to the genetic drift (Comes & Kadereit, 1998). Moreover, low levels of genetic diversity are expected in the populations at the wave front of expansion (Excoffier & Ray, 2008;Hewitt, 2000). Differently, we recorded low genetic diversity values in populations scattered over the CS group (Table   S7). The weak phylogeographic structure we detected in B. subacaulis may be explained mainly by a long-lasting period of connectivity among populations either in a scenario of extensive in situ survival and low degree of geographical isolation  or in a scenario of expansion during the glacial cycles followed by a short period of interglacial isolation, thus preventing the accumulation of among-populations genetic differentiation (Garcia et al., 2011;Gibbard & van Kolfschoten, 2004;Kropf et al., 2003;Stewart et al., 2010). This latter scenario is likely the one observed in B. subacaulis, as supported by its recent split between the two genetic groups (2.46 Kya; Table 2, Figure 4b), by SDMs and by ABC analysis.
However, in some species, the small population sizes in interglacial refugia may be more strongly affected by genetic drift, potentially overcoming temporal effects and generating deep phylogeographic structure (Bonatelli et al., 2014), but this does not seem to be the case for B. subacaulis, likely because of its very recent split between the two genetic groups and/or probably limited recent genetic drift.
Taken together, our results support the post-glacial range contraction hypothesis. These results contrast with the dominant pattern observed for endemics in the Alps, based on molecular analyses and SDMs, where it is hypothesized that cold-adapted species, both endemic (e.g. Androsace alpina) and not endemic (e.g. Phytheuma globulariifolium), primarily contracted their range, retreating into refugia during cold periods (Schönswetter et al., 2002(Schönswetter et al., , 2003(Schönswetter et al., , 2004(Schönswetter et al., , 2005Tribsch & Schönswetter, 2003). In fact, in most of the Alps, the massive ice extent during the LGM prevented cold-adapted endemics to reach or to survive in most of the potentially suitable areas, forcing them to retract in nunataks, when available, or to migrate to peripheral refugia (Schönswetter et al., 2005), occupied mainly by trees (Ravazzi, 2002) and steppe vegetation adapted to a cold and dry climate (Janská et al., 2017). Therefore, the reduced availability in suitable habitats may have resulted in range contractions for most of the cold-adapted endemics in the Alps.
The different pattern observed in B. subacaulis might be due to several factors. First, the SW Alps were characterized by a greater availability of ice-free terrain. In fact, here the Mediterranean Sea mitigated the climate during LGM, maintaining temperatures some degrees higher than in the rest of the Alps (Seguinot et al., 2018) and, consequently, the ice cover was less extensive (Federici et al., 2012).
Second, the SW Alps were characterized by relatively high precipitation (Janská et al., 2017), which, combined with the ice-free areas, might have allowed B. subacaulis to persist or even expand in most of the climatically suitable areas at high altitude during the LGM, while other cold-adapted endemics distributed in other parts of the Alps restricted their ranges because of the dryness of glacial steppes beyond the ice margins (Schmitt, 2007;Schönswetter et al., 2005).  (Lavergne et al., 2004) and stones coverage assuring the moisture necessary for plant growth (Körner, 2003;Pérez, 1998 (Birks, 2008;Ikeda & Setoguchi, 2007).
Our findings thus provide support to the idea that the response of cold-adapted endemics to past climate change depends not only on their ecological features, but also on where they occurred (Bennett & Provan, 2008;Nieto Feliner, 2014;Stewart et al., 2010), shedding light on the complexities of organisms' response to the Late Quaternary climatic oscillations.
In general, low levels of genetic diversity (Appendix 1, Table   S7), as observed in B. subacaulis populations, are expected in species with small geographic range (Hamrick & Godt, 1989) and are detected in other Tertiary species like Haberlea rhodopensis (Petrova et al., 2015) and Ramonda myconi (Dubreuil et al., 2008), endemic to other southern European mountains. Loss of genetic diversity in these species was hypothesized to result from severe genetic bottlenecks, likely as a consequence of range contractions during the cold periods (Dubreuil et al., 2008;Petrova et al., 2015). Contrary to the above-mentioned Tertiary species, B. subacaulis contracted its range after the LGM (Figure 4) and likely during other warm periods, such as the Last Interglacial (LIG, ~120 Kya; . As a consequence, its low genetic diversity may result from losses that occurred not only during the recent post-glacial contraction but also during the adverse warm periods. Moreover, because of the strong reduction in population sizes during bottlenecks (Nei et al., 1975), the mixed-mating system of B. subacaulis  might have favoured self-fertilization to ensure reproduction.
In conclusion, as opposed to the main response suggested for endemic species, attaining narrower distributions during glacial maxima (Pawłowski, 1970;Schönswetter & Schneeweiss, 2019;Schönswetter et al., 2005;Tribsch, 2004;Tribsch & Schönswetter, 2003), our results provide the evidence of post-glacial range contraction of an endemic plant in the European Alps through Late-Quaternary climatic oscillations. Future research in areas where the ice cover was less extensive, such as the SW Alps (the richest centre of endemism in the Alps;Aeschimann et al., 2011), will contribute to a more complete understanding of the role of climatic changes in shaping the highly endemic biota of the European Alps.

ACK N OWLED G EM ENTS
The authors thank the Ente di gestione Aree Protette Alpi Marittime for funding. We are grateful to CBNMed (

CO N FLI C T O F I NTE R E S T
The authors declare no conflicts of interest.