Survival in nunatak and peripheral glacial refugia of three alpine plant species is partly predicted by altitudinal segregation

Mountain biota survived the Quaternary cold stages most probably in peripheral refugia and/or ice‐free peaks within ice‐sheets (nunataks). While survival in peripheral refugia has been broadly demonstrated, evidence for nunatak refugia is still scarce. We generated RADseq data from three mountain plant species occurring at different elevations in the southeastern European Alps to investigate the role of different glacial refugia during the Last Glacial Maximum (LGM). We tested the following hypotheses. (i) The deep Piave Valley forms the deepest genetic split in the species distributed across it, delimiting two peripheral refugia. (ii) The montane to alpine species Campanula morettiana and Primula tyrolensis survived the LGM in peripheral refugia, while high‐alpine to subnival Saxifraga facchinii likely survived in several nunatak refugia. (iii) The lower elevation species suffered a strong population decline during the LGM. By contrast, the higher elevation species shows long‐term stability of population sizes due to survival on permanently ice‐free peaks and small population sizes at present. We found peripheral refugia on both sides of the Piave Valley, which acted as a major genetic barrier. Demographic modelling confirmed nunatak survival not only for S. facchinii but also for montane to alpine C. morettiana. Altitudinal segregation influenced the species' demographic fluctuations, with the lower elevation species showing a significant population increase at the end of the LGM, and the higher elevation species either showing decrease towards the present or stable population sizes with a short bottleneck. Our results highlight the role of nunatak survival and species ecology in the demographic history of mountain species.


| INTRODUC TI ON
The alternation of cold (glacial) and warm (interglacial) periods during the Quaternary caused large-scale expansions and contractions of glaciers and ice sheets in high latitudes and mountainous areas all over the world (Ehlers et al., 2011).During the Last Glacial Maximum (LGM, ∼20 ka; Clark et al., 2009), most of the European Alps were covered by ice (Ivy-Ochs et al., 2008), which forced the distribution ranges of mountain biota to contract to three main types of refugia (Holderegger & Thiel-Egenter, 2009).(i) Nunatak refugia were located on summits and ridges protruding from the ice sheet, potentially allowing the survival of extremely cold-hardy species (Schönswetter et al., 2005) in interior parts of the Alps.(ii) Peripheral glacial refugia were mountainous areas, which remained largely unglaciated and provided less challenging conditions; in the Alps, they were mainly situated at their southwestern, southern and eastern margin (Schönswetter et al., 2005;van Husen, 1987).Finally, (iii) lowland glacial refugia were located in the plains surrounding the Alps beyond the limits of the ice sheet (Holderegger & Thiel-Egenter, 2009).Lowland refugia are, however, expected to play a negligible role for high-elevation species of the Alps as they were covered by steppe, wetlands and coniferous forest even at the LGM (Ravazzi et al., 2014;Wesche et al., 2016).
Survival in nunatak and peripheral refugia can in principle be distinguished by differences in the spatial genetic structure and diversity of extant populations (Bettin et al., 2007;Lohse et al., 2011;Schönswetter & Schneeweiss, 2019;Stehlik et al., 2002;Stehlik et al., 2001;Westergaard et al., 2010), as long as no massive postglacial expansion replaced local lineages through genetic swamping (Todesco et al., 2016).Specifically, genetic groups exclusive to strongly glaciated areas, which do not overlap with potential peripheral refugia, are interpreted as descendants of nunatak survivors.
Such groups are additionally expected to exhibit reduced withinpopulation diversity as they likely underwent massive bottlenecks.
In contrast, a genetic group overlapping with an area with weak LGM glaciation exhibiting high-to-moderate levels of genetic diversity is likely derived from recent range expansion from a peripheral refugium (Excoffier et al., 2009;Mráz et al., 2007;Tonin et al., 2023).
However, different processes may result in convergent patterns of genetic structure and diversity, challenging a proper interpretation (Lawson et al., 2018).For instance, a founder event in the course of range expansion from a peripheral refugium might render a postglacially founded population strongly differentiated as a result of enhanced drift, leading to the incorrect inference of nunatak survival.The advent of genetic demographic modelling approaches in recent years provided a framework to discern different scenarios by testing models of divergence between populations, establishing estimates of divergence times and estimating the fluctuation of effective population size (Massatti & Knowles, 2016;Pan et al., 2020).Specifically, a better understanding of the temporal processes ultimately enables the discrimination between nunatak survival and postglacial colonization.
Mountain species adapted to different elevations have likely differentially responded to the Quaternary glaciations (e.g., Guerrina et al., 2022;Kropf et al., 2003;Surina et al., 2011;Theodoridis et al., 2017).Montane to low-alpine species have likely survived the glaciations in peripheral refugia and went extinct in interior ranges.
After the LGM, they likely expanded their distribution ranges when previously glaciated terrain became ice-free, with high potential for a massive postglacial expansion.In contrast, the most cold-hardy species currently occurring in the high alpine and subnival belts may have survived the glaciations on nunataks (Lohse et al., 2011;Pan et al., 2020;Schönswetter & Schneeweiss, 2019;Stehlik et al., 2002;Westergaard et al., 2010).After the LGM, the increased competition from upward migrating species restricted the occurrence of weak competitor cold-hardy species to the highest elevations, which limited the magnitude of postglacial population increase.
Our study area is situated in the Southern Limestone Alps, that is, the southern, mostly calcareous ranges of the Eastern Alps (Figure 1).It comprises large areas close to the range's periphery that remained weakly glaciated or unglaciated during the LGM (Seguinot et al., 2018;van Husen, 1987) and were consequently suggested as peripheral refugia (Escobar García et al., 2012;György et al., 2018;Ronikier et al., 2008;Schönswetter et al., 2005;Záveská et al., 2021).Specifically, we focus on the Dolomites and the southwestern Carnic PreAlps, an area of unique geomorphological and geological structures, best described as 'mainland archipelago' (Bosellini et al., 2003, Figure 1).During the LGM, nunatak 'islands' were isolated from each other by glaciers, which reached up to 2000 m above sea level in the central valleys of the Dolomites (van Husen, 1987).The highest peaks have likely provided snow-free patches, at least during summers (Seguinot et al., 2018), thus allowing nunatak survival, as suggested for the central Dolomites based on distribution patterns of endemic flowering plants (Pignatti & Pignatti, 2016;Prosser et al., 2019).The potential combination of peripheral as well as nunatak refugia makes the study area an ideal system to study the relevance of different types of glacial survival patterns in mountain plants.
The Dolomites and the southwestern Carnic Prealps are part of the UNESCO World Heritage List and are hotspots of biodiversity and endemism (Aeschimann et al., 2011;Pawlowski, 1970;Pitschmann & Reisigl, 1957;Tribsch, 2004).The border between both ranges is formed by the northeast-to-southwest running valley of the Piave river (Figure 1); it coincides with range limits of a suite of alpine species (Aeschimann et al., 2004) and is the eastern limit of a subunit (OAE5b, i.e., the Dolomites) within an area of endemism comprising the Southern Limestone Alps from the Dolomites to the Julian Alps (Tribsch, 2004).Hypothetical but not yet corroborated refugia to the west and to the east of the Piave river may have hosted genetically differentiated lineages of mountain species, whose postglacial range expansion has potentially been constrained by the deep river valley.Similarly, at the scale of the entire Alps, major valleys such as the Aosta Valley or the 'Brenner Line' were recurrently identified as important genetic break zones in various species (Gugerli et al., 2023;Thiel-Egenter et al., 2011).
Here, we studied the population genetic structure of three calcicolous plant species endemic to our study area based on RADseq data (Edwards et al., 2015).Subsequently, we tested demographic scenarios related to glacial survival using diffusion approximation | 3 of 17 ROTA et al. (Gutenkunst et al., 2009) and combined this with effective population size modelling (Liu & Fu, 2015).The study species include mon- increase is expected, as its occurrence is limited to small areas providing high-alpine to nival conditions.

