Survival in nunatak and peripheral glacial refugia of three alpine plant species is partly predicted by altitudinal segregation
Running title: Phylogeography of Southeastern Alps endemics
Authors: Francesco Rota1,2*, Pau Carnicero3*, Gabriele Casazza4, Juri Nascimbene5, Peter Schönswetter3, Camilla Wellstein1
*Francesco Rota and Pau Carnicero should be considered joint first author
Affiliations:
1 Faculty of Agricultural, Environmental and Food Sciences, Free University of Bozen‐Bolzano, Bolzano, Italy
2 Swiss Federal Institute for Forest, Snow and Landscape Research WSL, Zürcherstrasse 111, 8903 Birmensdorf, Switzerland
3 Department of Botany, University of Innsbruck, Innsbruck, Austria
4 Department of Earth, Environmental and Life Sciences (DISTAV), University of Genoa, Corso Europa 26, I-16132 Genova, Italy
5 BIOME Group, Department of Biological, Geological and Environmental Sciences, Alma Mater Studiorum - University of Bologna, Bologna, Italy
Correspondence:
Pau Carnicero, Department of Botany, University of Innsbruck, Innsbruck, Austria. Email:pau.carnicero@gmail.com
Camilla Wellstein, Faculty of Agricultural, Environmental and Food Sciences, Free University of Bozen‐Bolzano, Bolzano, Italy. Email:camilla.wellstein@unibz.it
Francesco Rota, Swiss Federal Institute for Forest, Snow and Landscape Research WSL, Birmensdorf, Switzerland. Email:francesco.rota@wsl.ch
Abstract
Mountain biota survived the Quaternary cold stages 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 andPrimula 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 modeling 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.
Keywords
endemic alpine plants, demographic modeling, Dolomites, Last Glacial Maximum, phylogeography, RADseq
Introduction
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 (van Husen, 1987; Schönswetter et al., 2005). 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., 2011), 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 within-population 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 modeling 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., 2011). After the LGM, the extreme conditions necessary for their survival were found only in the highest elevations, which limited the magnitude of postglacial population increase.
Our study area is situated in the Southern Limestone Alps, i.e., 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 (van Husen, 1987; Seguinot et al., 2018) 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 delimits an area of endemism (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 (Gutenkunst et al., 2009) and combined this with effective population size modeling (Liu & Fu, 2015). The study species include montane 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 focused 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 increase is expected, as its occurrence is limited to small areas providing high-alpine to nival conditions.
Material and Methods
Study species
We selected our study species C. morettiana , P. tyrolensisand 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 forC. 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 ten 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 and 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 andP. 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 single-end 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 analyzed 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/wtsinpg/illumina2bam) and the program process_radtags.pl in STACKS v. 2.3 (Catchen et al., 2011; 2013). RADseq loci were further assembled, and SNPs were called using thedenovo_map.pl pipeline 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.pl for values of the number of mismatches allowed between stacks to merge them into a putative locus (-M ) from 0–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 in the R package ADEGENET 2.1.3 (Jombart, 2008) 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 andpopulations was run again for the new dataset. The vcf files were further filtered with VCFtools 0.1.16 to 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 dataset (-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 optimalK 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 dataset 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.
Adegenet was used to calculate Nei’s distance matrices (Nei, 1972), basis to construct Neighbor-Nets (Bryant & Moulton, 2004) in SplitsTree 4.14.2 (Huson & Bryant, 2006). To infer phylogenetic relationships among individuals, we computed maximum likelihood (ML) phylogenetic analyses using RAXML v. 8.2.12 (Stamatakis, 2014). Invariant sites were removed from the original phylip format using the script “deleteAlignColumn.pl” (https://www.biostars.org/p/55555/) and Felsenstein’s ascertainment bias correction was used to account for missing invariant sites as recommended (Leaché et al. 2015). Tree searches were done under a General Time Reversible model with disabled rate heterogeneity among sites as recommended in the RAXML v. 8.2.X manual (ASC_GTRCAT; -V; Stamatakis, 2014). The best-scoring ML tree was bootstrapped using 1000 replicates and the frequency-based stopping criterion (Pattengale et al., 2010). Results were visualized with FIGTREE 1.4 (http://http://tree.bio.ed.ac.uk/software/figtree/).
Demographic modeling
To explore alternative glacial survival histories among genetic groups we used the diffusion approximation method δaδi (Gutenkunst et al., 2009) to analyze 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 analyzed 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 inC. 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 (Portik et al., 2017; Barratt et al., 2018; Firneno et al., 2020) 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 andP. tyrolensis , and three in S. facchinii ) originate from nunatak survival or recent founder events (Figure S1).
For all models, we performed consecutive rounds of optimizations to ensure convergence on similar log-likelihood scores (Portik et al., 2017). For each round, we ran multiple replicates using parameter estimates from the best-scoring replicate (highest log-likelihood) to seed searches in the following round. We used the default settings in the dadi_pipeline for each round (replicates = 10, 20, 30, 40; maxiter = 3, 5, 10, 15; fold = 3, 2, 2, 1), and optimized parameters using the Nelder-Mead method (optimize_log_fmin). 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 AIC (Akaike Information Criterion) 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 three-population models had to be excluded from these tests due to limited computational resources. Shortly, for each population pair modeled 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 modeled the effective population size (Ne ) 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/isaacovercast/easySFS) by downprojecting the datasets 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. MedianNe 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.
Results
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 #####). 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). Low admixture was 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.
The Neighbor-Net and the outgroup-rooted RAXML phylogenetic tree ofC. morettiana showed two main groups, one in the southwestern Carnic Prealps (population 1, BS 100%) and one in the Dolomites (2–18, Figures 2c, S4, BS 100%). Within the Dolomites, a southern (2–9), a northwestern (10–12) and a northeastern group (13–18) could be differentiated in the Neighbor-Net (Figure 2c). Similar to C. morettiana , the Neighbor-Net of P. tyrolensis showed two main groups, one in the southwestern Carnic Prealps (1–3) and one in the Dolomites (Figure 3c). In the RAXML tree these two groups had low support (BS 51–52%; Figure S5). The Neighbor-Net of S. facchinii showed a strong differentiation of the southern populations (1, 2), a northern group (5–13) with star-like pattern, and the central populations (3, 4) placed along the splits linking the northern and southern groups (Figure 4c). In the RAXML tree, the southern populations 1 and 2 were sister to the central and northern populations (3–13), which constituted a monophyletic group (Figure S6, BS 100%). The central populations were not monophyletic and formed the sister group of the northern clade (5–13, BS 100%). Populations 2, 3, 6, 10, 11, 12 and 13 constituted monophyletic groups with BS > 95%.
Demographic modeling
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, 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 Chi-squared 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, S7). Within the Dolomites, modeling 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 modeling was conducted at this level. In the case of S. facchinii , the second-best three-population model showed a much better fit between the modeled 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 theNe 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 Ne 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 Ne , 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).
Discussion
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 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 inC. morettiana (Figures 2, 3, S3–S5). In addition, high levels of genetic diversity (π, Table S1) and elevated numbers of private alleles (Figures 2, 3, Table S1) to the east as well as to the west of the Piave 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). The Piave Valley roughly coincides with major genetic breaks in other plants (Campanula barbata L., Hypochaeris unifloraVill. 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 modeling 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; Figure 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 (Schönswetter et al., 2005; Ronikier et al., 2008; Escobar García et al., 2012; György et al., 2018). The valley floor obviously acted as an efficient – albeit not impermeable (Figures 2, 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 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.
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 ofC. 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 ; 4–11 in P. tyrolensis ; Figures 2, 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, center and south of the species’ distribution area. Interestingly, the three northernmost massifs hosting populations 10–13 constituted divergent unique genetic groups at higher K s.
Indeed, a congruent pattern across the three study species is the divergence of the northernmost populations at high K s (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 neighboring groups. Alternatively, divergence may be caused by strong drift due to bottlenecks related to recent founder events (i.e., postglacial colonization). In population 18 ofC. morettiana and populations 11–13 of S. facchinii , the old divergence times allow concluding that these populations diverged before the LGM (Figures 2d, 4d, 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., 2001; Stehlik et al., 2002; Westergaard et al., 2011). To date, the most unequivocal support for nunatak survival came from range-wide studies unraveling unique divergent genotypes in areas with strong LGM glaciation, while widespread genotypes dominated in potential peripheral refugia (Stehlik et al., 2002; Escobar García et al., 2012). Recently, the application of spatially explicit demographic models has provided independent support for nunatak survival in different geographic settings (Carex chalciolepis Holm.: Massatti & Knowles, 2014; Pedicularis aspleniifolia Flörke ex Willd.: Pan et al., 2020).
Disentangling artifacts from the true demographic history
The inclusion of goodness-of-fit tests in our demographic modeling approach proved highly useful. Amongst 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, S8). While the exact reason for the inaccuracy cannot be ascertained, some factors may increase the probability of a misleading model selection (Adams & Huson, 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 & Huson, 2004; Liu & Fu, 2015; Robinson et al., 2014). While old and pronounced demographic events were accurately modeled 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 (Robinson et al., 2014; Lou et al., 2021) 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 ofC. morettiana , population 12 of P. tyrolensis and population 11 of S. facchinii (Figure 5) are interpreted as artifacts. 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.
Altitudinal segregation largely determines glaciation-induced demographic fluctuations
Differences in species’ ecological preferences can determine divergent demographic histories in codistributed species (Carnicero et al., 2022; Massatti & Knowles, 2014). In mountain regions, ecological differences are often reflected in altitudinal segregation of species, which may confer survival in different glacial refugia (‘displacement refugia model’; Kropf et al., 2003; Surina et al., 2011). 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 harboring these particular populations might have favored 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 ice-free areas above them (Van Husen, 1987) with harsh conditions, likely comparable to those of the nival belt at present (Körner, 2021). There, S. facchiniimight 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.
Conclusions
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).
References
Adams, A. M., & Hudson, R. R. (2004). Maximum-likelihood estimation of demographic parameters using the frequency spectrum of unlinked single-nucleotide polymorphisms. Genetics , 168 (3), 1699–1712.https://doi.org/10.1534/genetics.104.030171
Aeschimann, D., Lauber, K., Moser, D. M., & Theurillat, J. P. (2004). Flora alpina. Haupt Verlag.
Aeschimann, D., Rasolofo, N., & Theurillat, J. P. (2011). Analyse de la flore des Alpes. 1: historique et biodiversité. Candollea ,66 (1), 27-55.https://doi.org/10.15553/c2011v661a2
Antonelli, A., Nylander, J. A. A., Persson, C., & Sanmartín, I. (2009). Tracing the impact of the Andean uplift on Neotropical plant evolution.PNAS , 106 (24), 9749–9754.https://doi.org/10.1073/pnas.0811421106
Antonelli, A., & Sanmartín, I. (2011). Why are there so many plant species in the Neotropics? Taxon , 60 (2), 403–414.https://doi.org/10.1002/tax.602010
Barratt, C. D., Bwong, B. A., Jehle, R., Liedtke, H. C., Nagel, P., Onstein, R. E., Portik, D. M., Streicher, J. W., & Loader, S. P. (2018). Vanishing refuge? Testing the forest refuge hypothesis in coastal East Africa using genome-wide sequence data for seven amphibians. Molecular Ecology , 27 (21), 4289–4308.https://doi.org/10.1111/mec.14862
Bettin O, Cornejo C, Edwards PJ, Holderegger R. (2007). Phylogeography of the high alpine plant Senecio halleri (Asteraceae) in the European Alps: in situ glacial survival with postglacial stepwise dispersal into peripheral areas. Molecular Ecology ,16 (12), 2517–2524.https://doi.org/10.1111/j.1365-294X.2007.03273.x.
Bonato, L., Corbetta, A., Giovine, G., Romanazzi, E., Šunje, E., Vernesi, C., & Crestanello, B. (2018). Diversity among peripheral populations: genetic and evolutionary differentiation ofSalamandra atra at the southern edge of the Alps. Journal of Zoological Systematics and Evolutionary Research , 56 (4), 533–548.https://doi.org/10.1111/jzs.12224
Bosellini, A., Gianolla, P., & Stefani, M. (2003). Geology of the Dolomites. Episodes , 26 (3), 181–185.https://doi.org/10.18814/epiiugs/2003/v26i3/005
Bryant, D., & Moulton, V. (2004). Neighbor-Net: an agglomerative method for the construction of phylogenetic networks. Molecular Biology and Evolution , 21 (2), 255–265.https://doi.org/10.1093/molbev/msh018
Burnham, K. P., & Anderson, D. R. (Eds.). (2002). Model selection and multimodel inference (2nd ed.). Springer New York.https://doi.org/10.1007/b97636
Carnicero, P., Wessely, J., Moser, D., Font, X., Dullinger, S., & Schönswetter, P. (2022). Postglacial range expansion of high-elevation plants is restricted by dispersal ability and habitat specialization.Journal of Biogeography , 49 , 1739–1752.https://doi.org/10.1111/jbi.14390
Catchen, J. M., Amores, A., Hohenlohe, P., Cresko, W., & Postlethwait, J. H. (2011). Stacks: building and genotyping loci de novo from short-read sequences. G3: Genes, Genomes, Genetics , 1 (3), 171–182.https://doi.org/10.1534/g3.111.000240
Catchen, J., Hohenlohe, P. A., Bassham, S., Amores, A., & Cresko, W. A. (2013). Stacks: An analysis tool set for population genomics.Molecular Ecology , 22(11).https://doi.org/3124–3140.10.1111/mec.12354
Cavanaugh, J. E., & Neath, A. A. (2019). The Akaike information criterion: Background, derivation, properties, application, interpretation, and refinements. WIREs Computational Statistics ,11 (3), e1460.https://doi.org/10.1002/wics.1460
Charles, K. L., Bell, R. C., Blackburn, D. C., Burger, M., Fujita, M. K., Gvoždík, V., Jongsma, G. F. M., Kouete, M. T., Leaché, A. D., & Portik, D. M. (2018). Sky, sea, and forest islands: Diversification in the African leaf-folding frog Afrixalus paradorsalis (Anura: Hyperoliidae) of the Lower Guineo-Congolian rain forest. Journal of Biogeography , 45 (8), 1781–1794.https://doi.org/10.1111/jbi.13365
Clark, P. U., Dyke, A. S., Shakun, J. D., Carlson, A. E., Clark, J., Wohlfarth, B., Mitrovica, J. X., Hostetler, S. W., & Mccabe, A. M. (2009). The Last Glacial Maximum. Science , 325 , 710–714.https://doi.org/10.1126/science.1172873
Cupedo, F., & Doorenweerd, C. (2022). Mitochondrial DNA-based phylogeography of the large ringlet Erebia euryale (Esper, 1805) suggests recurrent Alpine-Carpathian disjunctions during Pleistocene (Nymphalidae, Satyrinae). Nota Lepidopterologica , 45 , 65–86.https://doi.org/10.3897/nl.45.68138
Doyle, J. J., & Doyle, J. L. (1987). A rapid DNA isolation procedure for small quantities of fresh leaf tissue. Phytochemistry Bulletin , 19, 11–15.
Dullinger, S., Gattringer, A., Thuiller, W., Moser, D., Zimmermann, N. E., Guisan, A., Willner, W., Plutzar, C., Leitner, M., Mang, T., Caccianiga, M., Dirnböck, T., Ertl, S., Fischer, A., Lenoir, J., Svenning, J.-C., Psomas, A., Schmatz, D. R., Silc, U., … Hülber, K. (2012b). Extinction debt of high-mountain plants under twenty-first-century climate change. Nature Climate Change ,2 (8), 619–622.https://doi.org/10.1038/nclimate1514
Edwards, S. V., Shultz, A. J., & Campbell-Staton, S. C. (2015). Next-generation sequencing and the expanding domain of phylogeography.Journal of Vertebrate Biology , 64(3), 187–206.https://doi.org/10.25225/fozo.v64.i3.a2.2015
Ehlers, J., Gibbard, P. L., & Hughes, P. D. (Eds.). (2011). Quaternary glaciations-extent and chronology: a closer look (Vol. 15). Elsevier.
Escobar García, P., Winkler, M., Flatscher, R., Sonnleitner, M., Krejčíková, J., Suda, J., Hülber, K., Schneeweiss, G. M., & Schönswetter, P. (2012). Extensive range persistence in peripheral and interior refugia characterizes Pleistocene range dynamics in a widespread Alpine plant species (Senecio carniolicus , Asteraceae). Molecular Ecology , 21 , 1255–1270.https://doi.org/10.1111/j.1365-294X.2012.05456.x
Evanno, G., Regnaut, S., & Goudet, J. (2005). Detecting the number of clusters of individuals using the software structure: a simulation study. Molecular Ecology , 14(8), 2611–2620.https://doi.org/10.1111/j.1365-294X.2005.02553.x
Excoffier, L., Foll, M., & Petit, R. J. (2009). Genetic consequences of range expansions. Annual Review of Ecology, Evolution, and Systematics , 40 , 481–501.https://doi.org/10.1146/annurev.ecolsys.39.110707.173414
Firneno, T. J., O’Neill, J. R., Portik, D. M., Emery, A. H., Townsend, J. H., & Fujita, M. K. (2020). Finding complexity in complexes: Assessing the causes of mitonuclear discordance in a problematic species complex of Mesoamerican toads. Molecular Ecology , 29 (18), 3543–3559.https://doi.org/10.1111/mec.15496
Gentili, R., Baroni, C., Caccianiga, M., Armiraglio, S., Ghiani, A., & Citterio, S. (2015). Potential warm-stage microrefugia for alpine plants: Feedback between geomorphological and biological processes.Ecological Complexity , 21 , 87–99.https://doi.org/10.1016/j.ecocom.2014.11.006
Guerrina, M., Theodoridis, S., Minuto, L., Conti, E., & Casazza, G. (2022). First evidence of post‐glacial contraction of Alpine endemics: Insights from Berardia subacaulis in the European Alps.Journal of Biogeography , 49 (1), 79–93.https://doi.org/10.1111/jbi.14282
Guillot, G., Estoup, A., Mortier, F., & Cosson, J. F. (2005). A spatial statistical model for landscape genetics. Genetics ,170 (3), 1261–1280.https://doi.org/10.1534/genetics.104.033803
Guisan, A., & Zimmermann, N. E. (2000). Predictive habitat distribution models in ecology. Ecological modelling , 135 (2-3), 147–186.https://doi.org/10.1016/S0304-3800(00)00354-9
Gugerli, F., Brodbeck, S., Lendvay, B., Dauphin, B., Bagnoli, F., van der Knaap, W. O., Tinner, W., Höhn, M., Vendramin, G. G., Morales-Molino, C., & Schwörer, C. (2023). A range-wide postglacial history of Swiss stone pine based on molecular markers and palaeoecological evidence. Journal of Biogeography , 00 , 1– 14.https://doi.org/10.1111/jbi.14586
Gutenkunst, R. N., Hernandez, R. D., Williamson, S. H., & Bustamante, C. D. (2009). Inferring the joint demographic history of multiple populations from multidimensional SNP frequency data. PLoS Genetics , 5 (10), e1000695.https://doi.org/10.1371/journal.pgen.1000695
György, Z., Tóth, E. G., Incze, N., Molnár, B., & Höhn, M. (2018). Intercontinental migration pattern and genetic differentiation of arctic-alpine Rhodiola rosea L.: A chloroplast DNA survey.Ecology and Evolution , 8 (23), 11508–11521.https://doi.org/10.1002/ece3.4589
Hoban, S., Archer, F. I., Bertola, L. D., Bragg, J. G., Breed, M. F., Bruford, M. W., Coleman, M. A., Ekblom, R., Funk, W. C., Grueber, C. E., Hand, B. K., Jaffé, R., Jensen, E., Johnson, J. S., Kershaw, F., Liggins, L., MacDonald, A. J., Mergeay, J., Miller, J. M., … Hunter, M. E. (2022). Global genetic diversity status and trends: towards a suite of Essential Biodiversity Variables (EBVs) for genetic composition. Biological Reviews , 97(4), 1511–1538.https://doi.org/10.1111/brv.12852
Holderegger, R., & Thiel‐Egenter, C. (2009). A discussion of different types of glacial refugia used in mountain biogeography and phylogeography. Journal of Biogeography , 36 (3), 476–480.https://doi.org/10.1111/j.1365-2699.2008.02027.x
Huson, D. H., & Bryant, D. (2006). Application of phylogenetic networks in evolutionary studies. Molecular Biology and Evolution ,23 (2), 254–267. https://doi.org/10.1093/molbev/msj030
Ivy-Ochs, S., Kerschner, H., Reuther, A., Preusser, F., Heine, K., Maisch, M., Kubik, P. W., & Schlüchter, C. (2008). Chronology of the last glacial cycle in the European Alps. Journal of Quaternary Science , 23 (6–7), 559–573.https://doi.org/10.1002/jqs.1202
Jombart T (2008). “adegenet: a R package for the multivariate analysis of genetic markers.” Bioinformatics , 24 , 1403–1405.https://doi.org/10.1093/bioinformatics/btn129
Kardos, M., Armstrong, E. E., Fitzpatrick, S. W., Hauser, S., Hedrick, P. W., Miller, J. M., Tallmon, D. A., & Funk, W. C. (2021). The crucial role of genome-wide genetic variation in conservation .118 , 2104642118.https://doi.org/10.1073/pnas.2104642118
Kopelman, N. M., Mayzel, J., Jakobsson, M., Rosenberg, N. A., & Mayrose, I. (2015). CLUMPAK: a program for identifying clustering modes and packaging population structure inferences across K. Molecular Ecology Resources , 15 (5): 1179–1191.https://doi.org/10.1111/1755-0998.12387
Körner, C. (2021). Alpine plant life: Functional plant ecology of high mountain ecosystems . Springer Cham.
Körner, C., & Hiltbrunner, E. (2021). Why is the alpine flora comparatively robust against climatic warming? Diversity ,13 (8).https://doi.org/10.3390/D13080383
Kropf, M., Kadereit, J. W., & Comes, H. P. (2003). Differential cycles of range contraction and expansion in European high mountain plants during the Late Quaternary: Insights from Pritzelago alpina (L.) O. Kuntze (Brassicaceae). Molecular Ecology , 12 (4), 931–949.https://doi.org/10.1046/j.1365-294X.2003.01781.x
Lawson, D. J., van Dorp, L., & Falush, D. (2018). A tutorial on how not to over-interpret STRUCTURE and ADMIXTURE bar plots. Nature Communications , 9 (1), 1–11.https://doi.org/10.1038/s41467-018-05257-7
Leaché, A. D., Banbury, B. L., Felsenstein, J., De Oca, A. N. M., & Stamatakis, A. (2015). Short tree, long tree, right tree, wrong tree: new acquisition bias corrections for inferring SNP phylogenies.Systematic Biology , 64 (6), 1032–1047.https://doi.org/10.1093/sysbio/syv053
Liu, X., & Fu, Y. X. (2015). Exploring population size changes using SNP frequency spectra. Nature Genetics , 47 (5), 555–559.https://doi.org/10.1038/ng.3254
Lohse, K., Nicholls, J. A., & Stone, G. N. (2011). Inferring the colonization of a mountain range – refugia vs. nunatak survival in high alpine ground beetles. Molecular Ecology , 20 , 394–408.https://doi.org/10.1111/j.1365-294X.2010.04929.x
Lou, R. N., Jacobs, A., Wilder, A. P., & Therkildsen, N. O. (2021). A beginner’s guide to low‐coverage whole genome sequencing for population genomics. Molecular Ecology , 00 , 1–28.https://doi.org/10.1111/mec.16077
Luqman, H., Wegmann, D., Fior, S., & Widmer, A. (2023). Climate-induced range shifts drive adaptive response via spatio-temporal sieving of alleles. Nature Communications , 14 , 1080.https://doi.org/10.1038/s41467-023-36631-9
Magri, D., Agrillo, E., di Rita, F., Furlanetto, G., Pini, R., Ravazzi, C., & Spada, F. (2015). Holocene dynamics of tree taxa populations in Italy. Review of Palaeobotany and Palynology , 218 , 267–284.https://doi.org/10.1016/j.revpalbo.2014.08.012
Malinsky, M., Trucchi, E., Lawson, D. J., & Falush, D. (2018). RADpainter and fineRADstructure: population inference from RADseq data.Molecular Biology and Evolution , 35 (5), 1284–1290.https://doi.org/10.1093/molbev/msy023
Massatti, R., & Knowles, L. L. (2014). Microhabitat differences impact phylogeographic concordance of codistributed species: Genomic evidence in montane sedges (Carex L.) from the Rocky Mountains.Evolution , 68 (10), 2833–2846.https://doi.org/10.1111/evo.12491
Massatti, R., & Knowles, L. L. (2016). Contrasting support for alternative models of genomic variation based on microhabitat preference: Species-specific effects of climate change in alpine sedges.Molecular Ecology , 25 , 3974–3986.https://doi.org/10.1111/mec.13735
Mráz, P., Gaudeul, M., Rioux, D., Gielly, L., Choler, P., & Taberlet, P. (2007). Genetic structure of Hypochaeris uniflora (Asteraceae) suggests vicariance in the Carpathians and rapid post-glacial colonization of the Alps from an eastern Alpine refugium. Journal of Biogeography , 34 (12), 2100–2114.https://doi.org/10.1111/j.1365-2699.2007.01765.x
Nei, M. (1972). Genetic distance between populations. The American Naturalist , 106 (949), 283–292.https://doi.org/10.1086/282771
Nieto Feliner, G. (2014). Patterns and processes in plant phylogeography in the Mediterranean Basin. A review. Perspectives in Plant Ecology, Evolution and Systematics , 16 (5), 265–278.https://doi.org/10.1016/j.ppees.2014.07.002
Orsenigo, S., Calderisi, G., Cogoni, D., Mondoni, A., Nascimbene, J., Rota, F., … & Fenu, G. (2023). Global and Regional IUCN Red List Assessments: 15. Italian Botanist , 15 , 177-191.https://doi.org/10.3897/italianbotanist.15.107040
Ossowski, S., Schneeberger, K., Lucas-Lledó, J. I., Warthmann, N., Clark, R. M., Shaw, R. G., Weigel, D., & Lynch, M. (2010). The rate and molecular spectrum of spontaneous mutations in Arabidopsis thaliana . Science , 327 (5961), 92–94.https://doi.org/10.1126/science.1180677
Pan, D., Hülber, K., Willner, W., & Schneeweiss, G. M. (2020). An explicit test of Pleistocene survival in peripheral versus nunatak refugia in two high mountain plant species. Molecular Ecology ,29 (1), 172–183.https://doi.org/10.1111/mec.15316
Paris, J. R., Stevens, J. R., & Catchen, J. M. (2017). Lost in parameter space: a road map for stacks. Methods in Ecology and Evolution , 8 (10), 1360–1373.https://doi.org/10.1111/2041-210X.12775
Pattengale, N. D., Alipour, M., Bininda-Emonds, O. R. P., Moret, B. M. E., & Stamatakis, A. (2010). How many bootstrap replicates are necessary? Journal of Computational Biology , 17 (3), 337–354.https://doi.org/10.1089/cmb.2009.0179
Paun, O., Turner, B., Trucchi, E., Munzinger, J., Chase, M. W., & Samuel, R. (2016). Processes driving the adaptive radiation of a tropical tree (Diospyros , Ebenaceae) in New Caledonia, a biodiversity hotspot. Systematic Biology, 65 (2), 212–227.https://doi.org/10.1093/sysbio/syv076
Pawlowski, B. (1970). Remarques sur l’endémisme dans la flore des Alpes et des Carpates. Vegetatio , 21 (4), 181–243.
Pignatti, E., & Pignatti, S. (2016). Plant life of the Dolomites: Atlas of flora. Springer Berlin.
Pitschmann, H., & Reisigl, H. (1957) Endemische Blutenpflanzen der Sudtiroler Dolomiten. Veröffentlichungen des Tiroler Landesmuseums Ferdinandeum, 37 , 5–17.
Pritchard, J. K., Stephens, M., & Donnelly, P. (2000). Inference of population structure using multilocus genotype data. Genetics ,155 , 945–959.https://doi.org/10.1093/genetics/155.2.945
Prosser, F., Bertolli, A., Festi, F., & Perazza, G. (2019). Flora del Trentino. Fondazione Museo civico di Rovereto, Rovereto.
Portik, D. M., Leaché, A. D., Rivera, D., Barej, M. F., Burger, M., Hirschfeld, M., Rödel, M. O., Blackburn, D. C., & Fujita, M. K. (2017). Evaluating mechanisms of diversification in a Guineo-Congolian tropical forest frog using demographic model selection. Molecular Ecology ,26 (19), 5245–5263.https://doi.org/10.1111/mec.14266
Qiu, Y. X., Fu, C. X., & Comes, H. P. (2011). Plant molecular phylogeography in China and adjacent regions: Tracing the genetic imprints of Quaternary climate and environmental change in the world’s most diverse temperate flora. Molecular Phylogenetics and Evolution , 59 (1), 225–244.https://doi.org/10.1016/j.ympev.2011.01.012
R Core Team (2020). R: A language and environment for statistical computing. R Foundation for Statistical Computing, Vienna, Austria.https://www.R-project.org/
Ravazzi, C., Pini, R., Badino, F., De Amicis, M., Londeix, L., & Reimer, P. J. (2014). The latest LGM culmination of the Garda Glacier (Italian Alps) and the onset of glacial termination. Age of glacial collapse and vegetation chronosequence. Quaternary Science Reviews , 105 , 26–47.https://doi.org/10.1016/j.quascirev.2014.09.014
Robinson, J. D., Coffman, A. J., Hickerson, M. J., & Gutenkunst, R. N. (2014). Sampling strategies for frequency spectrum-based population genomic inference. BMC Evolutionary Biology , 14 (1).https://doi.org/10.1186/s12862-014-0254-4
Rochette, N. C., & Catchen, J. M. (2017). Deriving genotypes from RAD-seq short-read data using Stacks. Nature Protocols ,12 (12), 2640–2659.https://doi.org/10.1038/nprot.2017.123
Ronikier, M., Costa, A., Aguilar, J. F., Feliner, G. N., Küpfer, P., & Mirek, Z. (2008). Phylogeography of Pulsatilla vernalis (L.) Mill. (Ranunculaceae): Chloroplast DNA reveals two evolutionary lineages across central Europe and Scandinavia. Journal of Biogeography ,35 (9), 1650–1664.https://doi.org/10.1111/j.1365-2699.2008.01907.x
Rota, F., Casazza, G., Genova, G. Midolo G., Prosser F., Bertolli A., Wilhalm T., Juri Nascimbene J. & Wellstein C., (2022) Topography of the Dolomites modulates range dynamics of narrow endemic plants under climate change. Scientific Reports, 12 , 1398.https://doi.org/10.1038/s41598-022-05440-3
Schönswetter, P., & Schneeweiss, G. M. (2019). Is the incidence of survival in interior Pleistocene refugia (nunataks) underestimated? Phylogeography of the high mountain plant Androsace alpina(Primulaceae) in the European Alps revisited. Ecology and Evolution , 9 (7), 4078–4086.https://doi.org/10.1002/ece3.5037
Schönswetter, P., Stehlik, I., Holderegger, R., & Tribsch, A. (2005). Molecular evidence for glacial refugia of mountain plants in the European Alps. Molecular Ecology , 14 (11), 3547–3555.https://doi.org/10.1111/j.1365-294X.2005.02683.x
Seguinot, J., Ivy-Ochs, S., Jouvet, G., Huss, M., Funk, M., & Preusser, F. (2018). Modelling last glacial cycle ice dynamics in the Alps.The Cryosphere , 12 (10), 3265–3285.https://doi.org/10.5194/tc-12-3265-2018
Španiel, S., & Rešetnik, I. (2022). Plant phylogeography of the Balkan Peninsula: spatiotemporal patterns and processes. Plant Systematics and Evolution , 308 (5), 38.https://doi.org/10.1007/s00606-022-01831-1
Stamatakis, A. (2014). RAxML version 8: a tool for phylogenetic analysis and post-analysis of large phylogenies. Bioinformatics ,30 (9), 1312–1313.