ARTICLE https://doi.org/10.1038/s41467-022-30394-5 OPEN Tempo and drivers of plant diversification in the European mountain system Jan Smyčka 1,2,3 ✉, Cristina Roquet 1,4, Martí Boleda1, Adriana Alberti5,6, Frédéric Boyer 1, Rolland Douzet7, Christophe Perrier 7, Maxime Rome7, Jean-Gabriel Valay7, France Denoeud 5, Kristýna Šemberová 3,8, Niklaus E. Zimmermann 9, Wilfried Thuiller1, Patrick Wincker 5, Inger G. Alsos 10, Eric Coissac1, the PhyloAlps consortium* & Sébastien Lavergne1 1234567890():,; There is still limited consensus on the evolutionary history of species-rich temperate alpine floras due to a lack of comparable and high-quality phylogenetic data covering multiple plant lineages. Here we reconstructed when and how European alpine plant lineages diversified, i.e., the tempo and drivers of speciation events. We performed full-plastome phylogenomics and used multi-clade comparative models applied to six representative angiosperm lineages that have diversified in European mountains (212 sampled species, 251 ingroup species total). Diversification rates remained surprisingly steady for most clades, even during the Pleisto- cene, with speciation events being mostly driven by geographic divergence and bedrock shifts. Interestingly, we inferred asymmetrical historical migration rates from siliceous to calcareous bedrocks, and from higher to lower elevations, likely due to repeated shrinkage and expansion of high elevation habitats during the Pleistocene. This may have buffered climate-related extinctions, but prevented speciation along elevation gradients as often documented for tropical alpine floras. 1 Univ. Grenoble Alpes, Univ. Savoie Mont Blanc, CNRS, LECA, FR-38000 Grenoble, France. 2 Center for Theoretical Study, Charles University and the Academy of Sciences of the Czech Republic, CZ-11000 Prague, Czech Republic. 3 Department of Botany, Faculty of Science, Charles University, CZ-12801 Prague, Czech Republic. 4 Systematics and Evolution of Vascular Plants (UAB) – Associated Unit to CSIC, Departament de Biologia Animal, Biologia Vegetal i Ecologia, Facultat de Biociències, Universitat Autònoma de Barcelona, ES-08193 Bellaterra, Spain. 5 Génomique Métabolique, Genoscope, Institut François Jacob, CEA, CNRS, Université Evry, Université Paris-Saclay, FR-91057 Evry, France. 6 Université Paris-Saclay, CEA, CNRS, Institute for Integrative Biology of the Cell (I2BC), FR-91190 Gif-sur-Yvette, France. 7 CNRS, Lautaret, Jardin du Lautaret, Université Grenoble Alpes, FR-38000 Grenoble, France. 8 Czech Academy of Sciences, Institute of Botany, CZ-25243 Průhonice, Czech Republic. 9 Swiss Federal Research Institute WSL, CH-8903 Birmensdorf, Switzerland. 10 UiT – The Arctic University of Norway, The Arctic University Museum of Norway, N-9037 Tromsø, Norway. *A list of authors and their affiliations appears at the end of the paper. ✉email: smyckaj@gmail.com NATURE COMMUNICATIONS | (2022)13:2750 | https://doi.org/10.1038/s41467-022-30394-5 | www.nature.com/naturecommunications 1 ARTICLE NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-022-30394-5 M ountain regions across the world are important biodi- The unique flora of alpine environments was likely assembled versity hotspots, owing to their high species richness, through a complex interaction between spatial and ecological endemism and faster pace of species diversification drivers that have influenced the divergence and migration of compared to lowlands1–4, with tropical mountain regions har- mountain plant lineages, but the relative influence of these drivers boring by far the highest concentrations of plant diversity1,5. The requires further investigation. It has long been considered that diversity of temperate mountain floras is less prominent at the allopatric divergence was the main speciation driver in European global scale, but still remarkably high compared to surrounding mountains11,15,36,37 However, this view is challenged by the evi- lowlands, specifically when considering the stressful climatic and dence that genetic structure and local adaptation across different edaphic conditions prevailing at high elevations and the dramatic elevation belts and bedrocks is pervasive in European mountain climate and glacier oscillations of the Pleistocene period6–8. The plants38,39, and that parapatric speciation along elevational gra- exceptional plant diversity and endemism of the European dient seems to be common in tropical and subtropical mountain mountains was early recognized, and its emergence has chal- floras1,5,40–44. In addition, the Pleistocene climatic oscillations lenged the understanding of botanists and biogeographers since may have unevenly impacted the rates of migration, speciation the pioneering work of von Haller 250 years ago9,10. Yet, how and and extinction of plant lineages across the major ecological gra- when plant lineages have diversified in temperate mountains still dients (e.g. bedrocks, elevation) and between mountain ranges, remains insufficiently documented. While the drivers (how) and depending on whether the particular areas were heavily glaciated the tempo (when) of species diversification have already been or remained ice-free during most of the Pleistocene—e.g. high explored for several European mountain plant lineages (e.g.11–13), elevation, siliceous areas of central Alps vs. low and mid-eleva- no consensus has yet emerged on the evolutionary history of the tion, calcareous areas of peripheral Alps, respectively45. It thus species-rich European alpine flora. This is mostly due to previous remains unclear whether species diversification in European clade-specific studies using both sparse and distinct genetic data, mountains was dominantly driven by geographic divergence, or diverse analytical methods, and as a result, sometimes even whether it was significantly complemented by speciation events yielding non-conclusive results11–15. To allow generalizations occurring along ecological gradients as known from tropical about the tempo and drivers of alpine plant diversification, we mountain systems (e.g.5,41–44). Moreover, it is not clear how these thus need to collect high-quality genomic data covering multiple evolutionary mechanisms interacted with glaciation dynamics, plant lineages, estimate reliable and comparable phylogenies, and particularly how climatic oscillations influenced diversification apply a multi-clade framework of modern phylogenetic com- and migration in the regions and parts of ecological gradients parative methods across study lineages. heavily impacted by glaciation. Disentangling these effects The tempo of species diversification within mountain biotas is requires simultaneously estimating past rates of migration and often described as a continuous, typically fast, process following cladogenesis within and between mountain ranges, elevation belts the uplift of mountain ranges16–18, which slowed down after and bedrock types, which has never been performed so far on any saturation of the whole available physical and ecological space temperate mountain flora. (e.g.19), eventually being modulated by climatic changes. A In our study, we aim at providing a window into the evolu- dominant climatic process that affected current mountain biotas tionary history of European mountain plants by inferring the worldwide11,17,18,20 was the gradual cooling of global temperatures tempo and drivers of speciation of six study plant lineages con- since the middle Miocene (15 Ma BP), that culminated with the sidered as representative cases of in situ diversification within onset of Pleistocene (2.6 Ma BP), a geological epoch characterized European mountains (Fig. 1). These clades are: Androsace sect. with low average temperatures and cyclical fluctuation of global Aretia (hereafter Androsace), Campanula sect. Heterophylla temperatures21. The effect of Pleistocene could have been positive (hereafter Campanula), Gentiana sections Gentiana, Ciminalis in terms of net diversification rates with speciation being stimu- and Calanthianae (hereafter Gentiana), Phyteuma, Primula sect. lated by climate-induced range dynamics (speciation pump Auriculata (hereafter Primula) and Saxifraga sect. Saxifraga mechanism,22,23), or negative with decreased speciation or (hereafter Saxifraga). We reconstruct their phylogenies using a increased extinction rates due to the impact of Pleistocene glacial full-plastome phylogenomic dataset obtained by low coverage periods. The latter perspective is motivated by the observation of shotgun sequencing. Using a likelihood-based multi-clade com- greater plant endemism in areas that have experienced relative parative model approach developed specifically for this purpose, climatic stability or less glaciation7,24–26. Contrary to subtropical we estimate past rates and tempo of species diversification across and tropical mountains such as the Northern Andes where mas- all study lineages, and also lineage-specific deviations from gen- sive plant diversification occurred during the Pleistocene eral trends. We show that overall, there is a weak or no effect of (e.g.27,28), it has long been considered that temperate mountains past climate cooling on diversification rates in our study lineages. such as the European Alps were so severely glaciated during Further, we use state-dependent diversification models to show Pleistocene cold periods that this posed limits to recent plant that the evolutionary assembly of these plant lineages across diversification15,29. However, the evidence for refugia at the per- different bedrocks and elevational belts was strongly driven by ipheries of glaciated European mountain ranges7,8,30,31 and also migration. Finally, we examine the importance of different modes on ice-free mountain tops protruding from glaciers (so-called of speciation based on sister species comparisons, and demon- nunataks32–35) suggests that unglaciated mountain habitats per- strate that the dominant mode of speciation was allopatry, sisted during glacial periods and may even have triggered plant complemented with bedrock-driven speciation. Our results sug- speciation. Whether Pleistocene climatic oscillations slowed down gest that past migrations across ecological gradients might have or spurred the diversification of European mountain plants, and buffered climate-related extinctions on one hand, but prevented what was the role of pre-Pleistocene gradual cooling, still remains ecological speciations on the other hand. unclear, with most recent reliable analyses concerning Primula- ceae only11. Documenting the tempo of alpine plant diversification in European mountains therefore requires estimating accurate Results divergence times for a set of distinct plant clades. This would allow Tempo of species diversification. The six study lineages started to carefully assess how diversification rates have varied through to diversify in Europe at variable dates within the last 40 Ma time and also between regions and environments that have been (Fig. 2a). The credibility intervals of crown age estimate largely differently impacted by Pleistocene climate and glacial oscillations. overlap between study lineages, but Saxifraga (95% highest 2 NATURE COMMUNICATIONS | (2022)13:2750 | https://doi.org/10.1038/s41467-022-30394-5 | www.nature.com/naturecommunications NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-022-30394-5 ARTICLE parameters and compared their combined fit to the multi-clade model. Our multi-clade approach is based on calculating joint likelihood function as a product of individual likelihood functions of models fitted on each of the six lineages (see Methods for details). In order to validate the performance of multi-clade approach and to address general issues with identifiability in diversification models46, we explored the models behavior with simulated data, showing that our approach can perform unbiased model selection and accurately estimate model parameters (Supplementary Methods 2). The best performing multi-clade model assumed constant speciation and extinction rates (Table 1). Time- or temperature- dependent models (median AICdiff between −1.16 and −1.71 [df=1] across 100 sets of Bayesian posterior trees) cannot in principle be rejected with confidence, but our sensitivity analysis showed that the used modeling approach and dataset have sufficient statistical power to detect temperature-dependent scenarios where Quaternary speciation rate dropped to 63% of Tertiary rate or less (Supplemen- tary Methods 3). The universal temperature-dependence of specia- tion across all clades in our dataset considered together was thus either absent or weaker than such drop. Importantly, the combination of lineage-specific models slightly outperformed the models with shared parameters (median AICdiff = 9.68 [df=10]), suggesting that the estimated diversification parameters in reality differed across lineages (Fig. 2b). In particular, the diversification dynamics of two lineages was better explained by models with non- constant rates (Table 1): Primula showed a slowdown of speciation either with time (median AICdiff = 2.04 [df=1]) or during colder periods (median AICdiff = 2.09 [df=1]); and Androsace showed support for a speciation slowdown in colder periods, although model performance was only slightly better than the constant rate model (median AICdiff = 0.40 [df=1]) in this clade. Evolutionary assembly: the relative influence of bedrock, ele- Fig. 1 Map of the European mountain system depicting the five major vation and geography. We used cladogenetic state-dependent geographic regions and the six mountain plant lineages that notably diversification models (ClaSSE,47,48) to analyze the evolutionary diversified in European mountains. 84% of all species from these lineages assembly of our study lineages across bedrock types (calcareous was sampled in this study and we sampled across all five mountain regions. vs. siliceous), elevation belts (high elevation vs. mid-elevation The colors used for the six different phylogenies are used accordingly habitats), and geographic regions (five major European mountain throughout following figures. The timescale unit of phylogenetic branch regions, see Fig. 1). Using AIC comparisons, we quantified the lengths is Ma. The background map is adapted from maps-for-free.com. importance of speciation associated with splits between bedrock types, elevation belts, or regions (which we term state-change probability density of crown age 12.5–38.1 Ma BP), Phyteuma speciation); speciation within the same bedrock type, elevation (13.7–24.5 Ma BP), Gentiana (6.4–28.5 Ma BP) and Campanula belt or region (constant-state speciation); and the importance of (10.1–16.5 Ma BP) were probably older than both Primulaceae change of bedrock type, elevation belt or region occurring along clades, namely Androsace (4.7–10.7 Ma BP) and Primula branches of the phylogeny, that is without been linked to spe- (3.2–7.9 Ma BP). All lineages accumulated considerable amounts ciation events (which we term migration). The ClaSSE models for of diversity both before and during the Pleistocene (2.6 Ma BP) so bedrocks and elevation belts were equivalent to a GeoSSE that none of them results from exclusively Tertiary nor exclu- model47, while the model for geographic regions represents a sively Quaternary diversification. Our estimates integrate sources generalization of GeoSSE for more than two regions (see Methods of uncertainty stemming from both molecular and fossil data in a and Supplementary Methods 4 for details). The best models conservative way and seem to be robust to alternative inter- inferred from model selection were then re-run in a Bayesian pretations of fossil record—for details about dating analyses and framework to obtain credibility intervals around parameter esti- handling of various sources of uncertainty, see Methods and mates. As in the previous analyses, all models were evaluated in Supplementary Methods 1. the multi-clade framework with a model sharing parameters To better understand the dynamics of clade diversification, we between clades, and subsequently as lineage-specific models with tested a range of hypotheses about diversification changes distinct parameter set for each clade. To validate that parameters through time by contrasting five models depicting different of multi-clade and lineage-specific models are identifiable, we diversification scenarios (Table 1). The model support was tested these approaches with simulated phylogenies reflecting the evaluated by difference in AIC between the focal model and a size and structure of our dataset (Supplementary Methods 4). nested null model not containing the focal parameter(s) (referred The state-change speciation between siliceous and calcareous to as AICdiff throughout the paper, see Methods for details). We bedrocks appeared to be an important driver of plant diversifica- first ran all the models in a multi-clade setup with parameters tion (AICdiff = 8.96 [df=1]). Model parameters suggest that a shared across the six lineages to seek for general patterns. We bedrock generalist lineage split into descendant lineages specializ- then separately re-ran the five models with lineage-specific ing on calcareous or siliceous habitats on average cca 0.5 times NATURE COMMUNICATIONS | (2022)13:2750 | https://doi.org/10.1038/s41467-022-30394-5 | www.nature.com/naturecommunications 3 ARTICLE NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-022-30394-5 Fig. 2 Tempo of species diversification for the six study lineages. a Lineage-through-time curves of each lineage. The thick lines represent maximum credibility phylogenetic reconstructions, while the semi-transparent lines represent 100 trees sampled from Bayesian posterior distributions. The dashed line marks the onset of the Pleistocene at 2.6 Ma before present. The number of lineages is plotted in logarithmic scale, i.e. the exponential growth expected under pure birth model would appear linear here. b Parameter estimates indicating the effect size of the temperature control on speciation rates in temperature-dependent models of species diversification (Δ speciation rate). The black dots represent the estimates from the multi-clade model with shared parameters among lineages. Higher values above the dashed line indicate lower speciation in colder geological periods. Thick dots correspond to estimates based on the maximum credibility phylogenetic trees, and semi-transparent dots on the 100 trees sampled from the posterior distribution. Source data are provided as a Source Data file. per Ma and such speciation events constituted 19% of all directional shifts along ecological gradients, as our analyses speciation events (Table 2). However, we found no evidence for inferred strong net directional migration between bedrock types state-change speciation between elevation belts (AICdiff = and elevation belts. The rate of migration from siliceous to −2 [df=1]) or between the five major mountain regions of Europe calcareous bedrock habitats was higher than in the opposite (AICdiff = −1.62 [df=1]). The latter result suggesting no allopatric direction (AICdiff = 5.59 [df=1], Fig. 3c), which better explains speciation between mountain regions is however likely biased due comparable present-day proportions of species growing on to the coarse scale of the considered mountain regions caused by silicate and calcareous bedrocks (Fig. 3e, Supplementary Fig. 2) computational limitations (see Methods for details). For this despite differing speciation rates. We found even more important reason we addressed the prevalence of allopatric speciation and its asymmetric migration rates between elevation belts, as we spatial scaling with complementary analysis of sister species (see inferred strong directional migration from high to mid- below). elevation habitats (AICdiff = 19.79 [df=1], Fig. 3d). This result Constant-state speciation rates differed between bedrock types captures events of secondary migration of high altitude ancestors (AICdiff = 3.51 [df=1]) and regions (AICdiff = 17.87 [df=4]), but toward lower elevations, leading to relatively high present-day not between elevation belts (AICdiff = −1.02 [df=1]). In parti- and equilibrium proportions of species inhabiting either mid- cular, we inferred higher speciation rates on siliceous than on elevations or occurring in both elevation belts (Fig. 3f, Supple- calcareous bedrock (mean estimate and credibility interval in mentary Fig. 2). Fig. 3a) and also higher speciation rates in the Alps and the Models of evolutionary assembly slightly differed among the Iberian mountains than in any other European mountain region six study lineages (AICdiff = 7.39 [df=30] for bedrock, (Supplementary Fig. 1), but no notable difference in speciation AICdiff = 21.07 [df=20] for elevation). The most notable outliers rates between high- and mid- elevation belt (Fig. 3b). Unlike concerning the assembly across bedrock types were Androsace speciation, the inferred extinction rates did not vary between and Phyteuma, which exhibit higher proportions of species bedrock types, elevation belts or regions (AICdiff < 0 in all the occurring on siliceous than on calcareous bedrock (Fig. 3e). In cases). Interestingly, the differences in diversification rates were Androsace, this was linked to marginally faster diversification not the major force affecting proportions of species across rates on siliceous than on calcareous bedrock and symmetrical bedrock types or elevation belts. In other words, higher speciation migration between the two bedrock types. Phyteuma showed rate in certain habitats did not necessarily result in higher species marginally faster diversification on calcareous bedrock, but strong richness in these habitats (Figs. 3e, f). Rather, the difference in directional migration towards siliceous bedrock, thus explaining contemporary species richness seems to originate from the higher species diversity of this genus on siliceous bedrock 4 NATURE COMMUNICATIONS | (2022)13:2750 | https://doi.org/10.1038/s41467-022-30394-5 | www.nature.com/naturecommunications NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-022-30394-5 ARTICLE Table 1 Median AIC differences (AICdiff) between multiple time- or temperature-dependent diversification models and the nested null model with constant speciation and constant extinction. Androsace Campanula Gentiana Phyteuma Primula Saxifraga shared Speciation constant, extinction constant 0.00 0.00 0.00 0.00 0.00 0.00 0.00 Speciation dependent on past temperature, extinction constant 0.40 −0.99 −1.42 -1.82 2.09 −1.70 −1.71 Speciation constant, extinction dependent on past temperature −2.00 −1.82 −0.25 −2.00 −2.00 −1.83 −1.16 Speciation exponentially dependent on time, extinction constant 0.00 −0.95 −0.45 −1.89 2.04 −1.81 −1.16 Speciation constant, extinction exponentially dependent on time −0.08 −1.53 −0.31 −2.00 1.62 −1.80 −1.25 Values with AICdiff > 2 are depicted in bold. Table 2 Comparison of ClaSSE models of evolutionary assembly containing parameters of state-change speciation between bedrocks (siliceous vs. calcareous), elevation belts (high elevation vs. mid-elevation) or five coarse geographic regions, with the nested models lacking state-change speciation terms. Rate of state-change speciation (Ma−1) Proportion of state-change to all speciation events AICdiff Bedrock types 0.453 0.19 8.96 Elevation belts <0.001 <0.001 −2.00 Geographic regions (coarse scale) 0.013 0.063 −1.62 Values with AICdiff > 2 are depicted in bold. (Fig. 3a, c). A notable exception concerning assembly across Discussion elevation belts was Androsace, which has a higher proportion of Here we explored the tempo and drivers of species diversification species occurring at high elevations (Fig. 3f) due to directional across a representative sample of six diverse plant lineages whose migration from mid- to high elevations (Fig. 3d). Based on evolutionary history is tied to the European mountain system. As ancestral state reconstructions (Supplementary Fig. 3b), Andro- we purposely narrowed our study down to clades that previously sace also contains two exceptional lineages (containing five received a systematic revision, it is clear that many of our studied species each) which most likely spent their evolutionary history species had already been sequenced, although always for a few exclusively in high elevation habitats. Gentiana also has a genes only and using variable methodologies. Instead, we relatively high proportion of species in high elevation habitats, employed a uniform next-generation sequencing approach to all but here the very broad credibility intervals of model parameters study clades and produced over two hundred whole plastid did not allow us to further distinguish scenarios of evolutionary genomes. Chloroplast genomes provide a common ground assembly. molecular information that can be applied across various groups and phylogenetic scales. While the resulting topologies may be prone to biases introduced by incomplete lineage sorting or Sister-species overlap in geographical and ecological space. Due hybridization events in some cases (see Methods for more to computational limitations, the ClaSSE models used above information and a brief discussion28), our chloroplast-based could only be run with a maximum of five geographic regions. approach avoids some of the pitfalls of nuclear phylogenomics Hence, we performed a complementary analysis to infer how (e.g. paralogy49). Our study therefore constitutes a significant step sister species vary in the dimension of overlap in terms of ele- forward toward documenting the evolutionary origins of alpine vation belt, bedrock type and geographic range, the latter being plant diversity in European mountains, because it is the first to estimated both with coarse-scale and fine-scale geographic harness next-generation sequencing on multiple plant clades in regions (five and 87 mountain regions, see Fig. 1 and Supple- parallel and to apply a multi-clade comparative framework to mentary Methods 5). Most sister pairs diverged over the fine explicitly model the tempo and drivers of mountain plant geographic scale (Fig. 4), with 37% sister species pairs showing no diversification. To make sure that our data collection strategy and geographic overlap at this scale. Similarly, as in ClaSSE models, analytical approaches can inform us about patterns of mountain sister species differed less across bedrock types (27% of sister plant evolution, we also extensively validated the identifiability pairs with no overlap) and over the coarse geographic scale (22% and sensitivity of our comparative models using simulated data- of sister pairs with no overlap). Only 9% of sister species pairs sets. Our results may thus allow certain generalizations across our showed no overlap in the occupied elevation belts. study lineages, confirming previous knowledge or botanists Although the order of speciation modes remained the same intuitions regarding plant diversification in European mountains, among all six study lineages, the degrees of overlap slightly documenting some quite surprising, and perhaps counter- differed between lineages (Fig. 4). Notably, the lineages with intuitive patterns, and finally identifying some interesting varia- lowest geographic overlap at the fine scale were Primula (63% of tions between the evolutionary histories of different study clades. pairs with no overlap) and Androsace (43% of pairs with no The onset of cold and variable Pleistocene climate has caused overlap), suggesting that allopatric speciation at finer geographic massive losses of plant diversity in European lowlands50 and was scale was particularly important for these two lineages. The long thought to have induced a diversification slowdown in lineage Campanula exhibited high variance in sister species European mountains29, which was also demonstrated with par- overlap distribution at fine geographic scale, on the one hand ticular plant lineages11,15. Surprisingly, we detected neither strong having a high proportion of sister pairs with no overlap (43%), nor homogeneous effects of past temperature on plant diversifi- but on the other hand having the highest proportion of sister cation rates across our six study lineages. Diversification species pairs showing complete geographic overlap (36%). dynamics in all study clades together were best explained by a NATURE COMMUNICATIONS | (2022)13:2750 | https://doi.org/10.1038/s41467-022-30394-5 | www.nature.com/naturecommunications 5 ARTICLE NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-022-30394-5 Fig. 3 Evolutionary assembly across bedrocks and elevational belts based on ClaSSE models. a, c The differences between siliceous and calcareous habitats in constant-state speciation and migration rates, respectively. b, d The differences between high and mid-elevation habitats in constant-state speciation and migration rates, respectively. The dots represent mean parameter estimates and the bars indicate 95% credibility intervals, based on 5000 MCMC samples from the ClaSSE model posterior. Black dots and bars represent the shared parameter estimates from the multi-clade model, where each of the 6 phylogenies is regarded as an independent realization of the same diversification process. e The proportions of species inhabiting siliceous, calcareous or both types of habitats (light, dark and middle gray, respectively), and f the proportion of species inhabiting high elevation, mid-elevation or both types of habitats (light, dark and middle gray, respectively). Source data are provided as a Source Data file. constant rates model, and showed no indication for overall dramatically colder Pleistocene climate. This result was slowdown of species diversification with decreasing global tem- further corroborated by a sensitivity analysis, suggesting that if perature. When analyzed at the level of individual lineages, only any temperature-driven Pleistocene reduction of diversification the two Primulaceae clades showed certain support for diversi- rates took place across all study lineages, it was most likely fication slowdown, as was also demonstrated in earlier weaker than a drop to 63% of pre-Pleistocene rates, given the studies11,15. Yet more interestingly, the other four lineages observed data. showed no individual sign of slowdown with decrease of tem- The weak and ambiguous effect of Pleistocene climate on perature. The past temperature thus did not universally diversification rates was also detectable in analyses across eleva- influence the diversification rates of our study lineages, which was tional gradient and geographical space: the high elevation habitats not only true for the gradual, and perhaps hardly detectable, the and also the Alps sensu stricto, which were most severely effect of cooling since middle Miocene, but also for the onset of impacted by Pleistocene glaciers, did not show lower rates of 6 NATURE COMMUNICATIONS | (2022)13:2750 | https://doi.org/10.1038/s41467-022-30394-5 | www.nature.com/naturecommunications NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-022-30394-5 ARTICLE European mountain biotas. In fact, it is interesting to note the similarity between the scenario of lineage migrations suggested by our results and the one reported by 58 on Mediterranean mountain birds. Why did diversification of some lineage slowdown while other remain steady? Given that we have six cases for comparison and we only observed a slowdown in the Primulaceae family (Primula and Androsace), we can only hypothesize the answer. A scenario of ecological opportunity driving species diversification following mountain uplift17,27 could provoke a diversification slowdown after saturation of the available ecological space. But such an effect should be observed in the oldest and most diverse lineages of our study clades, rather than in the youngest and relatively least diverse ones such as Primula and Androsace. The observed diversification slowdown could alternatively be explained by a decreasing rate of allopatric speciation, as the two lineages that slowed down towards the present (Primula, Androsace) also show the highest rates of sister species allopatry. The prevalence of allopatric speciation can theoretically lead to an intrinsic slow- down of speciation rates, particularly in lineages with low dis- persion capacity, once the spatial setting of species’ geographic ranges no longer enhances further genetic isolation59. Another explanation could be that Pleistocene climatic oscillations pro- moted migrations between populations in areas that would otherwise remain isolated, thus reconnecting populations and inhibiting allopatric diversification. Simulations of diversification across dynamically fragmented landscapes indeed showed that more connected landscapes or a faster pace of connection- disconnection events can impede species diversification under certain conditions23,60. Allopatric speciation has long been regarded as a prevalent mode of speciation in mountain environments11,15,36,37,61, and the high importance of allopatric speciation demonstrated by our results is generally consistent with this view. An additional interesting finding of our study is the relatively small spatial scale of allopatric speciation, which did not take place between major mountain ranges such as the Alps or Carpathians, but rather Fig. 4 Sister-species overlap in geographic and ecological space. Violin within these ranges. It would be interesting to further study the plots show distribution of Schoener’s D index of overlap of sister species same lineages with nuclear markers, as spatial structure may pairs across 87 operational geographic units (fine geo), five European sometimes emerge due to the usage of chloroplast-based maternal mountain regions (coarse geo), bedrock types (calcareous vs. siliceous) lineage phylogenies. (e.g.28). Although the small-scale allopatric and elevational belts (high vs. mid-elevations). The point within each plot speciation appeared to be the most important mode of speciation represents median overlap value. Low overlap between sister pairs is (37% of speciation events) in sister species of our study lineages, considered as an indication of frequent speciation along the respective our results further suggest that bedrock-driven divergence has geographic or ecological dimension. Source data are provided as a Source been another important driver (27%) of sister-species speciations. Data file. This contrasts with a previous study11 finding an almost exclusive mode of allopatric speciation and little bedrock divergence in speciation or higher rates of extinction. Instead, high elevation three European alpine Primulaceae genera. These previous lineages exhibited high rates of directional migration to lower results11 may however reflect a peculiar situation in Primulaceae - elevations. This may suggest that the habitats suitable for high a family containing Androsace and Primula - the two lineages in mountain plants were not severely reduced during glacial periods our dataset for which we found allopatry to be by far the most but rather shifted downwards. Such shifts instead of loss in important mode of speciation. habitats may thus have resulted in relatively weaker and more We estimated that bedrock shifts account for about one fourth lineage-specific declines in diversification rates. Strictly speaking, (27% in sister species analyses) or one fifth (19% in ClaSSE our findings are based on lineages that significantly diversified in analyses) of all speciation events across our study clades. Bedrock the European mountain system, and thus do not contradict the or edaphic species differentiation is well known from medi- paleontological evidence for massive Pleistocene extinctions in terranean floras62–65, while in European mountain species it has the European lowland flora51–53. The high rates of endemism and mostly been studied at the intraspecific level39, or in studies using genetic diversity in low elevation refugia at mountain non-quantitative comparisons to other speciation drivers66. To peripheries7,8,54–56, and the presence of mountain plants in our knowledge, our result thus provides the first quantitative lowland glacial palynological record (Androsace sp., Gentiana sp., evaluation of the importance of bedrock-driven speciation in Saxifraga sp., but also Dryas octopetala, Polygonum viviparum or European mountain plants. Still, it remains unclear to what extent Saussurea sp.57) however suggest that the above described these speciation events were truly parapatric, that is with reduced migration, diversification and survival processes may represent a gene flow between adjacent populations due to bedrock adapta- general model for other extant mountain plant lineages not tion, or whether bedrock specialization originated as a con- included in our study, and perhaps also for other components of sequence of allopatry. Experimental evidence suggests that new NATURE COMMUNICATIONS | (2022)13:2750 | https://doi.org/10.1038/s41467-022-30394-5 | www.nature.com/naturecommunications 7 ARTICLE NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-022-30394-5 species can indeed arise due to strong divergent selection of may have differentiated at intraspecific levels38,75, but the strong calcareous and siliceous ecotypes, even in the presence of gene migration from high to lower elevations in response to the severe flow67. Our findings have far reaching implications: it shows that Pleistocene climate oscillations has likely prevented the emer- contrary to long-standing expectations, many speciation events gence of reproductive barriers necessary for speciation events. may not or not only have been due to selectively neutral geo- Important historical migrations along elevation gradients may graphic divergence, but also due to divergent evolution across a thus provide an explanation why species diversification has not complex geological landscape requiring specific physiological been as explosive in the European mountain system as in tropical adaptations to the divergent proportions of important nutritional ones15,29. elements68. The lineage Androsace sect. Aretia constitutes a notable Another important finding of our study is the inference of exception regarding the general scenario of evolutionary assembly frequent migration events between bedrock types. Intriguingly, highlighted above, probably because this lineage shows strong we found that these migrations were asymmetrical, with silicate affinity to high elevation habitats. In comparison to the others, habitats being more often the source of limestone-dwelling this lineage shows a greater species richness in high elevations, lineages than the other way around. This finding seemingly likely due to directional migration from mid- to high elevation contrasts with previous literature suggesting that plant adaptation habitats. Moreover, according to the ancestral state reconstruc- to calcareous soils is more restrictive than adaptation to silicate tion it also contains two sub-lineages that likely evolved exclu- bedrocks67,69–71. The higher pleiotropic fitness costs of adapta- sively in high elevation habitats. In fact, several species of tion to calcareous bedrock may however prevent lineages from Androsace stretch to the highest limits of vascular plant life in the switchback to siliceous habitats, and calcareous habitats might European mountain system (above 4000 m a.s.l.,76), and in thus function as the evolutionary trap72,73. Moreover, the stron- Himalaya (6350 m a.s.l.,77). Unlike other study lineages, Andro- ger migration from silicates to limestones we detected here may sace has likely undergone specific adaptations facilitating the also be linked to geographic contingencies: the peripheral lower repeated entrance into the harsh adaptive zone of high elevation mountain systems in Europe may be both calcareous and silic- habitats and continued population persistence at high elevations eous, while central high elevation parts of European mountains throughout the Pleistocene. A clear manifestation of such adap- (such as Central Alps, Central Pyrenees or High Tatras) consist tation is that Androsace repeatedly developed a dense cushion life almost exclusively of silicate bedrocks45. Due to this triangular form, an architecture seemingly well adapted to plant life at high relationship, the elevation belts may constitute a confounding elevations not only in this genus78 but also in many other hidden factor (sensu74) in the models of evolutionary assembly angiosperm lineages79. Given the inferred long-term affinity of across bedrock types, where the inferred migrations from high to Androsace to high alpine environments, it can be reasonably low elevations and from siliceous to calcareous bedrocks may hypothesized that high-elevation sub-lineages of Androsace have partly reflect the same process of historical migration out of the survived glacial periods in situ in so-called nunatak refugia—i.e. central ranges of the mountain systems, possibly due to Pleisto- rocky outcrops at high elevations protruding the glaciers80— cene glaciations. The interaction between bedrock, elevational rather than in peripheral refugia, which may have been the pre- niche and possibly habitat type may in fact be rather complex, as valent scenario for the majority of the other study lineages. The demonstrated by Androsace and Phyteuma, which both deviate survival in the limited area of nunataks could also explain from the general trend and yielded greater migration rates Pleistocene diversification slowdown observed in Androsace, towards the silicate bedrocks. In Androsace this pattern is likely although the consequences of nunatak survival for species linked to lineages surviving in high elevation silicate rocky diversification remains unclear. The high elevation sub-lineages habitats, probably even during Pleistocene (discussed below). In of Androsace sect. Aretia identified in this study thus constitute a Phyteuma, silicate specialists are often species inhabiting sub- perfect system for further tests of the nunatak survival hypothesis alpine grasslands in lower silicate mountains with mild slopes, at smaller spatial and phylogenetic scales. such as Massif Central, Sudetes or eastern Carpathians, where In conclusion, our study provides an unprecedented window to these lineages may also have survived Pleistocene. More work the history of diversification of the temperate mountain flora in based on an extended population sampling is certainly needed to Europe. It shows that plant diversification within European further disentangle different historical scenarios of evolutionary mountains was a complex evolutionary process, with a strong assembly of alpine floras. interplay between altitudinal migration, allopatric speciation and The evolutionary assembly of our study clades across elevation bedrock adaptation of different lineages. Importantly, the onset of belts was apparently marked by massive directional migration Pleistocene climate did not cause a strong diversification slow- from high to mid-elevations, causing the mid-elevation habitats down, as was previously expected, but rather stimulated strong to host more species despite speciation rates being roughly migrations of mountain biota across elevation gradients, parti- equivalent in both elevation belts. This directional migration was cularly towards lower altitudes. We hypothesize that these mas- likely fostered by Pleistocene glaciations causing downward shifts sive altitudinal migration events on the one hand buffered of communities during glacial periods, which then may have left extinctions due to habitat loss during glacial periods, but on the relictual populations at lower elevations during interglacials. Such other hand impeded Pleistocene adaptive radiations across the climate-induced range shifts may have had two-fold con- elevation gradient, as is classically known from tropical moun- sequences on the dynamics of species diversification. On the one tains. We also found speciation events to be mostly driven by hand, it may have prevented species extinctions due to habitat geographic divergence but almost as frequently by bedrock shifts. loss during glacial periods, as discussed earlier. On the other Overall, the absence of obvious adaptive diversification across the hand, intense altitudinal migration may have hampered specia- elevation gradient and the prevalence of allopatric speciation tion across elevation belts, as we did not find evidence for spe- likely contributed to the lower richness and slower diversification ciation along elevational gradients. This finding makes an dynamics generally observed in temperate mountains compared important distinction with known mechanisms of species diver- to tropical ones. sification in the Andes or other tropical mountains, where spe- ciation across elevation gradients seems to have been a relatively Methods important driver of plant speciation5,41–44. More specifically, high Clade selection and sampling. The study lineages were selected as representative and mid-elevation populations in the European mountain system cases of plant species diversification in the European mountains by the following 8 NATURE COMMUNICATIONS | (2022)13:2750 | https://doi.org/10.1038/s41467-022-30394-5 | www.nature.com/naturecommunications NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-022-30394-5 ARTICLE criteria: they contain more than 20 species in total; they contain at least 10 species identified and extracted based on their positioning relative to the neighboring coding inhabiting alpine and nival elevational belts in the European Alps (based on ref. 81); regions. Coding and non-coding regions of ingroup and outgroup samples were they were recently subject to taxonomic or phylogenetic revision suggesting that further filtered in order to minimize missing data in our family-level alignment Europe is their center of diversity; they are eudicots; and they do not follow derived matrices (see below). The filtering resulted in 72 coding regions, among which 40 are life strategies as is myco-heterotrophy, parasitism or carnivory. Within the PhyloAlps shared across all study families, and 17 non-coding regions among which 5 are shared consortium, which is an extensive network of collaborating institutions (see Supple- across all the families, see Supplementary Data 3 for a table of regions and their usage mentary Note), we sampled 212 out of 251 known ingroup species and well- for each family. Samples with missing data for some of the regions were kept in established subspecies from these lineages (that is 84%), specifically 26 out of 29 specific cases: when they corresponded to the only representative of an ingroup known species and subspecies for Androsace (90%), 45 out of 50 for Campanula species or to an outgroup species defining the dating node and at the same time had (90%), 28 out of 35 for Gentiana (80%), all 27 for Phyteuma, all 24 for Primula, and so large portion of missing data that excluding missing regions would severely limit 62 out of 86 for Saxifraga (72%). For details on sample counts, see Supplementary the resolution of the whole family phylogeny. To make sure that inclusion/exclusion Methods 6. The well-established subspecies were treated and referred to as species had no impact on resulting topology or dating (with exception of dating-node throughout the study. The majority of ingroup (sub)species were covered by one outgroups), we reran phylogenetic reconstruction analyses without these samples and sample, but in several cases we included two or more samples for control purposes. compared the resulting trees with all-samples phylogenies. The samples with missing For all lineages, we followed the most up-to-date taxonomic treatments and took into data in the alignment matrix are indicated as such in the Supplementary Data 2. account phylogenetic studies providing group circumscription, see Supplementary For each of the four families, we aligned the coding regions gene-by-gene using Methods 6 for details. As outgroups for our analyses, we sampled an extensive MACSE90 acknowledging the triplet structure of codon alignment. All the positions collection of species from the families Campanulaceae, Primulaceae, Gentianaceae, were quality filtered by Gblocks 0.9191 (specifying that the data type corresponded Saxifragaceae and Grossulariaceae. The latter one was added to Saxifragaceae because to codon alignments), and concatenated together using FasConCat92. The non- of their close phylogenetic relationships82 and to compensate the lack of fossil cali- coding regions were aligned by Mafft 793, quality filtered by Gblocks (specifying bration points within Saxifragaceae. The majority of samples were collected in natural that the data corresponded to DNA alignments) and concatenated together. This habitats and immediately desiccated in silicagel, but in some cases we used collections sequence processing pipeline resulted in high-quality sequence matrices for each from botanical gardens or herbarium material, for details on sample identities see family ranging between 35471–47102 bp for coding and 2435–5112 bp for non- Supplementary Data 1 and 2. coding regions, with less than 3% missing data overall. DNA extraction and sequencing. DNA was extracted from collections of leaf tissues gathered during our sampling campaigns, using DNA extraction kits from Phylogenetic analyses and molecular dating calibrations. The resulting coding Macherey-Nagel (Düren, Germany) and Qiagen (Hilden, Germany). The shotgun and non-coding alignments for each of the four families were used for inferring libraries were prepared and sequenced with methodology depending on the dated phylogenies in BEAST 2.6.294. Non-coding alignments and each codon sequencing facility (see Supplementary Data 1 for sequencing facility used for each position in coding alignments were modeled with separate averaging-site models as sample): implemented in bModeltest 1.2.095. We used Yule tree prior and lognormal clock For the samples sequenced in Genoscope (Paris, France), the library preparation with uniform fossil dating priors on at least two nodes of each family. Note that we protocol applied before sequencing was chosen on the basis of the DNA extraction applied uniform priors between minimum and maximum bounds because the yield. When available, 250 ng of genomic DNA were sonicated using the E210 Covaris sparse fossil record of our study clades does not allow applying more informative instrument (Covaris, Inc., Woburn, MA, USA) and the NEBNext DNA Modules priors—this therefore constitutes a very conservative approach. Products (New England Biolabs, Ipswich, MA, USA) were used for end-repair, 3′- In some cases, the fossil record used for dating calibration allowed multiple adenylation and ligation of NextFlex DNA barcodes (Bioo Scientific Corporation, interpretations of age and positioning in the phylogenetic tree. For consistency with Austin, Texas, USA). After two consecutive 1x AMPure XP (Beckmann Coulter literature, we specifically used the interpretations mirroring previous influential Genomics, Danvers, Massachusetts, USA) clean ups, the ligated products were studies on the respective families14,16,82,86 for the downstream analyses. To test the amplified by 12 cycles PCR using Kapa Hifi Hotstart NGS library Amplification kit robustness of this choice, we however reran the phylogenetic analyses with (Kapa Biosystems, Wilmington, MA, USA), followed by 0.6x AMPure XP alternative dating calibration where necessary and explored the differences. For a purification. When the extraction yielded low DNA quantities, 10–50 ng of genomic detailed description of used fossil calibrations and their comparison with DNA were sonicated. Fragments were end-repaired, 3′-adenylated and NEXTflex alternative interpretations, see Supplementary Methods 1 and Supplementary DNA barcoded adapters were added by using NEBNext Ultra II DNA Library prep Software 1 and 2. kit for Illumina (New England Biolabs). After two consecutive 1x AMPure clean ups, The BEAST inference for each family was performed in four independent runs, the ligated products were PCR-amplified with NEBNext® Ultra II Q5 Master Mix each having 100 M generations and a burn-in of 30%. We evaluated convergence included in the kit, followed by 0.8x AMPure XP purification. All libraries were both visually and by checking whether the effective sample size (ESS) was >100 for subjected to size profile analysis conducted by Agilent 2100 Bioanalyzer (Agilent all the parameters. We used TreeAnnotator to generate maximum credibility trees. Technologies, Santa Clara, CA, USA) and qPCR quantification (MxPro, Agilent The maximum credibility trees and reduced posterior trees were pruned to only Technologies), then sequenced using 101 base-length read chemistry in a paired-end include ingroup species and one individual per species, for details on species tree flow cell on the Illumina HiSeq2000 sequencer (Illumina, San Diego, CA, USA). On inference see Supplementary Methods 6. We obtained well resolved species-level average, 5 billion useful paired-end reads were obtained. An Illumina filter was phylogenies with 87% of nodes receiving >0.95 posterior probability (see applied to remove the least reliable data from the analysis. The raw data were filtered Supplementary Software 1 for maximum credibility trees with node supports). to remove any clusters with too much intensity corresponding to bases other than the The phylogenies based on chloroplast genomes provide high resolution and called base. Adapters and primers were removed on the whole read and low-quality accurate dating due to the large number of orthologous genomic regions with nucleotides were trimmed from both ends (while the quality value is lower than 20). variable mutation rates, and the universality of the bioinformatic pipeline allowing Sequences between the second unknown nucleotide (N) and the end of the read were to obtain homogeneous phylogenetic information across virtually any angiosperm also removed. Reads shorter than 30 nucleotides after trimming were discarded. group49. The limitation of our approach is that chloroplast-based phylogenies are Finally, the reads and their mates that mapped onto run quality control sequences only tracking the evolution of maternal lineages, which might be problematic in (PhiX genome) were removed. These trimming steps were achieved using internal clades with frequent hybridization, introgression, or incomplete lineage sorting. software based on the FastX83. The shortcomings of plastome-based phylogeny might have impacted, even For the samples sequenced in Fasteris (Geneva, Switzerland) and BGI slightly, the relative branching times and topologies of reconstructed trees of our (Shenzhen, China), custom commercial protocol for preparation was applied. The lineages (see e.g.13 for Saxifraga), and hence some downstream analyses of species prepared libraries were sequenced on Illumina HiSeq2000 sequencer (Illumina, San diversification. It is important to point out that different downstream analyses we Diego, CA, USA), using 101 base-length read chemistry in a paired-end flow cell. present in our paper are differently prone to biases related to the usage of chloroplast markers: Time-dependent analyses only depend on distribution of branching events through time96, and are thus least likely to be systematically Chloroplast genome reconstruction, gene region selection and sequence influenced by the usage of maternal tree topologies. State-dependent analyses are processing. The resulting paired-end reads were used to reconstruct the chlor- phylogeny-explicit, but their overall results across multiple clades are also relatively oplast genomes using Org.Asm 1.0.384, a de novo organelle assembler based on De unlikely to be biased by the usage of chloroplast markers in a systematic manner. Bruijn graph implemented in Python 3. We retrieved complete circular plastomes Moreover, the main result of high mobility of lineages across elevational gradient except for samples belonging to Campanulaceae (Campanula and Phyteuma and bedrocks would likely be underestimated rather than overestimated by usage of lineages) where the assembly typically resulted in several discontinuous contigs that presumably less mobile maternal lineages. Perhaps the most bias-prone analyses however covered the majority of chloroplast coding regions. are sister species comparisons, because they directly depend on terminal tree We aimed to work with all chloroplast coding regions that were present in our topologies. The absence of polyphyletic species and the low number of identified chloroplast genome assembled sequences. To do so, we detected all open reading paraphyletic species (see Supplementary Methods 6 for details) roughly suggest frames in our circular and fragmented reconstructions and automatically compared that neither hybridization and introgression nor incomplete lineage sorting are very them with a curated database of annotated genes from GenBank, following the common in the terminal parts of our phylogenies. More exact prevalence of Org.Annot procedure implemented in Org.Asm84. The non-coding regions used in hybridization, introgression and incomplete lineage sorting in our study clades our study were selected on the basis of their universal phylogenetic informativeness85 should however be investigated in future phylogenomic studies using nuclear or previous use in one of the focal lineages12–14,19,86–89. The non-coding regions were markers and denser population sampling, especially if results of interest are specific NATURE COMMUNICATIONS | (2022)13:2750 | https://doi.org/10.1038/s41467-022-30394-5 | www.nature.com/naturecommunications 9 ARTICLE NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-022-30394-5 tree topologies and fine-scale processes, rather than overall diversification patterns even if there was one completely non-informative parameter added on top of the and dynamics we present in this study. focal one. We used equivalent interpretation of AICdiff values also in other AIC comparisons throughout the paper. To test whether the dynamics of diversification are indeed shared across the six Ecological and geographic characteristics of species. The information about lineages as assumed by the multi-clade models, or whether they quantitatively bedrock, elevational niche and geographical occurrences for each ingroup species differ among them, we compared the AIC of the multi-clade model with the was obtained from the following regional floristic literature: Flora Alpina97, Flora respective sum of AIC values across the six lineage-specific models. Such a Iberica98, Flóra Slovenska99, Wildpflanzen Siebenbürgens100, and Flora Srbije101. comparison is meaningful, because the sum of AIC values of the six lineage-specific We used published ecological and biogeographical studies focused on certain models is equal to the AIC calculated from a model with joint likelihood function, lineages to obtain additional information about genera Androsace102–104, but with each model parameter kept lineage-specific. Phyteuma12 and Saxifraga105. The information about Campanula was in part The time-dependent diversification models were recently criticized due to compiled also based on information available in herbarium specimens, local identifiability problems46. In our analyses we address this issue in several ways: taxonomic literature and from field observations. In several cases, we contacted First, we use diversification models corresponding to explicit hypotheses of past local taxonomists and made categorization based on provided information. Where diversification dynamics, rather than hypothesis free approaches that were the possible, the information from floristic literature was compared and supplemented main subject of criticism. Second, following the recommendations in ref. 46, we use with previously published11 or publicly accessible106 sources of point parametrization where only speciation or only extinction is variable in time, and we occurrence data. interpret their results acknowledging that speciation or extinction variability may The calcareous bedrock niche was defined by the regular presence of species on fall in the same congruence classes. Finally, we perform validation tests and a calcareous, dolomitic or ultrabasic bedrocks, the siliceous bedrock niche was sensitivity analysis with values realistic for our dataset to show that both single- defined by regular presence of species on any bedrocks that do not fall in the lineage and multi-clade models correctly identify parameter values from the calcareous category. For instance, in the case of species covered by Flora Alpina97, simulated data (see Supplementary Methods 2 and 3). species is considered present to a calcareous niche if regular presence is indicated either from limestones (ca) or serpentinites (ser), and to siliceous niche if regular presence is indicated from silicates (si), intermediate substrates (ca/si) or volcanic Evolutionary assembly across bedrock types, elevation belts and geographic rocks (bas). For details on the attribution of bedrock niches to individual species, regions. We used cladogenetic state-dependent speciation-extinction models47,48 see Supplementary Data 4. (ClaSSE), to study separately how diversification and migration occurred across The mid-elevation niche was defined by regular presence of species in habitats bedrock types, elevation belts and regions. In the models for bedrock types and below timberline, i.e. up to subalpine elevational zone sensu Flora Alpina97. High- elevational belts we used parametrization equivalent to GeoSSE47 that attributes elevation niche was defined by regular presence above timberline, i.e. in alpine and one of 3 states to each species—exclusive for one bedrock type or elevation belt; nival zone sensu Flora Alpina97. In the cases in which floristic information did not exclusive for another bedrock type or elevation belt; present in both bedrock types refer to elevational belts, specifically Flora Iberica98, we marked species presence in or elevation belts. For the model for geographic regions we developed a general- mid-elevation niche if the species was inhabiting habitats lower that 200 m below ization of GeoSSE for more than two regions (see Supplementary Methods 4 for regional timberline and in high elevation niche if the species was inhabiting script). In this generalization, the number of model states is growing exponentially habitats 200 m above regional timberline. No such treated species was restricted to with number of regions, which prevented us from using a more detailed division of the range ±200 m around timberline. For details on the attribution of elevational European mountain system than 5 regions (and thus 31 states). All analyses were niches to individual species, see Supplementary Data 4. run with the maximum credibility phylogenies for each of the six focal lineages. In order to compile information about species geographical distribution, we Similarly to time-dependent diversification models, we ran all models for each defined the smallest operational geographical units across Europe for which we lineage separately, and also a multi-clade model with parameters shared for all the could get credible presence/absence information of every species in our dataset. lineages, using the same procedure of likelihood multiplication as described for Specifically, we used operational geographic units in Flora Alpina97 for the Alps, time-dependent models and equivalent to one in refs. 47 or108. and mountain regions based on Körner et al.107 for other European mountains, For each evolutionary assembly model (on bedrock types, elevation belts, with subsequent modifications that better reflect structuring of biogeographic regions), we performed a series of AICdiff comparisons of maximum likelihood information in local floristic literature (see Supplementary Methods 5). We used model fits to test the importance of different model features: (i) presence state- this approach rather than grid-based processing of point occurrence data (as e.g. in change speciation (following the terminology of ClaSSE models48, that is, ref. 11), because of imbalanced point data quality depending on the country and speciation associated with splits of ancestor species occupying both bedrock or broader European region. For the purpose of ClaSSE models, operational elevation types, or multiple regions), (ii) difference of constant-state speciation geographic units were merged into five major mountain regions of Europe and rates for each bedrock type, elevation belt or region, (iii) difference of extinction surrounding areas (Fig. 1). For details on operational geographic units and how rates for each bedrock type, elevation belt or region, (iv) directionality of migration they were merged into five major regions, please refer to Supplementary Methods 5 between bedrock types, elevation belts and regions. The inferred best model was and Supplementary Data 4. subsequently rerun in a Bayesian setup in order to obtain probability envelopes of parameter estimates, using slice MCMC sampler run for 10,000 iterations and a burn-in of 5000 iterations. The proportions of state-change to all speciation events Tempo of species diversification. Based on the inferred phylogenies of each were obtained by multiplying the present-day numbers of species in different lineage, we fitted five models depicting different temporal dynamics of species model states (e.g. calcareous specialist, siliceous specialist, bedrock generalist) by diversification (see Table 1). For the models with temperature-dependent specia- their respective speciation rates. The analyses were performed using the R package tion or extinction, we used exponential dependence on 18O isotope ratio time- diversitree 0.9-11111. series from Greenland ice cores21. To accommodate for phylogenetic uncertainty, The SSE models are known to suffer from elevated Type I errors when testing we fitted every model on 100 trees randomly selected from the Bayesian posterior constant-state diversification rates differences between states, and various distribution of phylogenetic trees, and report either the whole distribution (para- correction strategies were proposed to address this issue74,112–116. However, this meter estimates) or median of values (AIC comparisons). issue is relevant mostly when the evolutionary states are stable74,112, which is To explore the overall diversification dynamics across all study lineages, we clearly not the case of our study plant lineages with the important role of migration developed a multi-clade framework to fit the above-described diversification and state-change speciation across bedrocks, elevational belts and regions. models across multiple evolutionary lineages simultaneously. We assumed that Moreover, all constant-state speciation and extinction rate differences presented in each of our lineages is an independent realization of shared diversification this paper turned unsupported or only weakly supported even without these dynamics, which allowed us to construct the joint likelihood functions of temporal corrections, and are interpreted as such. Instead, we used the SSE framework for diversification models as a product of the likelihood functions of each of the six testing and interpreting rates of migration and state-change speciation, and lineages, and harnessing the lineage-specific parameters into shared values. A subsequent ancestral state reconstructions, where the SSE methodology is similar approach with joint likelihood was previously used for state-dependent adequate117. For estimating state-change dynamics and ancestral state diversification models47,108. We optimized model parameters using the joint reconstructions, the SSE methodology is known to statistically outperform the likelihood function with a simplex routine, equivalent to the default alternative approaches not explicitly accounting for diversification dynamics, as is implementation of single-lineage models in the R package RPANDA 1.5109 (see BioGeoBEARS118 or Mkn model117. As recommended74,115, we also tested the Supplementary Methods 2 for implementation details). All calculations using R ability of used single lineage and multi-clade SSE to correctly recover parameters, packages were performed in R 3.5.2. using simulated datasets resembling our data (Supplementary Methods 4). To test for time- and temperature-dependence of diversification in both multi- clade and single-lineage analyses, we calculated AIC difference (AICdiff) between constant speciation and constant extinction model vs. the other respective models Sister-species overlap in geographical and ecological space. We performed an (Table 1). The constant speciation and constant extinction model is nested in the overlap analysis between sister species in order to address the importance and other respective models, and has one less degree of freedom. AICdiff = −2 thus frequency of allopatric speciation at two different spatial scales, along with spe- suggests that the focal model parameter does not improve likelihood at all, ciation between bedrock types and elevation belts. To do this, we identified all sister AICdiff = 0 suggests that both models are equally valuable from information- species pairs in maximum credibility phylogenies. For each of the pairs we cal- theoretic point of view, whereas AICdiff>2 suggests substantial support for the culated the Schoener’s D niche overlap index119 for 87 operational geographic focal model110, i.e. that the focal model would be outperforming the null model units (fine scale geography), five European mountain regions (coarse scale 10 NATURE COMMUNICATIONS | (2022)13:2750 | https://doi.org/10.1038/s41467-022-30394-5 | www.nature.com/naturecommunications NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-022-30394-5 ARTICLE geography), elevational belts (high elevation vs. mid-elevation) and bedrocks 18. Ding, W. N., Ree, R. H., Spicer, R. A. & Xing, Y. W. Ancient orogenic and (calcareous vs. siliceous). The overlap estimates were calculated using the R monsoon-driven assembly of the world’s richest temperate alpine flora. package spaa 0.2.2120. Science 369, 578–581 (2020). 19. Roquet, C., Boucher, F. C., Thuiller, W. & Lavergne, S. Replicated radiations of Reporting summary. Further information on research design is available in the Nature the alpine genus Androsace (Primulaceae) driven by range expansion and Research Reporting Summary linked to this article. convergent key innovations. J. Biogeogr. 40, 1874–1886 (2013). 20. Luebert, F. & Muller, L. A. H. Biodiversity from mountain building. Front. Genet. 6, (2015). 21. Zachos, J. C., Dickens, G. R. & Zeebe, R. E. An early Cenozoic perspective on Data availability greenhouse warming and carbon-cycle dynamics. Nature 451, 279–283 Raw genomic data generated in this study are available in European Nucleotide Archive (2008). under project accession codes https://www.ebi.ac.uk/ena/browser/view/PRJEB43865, 22. Haffer, J. Speciation in Colombian forest birds west of the Andes. Am. https://www.ebi.ac.uk/ena/browser/view/PRJEB48693, https://www.ebi.ac.uk/ena/ Museum Novit. 2294, 1–58 (1967). browser/view/PRJEB48874 and https://www.ebi.ac.uk/ena/browser/view/PRJEB50489. 23. Aguilée, R., Claessen, D. & Lambert, A. Adaptive radiation driven by the Accession codes for individual samples are provided in Supplementary Data 1. Species- interplay of eco-evolutionary and landscape dynamics. Evolution 67, level phylogenies, both maximum credibility and posterior samples, are available in 1291–1306 (2013). Supplementary Software 1. Geographic and ecological characteristics of different species 24. Feng, G., Mao, L., Sandel, B., Swenson, N. G. & Svenning, J. C. High plant are available in Supplementary Data 4. Source data are provided with the paper. endemism in China is partially linked to reduced glacial-interglacial climate change. J. Biogeogr. 43, 145–154 (2016). 25. Molina-Venegas, R., Aparicio, A., Lavergne, S. & Arroyo, J. Climatic and topographical correlates of plant palaeo- and neoendemism in a Code availability Mediterranean biodiversity hotspot. Ann. Bot. 119, 229–238 (2017). Custom R functions used for fitting multi-clade comparative models are available as a 26. Saladin, B. et al. Rapid climate change results in long-lasting spatial GitHub repository containing a R package (https://github.com/smyckaj/multidiv/tree/ homogenization of phylogenetic diversity. Nat. Commun. 11, 1–8 (2020). 6clades, https://doi.org/10.5281/zenodo.6341727 121). Examples with simulated data are 27. Hughes, C. & Eastwood, R. Island radiation on a continental scale: exceptional available in Supplementary Methods 2 and 4. rates of plant diversification after uplift of the Andes. Proc. Natl Acad. Sci. 103, 10334–10339 (2006). 28. Pouchon, C. et al. Phylogenomic analysis of the explosive adaptive radiation of Received: 7 October 2021; Accepted: 26 April 2022; the Espeletia complex (Asteraceae) in the tropical Andes. Syst. Biol. 67, 1041–1060 (2018). 29. Kadereit, J. W. The role of in situ species diversification for the evolution of high vascular plant species diversity in the European Alps—a review and interpretation of phylogenetic studies of the endemic flora of the Alps. Perspect. Plant Ecol. Evol. Syst. 26, 28–38 (2017). References 30. Escobar García, P. et al. Extensive range persistence in peripheral and 1. Hughes, C. E. & Atchinson, G. W. The ubiquity of alpine plant radiations: interior refugia characterizes Pleistocene range dynamics in a widespread from the Andes to the Hengduan Mountains. N. Phytol. 207, 275–282 (2015). Alpine plant species (Senecio carniolicus, Asteraceae). Mol. Ecol. 21, 2. Rahbek, C. et al. Humboldt’s enigma: what causes global patterns of mountain 1255–1270 (2012). biodiversity? Science 365, 1108–1113 (2019). 31. Lohse, K., Nicholls, J. A. & Stone, G. N. Inferring the colonization of a 3. Antonelli, A. et al. Geological and climatic influences on mountain mountain range-refugia vs. nunatak survival in high alpine ground beetles. biodiversity. Nat. Geosci. 11, 718–725 (2018). Mol. Ecol. 20, 394–408 (2011). 4. Quintero, I. & Jetz, W. Global elevational diversity and diversification of birds. 32. Stehlik, I. Resistance or emigration? Response of alpine plants to the ice ages. Nature 555, 246–250 (2018). Taxon 52, 499–510 (2003). 5. Merckx, V. S. F. T. et al. Evolution of endemism on a young tropical 33. Schneeweiss, G. M. & Schönswetter, P. A re-appraisal of nunatak survival in mountain. Nature 524, 347–350 (2015). arctic-alpine phylogeography. Mol. Ecol. 20, 190–192 (2011). 6. Körner, C. Alpine Plant Life (Springer, 1999). 34. Westergaard, K. B. et al. Glacial survival may matter after all: Nunatak 7. Smyčka, J. et al. Reprint of: Disentangling drivers of plant endemism and signatures in the rare European populations of two west-arctic species. Mol. diversification in the European Alps - a phylogenetic and spatially explicit Ecol. 20, 376–393 (2011). approach. Perspect. Plant Ecol. Evol. Syst. 30, 31–40 (2018). 35. Bettin, O., Cornejo, C., Edwards, P. J. & Holderegger, R. Phylogeography of 8. Schönswetter, P., Stehlik, I., Holderegger, R. & Tribsch, A. Molecular evidence the high alpine plant Senecio halleri (Asteraceae) in the European Alps: In situ for glacial refugia of mountain plants in the European Alps. Mol. Ecol. 14, glacial survival with postglacial stepwise dispersal into peripheral areas. Mol. 3547–3555 (2005). Ecol. 16, 2517–2524 (2007). 9. Haller, A. von. Enumeratio Methodica Stirpium Helvetiae indigenarum. 36. Tomasello, S., Karbstein, K., Hodač, L., Paetzold, C. & Hörandl, E. (Officina Academica Abrami Vandenhoek, 1742). Phylogenomics unravels Quaternary vicariance and allopatric speciation 10. de Candolle, A. Sur les causes de l’inégale distribution des plantes rares dans la patterns in temperate-montane plant species: a case study on the Ranunculus chaîne des Alpes. Atti del Congr. Internazionale Bot. Tenuto Firenze. 92–104 auricomus species complex. Mol. Ecol. 29, 2031–2049 (2020). (1875). 37. Ozenda, P. L’endémisme au niveau de l’ensemble du Système alpin. Acta Bot. 11. Boucher, F. C., Zimmermann, N. E. & Conti, E. Allopatric speciation with Gall. 142, 753–762 (1995). little niche divergence is common among alpine Primulaceae. J. Biogeogr. 43, 38. Rolland, J., Lavergne, S. & Manel, S. Combining niche modelling and 591–602 (2016). landscape genetics to study local adaptation: A novel approach illustrated 12. Schneeweiss, G. M. et al. Molecular phylogenetic analyses identify Alpine using alpine plants. Perspect. Plant Ecol. Evol. Syst. 17, 491–499 (2015). differentiation and dysploid chromosome number changes as major forces for 39. Alvarez, N. et al. History or ecology? Substrate type as a major driver of spatial the evolution of the European endemic Phyteuma (Campanulaceae). Mol. genetic structure in Alpine plants. Ecol. Lett. 12, 632–640 (2009). Phylogenet. Evol. 69, 634–652 (2013). 40. Gao, Y.-D., Gao, X.-F. & Harris, A. Species boundaries and parapatric 13. Tkach, N. et al. Molecular phylogenetics, morphology and a revised speciation in the complex of alpine shrubs, Rosa sericea (Rosaceae), based on classification of the complex genus Saxifraga (Saxifragaceae). Taxon 64, population genetics and ecological tolerances. Front. Plant Sci. 10, 1–16 1159–1187 (2015). (2019). 14. Favre, A. et al. Out-of-Tibet: the spatio-temporal evolution of Gentiana 41. Knox, E. B. Adaptive radiation of African montane plants. In Adaptive (Gentianaceae). J. Biogeogr. 43, 1967–1978 (2016). Speciation (eds. Dieckmann, U., Doebeli, M., Metz, J. A. J. & Tautz, D.) 15. Kadereit, J. W., Griebeler, E. M. & Comes, H. Quaternary diversification in 345–361 (Cambridge University Press, 2004). European alpine plants: pattern and process. Philos. Trans. R. Soc. Lond. B. 42. Segar, S. T. et al. Speciation in a keystone plant genus is driven by elevation: a Biol. Sci. 359, 265–274 (2004). case study in New Guinean Ficus. J. Evol. Biol. 30, 512–523 (2017). 16. Xing, Y. & Ree, R. H. Uplift-driven diversification in the Hengduan 43. Pouchon, C. et al. Phylogenetic signatures of ecological divergence and Mountains, a temperate biodiversity hotspot. Proc. Natl Acad. Sci. 114, leapfrog adaptive radiation in Espeletia. Am. J. Bot. 108, 113–128 (2021). E3444–E3451 (2017). 44. Luebert, F. & Weigend, M. Phylogenetic insights into Andean plant 17. Lagomarsino, L. P., Condamine, F. L., Antonelli, A., Mulch, A. & Davis, C. C. diversification. Front. Ecol. Evol. 2, 1–17 (2014). The abiotic and biotic drivers of rapid diversification in Andean bellflowers 45. Nagy, L. & Grabherr, G. The Biology of Alpine Habitats (Oxford University (Campanulaceae). N. Phytol. 210, 1430–1442 (2016). Press, 2009). NATURE COMMUNICATIONS | (2022)13:2750 | https://doi.org/10.1038/s41467-022-30394-5 | www.nature.com/naturecommunications 11 ARTICLE NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-022-30394-5 46. Louca, S. & Pennell, M. W. Extant timetrees are consistent with a myriad of 76. Dentant, C. & Lavergne, S. Plantes de haute montagne: état des lieux, diversification histories. Nature 580, 502–505 (2020). évolution et analyse diachronique dans le massif des Écrins (France). Bull. Soc. 47. Goldberg, E. E., Lancaster, L. T. & Ree, R. H. Phylogenetic inference of linn. Provence 64, 83–98 (2013). reciprocal effects between geographic range evolution and diversification. Syst. 77. Dentant, C. The highest vascular plants on Earth. Alp. Bot. 128, 97–106 (2018). Biol. 60, 451–465 (2011). 78. Boucher, F. C. et al. Reconstructing the origins of high‐alpine niches and 48. Goldberg, E. E. & Igić, B. Tempo and mode in plant breeding system cushion life form in the genus Androsace sl (Primulaceae). Evolution 66, evolution. Evolution 66, 3701–3709 (2012). 1255–1268 (2012). 49. Gitzendanner, M., Soltis, P., Yi, T.-S., Li, D.-Z. & Soltis, D. Plastome 79. Boucher, F. C., Lavergne, S., Basile, M., Choler, P. & Aubert, S. Evolution and Phylogenetics: 30 years of inferences into plant evolution. In Advances in biogeography of the cushion life form in angiosperms. Perspect. Plant Ecol. Botanical Research 293–313 (Elsevier, 2018). Evol. Syst. 20, 22–31 (2016). 50. Birks, H. H. The late-quaternary history of arctic and alpine plants. Plant Ecol. 80. Schönswetter, P. & Schneeweiss, G. M. Is the incidence of survival in interior Divers. 1, 135–146 (2008). Pleistocene refugia (nunataks) underestimated? Phylogeography of the high 51. Mai, D. Tertiäre Vegetationsgeschichte Europas—Metoden und Ergebnisse. mountain plant Androsace alpina (Primulaceae) in the European Alps (Gustav Fischer Verlag, 1995). revisited. Ecol. Evol. 9, 4078–4086 (2019). 52. Svenning, J. C. Deterministic Plio-Pleistocene extinctions in the European 81. Aeschimann, D., Rasolofo, N. & Theurillat, J. P. Analyse de la flore des Alpes. cool-temperate tree flora. Ecol. Lett. 6, 646–653 (2003). 2: Diversité et chorologie. Candollea 66, 225–253 (2011). 53. Fauquette, S. et al. The Alps: a geological, climatic and human perspective 82. Ebersbach, J. et al. In and out of the Qinghai-Tibet Plateau: divergence time on vegetation history and modern plant diversity. In Mountains, Climate estimation and historical biogeography of the large arctic-alpine genus and Biodiversity (eds. Hoorn, C., Perrigo, A. & Antonelli, A.) 413 (Wiley- Saxifraga L. J. Biogeogr. 44, 900–910 (2017). Blackwell, 2018). 83. Hannon, G. FASTX. http://hannonlab.cshl.edu/fastx_toolkit/ (2014). 54. Mráz, P. et al. Vascular plant endemism in the Western Carpathians: spatial 84. Coissac, E. The ORGanelle ASseMbler 1.0.3. https://git.metabarcoding.org/ patterns, environmental correlates and taxon traits. Biol. J. Linn. Soc. 119, org-asm/org-asm/wikis/home (2016). 630–648 (2016). 85. Shaw, J. et al. Chloroplast DNA sequence utility for the lowest phylogenetic 55. Puşcaş, M. et al. Post-glacial history of the dominant alpine sedge Carex and phylogeographic inferences in angiosperms: the tortoise and the hare IV. curvula in the European Alpine System inferred from nuclear and chloroplast Am. J. Bot. 101, 1987–2004 (2014). markers. Mol. Ecol. 17, 2417–2429 (2008). 86. Mansion, G. et al. How to handle speciose clades? Mass taxon-sampling as a 56. Puşcaş, M., Taberlet, P. & Choler, P. No positive correlation between species strategy towards illuminating the natural history of Campanula and genetic diversity in European alpine grasslands dominated by Carex (Campanuloideae). PLoS ONE 7, e50076 (2012). curvula. Divers. Distrib. 14, 852–861 (2008). 87. Rossi, M. Taxonomy, phylogeny and reproductive ecology of Gentiana lutea L 57. Magyari, E. K. et al. Late Pleniglacial vegetation in eastern-central Europe: are (University in Bologna, 2011). there modern analogues in Siberia? Quat. Sci. Rev. 95, 60–79 (2014). 88. Hämmerli, M. Molecular Aspects in Systematics of Gentiana Sect. Calathianae 58. Prodon, R., Thibault, J. C. & Dejaifve, P. A. Expansion vs. compression of bird Froel (Université de Neuchâtel, 2007). altitudinal ranges on a Mediterranean island. Ecology 83, 1294–1306 (2002). 89. Hungerer, K. B. & Kadereit, J. W. The phylogeny and biogeography of 59. Moen, D. & Morlon, H. Why does diversification slow down? Trends Ecol. Gentiana L. sect. Ciminalis (Adans.) Dumort.: A historical interpretation of Evol. 29, 190–197 (2014). distribution ranges in the European high mountains. Perspect. Plant Ecol. Evol. 60. Aguilée, R., Gascuel, F., Lambert, A. & Ferriere, R. Clade diversification Syst. 1, 121–135 (1998). dynamics and the biotic and abiotic controls of speciation and extinction rates. 90. Ranwez, V., Harispe, S., Delsuc, F. & Douzery, E. J. P. MACSE: Multiple Nat. Commun. 9, 1–13 (2018). alignment of coding SEquences accounting for frameshifts and stop codons. 61. Vargas, P. Molecular evidence for multiple diversification patterns of alpine PLoS One 6, e22594 (2011). plants in Mediterranean Europe. Taxon 52, 463–476 (2003). 91. Castresana, J. Selection of conserved blocks from multiple alignments for their 62. Kruckeberg, A. R. An essay: the stimulus of unusual geologies for plant use in phylogenetic analysis. Mol. Biol. Evol. 17, 540–552 (2000). speciation. Syst. Bot. 11, 455–463 (1986). 92. Kück, P. & Meusemann, K. FASconCAT: convenient handling of data 63. Cowling, R. M. & Holmes, P. M. Endemism and speciation in a lowland flora matrices. Mol. Phylogenet. Evol. 56, 1115–1118 (2010). from the Cape Floristic Region. Biol. J. Linn. Soc. 47, 367–383 (1992). 93. Katoh, K., Kuma, K. I., Toh, H. & Miyata, T. MAFFT version 5: Improvement in 64. Lexer, C. et al. Genomics of the divergence continuum in an African plant accuracy of multiple sequence alignment. Nucleic Acids Res. 33, 511–518 (2005). biodiversity hotspot, I: drivers of population divergence in Restio capensis 94. Bouckaert, R. et al. BEAST 2: a software platform for Bayesian evolutionary (Restionaceae). Mol. Ecol. 23, 4373–4386 (2014). analysis. PLoS Comput. Biol. 10, 1–6 (2014). 65. Anacker, B. L. & Strauss, S. Y. The geography and ecology of plant speciation: 95. Bouckaert, R. R. & Drummond, A. J. bModelTest: Bayesian phylogenetic site range overlap and niche divergence in sister species. Proc. R. Soc. B Biol. Sci. model averaging and model comparison. BMC Evol. Biol. 17, 1–11 (2017). 281, 20132980 (2014). 96. Morlon, H. Phylogenetic approaches for studying diversification. Ecol. Lett. 66. Moore, A. J. & Kadereit, J. W. The evolution of substrate differentiation in 17, 508–525 (2014). Minuartia series Laricifoliae (Caryophyllaceae) in the European Alps: in situ 97. Aeschimann, D., Lauber, K., Moser, D. M. & Theurillat, J. P. Flora Alpina origin or repeated colonization? Am. J. Bot. 100, 2412–2425 (2013). (Editions Belin, 2004). 67. Guggisberg, A. et al. The genomic basis of adaptation to calcareous and 98. Castroviejo, S. Flora Iberica (Real Jardin Botanico CSIC, 2012). siliceous soils in Arabidopsis lyrata. Mol. Ecol. 27, 5088–5103 (2018). 99. Goliášová, K. & Michalková, E. Flóra Slovenska (Vydavateľstvo Slovenskej 68. Gigon, A. Vergleich alpiner Rasen auf Silikat- und auf Karbonatboden— akadémie vied, 2012). Konkurrenz—und Stickstofformenversuche sowie standortskundliche 100. Speta, E. & Rákosy, L. Wildpflanzen Siebenbürgen (Naturhistorisches Museum Untersuchungen im Nardetum und im Seslerietum bei Davos. (ETH Zuerich, Wien, 2010). 1971). 101. Sarić, M. Flora Srbije (Srpska akademija nauka i umetnosti, 1992). 69. Davies, M. S. & Snaydon, R. W. Physiological differences among 102. Schönswetter, P. & Schneeweiss, G. M. Androsace komovensis sp. nov., a long populations of Anthoxanthum odoratum L. collected from the park grass mistaken local endemic from the southern Balkan Peninsula with experiment, Rothamsted. I. Response to calcium. J. Appl. Ecol. 10, 33–45 biogeographic links to the Eastern Alps. Taxon 58, 544–549 (2009). (1973). 103. Schönswetter, P., Magauer, M. & Schneeweiss, G. M. Androsace halleri subsp. 70. Snaydon, R. W. Rapid population differentiation in mosaic environment. I. nuria Schönsw. & Schneew. (Primulaceae), a new taxon from the eastern The response of Anthoxantum odoratum populations to soils. Evolution 24, Pyrenees (Spain, France). Phytotaxa 201, 227–232 (2015). 257–269 (1970). 104. Schneeweiss, G. M. & Schonswetter, P. The wide but disjunct range of the 71. Zohlen, A. & Tyler, G. Soluble inorganic tissue phosphorus and calcicole- European mountain plant Androsace lactea L. (Primulaceae) reflects Late calcifuge behaviour of plants. Ann. Bot. 94, 427–432 (2004). Pleistocene range fragmentation and post-glacial distributional stasis. J. 72. Kassen, R., Llewellyn, M. & Rainey, P. B. Ecological contraints on Biogeogr. 37, 2016–2025 (2010). diversification in a model adaptive radiation. Nature 431, 984–988 (2004). 105. Webb, D. A. & Gornall, R. J. Saxifrages of Europe (Timber Press, 1989). 73. MacLean, R. C., Bell, G. & Rainey, P. B. The evolution of a pleiotropic fitness 106. GBIF. https://www.gbif.org/ (2018). tradeoff in Pseudomonas fluorescens. Proc. Natl Acad. Sci. USA 101, 107. Körner, C. et al. A global inventory of mountains for bio-geographical 8072–8077 (2004). applications. Alp. Bot. 127, 1–15 (2017). 74. Rabosky, D. L. & Goldberg, E. E. Model inadequacy and mistaken inferences 108. Anacker, B. L., Whittall, J. B., Goldberg, E. E. & Harrison, S. P. Origins and of trait-dependent speciation. Syst. Biol. 64, 340–355 (2015). consequences of serpentine endemism in the California flora. Evolution 65, 75. Kolář, F. et al. Northern glacial refugia and altitudinal niche divergence shape 365–376 (2011). genome-wide differentiation in the emerging plant model Arabidopsis 109. Morlon, H. et al. RPANDA: An R package for macroevolutionary analyses on arenosa. Mol. Ecol. 25, 3929–3949 (2016). phylogenetic trees. Methods Ecol. Evol. 7, 589–597 (2016). 12 NATURE COMMUNICATIONS | (2022)13:2750 | https://doi.org/10.1038/s41467-022-30394-5 | www.nature.com/naturecommunications NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-022-30394-5 ARTICLE 110. Burnham, K. & Anderson, D. Model Selection and Multimodel Inference Framework Programme FP7/2007-2013 grant agreement 281422 (TEEMBIO, W.T.) and (Springer, 2002). by the SNF grant (31003A_149508, N.E.Z.), the sampling campaigns in the Balkans and 111. Fitzjohn, R. G., Maddison, W. P. & Otto, S. P. Estimating trait-dependent in the Iberian Peninsula were funded by the French Ecological Society (SFE²) awarded to speciation and extinction rates from incompletely resolved phylogenies. Syst. J.S., and by the Systematics Research Fund 2016 of the Linnean Society of London and Biol. 58, 595–611 (2009). the Systematics Association awarded to C.R. J.S. was supported by a grant from the 112. O’Meara, B. C. & Beaulieu, J. M. Past, future, and present of state-dependent Doctoral school of Chemistry and Life Sciences within the Univ. Grenoble Alpes, by the models of diversification. Am. J. Bot. 103, 792–795 (2016). ANR project Sphinx (ANR-16-CE02-0011, S.L.) and by Czech Science Foundation 113. Beaulieu, J. M. & O’Meara, B. C. Detecting hidden diversification shifts in (GAČR 20-29554X, J.S.). K.Š. was supported by Charles University project GAUK models of trait-dependent speciation and extinction. Syst. Biol. 65, 583–601 (815516, K.Š.), the Mobility funds of Charles University and a long-term research project (2016). of the Czech Academy of Sciences (RVO 67985939, K.Š.). 114. Herrera-Alsina, L., Van Els, P. & Etienne, R. S. Detecting the dependence of diversification on multiple traits from phylogenetic trees and trait data. Syst. Biol. 68, 317–328 (2019). Author contributions J.S., C.R., and S.L. designed the study, J.S. performed the analyses, and J.S., C.R., and S.L. 115. Onstein, R. E. et al. To adapt or go extinct? The fate of megafaunal palm fruits wrote the first draft of the manuscript. J.S., C.R., M.B., R.D., K.Š., N.E.Z. W.T., I.G.A., S.L. under past global change. Proc. R. Soc. B Biol. Sci. 285, (2018). and other members of PhyloAlps consortium contributed to the sample collection and its 116. Rabosky, D. L. & Goldberg, E. E. FiSSE: a simple nonparametric test for the coordination. C.P., M.R., R.D., and JG.V. coordinated and mounted the reference her- effects of a binary character on lineage diversification rates. Evolution 71, barium collection. M.B., A.A., F.D, P.W., F.B., and E.C. designed and performed the wet 1432–1442 (2017). lab procedures, sequencing and bioinformatic pipelines. The project as a whole was 117. Holland, B. R., Ketelaar-Jones, S., O’Mara, A. R., Woodhams, M. D. & coordinated by S.L. All authors contributed to the final version of the manuscript. Jordan, G. J. Accuracy of ancestral state reconstruction for non-neutral traits. Sci. Rep. 10, 1–10 (2020). 118. Ree, R. H. & Sanmartín, I. Conceptual and statistical problems with the DEC+J Competing interests model of founder-event speciation and its comparison with DEC via model The authors declare no competing interests. selection. J. Biogeogr. 45, 741–749 (2018). 119. Schoener, T. W. Nonsynchronous spatial overlap of lizards in patchy habitats. Ecology 51, 408–418 (1970). Additional information 120. Zhang, J. spaa: SPecies Association Analysis 0.2.2. https://cran.r-project.org/ Supplementary information The online version contains supplementary material package=spaa (2016). available at https://doi.org/10.1038/s41467-022-30394-5. 121. Smyčka, J. Tempo and drivers of plant diversification in the European mountain system. multidiv, https://doi.org/10.5281/zenodo.6341727 (2022). Correspondence and requests for materials should be addressed to Jan Smyčka. Peer review information Nature Communications thanks Gabriele Casazza, Graham Acknowledgements Wallis and the other, anonymous, reviewer(s) for their contribution to the peer review of This paper is dedicated to the memory of Prof. Serge Aubert (1966‐2015) who devoted this work. Peer reviewer reports are available. his career to the study of alpine plants and was an invaluable source of motivation for the PhyloAlps project. We thank P. Taberlet, F. Condamine, C. Graham, G. Schneeweiss, Reprints and permission information is available at http://www.nature.com/reprints P. Schönswetter, N. Alvarez, F. Kolář, M. Slovák, T. Herben, D. Storch, M. Elias and L. Pollock for helpful comments, guidance and insights at different stages of the project, Publisher’s note Springer Nature remains neutral with regard to jurisdictional claims in M. Smyčková, D. Požárová, M. Kolářová, J. Pilátová, T. Figura, L. Boulangeat, A. Frat- published maps and institutional affiliations. taroli, P. Koutecký, C. Lagaye, J. Šibík, P. Anastasiu, L. Filipas, R. Giordano, P. D. Turtureanu, A. Guttová, J. Kučera, M Valachovič, I. Hodalová, J. Kochjarová, E. Štub- ňová, A. Pleceníková, S. Španiel, G. Casazza, S. Wipf, R. Geremia, P. Choler, N. Passa- Open Access This article is licensed under a Creative Commons lacqua, L. Hugot, T. Engelskjon, T. Alm, S. Eie, H. Johansen, H. Edvardsen, E. W. Attribution 4.0 International License, which permits use, sharing, Hanssen, C. Bay, B.-G. Osterkloft and L. Sáez for help with sampling, and M. K. F. adaptation, distribution and reproduction in any medium or format, as long as you give Merkel and Y. Lammers for help with molecular and bioinformatic analyses. The appropriate credit to the original author(s) and the source, provide a link to the Creative National park of High Tatra, the DNA Bank of the Natural History Museum Oslo, Commons license, and indicate if changes were made. The images or other third party Tromsø Museum, University in Innsbruck, St Andrews Botanic Garden and the Bota- material in this article are included in the article’s Creative Commons license, unless nical Garden of the Charles University in Prague kindly provided us plant materials. The indicated otherwise in a credit line to the material. If material is not included in the research was funded by the joint ANR-SNF project Origin-Alps (ANR-16-CE93-0004, article’s Creative Commons license and your intended use is not permitted by statutory SNF-310030L_170059, S.L. and N.E.Z.) and the OSUG@2020 labex (ANR10 LABX56, regulation or exceeds the permitted use, you will need to obtain permission directly from S.L.). The sequencing was performed within the framework of the PhyloAlps project, the copyright holder. To view a copy of this license, visit http://creativecommons.org/ funded by France Génomique (ANR-10-INBS-09-08, P.W.), and the PhyloNorway licenses/by/4.0/. project funded by the Research Council of Norway (226134/F50, I.G.A.) and the Nor- wegian Biodiversity Information Centre (14-14, 70184209, I.G.A.). Bioinformatics and statistical analyses were carried out with the GRICAD infrastructure (https://gricad.univ- © The Author(s) 2022 grenoble-alpes.fr). The sampling campaign and preliminary genomic analyses were partly funded by the European Research Council under the European Community’s Seventh the PhyloAlps consortium Cristina Roquet 1,4, Martí Boleda1, Adriana Alberti5,6, Frédéric Boyer 1, Rolland Douzet7, Christophe Perrier 7, Maxime Rome7, Jean-Gabriel Valay7, France Denoeud 5, Niklaus E. Zimmermann 9, Wilfried Thuiller1, Patrick Wincker 5, Inger G. Alsos 10, Eric Coissac1 & Sébastien Lavergne1 A full list of members and their affiliations appears in the Supplementary Information. NATURE COMMUNICATIONS | (2022)13:2750 | https://doi.org/10.1038/s41467-022-30394-5 | www.nature.com/naturecommunications 13