Phylogeographic Insights into a Peripheral Refugium: The Importance of Cumulative Effect of Glaciation on the Genetic Structure of Two Endemic Plants

Quaternary glaciations and mostly last glacial maximum have shaped the contemporary distribution of many species in the Alps. However, in the Maritime and Ligurian Alps a more complex picture is suggested by the presence of many Tertiary paleoendemisms and by the divergence time between lineages in one endemic species predating the Late Pleistocene glaciation. The low number of endemic species studied limits the understanding of the processes that took place within this region. We used species distribution models and phylogeographical methods to infer glacial refugia and to reconstruct the phylogeographical pattern of Silene cordifolia All. and Viola argenteria Moraldo & Forneris. The predicted suitable area for last glacial maximum roughly fitted current known distribution. Our results suggest that separation of the major clades predates the last glacial maximum and the following repeated glacial and interglacial periods probably drove differentiations. The complex phylogeographical pattern observed in the study species suggests that both populations and genotypes extinction was minimal during the last glacial maximum, probably due to the low impact of glaciations and to topographic complexity in this area. This study underlines the importance of cumulative effect of previous glacial cycles in shaping the genetic structure of plant species in Maritime and Ligurian Alps, as expected for a Mediterranean mountain region more than for an Alpine region.


Introduction
Quaternary glaciations and mostly late Pleistocene glaciation (0.120-0.018 Ma) have shaped the contemporary distribution of many species [1]. In general, species survived cold adverse periods in so-called glacial refugia and later recolonised areas newly available during warmer post-glacial periods [2,3]. Populations that survived glaciations in refugia could have accumulated large amounts of genetic diversity, whereas those now found in post-glacially recolonised regions may have less diversity as a result of bottlenecks during their recent expansions [3,4].
In the last few years several phylogeographical studies have supported this general view concerning the importance of the last glacial maximum (LGM; about 0.020 Ma) in shaping species distribution [5] and contemporary intraspecific genetic diversity [6,4]. However, it has also been recently demonstrated that Mediterranean refugia represent climatically stable areas for which it is necessary to consider the cumulative effects of changes, rather than just the changes that occurred during the last glacial period [7]. The high genetic diversity detected in the Mediterranean refugia is due to the preservation of genotypes that went extinct in other places and to the intensity and accumulation of a number of processes in a patchy landscape across a varied topography [8]. A detailed understanding of the effects of past climate changes on the distribution and genetic pattern of organisms may help us to better predict the effects of ongoing climate change [9,10]. It has been recently demonstrated that the combination of paleo-distribution modelling with phylogeographical approaches may led to new interpretations of population genetic patterns and to new hypotheses about glacial survival and postglacial colonization [11,12,13].
Situated at the crossroads of the Mediterranean Basin and the Alps, the Maritime and Ligurian Alps are a hotspot for plant biodiversity [14,15] and well-known peripheral refugia [16]. Molecular investigations on the endemic plants of this region have underlined the role of vicariance events in modelling the genetic pattern of species during the last glaciation [17,18,19]. Nevertheless, Casazza et al. [13] recently demonstrated that the narrow endemic Primula allionii Loisel. was separated into two geographically disjoined groups during the Early/Mid Pleistocene border (0.781 Ma) and the region probably contains other paleoendemic species representing living examples of an ancient Tertiary flora present prior to the onset of the Mediterranean climate regime and subsequent Quaternary glaciations. However, the low number of endemic species that have been studied in this region limits our capacity to understand the processes that took place within the Maritime and Ligurian Alps [13,17,18,19,20].
In the present study we used species distribution models (SDMs) and phylogeographical methods to infer glacial refugia and to reconstruct the phylogeographical pattern of two endemic species of the area, Silene cordifolia All. and Viola argenteria Moraldo & Forneris (Fig  1). Silene cordifolia is a strict endemic species restricted to the Maritime Alps that represents an old Tertiary lineage [21] and currently occurs on siliceous cliffs and erratic rocks between 700 and 3000 m of altitude. Viola argenteria is a disjoint endemic species (sensu [22]) that has a small number of populations on Corsica (Fig 1). This species grows on scree slopes from 2000 to 3000 m altitude and may have originated more recently, as suggested by morphological, cytogenetic and molecular evidences [23]. More specifically, we asked the following questions: (1) Is there any evidence of haplotype divergence within S. cordifolia and V. argenteria before the last Pleistocene glaciation? (2) Is there any evidence that populations of S. cordifolia and V. argenteria were influenced by Pleistocene climate fluctuations, and, if so, how?