| Study species
We selected our study species C. morettiana, P. tyrolensis and S. facchinii (Figure 1) among the species with occurrence databases from Rota et al. (2022).The selected species had to match the following criteria.(i) Their distribution should be restricted to the southwestern Carnic Prealps and/or the Dolomites.(ii) Their altitudinal distribution should differ, but in principle span the montane to nival belts between 1400 and 3300 m a.s.l.(iii) They should not be closely related and (iv) grow in similar rocky habitats on carbonate bedrock.

| Plant material
Leaf samples were collected in 43 sites in 2019 and 2020 covering the entire distribution range of each species, totalling 18 sites for C. morettiana, 12 for P. tyrolensis and 13 for S. facchinii plus suitable outgroups (Table S1).Young and healthy leaves of three to five individuals per population were sampled and stored in silica gel for the extraction of DNA; in population number 1 of C. morettiana 10 individuals were collected.In some small populations, the number of sampled individuals was reduced (Table S1).Herbarium specimens were not collected given the rarity and conservation status of the species.Collection permits were obtained from the administrative bodies when needed.Flow cytometry was used to establish the ploidy level of each investigated population (Supplementary Text 1).

| DNA extraction
Total genomic DNA was extracted from ca. 10-20 mg dried leaf material following a CTAB protocol (Doyle & Doyle, 1987) with modifications (Tel-Zur et al., 1999).The ground leaf material was washed with a wash buffer containing sorbitol, three times for C. morettiana and P. tyrolensis, and twice for S. facchinii.The quality of the extracts was examined with a Nanodrop spectrometer (NanoDrop ND-2000, Thermo Scientific); subsequently, they were purified with the Nucleospin gDNA clean-up kit (Macherey-Nagel).DNA concentration was estimated using a Qubit 4 fluorometer (ThermoFisher Scientific).

