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.