Species distribution modelling
We collected all occurrence data from the SILENE database (www.silene.eu) of the Conservatoire Botanique National Méditerranéen de Porquerolles (CBNMed), from Florence, Turin and Genoa Herbaria (FI, TO, GE) and from our own field surveys. We thus gathered a total of 369 points (266 from SILENE, 49 from Herbaria and 54 from field surveys) for S. cordifolia and 224 points (121 from SILENE, 54 from Herbaria and 49 from field surveys) for V. argenteria.
At the local scale, species distribution models provide excellent spatial projections by additionally using non-climatic predictors [24]. For this reason, a layer reporting the presence/ absence of suitable substrate was elaborated, starting from the global lithological map dataset, GLiM [25]. The bioclimatic layers for the current (30 arcsec resolution) and LGM (2.5 arcmin resolution) conditions, data were downloaded from the WorldClim 1.4 database website (http://www.worldclim.org). We used data from the Community Climate System Model (CCSM), the Model for Interdisciplinary Research on Climate (MIROC) and from the Max Planck Institute for Meteorology (MPI) to hindcast LGM (http://www.worldclim.org/paleo-climate1#com). These data were statistically downscaled from the original resolution to 30 arcsec with the Delta method [26]. This downscaling method is based on the interpolation of anomalies (differences between LGM and current conditions) at 30 arcsec resolution through the thin-plate smoothing spline algorithm. The interpolated anomalies were finally added to the layers of current condition. To reduce the multicollinearity between predictors and to minimize model overfitting, after normalizing the 19 bioclimatic predictors we performed a pairwise Pearson correlation between them (S1 Table) using R v. 3.1.1 [27]. We used the predictors that are considered physiologically important for plants [28] and that were not strongly correlated with each other (Pearson correlation coefficient, r 2 <|0.80|). To model species distributions, we used the following six variables: BIO3, isothermality; BIO4, temperature seasonality; BIO14, precipitation of driest month; BIO19, precipitation of coldest quarter; the presence/absence of suitable substrate. We removed duplicates using ENMTools [29]. Finally, 360 for S. cordifolia and 223 presence records for V. argenteria were applied in SDM analysis.
To account for model-based uncertainties, we applied six SDM techniques: multivariate adaptive regression splines (MARS- [30]), generalized linear models (GLM- [31]), classification tree analysis (CTA- [32]), flexible discriminant analysis (FDA- [33]), random forest (RF- [34]), and Maximum entropy modelling (MaxEnt- [35]). These techniques belong to three different categories of models (i.e. regression methods-MARS and GLM; classification methods-CTA and FDA; and machine learning algorithms-RF and MaxEnt). Analyses were implemented with BIOMOD2 v. 1.3.5 package [36] for R v. 3.1.1 [27]. The number of replicate sets of pseudo-absences may influence the predictive accuracy of the models. For this reason, we chose the best strategy in pseudo-absences selection according to Barbet-Massin et al. [37] tips (S2 Table). The predictive performance of the models was evaluated for each pseudo-absence run using a random subset (70%) of the initial dataset each time to calibrate the models, while the remaining 30% were used to evaluate the models. The predictive performance of the models was evaluated based on the true skill statistic (TSS), which takes into account both omission and commission errors and it is not affected by prevalence [38]. To transform the inferred continuous probability values to binary presence-absence form we used the TSS. Consensus projection under current conditions was created summing binary predictions of the six modelling techniques. We also created consensus projection under LGM conditions summing the binary predictions of the six modelling techniques for the three LGM climates, thus outlining those areas predicted by most modelling techniques under the different climates.

DNA sampling and sequencing
A total of 94 individuals from 19 populations for S. cordifolia and a total of 80 individuals from 21 populations for V. argenteria were sampled, covering the entire geographical range of species (Table 1). Permissions to collect samples in field were obtained from Parc National du Mercantour and form Parco Naturale Alpi Marittime. Total DNA was extracted from leaf tissue using the DNeasy Plant Mini Kit (Qiagen AG, Hombrechtikon, Switzerland), following the manufacturer's instructions. Seven cpDNA regions in S. cordifolia (rps16 gene, trnQ-rps16, trnG2G-trnG, trnH-psbA, trnC-ycf6R, trnT a -trnL b , trnL-trnF) and seven cpDNA regions in V. argenteria (atpF-atpH, rpoC1 gene, rps12-rpL20, trnG2G-trnG, trnH-psbA, trnL-trnF, trnT a -trnL b ) were chosen (S3 Table). PCR reactions were performed using Eppendorf thermocycler, with 20 μL total volume reactions containing 5-30 ng DNA, 10 μL GoTaq Green Master Mix (Promega, Milan, Italy), 10 pmol of each primer and purified water. We used the following PCR conditions: 2 min pre-treatment at 94˚C, followed by 35 cycles of 120 s denaturation (94˚C), 45 s annealing (48-63˚C), 75-120 s extension (72˚C). After the last cycle the temperature was kept at 72˚C for the last 7 min of extension and then lowered to 4˚C (S3 Table). The PCR products were checked on 2% agarose gels in biotium staining. The amplification products were then purified using the commercial kit QIAquick PCR Purification Kit (Qiagen, USA). All DNA sequencing was performed by Macrogen company (Korea). Sequences were aligned using MAFFT v. 6.903 [39] and adjusted by hand. The cpDNA regions were concatenated in a single matrix. We considered indels as character because previous studies exploring the utility of cpDNA indels have concluded that they are informative and that they might increase resolution between geographical regions of high diversity [40]. All sequences have been deposited in GenBank database (S4 Table).

Genetic diversity and differentiation
Haplotype number (H), haplotype diversity (H d ), number of segregating sites (S) and nucleotide diversity (π) [41] were calculated with DNASP 4.10 [42]. The total gene diversity (H T ), population differentiation (G ST ) and an estimation of population subdivision for phylogenetically ordered alleles (N ST ) were calculated using the program PERMUT [43]. We further tested the presence of phylogeographical structure (sensu [43]), comparing N ST to G ST using 1000 permutations to obtain statistical significance with the program PERMUT. To evaluate the genetic structure, we conducted an analysis of molecular variance (AMOVA) in ARLEQUIN 3.1 [44], partitioning the genetic diversity into the diversity among groups, the diversity among populations within groups, and the diversity within populations. The correlation between the genetic (F ST ) and geographic pairwise population distance matrices was evaluated using a Mantel test [45] with 10,000 permutations using ARLEQUIN 3.1. The Mantel test was Table 1. Information about the populations of S. cordifolia and V. argenteria sampled for the analysis. Population code (C), Latitude N (Lat), Longitude E (Long) and elevation (Elev.) as metres above sea level, are reported. Populations from Corsica were sampled in Florence Herbarium (FI) and the coordinates are approximate.

Taxon
Population C Lat Long Elev.

Analysis of haplotypes
To determine the genetic relationship among populations, we followed the method implemented by Muñoz-Pajares [46] in the SIDIER package. The method gathers the evolutionary information contained in insertions and deletions (indels) and combines such information with the p-distance matrix inferred from substitutions. This method improves the resolution of phylogenetic inference of low-divergence organisms, yielding results compatible with other evolutionary procedures [46]. Based on recommendations by Simmons et al. [47], indels in the concatenated dataset were coded as single mutational event using a variation of modified complex indel coding (MCIC-sensu [48]). Spatial analysis of molecular variance (SAMOVA- [49], http://web.unife.it/progetti/ genetica/Isabelle/samova.html) was used to identify groups of populations that are geographically homogeneous and maximally differentiated from each other without assuming any prior association between populations. The algorithm identifies the optimal number of groups of populations (k) by maximizing F CT (the proportion of total genetic variance due to differences among groups of populations) and minimizing F SC (genetic differentiation among populations within groups). SAMOVA was run for 10,000 iterations from each of 100 random initial conditions, and testing the predefined number of groups (k) from 2 to 10.
Best-fit models of nucleotide substitution were determined to be a TIM1 + I model for all S. cordifolia sequences and F81 model for V. argenteria, by the Akaike Information Criterion using JMODELTEST 2.1.1 [50]. Divergence time among haplotypes was estimated using BEAST v.1.8.0 [51]. We used a strict molecular clock and a coalescent Gaussian Markov random field (GMRF) Bayesian skyride tree prior. We sampled all parameters once every 2000 steps from 1 x 10 9 MCMC steps after a burn-in of 2.5 x 10 8 steps. The program TRACER [52] was used to examine convergence of chains to the stationary distribution [effective sample size (ESS) > 200]. In this study, we used normal distribution priors with a mean of 2 x 10 −9 and a standard deviation of 5 x 10 −10 substitutions per site per year to cover the cpDNA substitution rates for angiosperm (1-3 x 10 −9 substitutions per site per year; [53,54,55]).

Species distribution models
Under current climate conditions, following Swets' scale modified by Araújo et al. [56], TSS ( Table 2) indicated excellent model performance for both species. The SDMs projections for current conditions (Fig 2A and 2C) roughly fitted the current species distribution. The predicted potential distribution during the LGM was different among climates and algorithms (S1 Fig). In both species, LGM models based on CCSM and MIROC climate identified regions similar to present models, but reduced in range. On the contrary, MPI climate predicted a range wider than present in S. cordifolia while in V. argenteria it predicted a range wider than present in Corsica and narrower in the Maritime and Ligurian Alps. Despite this variation, the ensemble maps of all climates and algorithms suggested the presence of suitable areas in the current distributional range (Fig 2B and 2D).

Sequence data
The total combined length of the aligned sequences of the seven cpDNA was 4461 bp in S. cordifolia and 3730 bp in V. argenteria. We detected variations in four out of seven cpDNA regions both in S. cordifolia and in V. argenteria (Table 3). In particular, we recorded 12 variable and 9 parsimony-informative sites in S. cordifolia and 15 variable and parsimony-informative sites in V. argenteria. Alignments of the sequences of S. cordifolia and V. argenteria showed 1 and 17 insertion/deletion sites (indels) respectively. In all, there were 13 variable sites among the sampled individuals of S. cordifolia and 65 in V. argenteria (S5 Table). We recovered 12 haplotypes within S. cordifolia and 16 haplotypes within V. argenteria (Table 3).    Phylogeographical relationships and divergence times of haplotypes The dendrogram shows a clear divergence between two groups of haplotypes in both species. In S. cordifolia the first group (grey haplotypes in Fig 3A) is present in the north-western populations while the second group (white haplotypes in Fig 3C) is present throughout the entire range of distribution. In V. argenteria the first group (grey haplotypes in Fig 3B) is present in the north-western populations while the second group (white haplotypes in Fig 3D) is only present in the south-eastern populations. Diversification of S. cordifolia and V. argenteria haplotype lineages occurred during the Pleistocene, with the main lineages diverging during the middle Pleistocene. Phylogenetic analysis of cpDNA haplotypes showed that S. cordifolia differentiated into two main lineages at 0.49 (0.14-1.12) Ma. The first (H01-H05, H11) is mainly present in south-eastern populations (H01-H04) but also in some north-western populations (H05 and H11) while the second lineage (H06-H10, H12) is mainly present in north-western populations (H06-H07, H09-H10, H12) but also in some south-eastern populations (H08). Similarly, V. argenteria differentiated into two main lineages at 0.37 (0.12-0.79) Ma. The first lineage (H01 and H02) is exclusive of the north-western populations while the second lineage (H03-H16) is exclusive of the southeastern and Corsican populations. However, the node age estimates should be interpreted with caution given the width of the 95% confidence intervals for each node (Fig 4).  within S. cordifolia. Allelic differentiation (N ST = 0.749, SE = 0.0835) is greater than but not significantly different from G ST , indicating a weak genetic differentiation among populations. In V. argenteria 14 populations contained just one haplotype, five populations contained two haplotypes and the remnant two populations contained three and four haplotypes, respectively. Ten populations, mainly located in the southern part of Maritime Alps and in Corsica, showed private haplotype (4, 7, 11, 12, 16, 17, 18, 19, 20 and 21).  Fig 5) and V. argenteria (A V , B V and C V in Fig 5). Although the F CT continued to increase slightly and the F SC continued to decrease slightly as the number of groups was increased for each species, the F CT largely plateaued and the F SC markedly decreased at these respective groupings (S2 Fig). Basing on population distance, SIDIER analysis clearly separates populations into two discrete clusters: north-western and south-eastern populations (Fig 3A and 3B) in both species. In particular, in S. cordifolia the north-western cluster corresponds to SAMOVA group A S while southeastern cluster corresponds to B S and C S groups. Similarly, in V. argenteria the north-western cluster corresponds to SAMOVA group A V while south-eastern cluster corresponds to B V and C V groups of SAMOVA.

Palaeoclimatic models and haplotype divergence
In general, Alpine species are thought to have persisted LGM (0.12-0.01 Ma- [16,57]) in glacial refugia and later to have recolonised newly available areas [4,6]. As a consequence, the LGM had a strong influence on the genetic structure, distribution and evolution of plant species [2,6,58]. In fact, the LGM is recognized as a primary factor causing distributional patterns of endemic plants in Alps [5] while, the current ranges of many Alpine plants may also have been shaped by delayed Holocene recolonisation of suitable sites from refugial areas [59].
The palaeo-distribution models inferred in this work together with the low and patchy ice sheet coverage of LGM in the Maritime and Ligurian Alps [60,61] suggest that the two study species, S. cordifolia and V. argenteria, probably survived the LGM in the sites where they were already present (in situ survival sensu [57]). The weak impact of the LGM in Maritime and Ligurian Alps may have been due to the climatic mitigation effect of the nearby Mediterranean Sea and to the relatively steep relief of this sector which was probably unfavourable to the formation of major glacier basins [61]. Recently, Patsiou et al. [62] showed that Saxifraga florulenta Moretti, endemic to Maritime Alps, survived the last glaciation in several microrefugia close to sites of current occurrence. Our finding supports the idea of a weak impact of the LGM on the distribution pattern of endemic species in the Maritime and Ligurian Alps, in marked contrast to its probable effect in other parts of the Alps [2,4]. In line with the possibility of in situ survival during the LGM detected by SDM and previously debated, our analyses (Fig 4) suggest a divergence time between lineages predating the Late Pleistocene glaciation (0.115-0.018 Ma) both in S. cordifolia (0.490 Ma; 95% HPD:0.140-1.120 Ma) and V. argenteria (0.370 Ma; 95% HPD: 0.120-0.790 Ma). We thus suggest that these two species illustrate the persistence of an ancient genetic structure through several Pleistocene climatic cycles. A similar time frame has been previously suggested for other animal and plant species from other biogeographical regions [63,64]. Our results are also congruent with the palaeo-climatic reconstructions. In fact, the most extensive glaciation in Mediterranean mountains occurred during the Middle Pleistocene (0.780-0.120 Ma), in particular a major phase probably occurred during Mindel glaciation (0.470-0.420 Ma- [65,66]). Moreover, pollen-inferred historical climate indicates that maxima temperatures for interglacial stages in the Middle Pleistocene (0.780-0.120 Ma) were similar to the lowest values observed during the Late Pleistocene glacial period (0.115-0.018 Ma- [67]). In contrast, in the central Alps there is evidence of extensive glaciation during Riss (0.190-0.115 Ma- [68]) and Würm glaciation (0.115-0.018 Ma- [60]). Although the differentiation between main lineages is most likely to have begun in the Middle Pleistocene, haplotype divergence within main lineages probably occurred successively (Fig 4). This is congruent with the observation that in some Alpine animal species the separation of the major clades occurred long before the LGM and the following repeated glacial and interglacial periods may have been the motor for marked differentiation that is currently observed [69,70]. A similar time frame was detected in another endemic species in the calcareous cliffs of the Maritime and Ligurian Alps, P. allionii where the separation into two groups dates back to the Middle Pleistocene [13]. In contrast, in S. florulenta the differentiation of two ancestral gene pools seems to have been influenced by range contractions that were more recent than the glacial oscillations [20], an interpretation that is supported by SDMs [62]. Nevertheless, the difference in differentiation time for S. florulenta compared to those of our study species agrees with the observation that, although recurrent patterns exist, some species may exhibit distinct patterns reflecting the unique, rather than the shared, aspects of species' histories [71,72].
In general, our results for S. cordifolia and V. argenteria point out the importance of not just the LGM but also previous glacial cycles in shaping the contemporary genetic structure of plant species in the Maritime and Ligurian Alps. The results thus corroborate the notes inferring the necessity to go beyond interpretations that relate only to the role of a small number of major historical events, such as the LGM and the Messinian Salinity Crisis [7].

The distinctiveness of current populations
In general, current patterns of genetic variation are explained by contraction/expansion model in association with climatic oscillations. In this model refuge populations are associated with multiple genetic lineages and higher genetic diversity compared to more recently established populations [3]. Due to the recurrent founder effect affecting colonizing populations, the genetic diversity is supposed to decrease by increasing the distance from refugia, producing a genetic signature [73]. Nevertheless, this pattern can be obscured by secondary contact of vicariant lineages which results in admixing populations with high heterozygosity and genetic diversity [74].
The total gene diversity and haplotype diversity (  [77]). Although both species show relatively high genetic differentiation among populations, the within-population values of haplotype and nucleotide diversity are lower and the number of private haplotypes higher in V. argenteria than in S. cordifolia (Table 3) Table 4). Low within population diversity we observed suggests bottleneck or founder effects during glacial population contraction and post-glacial expansions, resulting in few haplotype being fixed in most populations. This result, together with the weak phylogeographical structure, also suggests that a repeated pattern of contraction/expansion, throughout the glacial cycles, may have shaped the contemporary genetic structure of the studied species. This result is in line with our SDM and molecular dating findings. The level of genetic differentiation detected among populations of both species even indicates low levels of recurrent gene flow among populations, as also confirmed by IBD results. The high haplotype and gene diversity, taken together with the patchy distribution and the numerous private haplotype, may be explained by a long history of differentiation and in situ glacial survival. Two factors may contribute to the higher cpDNA-based population subdivision in V. argenteria compared to S. cordifolia. First, vicariance events due to glaciation cycles may have had a stronger impact on V. argenteria, growing at higher altitude. Second, a more efficient seed dispersal mechanism in S. cordifolia may contribute to the lack of correspondence between population groups and haplotype clades in this species. This last explanation seems unlikely because V. argenteria may be dispersed by bird (see discussion below). However, the IBD result suggests a similar immigration by seed in Maritime and Ligurian Alps populations of both species (S. cordifolia: r = 0.417, P < 0.0001; V. argenteria: r = 0.586, P < 0.0001). In addition, the high proportion of variation residing among the population groups identified by SAMOVA indicates that past climatic oscillations would have stimulated a stronger intraspecific diversification and population genetic differentiation in V. argenteria. According to this result, the higher within-population diversity detected in S. cordifolia may be due to its larger altitudinal range. This may have facilitated short-distance altitudinal migrations likely decreasing the impact of founder effects and the genetic signature of glaciations. A similar pattern of short-distance altitudinal shift and recurrent re-colonizations of neighbouring populations was reported in Armeria genus growing at high altitude [78]. This is congruent with the higher haplotype diversity observed in central populations for S. cordifolia and in southern populations for V. argenteria, and with the stronger effect of glaciation on the latest. V. argenteria growing only in the alpine belt probably survived glaciation in a southern peripheral refugium while S. cordifolia may have survived to last glaciation in several refugia at lower altitudes, showing so its higher diversity in the contact zone.
The recent origin of haplotypes detected in Corsica suggests that the disjoint Corsican populations of V. argenteria may have originated by a rare long-distance dispersal of seed by birds rather than by vicariance due to Corsica-Sardinia microplate disjunction (15 Ma). In fact, seeds of violets have been shown to be capable of bird-mediated endozoochorous dispersal [79]. Moreover, Corsica is located along one of the main present-day migration routes of birds from European continental regions to the African areas [80]. This random dispersion is congruent with the weak IBD detected among all populations of V. argenteria (r = 0.296, P < 0.0001).
The complex phylogeographical pattern observed in the studied species would indicates that, even if some species-specific differences are detected, both population and genotype extinction was minimal during the LGM, probably due to the low impact of glaciations and to topographic complexity in this area [81]. This result is in line with evidences from SDM and divergence time among haplotypes. In the south western Alps, survival in peripheral refuges and nunataks during the LGM has recently been demonstrated in Primula spp. using both SDM and phylogeographical approaches [12]; however, our results are congruent with previous findings in the Maritime and Ligurian Alps where glaciations seem to have had a low influence on plant distribution and their effect seems to be weakened by high level of postglacial migrations [81]. These findings are more similar to that expected for a Mediterranean mountain region [7,8] rather than for an Alpine region [16,58], where patterns were simplified to a great extent by major losses of diversity during glacial periods.
Supporting Information S1