| RADseq: Library preparation, identification of RADseq loci and SNP calling
Single-digest restriction-site associated (RADseq) libraries were prepared using the restriction enzyme PstI (New England Biolabs) and a protocol adapted from Paun et al. (2016).Briefly, we started with 95-140 ng DNA per individual and ligated 100 mM P1 adapters to the restricted samples.Shearing by sonication was performed with a M220 Focused-ultrasonicator (Covaris) with settings targeting a size range of 200-800 bp and a mode at 400 bp (peak in power: 50, duty factor 10%, 200 cycles per burst and treatment time 90 s at 20°C).Libraries were sequenced on Illumina HiSeq as 100 bp singleend reads or NovaSeq 6000 SP as 150 bp paired-end reads at VBCF NGS Unit (http:// www.vbcf.ac.at/ ngs/ ).Sequences were further processed and analysed using the LEO HPC infrastructure of the University of Innsbruck.
Raw reads were quality-filtered and demultiplexed based on individual-specific barcodes using Picard BamIndexDecoder included in the Picard Illumina2bam package (https:// github.com/ wtsin pg/ illum ina2bam) and the program process_radtags.pl in STACKS v. 2.3 (Catchen et al., 2011(Catchen et al., , 2013)).RADseq loci were further assembled, and SNPs were called using the denovo_map.plpipeline as implemented in STACKS.To select the optimal parameters, a preliminary optimization step following the 80% rule (Paris et al., 2017) was conducted.First, we selected 12-15 samples of each species and ran denovomap.plfor values of the number of mismatches allowed between stacks to merge them into a putative locus (-M) from 0 to 8 and a percentage of individuals that must possess a particular locus for it to be included in calculation of population-level statistics (-r) of 80%.The minimum number of raw reads required to form a stack (-m) and the maximum number of differences among loci to be considered as orthologous across multiple samples (-n) were given values equal to M. After optimization, M = 5 was selected as the value optimizing the number of de novo assembled loci for C. morettiana and P. tyrolensis; M = 4 was selected for S. facchinii.These settings were used for subsequent runs of the pipeline with all individuals.
The program populations implemented in STACKS was used to export the selected loci and generate population genetics statistics, applying different filters for the subsequent analyses.Preliminary exploratory analyses of files generated under different filtering parameters were conducted with the R package ADEGENET 2.1.3(Jombart, 2008) in R 3.6.3 (R Core Team, 2020).This allowed us to select a filtering scheme with a good balance of missing data and amount of informative characters.For phylogenetic tree construction and phylogenetic networks, a vcf file was exported including the outgroup and using a minor allele frequency filter (--min_maf) of three individuals and a maximum heterozygosity of 65% (--max_ obs_het), following Rochette and Catchen (2017).The vcf files were then explored in ADEGENET; samples with a content of missing data higher than a certain threshold (50%) were excluded and populations was run again for the new data set.The vcf files were further filtered with VCFtools 0.1.16to keep loci with a minimum coverage of 10× and loci present in a minimum of 50% of the samples.For STRUCTURE, a minimum of 80% of individuals in the data set (-R) had to contain the locus for it to be processed; only the first SNP per locus was kept (--write_single_snp) to meet the assumption of unlinked SNPs, and the same minor allele frequency and maximum heterozygosity filters as for the previous analyses were applied.

| Analyses of SNP data
The optimal grouping of individuals was determined using Bayesian clustering in STRUCTURE 2.3.4,using the admixture model with uncorrelated allele frequencies (Pritchard et al., 2000).Ten replicate runs for K (number of groups) ranging from 1 to 10 were carried out using a burn-in of 50,000 iterations followed by 500,000 additional MCMC iterations.CLUMPAK (Kopelman et al., 2015) was used to summarize the results across different Ks and to produce plots.The optimal K was where the increase in likelihood started to flatten out, the results of replicate runs were similar, and there were no 'ghost clusters' (clusters which are not modal in any individual, Guillot et al., 2005).Additionally, the deltaK criterion was employed, reflecting an abrupt change in likelihood of runs across different Ks (Evanno et al., 2005).FineRADstructure (Malinsky et al., 2018) was used to infer the coancestry matrix of the same data set but keeping all SNPs per RAD-locus.We used the program populations in STACKS to estimate summary statistics including number of private alleles and nucleotide diversity (π) per population.

| Demographic modelling
To explore alternative glacial survival histories among genetic groups, we used the diffusion approximation method δaδi (Gutenkunst et al., 2009) to analyse joint site frequency spectra (JSFS).We fitted two-and three-population demographic models using dadi_pipeline v3.1.4(Portik et al., 2017; Figure S1).The genetic groups analysed were chosen according to genetic structure and phylogenetic results.For population pairs with strong differences in sampling size, the JSFSs were computed for a subsample of the individuals of the group with the highest sample size.We fit different models corresponding with three hypotheses.First, related to hypothesis 1, we fitted 23 two-population models comprising vicariance and founder events in alternative directions (Portik et al., 2017;Charles et al., 2018;Záveská et al., 2021;Figure S1) for the divergence between Dolomites and southwestern Carnic Alps in C. morettiana and P. tyrolensis.To test different refugial hypotheses within the Dolomites related to Hypothesis 2, nine two-population models comprising vicariance and founder events from the south to the north(west) were fit in the case of C. morettiana (Figure S1).For S. facchinii, five three-population models (Barratt et al., 2018;Firneno et al., 2020;Portik et al., 2017) were fit to explore whether the observed admixture in the central group was caused by gene flow or by hybrid origin (Firneno et al., 2020).Finally, seven two-population models were fit to test whether the strongly differentiated northernmost populations (one in C. morettiana and P. tyrolensis, and three in S. facchinii) originate from nunatak survival or recent founder events (Figure S1).
Across all analyses, we used the optimized parameter sets of each replicate to simulate the JSFS; the multinomial approach was used to estimate the log-likelihood of the JSFS given the model.For each analysis, the Akaike Information Criterion (AIC) and log likelihoods were inspected, and ΔAIC scores were used to calculate relative likelihoods and Akaike weights (ω i ).The model with highest Akaike weight was selected as the most likely for each divergence event (Burnham & Anderson, 2002).For the two best models of each analysis, the departure of simulated and empirical JSFSs and the progression of the log-likelihood along optimization rounds were inspected.
To verify that our models were reasonable explanations of the JSFS and to estimate the divergence time between genetic groups, we performed goodness-of-fit tests of the above-mentioned two-population models following Barratt et al. (2018).The threepopulation models had to be excluded from these tests due to limited computational resources.Shortly, for each population pair modelled with δaδi, we fit the top-ranked model using our optimized parameters, scaled the resulting model JSFS by the inferred theta value and used the model JSFS to generate 100 Poisson-sampled JSFS.We then optimized each simulated JSFS to obtain a distribution of log-likelihood scores and Pearson's chi-squared test statistic and subsequently determined whether our empirical values were within these distributions.If the estimated statistics were not within the posterior distributions, we concluded that estimation of parameters under the selected model was inaccurate (Wegmann et al., 2010;Barratt et al., 2019).Estimated divergence times for each population pair were plotted and confidence intervals calculated and transformed to conventional units, using a generation time of 10 years for all study species (Gutenkunst et al., 2009).In the absence of reliable data for close relatives, we used the mutation rate of Arabidopsis thaliana (7 × 10 −9 , Ossowski et al., 2010).If the distribution of the divergence times was not Gaussian and showed extremely broad confidence intervals, we concluded that parameter estimation failed.Acknowledging the bias of AIC towards more complex models (Cavanaugh & Neath, 2019), the goodness-of-fit protocol was applied to the second-best model if the first failed, but only in the case that it had fewer parameters and the error between the empirical and simulated JSFS was not clearly worse than the best model.The same criteria as above were applied to the second-best model to estimate divergence times.
To explore population size changes associated with the last glaciation, we modelled the effective population size (N e ) using the software Stairway plot (Liu & Fu, 2015).First, we computed the folded SFS of each species and genetic group using easySFS (https:// github.com/ isaac overc ast/ easySFS) by downprojecting the data sets to a minimum sample size, maximizing the number of SNPs kept (Gutenkunst et al., 2009).We ran Stairway plot on 200 bootstrap replicates drawn from the calculated SFSs.Median N e and confidence intervals were obtained based on 200 estimations.To obtain time estimates, we used the mutation rate of Arabidopsis thaliana (7 × 10 −9 , Ossowski et al., 2010) and a generation time of 10 years, as above.

| Genetic structure and phylogenetic relationships
No intraspecific ploidy variation was encountered (Table S1).The average number of high-quality reads per sample retained after demultiplexing and quality filtering was 1.44 million (SD = 0.85).The data have been deposited in the NCBI sequence read archive (BioProject PRJNA1089485).The resulting number of SNPs per species and filtering scheme is presented in Table S2.
In the STRUCTURE analysis of C. morettiana (Figure S2a), the southwestern Carnic Prealps population 1 was split from the populations in the Dolomites at K = 2 (Figure 2a,b).At K = 5, the same southern, northwestern and northeastern groups detected in the Neighbor-Net for populations in the Dolomites were identified, with the northernmost population 18 constituting a single-population group.A certain degree of admixture was observed across all groups, particularly in the northeastern group.The fineRADstructure analysis showed similar results, with well-defined southwestern Carnic Prealps and Dolomites groups, and a high degree of substructure within the Dolomites (Figure S3a).For P. tyrolensis (Figure S2b), at K = 2, STRUCTURE resolved the southwestern Carnic Prealps populations 1-3 and populations 4-12 from the Dolomites (Figure 3a,b).
Signs of weak admixture with the southwestern Carnic Prealps populations were detected in populations 4-7.At K ≥ 4, the northernmost population 12 was differentiated.FineRADstructure showed similar results (Figure S3b).In S. facchinii (Figure S2c), the best solution at K = 2 identified a southern (populations 1, 2) and a northern group (populations 5-13); the two central populations 3 and 4 were strongly admixed.At K = 4, the northern populations 10 and 11 constituted two non-admixed single-population groups, while the geographically close populations 12 and 13 formed a separate cluster (Figure 4a,b).The northern populations 5-9 were strongly admixed, with contributions of the three northernmost clusters, while populations 3 and 4 from the center of the distribution area were strongly admixed between the northern and southern genetic groups.The fineRADstructure analysis showed similar results, with well-defined southern and northern groups, which included clearly defined populations (Figure S3c).
Similarly to the STRUCTURE analysis at K = 4, the two central populations 3 and 4 were admixed.Population genetics summary statistics for all three species are presented in Table S1 and Figures 2-4.
F I G U R E 2 Genetic structure and demographic modelling of Campanula morettiana.(a) Sampled populations and their admixture proportions at K = 5 derived from STRUCTURE analyses.Sizes of pie charts reflect the number of private alleles (PA).Numbers correspond to population identifiers.Coloured curves indicate groups of populations used for demographic modelling.Light-blue areas indicate the maximum extent of glaciers during the LGM 21 ka (based on Ehlers et al., 2011, Van Husen 1987).The map was generated with UTM 32 N projection, with datum WGS84.(b) STRUCTURE bar plots depicting the cluster memberships of each individual at the two best Ks (2 and 5).(c) Neighbor-Nets with populations colour-coded according to the most likely clustering solution in STRUCTURE analyses at K = 5.Populations with <70% assignment to a group are indicated with a black dotted line.(d) Two-population demographic modelling results; colour coding of population groups match those presented in panel (a).The left column shows a graphical representation of the selected two-population models in δaδi; "2nd best model" indicates that the depicted model does not correspond to the best model according to the Akaike Information Criterion (see the main text for an explanation).The right column shows the distribution of the estimated divergence time after the goodness-of-fit test with 100 replicates.The histogram represents the distribution, the blue curve shows the density distribution, and the red line represents the LGM (21 ka).

| Demographic modelling
The number of SNPs used for each two-population scenario is reported in Figure S7; for the three-population model in S. facchinii, 30,450 SNPs were used.The log-likelihood increased along optimization rounds as expected for all two-and three-population models simulated with δaδi.The departure of simulated and empirical JSFSs was in some cases strong, indicating potential model fit problems (Figures S7 and S8).The goodness-of-fit tests indicated in some cases that the best model according to the AIC was not able to accurately estimate the parameters, indicated by broad distributions of the time estimation and/or deviations of the empirical log-likelihood or chisquared from the estimated distribution (Figure S7).In some of these cases, the second-best model provided a more reliable time estimate.
The divergence between populations from the Dolomites and the southwestern Carnic Prealps was best described by founder event models.However, in C. morettiana, the estimated proportion of founder individuals from the original population was 0.5, and we therefore considered it vicariance.Moreover, founder events are per default selected when a population experienced exponential growth after divergence, which could as well occur after vicariance.Therefore, we are cautious in the discrimination between vicariance and founder events.In both cases, the divergence largely predated the LGM (Figures 2d, 3d and S7).Within the Dolomites, modelling the split between the main groups of C. morettiana resulted in poor estimates using the best model according to the AIC, indicating lack of ability of the model to estimate the parameters, and therefore, the second-best model was used.The divergence times between the northwestern and southern, as well as between the northern and southern populations predated the LGM.Since no clear grouping of P. tyrolensis populations within the Dolomites arose from the genetic structure analyses, no demographic modelling was conducted at this level.In the case of S. facchinii, the second-best three-population model showed a much better fit between the modelled and the empirical JSFS (Figure S8).Both models indicated secondary contact as the cause for the admixed pattern observed in the two central populations 5 and 6.
The divergence of the northernmost population 18 from the northern group in C. morettiana occurred by vicariance long before the LGM (Figure 2d).The poor performance of the goodness-of-fit test (wide confidence intervals, deviation of the log-likelihood from its distribution) of the analysis with the southern group as potential sister indicates that the population pair probably does not fit true divergence (i.e., the true divergence occurred between the northern group and population 18, Figure S7).Similarly, inconclusive results were obtained for the northernmost population 12 of P. tyrolensis and for the northeasternmost population 10 of S. facchinii (Figure S7).Populations 11 and 12-13 of S. facchinii diverged via vicariance events in strict isolation from the northern group before the LGM (Figure 4d).
The number of SNPs used for the Stairway plot simulations of the N e are reported in Table S3.The results revealed increases towards the present in all studied populations of C. morettiana and P. tyrolensis, with exception of the northernmost populations 18 of C. morettiana and 12 of P. tyrolensis, which showed no significant changes and a narrow confidence interval (Figure 5a,b).These populations had particularly small sample sizes, which could negatively impact the estimations, and the results are therefore considered spurious (see Discussion, Table S1, Liu & Fu, 2015).These N e increases occurred at different times before the LGM, in most cases early enough to accommodate for the uncertainty introduced by imprecise mutation rates and generation times.Saxifraga facchinii showed an ancestral N e , which either suffered a temporal drop followed by recovery (southern and northern group), or gradually decreased towards the present (northernmost populations 10-13 and central populations 3-4; Figure 5c).

| The Piave Valley separates refugia in the Dolomites and the Carnic Prealps
Biogeographers have, since the origin of the discipline, identified barriers all over the world, which define distribution limits across different organisms (Wallace, 1863).Barriers coinciding with deep valleys have been identified in many mountainous regions of the world; at the intraspecific level, these valleys define genetic splits, which are often congruent across codistributed species (Antonelli et al., 2009;Qiu et al., 2011;Španiel & Rešetnik, 2022).The role of F I G U R E 4 Genetic structure and demographic modelling of Saxifraga facchinii.(a) Sampled populations and their admixture proportions at K = 4 for derived from STRUCTURE analyses.Sizes of pie charts reflect the number of private alleles (PA).Numbers correspond to population identifiers.Coloured curves indicate groups of populations used for demographic modelling.Light-blue areas indicate the maximum extent of glaciers during the LGM 21 ka (based on Ehlers et al., 2011;van Husen, 1987).The map was generated with UTM 32 N projection, with datum WGS84.(b) STRUCTURE bar plots depicting the cluster memberships of each individual at the two best Ks (2 and 5).(c) Neighbour-Nets with populations colour-coded according to the most likely clustering solution in STRUCTURE analyses at K = 4. Populations with <70% assignment to a group are indicated with a black dotted line.(d) Two-and three population demographic modelling results; colour coding of population groups match those presented in panel (a).The left column shows a graphical representation of the selected two-population models in δaδi; "2nd best model" indicates that the depicted model does not correspond to the best model according to the Akaike Information Criterion (see the main text for an explanation).The right column shows the distribution of the estimated divergence time after the goodness-of-fit test with 100 replicates.The histogram represents the distribution, the blue curve shows the density distribution, and the red line represents the LGM (20 ka).In case of three-population models, the goodness-of-fit test could not be performed, instead a graphical representation of the second-best three-population model is shown.major valleys as barriers has also been highlighted in the European Alps (Gugerli et al., 2023;Thiel-Egenter et al., 2011;Tribsch, 2004), with the Aosta Valley and the "Brenner Line" emerging as the most prominent barriers.
In this line, our data identified the Piave Valley in the Southern Limestone Alps, the geographic border between the Dolomites and the Carnic Prealps (Figure 1), as the main genetic break in C. morettiana and P. tyrolensis.This strong differentiation is reflected in clusters in the Structure analyses, main groups in the Neighbor-Nets and fineRADstructure, and the main phylogenetic split in C. morettiana ( Figures 2, 3 and S3-S5).In addition, high levels of genetic diversity to the east as well as to the west of the Piave (π, Table S1) and the strongly elevated numbers of private alleles in population 1 of C. morettiana (Figure 2, Table S1) are indicative for separate refugia (e.g., Carnicero et al., 2022;Schönswetter et al., 2005;Westergaard et al., 2019;Záveská et al., 2021).We emphasize that population 1 of C. morettiana harbours so many private alleles because it is the only population constituting the divergent orange-coded gene pool in Figure 2. The Piave Valley roughly coincides with major genetic breaks in other plants (Campanula barbata L., Hypochaeris uniflora Vill.and Phyteuma betonicifolium Vill., Thiel-Egenter et al., 2011; Dianthus sylvestris Wulfen, Luqman et al., 2023), as well as in the rupicolous snail Charpentiera itala (Xu & Hausdorf, 2021) or the butterfly Erebia euryale (Cupedo & Doorenweerd, 2022).
Our demographic modelling approach confirmed the existence of glacial refugia east and west of the Piave Valley, by supporting divergence times between the groups largely predating the LGM (Figures 2, 3 and S1).These results are the first molecular-based evidence of separate refugia for calcicolous species in the Dolomites and the southwestern Carnic Prealps.To date, only the siliceous parts of the southern Dolomites were identified as glacial refugium for silicicolous mountain plants (Escobar García et al., 2012;György et al., 2018;Ronikier et al., 2008;Schönswetter et al., 2005).The valley floor obviously acted as an efficient-albeit not impermeable (Figures 2 and 3) barrier, promoting the observed differentiation.
Most of the time, it was likely unsuitable for alpine species, either due to the presence of a glacier tongue during glacials, or due to expansion of forests during interglacials (Magri et al., 2015;Seguinot et al., 2018).For lower-elevation species linked to forests such as Aposeris foetida and Helleborus niger or the amphibian Salamandra atra, the Piave Valley delimited no major intraspecific genetic groups (Bonato et al., 2018;Záveská et al., 2021;Voisin et al., in prep.).In 4-11 in P. tyrolensis; Figures 2 and 3) suggests a single peripheral refugium in the southern Dolomites.In C. morettiana, the western (10-12) and the northern populations (13-18) are clearly no postglacial descendants of this refugium as their divergence largely predates the LGM (Figure 2d).These findings reject the hypothesized postglacial northward migration from a peripheral refugium in this species, and rather suggest survival on nunataks in addition to the peripheral refugium.The marked bottleneck observed in the northern populations of C. morettiana (Figure 5a) may be indicative of the small spatial extent of this refugium, presumably caused by the high elevation of the LGM glacier surface in the central Dolomites (van Husen, 1987).In P. tyrolensis, the differentiation of the northernmost population 12 as a unique genetic group at K ≥ 4 is likely postglacial (discussed below).
For subnival S. facchinii, restricted to interior parts of the Dolomites, several glacial refugia on nunataks were expected.This hypothesis was confirmed by the presence of four, partially geographically correlated genetic clusters (Figure 4a,b), as well as by the divergence times among population groups preceding the LGM (Figure 4d).Further, our models indicate that the strong admixture in the central populations 3 and 4 results from secondary contact of the southern and northern groups with the central populations, and not from a recent admixed origin (Figure 4d).Altogether, our models support the presence of several nunatak refugia for S. facchinii in the north, centre and south of the species' distribution area.
Interestingly, the three northernmost massifs hosting populations 10-13 constituted divergent unique genetic groups at higher Ks.
Indeed, a congruent pattern across the three study species is the divergence of the northernmost populations at high Ks (population 18 of C. morettiana; population 12 of P. tyrolensis; populations 10-13 of S. facchinii; Figures 2-4).Divergence of populations is often caused by strong genetic drift, which may result from different demographic histories (Lawson et al., 2018).Specifically, divergent populations may be old populations (i.e., predating the LGM), which strongly drifted due to small population size and lack of gene flow with neighbouring groups.Alternatively, divergence may be caused by strong drift due to bottlenecks related to recent founder events (i.e., postglacial colonization).In population 18 of C. morettiana and populations 11-13 of S. facchinii, the old divergence times allow concluding that these populations diverged before the LGM (Figures 2d, 4d and S7), and therefore survived the LGM in nunatak refugia (Holderegger & Thiel-Egenter, 2009).In contrast, models were inconclusive for the northernmost population of P. tyrolensis and for population 10 of S. facchinii (Figure S1).
While previous studies have demonstrated a prevalence of peripheral glacial refugia (Schönswetter et al., 2005), nunatak survival has been invoked in several studies on alpine or arctic species as the most parsimonious hypothesis to explain observed genetic structure patterns (Bettin et al., 2007;Lohse et al., 2011;Schönswetter & Schneeweiss, 2019;Stehlik et al., 2001Stehlik et al., , 2002;;Westergaard et al., 2010).To date, the most unequivocal support for nunatak survival came from range-wide studies unravelling unique divergent genotypes in areas with strong LGM glaciation, while widespread genotypes dominated in potential peripheral refugia (Escobar García et al., 2012;Stehlik et al., 2002).Recently, the application of spatially explicit demographic models has provided independent support for nunatak survival in different geographic settings (Carex chalciolepis

| Disentangling artefacts from the true demographic history
The inclusion of goodness-of-fit tests in our demographic modelling approach proved highly useful.Among others, they allowed us to approximate confidence intervals around divergence times estimates.
Most importantly, goodness-of-fit tests provided a double-check of model accuracy (Barratt et al., 2018); often, the best model selected based on the AIC showed a deficient fit and inability to estimate the parameters of the model (Figures S1,S7 and S8).While the exact reason for the inaccuracy cannot be ascertained, some factors may increase the probability of a misleading model selection (Adams & Hudson, 2004;Robinson et al., 2014).A central one is lack of statistical power, understood here as the ability to detect the true model.
It results from an insufficient amount of data to inform model selection, either due to a low number of sampled individuals per population, insufficient genomic information (i.e., number of SNPs), or both.
Previous simulation studies have highlighted the importance of sufficient sample size for selecting the correct model and estimating the parameters (Adams & Hudson, 2004;Liu & Fu, 2015;Robinson et al., 2014).While old and pronounced demographic events were accurately modelled in these studies based on a sampling of three to five individuals per genetic group, larger sample sizes were necessary for correctly describing recent and more subtle events.
Moreover, these simulations were performed based on a number of SNPs greater than 10,000, much higher than here and in most RADseq-based publications.In these cases, a higher sample size (Lou et al., 2021;Robinson et al., 2014) should be considered to achieve satisfactory results.In our case, the demographic analyses of single northernmost populations do clearly not meet these thresholds as sample sizes range between four and six individuals.Thus, the stable effective population sizes with a narrow confidence interval observed in the northernmost population 18 of C. morettiana, population 12 of P. tyrolensis and population 11 of S. facchinii (Figure 5) are interpreted as artefacts.Additionally, the bias of AIC towards more complex models is well-known (Cavanaugh & Neath, 2019), and it is therefore recommended to compare the best model with models with fewer parameters.In cases braving the edge of recommended sampling size and number of SNPs, we strongly advocate a careful and conservative interpretation of the results, application of additional tests such as the goodness-of-fit, and visual inspection of simulated and empirical SFSs.
Here, we suggest that altitudinal segregation of our study species is causing the major differences in the temporal fluctuation of effective population size across species.On the one hand, C. morettiana and P. tyrolensis showed reduced effective population sizes in the past, followed by an increase towards the present (Figure 5a,b), which we interpret as a recovery of population size after unsuitable conditions during glacial periods.
In contrast, S. facchinii exhibited similar or higher effective population sizes in the past than at present for all main genetic groups (Figure 5c).This pattern suggests the availability of large climatically suitable areas for the species' survival during cold periods, in some cases even allowing for larger effective population sizes than at present (populations 3-4, 10, and 12-13; Figure 5c).However, the marked, short-time bottlenecks for the northern (5-9) and southern (1, 2) populations indicate a short period of reduced climatic suitability, probably corresponding to particularly harsh conditions at the LGM.We suggest that the plateau morphology of the mountains harbouring these particular populations might have favoured formation of permanent snowfields and ice caps during the LGM, which strongly reduced the available area for S. facchinii during a relatively short period.
Taken together, the different demographic trajectories reconstructed for the three study species can be largely explained by their altitudinal segregation and the landscape of the Dolomites during glacial periods.Then, glaciers covered all main valleys, leaving icefree areas above them (van Husen, 1987) with harsh conditions, likely comparable to those of the nival belt at present (Körner, 2021).
There, S. facchinii might have found enough suitable habitat to maintain populations, which might not have been much smaller than at present.Instead, C. morettiana and P. tyrolensis likely survived in patches of suitable habitat much smaller than the current ones, resulting in a significant postglacial increase of the effective population size.

| CON CLUS IONS
We provide the first evidence for separate refugia for calcicolous species in the Dolomites and the southwestern Carnic Prealps, which reflect the area's geography as well as previously identified areas of endemism (Tribsch, 2004).We additionally highlight the importance of nunatak survival in the Dolomites in addition to peripheral refugia (Holderegger & Thiel-Egenter, 2009;Pan et al., 2020).
Glacial survival patterns strongly depend on ecological preferences of the study species with nunatak survival in S. facchinii, nunatak and peripheral refugia for C. morettiana, and no evidence for nunatak survival for P. tyrolensis.Altitudinal segregation strongly influenced the species' demographic fluctuations, with montane to alpine C. morettiana and P. tyrolensis showing higher effective population sizes at present than during glacial periods, and high-alpine to subnival S. facchinii showing the opposite trend, as expected for cold-hardy species (Gentili et al., 2015).Whereas it has been argued that cold-hardy species can respond to warming through shifts into topographical micro-refugia maintaining cold conditions (Körner & Hiltbrunner, 2021), we show how warming since the LGM had a negative impact on effective population size of several populations of S. facchinii.In this line, the IUCN threat category of the species has been recently raised to EN, reflecting its range fragmentation and projected decline (Orsenigo et al., 2023).Even if cold-hardy species can survive the next few decades in micro-refugia, the decrease of effective population size could reach a critical point, eventually leading to extinction (Dullinger et al., 2012;Kardos et al., 2021), further increasing the need for monitoring and conservation actions (Hoban et al., 2022).
tane to alpine Campanula morettiana Rchb.(Campanulaceae) and Primula tyrolensis Schott ex Rchb.(Primulaceae), and high alpine to nival Saxifraga facchinii Koch (Saxifragaceae, Figure 1).We focussed on the following hypotheses.(1) The deep, northeast-to-southwest running Piave Valley does not only form a major break in the distribution ranges of C. morettiana and P. tyrolensis (Figure 1), but also coincides with the deepest genetic split in these species.As unglaciated or weakly glaciated limestone areas existed on both sides of the valley, we accordingly hypothesize one peripheral refugium in the southern Dolomites to the west of the Piave river and one in the southwestern Carnic Prealps to the east of it.Consequently, we hypothesize divergence predating the LGM.(2) According to the observed differences in their altitudinal ranges, the study species are expected to show differential patterns of genetic structure and glacial survival within the Dolomites.Specifically, (2a), we hypothesize that montane to alpine C. morettiana and P. tyrolensis exhibit a shallow genetic structure as they likely survived in a peripheral refugium and recolonized northerly adjacent areas after the LGM through founder events.In contrast (2b), high-alpine to subnival S. facchinii is expected to show a marked genetic structure as it likely survived in several separate nunatak refugia in central and northern massifs of the Dolomites.(3) Finally, we hypothesize that C. morettiana and P. tyrolensis suffered a massive glaciation-induced population decline followed by a strong population increase because much of their present-day distribution ranges were covered by ice during the LGM.Cold-hardy S. facchinii is expected to show bottlenecks related to survival in nunatak refugia, but in contrast to the other two study species, a less pronounced postglacial population F I G U R E 1 Study area and investigated species.(a) View of the Dolomites from the north (Conturines massif) towards the southwest, showing the typical geomorphology which results in a sky-island system of alpine habitats.(b) Density of occurrence along the elevation gradient for all known populations of each species.(c) Distribution ranges of Campanula morettiana (blue dots and solid line), Primula tyrolensis (pink dots and short-dashed line) and Saxifraga facchinii (red dots and long-dashed line).White triangles highlight the Piave Valley.Green-coloured areas represent carbonate bedrock, which is suitable for the study species.The grey scale in the background represents the elevation from 100 m to 3343 m a.s.l.Maps were generated with UTM 32 N projection, with datum WGS84.Photo of landscape, C. morettiana and S. facchinii by F. Rota, photograph of P. tyrolensis by J. Nascimbene.

F
Genetic structure and demographic modelling of Primula tyrolensis.(a) Sampled populations and their admixture proportions at K = 4 derived from STRUCTURE analyses.Sizes of pie charts reflect the number of private alleles (PA).Numbers correspond to population identifiers.Coloured curves indicate groups of populations used for demographic modelling.Light-blue areas indicate the maximum extent of glaciers during the LGM 21 ka (based on Ehlers et al., 2011, Van Husen 1987).The map was generated with UTM 32 N projection, with datum WGS84.(b) STRUCTURE bar plots depicting the cluster memberships of each individual at the two best Ks (2 and 4).(c) Neighbor-Nets with populations colour-coded according to the most likely clustering solution in STRUCTURE analyses at K = 4.(d) Two-population demographic modelling results; colour coding of population groups match those presented in panel (a).The left column shows a graphical representation of the selected two-population models in δaδi.The right column shows the distribution of the estimated divergence time after the goodness-of-fit test with 100 replicates.The histogram represents the distribution, the blue curve shows the density distribution, and the red line represents the LGM (21 ka).

FIGURE 4
FIGURE 4 Legend on next page

F
I G U R E 5 Demographic modelling analyses of (a) Campanula morettiana, (b) Primula tyrolensis and (c) Saxifraga facchinii.The upper panels illustrate the geographic configuration of the modelled population groups, indicated by a coloured outline.Black dots represent sampled populations.The base map is a digital elevation model overlaid with the extension of glaciers at the Last Glacial Maximum (LGM).The lower panels show the change of the effective population size (N e ) over time modelled with the software Stairway plot.Colour coding of population groups match those presented in the maps.The coloured curves indicate the fluctuation of N e of the analysed population groups, 95% confidence intervals are shown as a shaded background.
line with this finding, previous studies (e.g., Strait of Gibraltar: Nieto Feliner, 2014; the Andes: Antonelli & Sanmartín, 2011) have demonstrated that the same geographic feature can act as either barrier or corridor for species with different ecological preferences.4.2 | Evidence for peripheral and nunatak refugia within the Dolomites Within the Dolomites, the retrieved genetic structure supported species-specific glacial histories.Whereas survival in a peripheral refugium was demonstrated both for C. morettiana and for P. tyrolensis, nunatak survival was supported for northern populations of C. morettiana and several population groups within S. facchinii.Specifically, for C. morettiana and P. tyrolensis, the lack of differentiation among the southern populations (2-9 in C. morettiana;
., C.W., P.C. and P.S. were involved in conceptualization and sampling design.F.R. (lead), J.N. and C.W. were involved in fieldwork.F.R. was involved in labwork.P.C., F.R. and G.C. were involved in data analysis.F.R., P.C., P.S., G.C. and C.W. were involved in investigation.F.R., P.C., P.S., G.C., J.N. and C.W. were involved in writing, review and editing.C.W. was involved in funding.All authors read and approved the final manuscript.