Molecular Phylogenetics and Evolution 189 (2023) 107928 Contents lists available at ScienceDirect Molecular Phylogenetics and Evolution journal homepage: www.elsevier.com/locate/ympev Impact of the climatic changes in the Pliocene-Pleistocene transition on Irano-Turanian species. The radiation of genus Jurinea (Compositae) Sonia Herrando-Moraira a, Cristina Roquet b, *, The Cardueae Radiations Group (in alphabetical order): a Botanic Institute of Barcelona (IBB, CSIC-Ajuntament de Barcelona), Pg. del Migdia, s.n., 08038 Barcelona, Spain b 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, 08193 Bellaterra, Spain A R T I C L E I N F O A B S T R A C T Keywords: The Irano-Turanian region is one of the world’s richest floristic regions and the centre of diversity for numerous Asteraceae xerophytic plant lineages. However, we still have limited knowledge on the timing of evolution and biogeo­ Biogeography graphic history of its flora, and potential drivers of diversification remain underexplored. To fill this knowledge Diversification gap, we focus on the Eurasian genus Jurinea (ca. 200 species), one of the largest plant radiations that diversified Hyb-seq in the region. We applied a macroevolutionary integrative approach to explicitly test diversification hypotheses Irano-Turanian region Phylogeny and investigate the relative roles of geography vs. ecology and niche conservatism vs. niche lability in speciation processes. To do so, we gathered a sample comprising 77% of total genus richness and obtained data about (1) its phylogenetic history, recovering 502 nuclear loci sequences; (2) growth forms; (3) ecological niche, compiling data of 21 variables for more than 2500 occurrences; and (4) paleoclimatic conditions, to estimate climatic stability. Our results revealed that climate was a key factor in the evolutionary dynamics of Jurinea. The main diversification and biogeographic events that occurred during past climate changes, which led to colder and drier conditions, are the following: (1) the origin of the genus (10.7 Ma); (2) long-distance dispersals from the Iranian Plateau to adjacent regions (~7–4 Ma); and (3) the diversification shift during Pliocene-Pleistocene Transition (ca. 3 Ma), when net diversification rate almost doubled. Our results supported the pre-adaptation hypothesis, i. e., the evolutionary success of Jurinea was linked to the retention of the ancestral niche adapted to aridity. Interestingly, the paleoclimatic analyses revealed that in the Iranian Plateau long-term climatic stability favoured old-lineage persistence, resulting in current high species richness of semi-arid and cold adapted clades; whereas moderate climate oscillations stimulated allopatric diversification in the lineages distributed in the Circumboreal region. In contrast, growth form lability and high niche disparity among closely related species in the Central Asian clade suggest adaptive radiation to mountain habitats. In sum, the radiation of Jurinea is the result of both adaptive and non-adaptive processes influenced by climatic, orogenic and ecological factors. 1. Introduction (Hughes and Atchison, 2015). This type of research has allowed in­ vestigators to identify the main evolutionary and environmental factors The study of evolutionary radiations is a hot topic in plant science shaping biodiversity patterns in these hotspots and assess their relative research (Hughes et al., 2015), with the number of published articles on contributions. However, there are still species-rich regions where con­ plant radiations increasing dramatically since 2000s (Stroud and Losos, tinental radiations have taken place that remain understudied from an 2016). Most studies to date have focused on lineages with outstanding evolutionary and biogeographic perspective. One of such knowledge species richness in plant diversity hotspots such as the Cape Floristic gaps concerns the Irano-Turanian floristic region (hereafter IT; Takh­ Region, the Mediterranean Basin, oceanic islands and archipelagos such tajan, 1986), whose biodiversity origins and biogeographic relationships as Macaronesia, or mountain ranges such as the Andes or the Himalayas with other regions remain poorly studied despite its relevance in terms * Corresponding author at: 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, 08193 Bellaterra, Spain. E-mail address: cristina.roquet@uab.cat (C. Roquet). https://doi.org/10.1016/j.ympev.2023.107928 Received 30 March 2023; Received in revised form 7 September 2023; Accepted 12 September 2023 Available online 14 September 2023 1055-7903/© 2023 The Author(s). Published by Elsevier Inc. This is an open access article under the CC BY license (http://creativecommons.org/licenses/by/4.0/). S. Herrando-Moraira and C. Roquet Molecular Phylogenetics and Evolution 189 (2023) 107928 of plant diversity (32,000 species; Takhtajan, 1986; Sales and Hedge, processes, and the ecological niche has been preserved among sister taxa 2013) and endemism (25–40%; Zohary, 1981; Takhtajan, 1986). Indeed, (niche conservatism; Wiens et al., 2010). In the second, lineages this region harbours two global biodiversity hotspots: the Irano- diverged by adaptation to new ecological zones (niche divergence). Anatolian region and the Mountains of Central Asia (Mittermeier Concerning the biogeographic interactions of IT with other regions, et al., 2011). previous studies have suggested that IT constituted an important source The IT region is located in Western-Central Asia, surrounded by the of taxa to adjacent regions, especially contributing to the assemblage of Circumboreal floristic region in the north, the eastern Mediterranean in the Mediterranean flora (Jabbour and Renner, 2011; Salvo et al., 2011; the west, Arabia in the south, and East Asia in the east (see Manafzadeh Manafzadeh et al., 2014; Lauterbach et al., 2019; Moharrek et al., 2019; et al., 2017 for a review of its circumscription). One of most remarkable Peterson et al., 2019; Mahmoudi Shamsabad et al., 2021). Two main characteristics of IT is its continental climate, with cold winters, warm- hypotheses have been proposed regarding the temporal framework of dry summers, and a high precipitation seasonality (Djamali et al., 2012). dispersals from IT to the Mediterranean region: (1) during the early to Topographically, IT presents high landscape heterogeneity and wide late Miocene (Zohary, 1973; Manafzadeh et al., 2014); and (2) during altitudinal gradients, with large mountains ranges (Taurus, Zagros, Late Quaternary (Magyari et al., 2008). On the other hand, it is unknown Alborz, Kopet-Dagh, Pamir, Tian Shan, Western Himalayas) and broad whether IT also constituted an important source of xerophytic lineages plateaus (Anatolian, Iranian, and Qinghai-Tibetan plateaus). The geo- for other hotspots of the rest of Eurasia. climatic configuration of IT has likely favoured the diversification of Integrative studies on mega-genera based on well-resolved phylog­ xerophytes, which are adapted to arid or semi-arid environments under enies are needed to significantly advance in our understanding of the low water availability conditions, such as steppe elements (e.g., mega- relative roles of geological, climatic, biogeographic and ecological fac­ diverse genera like Acantholimon, Acanthophyllum, Astragalus, Cousinia, tors in the diversification of the IT flora (Manafzadeh et al., 2017). Here, Eremostachys and Eremurus; Manafzadeh et al., 2017) or C4 metabolism we studied one of the largest plant radiations of the IT: the genus Jurinea plants (e.g. Eragrostis, Stipagrostis; Rudov et al., 2020). Diversification is Cass. (tribe Cardueae, Compositae), which includes ca. 200 species, of especially noticeable in mountain habitats: for instance, endemic rich­ which ca. 70% occur in the IT (Susanna and Garcia-Jacas, 2007). Jurinea ness of Compositae in Iran reaches maximum values at mid-elevations in spans throughout Eurasia from latitude 20◦ N to 60◦ N, and from longi­ mountain ranges (1400–2100 m a.s.l.; Noroozi et al., 2019). tude 10◦ W to 120◦ E, and it is found in species-rich areas of the IT, However, limited information is available about the origins and especially in Iran, Afghanistan, and Central Asia Mountains (Rechinger drivers of diversification of the IT flora (see a synthesis in Manafzadeh and Wagenitz, 1979; Szukala et al., 2019). The genus is also present et al., 2017). Regarding the timing of the origin of IT elements, several along the Mediterranean Basin, central-eastern Europe, European hypotheses have been proposed. The oldest scenario points to a Creta­ Russia, Himalayas, and East Asia. The predominant altitudinal range is ceous origin (145–66 Ma; Takhtajan, 1986). Other studies suggest the around 1000–2000 m, although some species grow near the sea-level early Eocene as a period of initial divergence of lineages (55 Ma; Man­ (like J. kilaea Azn.; Danin and Davis, 1975) and others up to the afzadeh et al., 2014), or during the drying of the Tethys Sea in the alpine belt (4500 m like J. gilesii (Hemsl.) N.Garcia, Herrando & Sus­ Neogene (34 Ma; Zohary, 1973). Nevertheless, recent time-calibrated anna; Rechinger and Wagenitz, 1979). Ecologically, Jurinea species phylogenies of diverse groups suggest that IT elements originated mainly grow on mountain and rocky slopes, grasslands, steppes, pine- much more recently, from the Miocene onwards (Wu et al., 2015; Lau­ forest stands, forest margins, lake and riverbanks, and subalpine- terbach et al., 2019), with special emphasis on the Pliocene-Pleistocene alpine meadows (Iljin, 1962; Danin and Davis, 1975; Rechinger and (Moharrek et al., 2019; Mahmoudi Shamsabad et al., 2021). Regarding Wagenitz, 1979). Several growth forms are found within the genus, each potential drivers of diversification, the orogenic activity during the associated to particular habitats and reflecting different ecological Miocene deeply shaped the IT landscape (Manafzadeh et al., 2017) and strategies: (1) the only annual species, J. modesta Boiss., shows a pecu­ likely impacted the evolution of its flora: the second collision of the Afro- liar growth form intermediate between the acaulescent rosette and the Arabian plate against the Eurasian plate originated the uplift of the scapose form, and is adapted to the extreme arid habitat found in the Iranian Plateau ca. 13 Ma, and subsequent uplift and deformations foothills of the mountains of Afghanistan; (2) biennial plants charac­ occurred from 15 to 5 Ma along the Zagros, Alborz, Kopeth Dagh, and terized by long-foliose stems with soft and wide leaves (very different Caucasus mountains (Mouthereau, 2011). Such orogenic activity may from the usually narrow, scariose or leathery leaves of all other Jurinea have been an important driver of species diversification increasing species) are found in humid valleys south of the Himalayas; (3) rosetted habitat heterogeneity and thus ecological opportunities and favouring acaulescent perennials are always mountain-adapted plants; (4) scapose allopatric speciation through vicariance. On the other hand, the distri­ rosetted perennials usually represent an adaptation to Mediterranean bution and diversity of organisms have also likely been impacted by past climate or steppe with summer drought; and (5) foliose subshrubs are climatic changes such as: (1) the aridification of central Asia during mid- mostly adapted to semideserts. late Miocene (17–5 Ma; Miao et al., 2012); (2) the global climate cooling The complex taxonomy of Jurinea has contributed to its understudied between ~3.2–2.7 Ma, a period known as the Pliocene-Pleistocene status. Historically, up to nine satellite genera have been recognized transition (PPT; Lisiecki and Raymo, 2005); and (3) the drastic cli­ within it, though a recent phylogeny of subtribe Saussureinae based on matic oscillations during the Quaternary glacial-interglacial cycles Hyb-Seq data has shown that all of them are more properly accommo­ (Zachos et al., 2001). dated as part of the large genus Jurinea (Herrando-Moraira et al., 2020, Recent methodological advances in comparative phylogenetic ana­ which included ca. 40 Jurinea spp.). The most complete phylogeny of lyses (Borges et al., 2019), biogeographic history inference (Matzke, Jurinea published to date includes only 81 spp., which represents one 2013), and models of diversification (e.g. Etienne et al., 2012; Morlon third of the genus diversity (Szukala et al., 2019). Although it provides a et al., 2016), together with ecological niche modelling and high- substantial contribution on the phylogenetic relationships, it is still resolution paleoclimatic datasets (PaleoClim; Brown et al., 2018), notably incomplete in terms of taxonomic and biogeographic sampling, make it now possible to investigate the relative importance of abiotic as it lacks many representatives from Central Asia, and it does not and biotic factors on the diversification of IT taxa, and allow testing include a diversification analysis with explicit temporal and spatial hypotheses about macroevolutionary and biogeographic dynamics of frameworks. Therefore, the evolutionary history of the group is still to be the IT flora proposed by previous authors (reviewed in Manafzadeh unraveled. et al., 2017). Two main scenarios, not mutually exclusive, may explain In sum, our study aims to decipher the diversification history of IT diversification: (1) vicariance through uplifts of mountain ranges Jurinea with the central goal to fill the knowledge gap concerning the (geographic speciation); and (2) ecological speciation. In the first sce­ evolution of the IT flora, investigating the three classical key issues nario, species diverged by allopatry driven by geologic or climatic addressed in radiation studies: tempo, mode, and diversification drivers. 2 S. Herrando-Moraira and C. Roquet Molecular Phylogenetics and Evolution 189 (2023) 107928 We did so by integrating phylogenetic data with biogeographic history initial supermatrix (159,512 bp), the same unsupported nodes remained inference, diversification modelling, ancestral state reconstruction, and (see comparison in Supplementary Fig. 1). Because of that, we used the ecological niche modelling. Specifically, we first obtained a highly initial unfiltered dataset in downstream analyses. resolved time-calibrated tree based on high throughput molecular methods yielding 502 family-specific loci and including a wide taxon 2.4. Phylogenetic inference and dating analyses sampling (77%). Then, we conducted an analysis of biogeographic his­ tory inference to elucidate the relative contributions of in situ speciation Phylogenetic inference of Jurinea was conducted twice, using versus vicariance to the IT flora, and to explore the potential impact of concatenation and coalescence methods. For the first method, we per­ range expansions and contractions in the diversification dynamics of formed a maximum likelihood (ML) analysis based on the supermatrix Jurinea. Next, we explicitly tested whether species diversification varied using RAxML-NG (Kozlov et al., 2019) from the CIPRES Science with time, climate or species-diversity; and whether the major climatic Gateway v.3.1 (Miller et al., 2010). Specifically, we ran a search for the shifts that occurred during the PPT significantly impacted diversifica­ best-scoring ML tree and a slow bootstrap (BS) analysis, with ten ran­ tion dynamics. Furthermore, we investigated the relative importance of domized and parsimony starting trees and setting each locus as a niche conservatism vs. niche divergence in speciation by measuring the different partition under the GTR + G nucleotide model. For the BS phylogenetic signal of the environmental niche and growth form lability analysis, we applied the auto-stopping criterion “autoMRE”. Branches (i.e. the propensity of lineages to change growth form). Last, we iden­ with BS values >70% were considered as statistically supported (Hillis tified stable climatic areas since the Pliocene to the present to assess how and Bull, 1993). The resulting tree was visualized with FigTree v.1.4.3 regional climatic fluctuations may have impacted Jurinea. (https://tree.bio.ed.ac.uk/software/figtree/). The coalescence infer­ ence was performed with ASTRAL-III v.5.5.3 (Zhang et al., 2018) based 2. Materials and methods on the individual gene trees for each locus previously obtained with a RAxML analysis applying the same settings as for the supermatrix, but Methods are described in greater detail with full references in the with 200 resampling replicates. The local posterior probability (LPP) Supplementary Information. was used as a branch support metric, considering strongly supported branches with a LPP > 0.95 (Sayyari and Mirarab, 2016). 2.1. Sampling strategy Due to the large size of the dataset, we used the penalized-likelihood approach implemented in treePL (Smith and O’Meara, 2012) to time- To reconstruct the evolutionary history of Jurinea, we sampled 187 calibrate the best-scoring ML tree obtained with the concatenation species representing ca. 77% of the genus diversity covering all distri­ approach. We applied the same procedure and calibration points as in bution areas (see Appendix A for a complete list of species). Plant ma­ Herrando-Moraira et al. (2019a). As a brief description, we obtained the terial was extracted from dried leaves of herbarium collections or field confidence intervals (95% CI) in the estimated node ages running 100 expeditions. Based on previous phylogenetic studies (Susanna et al., independent treePL analyses with calibration points constrained to a set 2006; Barres et al., 2013; Herrando-Moraira et al., 2019a), we added 58 of random values generated under normal and lognormal distributions, outgroup species to estimate the time divergence of Jurinea lineages. In depending on each point (see Herrando-Moraira et al., 2019a). The total, 245 species were analyzed, from which 131 have been sequenced resultant time-calibrated 100 trees were introduced in TreeAnnotator for the first time in this study and 114 taken from previous ones (Her­ v.1.7.5 (Drummond et al., 2012) to obtain the maximum clade credi­ rando-Moraira et al., 2018, 2019a, 2020; Jones et al., 2019;; Mandel bility (MCC) tree chronogram with median node heights and corre­ et al., 2014, 2017, 2019). sponding CIs. This final dated MCC tree is provided in Supplementary Fig. 2. 2.2. Laboratory workflow 2.5. Ancestral area reconstructions We followed the library prep protocol and sequence capture work­ flow detailed in Herrando-Moraira et al. (2018, 2019a, 2020) based on To infer the most probable ancestral ranges of Jurinea lineages, we Compositae COS 1061 loci kit (Mandel et al., 2014) except for the performed a biogeographical analysis using the software RASP 4.02 (Yu pooling step: here, groups up to eight libraries were arranged (500 ng in et al., 2015), which implements the R package BioGeoBEARS (Matzke, total per pool). The final enriched libraries were sequenced using 100 bp 2013). As input we used the chronogram tree generated from treePL in paired-end reads in the DNA Sequencing Core CGRC/ICBR of the Uni­ which all species other than Jurinea were pruned. We defined nine versity of Florida using an Illumina HiSeq 3000 (Illumina, USA) or in biogeographic areas (see Fig. 1; Supplementary Table 1) according to Macrogen Co. (Seoul, South Korea) in an Illumina HiSeq 4000 (Illumina, data on floristic regions (Takhtajan, 1986), biogeographic syntheses of USA). Raw DNA sequence reads were deposited in the NCBI Short Read the IT (White and Léonard, 1991; Djamali et al., 2012; Manafzadeh Archive database (SRA, https://www.ncbi.nlm.nih.gov/sra) under the et al., 2017), and global climate delineations (Köppen and Geiger, 1936; BioProject accession number PRJNA1013784. Beck et al., 2018). The maximum number of ancestral areas was set to nine. 2.3. Target sequence extraction A first model comparison step was run for the following models: Dispersal–Extinction–Cladogenesis (DEC; Ree et al., 2005), and the A first read quality assessment was conducted with FastQC v.0.10.1 likelihood versions of Dispersal–Vicariance (DIVA-like), and BayArea (https://www.bioinformatics.babraham.ac.uk/projects/fastqc/). Trim­ (BAYAREALIKE). The addition of the parameter jump dispersal or momatic v.0.36 (Bolger et al., 2014) was used to remove the adapters founder event speciation “J” (Matzke, 2013) to the models was dis­ and clean the reads, and HybPiper v.1.3.1 (Johnson et al., 2016) to carded due to the recently reported possible statistical and conceptual extract the target loci. We then removed in downstream analyses all loci problems derived from its inclusion (Ree and Sanmartín, 2018). The flagged as potentially paralogous by HybPiper. Finally, MAFFT v.7.266 DEC model resulted as the best-fit one according to AIC values (Sup­ (Katoh and Standley, 2013), trimAl v.14 (Capella-Gutiérrez et al., 2009), plementary Table 2), and thus it was used to compute the final analysis. and FASconCAT-G v.1.02 (Kück and Longo, 2014) were used sequen­ All models were run under a non-time-stratified analysis and without tially to obtain the loci-separated alignments and the concatenated dispersal multipliers matrix. supermatrix. In order to reduce phylogenetic noise, we excluded posi­ tions with high substitution rates, following Fragoso-Martínez et al. (2017). However, after removing a total of 6109 positions from the 3 S. Herrando-Moraira and C. Roquet Molecular Phylogenetics and Evolution 189 (2023) 107928 Fig. 1. Distribution of Jurinea species di­ versity (taxa included in the study) across the nine biogeographic areas (justification in supplementary text and descriptions in Supplementary Table S1) and buffer zones, which represent areas around 20 km of each occurrence record (2691 points; Ap­ pendix B). Phylogenetic reconstruction of Jurinea diversification inferred from the supermatrix of 502 COS loci by maximum likelihood analysis (in RaxML), time- calibrated (in treePL), and reconstructed ancestral biogeographic areas (in RASP, BioGeoBEARS). Colours of tree nodes indicate the most likely estimation of ancestral areas. Tips are also coloured ac­ cording to the current geographic distri­ bution of each species. Abbreviations of areas for inferred ancestral states are only specified for the main tree backbone nodes, for information about the rest of nodes and all colour correspondence see Supplemen­ tary Table S6 and Supplementary Figs. S10 and S11. Gray bars on nodes show the 95% of confidence intervals (CI). The graphs below phylogeny show: the upper one, a lineage-through-time plot (LTT) that in­ dicates the number of lineages emerged (log-transformed) in a temporal line. The black line represents the maximum clade credibility (MCC) tree chronogram and the grey ones correspond to 100 time- calibrated trees. The lower graphic shows the evolution of global mean temperature (◦ C) since Miocene to present. Abbrevia­ tions of areas: MED = Mediterranean; CB = Circumboreal; IT = Irano-Turanian; HIM = Himalayas; EA = Eastern Asia. Abbreviations of time stages: C = Chiba­ nian; U = “Upper”. 4 S. Herrando-Moraira and C. Roquet Molecular Phylogenetics and Evolution 189 (2023) 107928 2.6. Ancestral state reconstruction of growth form obtained after dropping the EA subclade plus J. mesopotamica (hereafter we refer to this subset tree as CB + CA), with the main goal to further Reconstruction of ancestral growth form was performed on the investigate the drivers of diversification of the two main subclades that Jurinea chronogram. We classified Jurinea species into five growth form constitute the crown radiation of Jurinea (Circumboreal and Central- categories based on our own observations and descriptions from local Asian ones). floras: (A) annuals with dwarf semiscapose basal rosette; (B) biennials presenting leafy stems with soft and wide leaves; (C) acaulescent pe­ 2.8. Environmental niche: Data assembly and analyses rennials; (D) caulescent perennials formed by basal rosettes with leafless flowering stems; and (E) caulescent perennials with leafy flowering To explore the impact of environmental factors in the evolution of stems (see dichotomic diagram in Fig. 2). Inferences about growth form main Jurinea lineages, we conducted comparative niche analyses. In evolution were done using Markov models (Mk) in the R package total, we gathered 2691 records for 158 species, representing 84.5% of diversitree (FitzJohn, 2012). Four possible discrete Mk models were our taxon sampling (see Appendix B for the complete list of occurrences fitted and compared using AIC values, a null model with different and sources). The gathered localities were extracted from several sour­ transition rates (ARD) plus three constrained models: a symmetrical ces: GBIF portal (https://www.gbif.org/), research papers, online floras model in which transition rate between any two states do not differ and atlases, herbarium collections, and virtual herbaria. As environ­ (SYM), a symmetrical model with equal rates for the three types of mental predictors, we collected an initial set of 45 variables (at 2.5 arc perennial lifeforms (SYM-ERp), and a model with all rates equal (ER). min or coarser resolutions) related with climate, topography, vegeta­ tion, soil properties, habitat heterogeneity, and solar radiation (Sup­ 2.7. Diversification analyses plementary Table 3). We extracted the values of variables for each locality with ArcGIS v.10.2.2 (Esri, Redlands, California, USA 2014). We As a first exploratory analysis, we outlined the accumulation of lin­ reduced the initial 45 variables to a set of 21 uncorrelated ones based on eages over time with a lineage-through-time plot (LTT; Nee et al., 1992) the results of a Pearson correlation analysis (Supplementary Fig. 3) and a with the R package ape and function mltt.plot (Paradis and Schliep, 2019) principal component analysis (PCA) (Supplementary Table 4) to avoid setting in log scale axis “y” (number of lineages). As input, we used the model overfitting. See details of variable selection in Supplementary maximum clade credibility tree chronogram and the 100 independent Information. Two datasets were tested, one with all occurrences and dated trees extracted from treePL, which were used to account for the values, and another with the mean values for each species. Given that uncertainty of the node age estimations (see details in Herrando-Moraira similar results were obtained (Supplementary Fig. 4C and 4D), we used et al., 2019a). the mean values as a summary of the species average potential niche. We investigated diversification dynamics of Jurinea following a We identified the variables responsible for niche variation of the hypothesis-driven approach to avoid identifiability problems (Louca and main phylogenetic lineages of Jurinea through a standard PCA analysis Pennell, 2020; Morlon et al. 2022), i.e. we defined a set of diversification conducted in R v.3.6.3 (R Development Core Team, 2019) using Facto­ models to be fitted and compared in order to test explicit evolutionary MineR, factoextra (for PCAs) and ggplot2. To test the hypothesis of niche hypotheses. Specifically, we tested whether species diversification rates conservatism, expansion, and divergence among lineages, we also con­ varied (linearly or exponentially) with time, climate, or species- ducted the PCA-env analysis designed by Broennimann et al. (2012) diversity, and included as null models constant-rate ones. Diversity- with the R script reported in Herrando-Moraira et al. (2019b). dependent models include an additional parameter: the clade carrying We reconstructed the most probable ancestral niche states with the R capacity (K), which corresponds to the number of species that can exist function fastAnc from phytools package (Revell, 2012). As input vari­ at the same time (Etienne et al., 2012). Concerning diversity-dependent ables for each species, we used the PC1 and PC2 scores from the standard models, we fitted three main types of models: (1) models in which PCA analysis, and the most variable environmental factors at ancestral speciation rate was dependent on species-diversity, and a single clade nodes or at interspecific level. The changes in environmental values carrying capacity was estimated for the entire study-group; (2) models were visualized with the contMap function (also from phytools). in which speciation rate also varied with species-diversity, but all esti­ mated parameters were allowed to shift at 3 Ma, with the aim to test for 2.9. Phylogenetic signal in phenotype and environmental variables the potential effect of the Pliocene-Pleistocene Transition (PPT) on the net diversification rate and/or the clade carrying capacity; and (3) We measured the phylogenetic signal of the categorical variable models in which we allow one of the three main subclades to undergo its “growth form” with the δ statistic (Borges et al., 2019) computed in R, own diversification dynamics by enabling a decoupling of diversifica­ and of the quantitative variables “altitude” and “environmental niche” tion rates and clade carrying capacities (Etienne et al., 2012), i.e., the (PC1 and PC2 values) with Blomberg’s K (Blomberg et al., 2003) estimated parameters (speciation and extinction rates plus the clade computed in RASP. These analyses were performed to discern between carrying capacity) were separately estimated for the main clade and the non-adaptive radiation (high phylogenetic signal expected, i.e. closely decoupled subclade. The latter type of models was fitted to account for related species present similar phenotype and niche traits), and adaptive the fact that Jurinea is composed of three main subclades with different radiation (low phylogenetic signal expected, i.e. close relatives exhibit ages and geographical contexts (and thus each subclade may be at a notable differences in trait values). different stage of diversification). We only fitted models allowing vari­ ations in speciation rate to avoid potential flaws as indicated by Burin 2.10. Climatic stability inference et al. (2019) on models assuming variable extinction rates. To account for incomplete taxon sampling, we applied in each model the analytical We implemented a recently published approach (Herrando-Moraira correction corresponding to the sampling fraction. The best fitting- et al., 2022) to assess the impact of climatic fluctuations on Jurinea model was selected based on AICc. All analyses were run with R pack­ diversification and to identify stable refugia or unstable climatic areas ages RPANDA (Morlon et al., 2016) and DDD (Etienne et al., 2012). where the species occur (see schematic overview in Supplementary Given that the best-fitting model according to AICc was a diversity- Fig. 5). In this approach, a Climate Stability Index (CSI) is estimated at a dependent model in which the diversification dynamics of the Eastern- spatial resolution of 2.5 arc min for the period comprised between Asian (EA) subclade were decoupled from the main clade, we per­ middle-late Pliocene (3.3 Ma) and the present. To do this, CSI uses the formed a second round of diversification analyses: we fitted the whole reconstructed paleoclimatic data of the 12 time periods (Supplementary set of models (except those models in which the diversification param­ Fig. 5A) included in the PaleoClim database (Brown et al., 2018) rep­ eters were decoupled for a particular subclade) to the phylogenetic tree resenting both warm and cold cycles. The CSI values for each 2.5 arc min 5 S. Herrando-Moraira and C. Roquet Molecular Phylogenetics and Evolution 189 (2023) 107928 Fig. 2. Results of analyses on growth form evolution. Dichotomic scheme of morphological classification of Jurinea according to five main morphological syndromes (coded as A, B, C, D, and E), followed by an example species of each category. For each morphotype, boxplots representing variance of several climate variables based on their current distribution areas (aridity index, bio 9, and bio 15), altitude, and beanplot showing the magnitude of paleoclimatic variation (climatic stability index). Phylogenetic reconstruction of Jurinea inferred from the supermatrix of 502 COS loci by maximum likelihood analysis (with RaxML), time-calibrated (with treePL), and reconstructed ancestral morphological syndromes (with diversitree). Colours of tree nodes indicate the most likely estimation of ancestral morphologies. Tips are also coloured according to the current morphotype of each species. Abbreviations of growth forms for inferred ancestral states are only specified for the main tree backbone nodes. Abbreviations of time stages: C = Chibanian; U = “Upper”. 6 S. Herrando-Moraira and C. Roquet Molecular Phylogenetics and Evolution 189 (2023) 107928 cell are based on estimating the amount of variation or dispersion along harboured species with four of the five distinct growth forms (all except time series by using the standard deviation (Supplementary Figs. 6 and annuals) (Fig. 2). Conversely, the CB clade was the most uniform, with 7; see more details in Herrando-Moraira et al., 2022). The procedure is type “D” (perennials with leaf rosettes and leafless flowering stems) schematized in Supplementary Fig. 5C. For each species occurrence re­ predominating and with few transitions to types “C” (acaulescent pe­ cord, we extracted its corresponding climatic stability index score with rennials with leaf rosettes) and “E” (subshrubs with leafy flowering the ArcGIS tool “Extract Multi Values to Points”. Then, the average was stems). The CA clade showed a wide heterogeneity of growth forms calculated for the species with more than one occurrence to obtain a (syndromes “C”, “D”, and “E”, being the most abundant the last two), final dataset with a unique index value per species. Climatic stability besides exhibiting large disparities among closely related species. data plots were generated with boxplot from graphics R package and bean plots from portal https://shiny.chemgrid.org/boxplotr/, a web-based R 3.4. Diversification analyses application. The R scripts used for the paleoclimatic analyses are included in Appendix C. The maximum slope in LTT plot occurred around 3–1.5 Ma (Fig. 1). We compared the fit of 30 models (Supplementary Table 8) and, ac­ 3. Results cording to AICc, the best-fitting model was the diversity-dependent one in which the diversification dynamics of the subclade EA was decoupled 3.1. Phylogenetic analyses from those of the main clade. This model indicated that the speciation rate of EA was nearly half of the main clade (λEA = 0.85; λmain clade = From the total 1061 target loci, 502 were finally retained after dis­ 1.84) whereas the extinction rate was one order of magnitude higher in carding potential paralogs. The supermatrix resulted in an alignment of EA (μEA = 0.00023; μmain clade = 0.00001). Following this model, EA had 161,485 bp. Phylogenetic reconstruction using concatenation and reached its clade carrying capacity (KEA = 25, which is indeed the coalescent-based approaches returned similar topologies at deep tree estimated current species diversity), whereas the main clade had not yet nodes with strong branch supports (Supplementary Figs. 8–9), except for reached it (Kmain clade = 342, while Jurinea diversity estimations are of the position of J. mesopotamica Hand.-Mazz. Three main clades were ca. 200 species). recovered, which mostly corresponded to the three main geographic Concerning the diversification models fitted for the CB + CA subset groups (Fig. 1). Hereafter, we refer to these clades as “Eastern Asia (Supplementary Table 9), the best-fitting model according to AICc was clade” (EA clade), “Circumboreal clade” (CB clade), and “Central Asia the diversity-dependent model with a major shift of diversification dy­ clade” (CA clade). namics at the Pliocene-Pleistocene transition (PPT, ca. 3 Ma) and null extinction. This model indicated that the speciation rate nearly doubled 3.2. Divergence times of Jurinea and biogeographic history inference after the shift (λ1 = 0.93, λ2 = 1.77). The estimated values of clade carrying capacity (K) before and after the shift (K1 = 6163 and K2 = 331, All numerical results from divergence time and ancestral area esti­ respectively) suggested that diversity-dependence developed mainly mates are compiled in Supplementary Material: node age values (Sup­ after PPT. The ΔAIC between the best and second-best model was above plementary Table 5), corresponding node identifiers (hereafter IDs; the threshold of 2, which is typically chosen to discern models with Supplementary Fig. 10), probabilities of ancestral areas for each node confidence for both datasets (all Jurinea and CB + CA subset). (Supplementary Table 6; Supplementary Fig. 11) and inferred areas appended to nodes on the phylogenetic tree (Supplementary Fig. 11). 3.5. Phylogenetic patterns in environmental variables The estimated stem and crown ages for Jurinea were 12.8 Ma (11.5–14.5 95% CI) and 10.7 Ma (9.5–12.0), respectively. The most Environmental variables that mainly contributed to the separation of probable ancestral area for the root node of Jurinea was IT1–IT2 region the phylogenetic clades were aridity (12.8% of contribution to PC1 axis (Mesopotamia-Iranian Plateau; Fig. 1). The first main divergence event of standard PCA) and altitude (19.2% in PC2 axis; Fig. 3A–3B; Supple­ occurred at 8.9 Ma (7.8–10.0 95% CI) when the EA segregated. In the mentary Table 10). In the present study case, altitude was highly Messinian period, the EA clade started to diverge at 5.9 Ma (5.2–6.5 95% correlated (Pearson’s r > 0.4) with soil quality and other topographic CI) from an IT2 ancestor, giving rise mainly to Himalayan and East Asian variables (Supplementary Fig. 2). Standard PCA (Fig. 3B) and PCA-env lineages. Subsequent tree backbone nodes showed that the most likely (Supplementary Fig. 3B) analyses revealed that species environmental ancestral area was IT2, pointing out that Iranian Plateau was the focal niches of the three main clades (EA, CB, and CA) overlap considerably zone of subsequent dispersals to adjacent regions (Fig. 1). In the (D = ca. 30%; Supplementary Table 11). The clade EA showed the Miocene-Pliocene transition a second main dispersal occurred to the CB highest niche breadth and turned out to be equivalent to niches found in region, which started to diverge around 4.5 Ma (5.1–4.0 95% CI) in a CB and CA (p-value < 0.05; Supplementary Table 11). Comparing the IT2–CB transition, and spread later towards the following regions: IT2, three main clades, the CB showed high niche expansion values with IT3, IT4, IT5, CB, and MED (Fig. 1). Similarly, in the Miocene-Pliocene respect to the other two clades (33% with respect to EA and 41% with transition the other main range expansion took place, dated to around respect to CA; Supplementary Table 11). 3.8 Ma (3.4–4.6 95% CI) with an ancestral area in IT2–IT4, triggering Ancestral reconstruction of the first two PC axes (Fig. 3C) showed the origin of the CA clade, which contains lineages that colonized all IT that CB clade was the one including species with most diverging niches. regions, CB, and MED (Fig. 1). Our results indicated that the majority of Jurinea species are present in territories with an aridity index < 0.5, which corresponds to drylands 3.3. Growth form evolution according to the criterion of Middleton & Thomas (1997); the species found in humid zones (aridity index > 0.65) were mostly clustered in the According to AIC values, the best fitting model of growth form CB and EA clades (Fig. 4). All the species from Caucasus (which evolution was SYM, followed by SYM-Erp (ΔAIC = 0.47; Supplementary constitute a single subclade of 14 species within the CB clade) were Table 7). According to SYM, transitions were most frequent between found in humid or sub-humid areas (aridity index ≥ 0.5; Fig. 4). different types of perennials (q54 = 0.104, q43 = 0.047, q53 = 0.037) Regarding altitude, the majority of species (57%) of the CB clade are than between perennials and biennials/annuals (q values < 0.02). found at low elevation (<1000 m). In contrast, EA and CA clades present Congruently, according to the second best model, transitions between a much lower proportion of low altitude species (14 and 33%, perennial growth forms were also the most frequent ones (qperennial = respectively). 0.073). Indeed, most of Jurinea species show perennial biotypes. The Statistically significant levels of phylogenetic signal (Blomberg’s K most heterogeneous group in growth forms was the EA clade, which > 1; p-value = 0.001) were recovered only in CB clade for environmental 7 S. Herrando-Moraira and C. Roquet Molecular Phylogenetics and Evolution 189 (2023) 107928 Fig. 3. (A) Classification of clades (see Fig. 1). (B) Graphical output of Principal Component Analysis (PCA) analysis of the values of 21 climatic variables for each species, representing the global environmental space and the distribution of each phylogenetic group. (C) Ancestral state reconstructions (obtained with the fasAnc function of R package phytools) of PC1 and PC2 scores for each species, representing the ecological niche evolution along the phylogenetic history of Jurinea. Loading values for each species of both PCs are represented in vertical barplots besides the time-calibrated phylogeny. niche and altitude. The δ-value related to growth form conservatism was Asia. As a general pattern, climatic fluctuations affected the northern­ also three-fold higher in CB than for the other clades (Table 1). most areas (especially for temperature oscillations) and East Asia (especially for precipitation; Fig. 6). Accordingly, the CB clade, which has more northerly distributed species, showed higher median and range 3.6. Paleoclimatic oscillations and long-term stable regions values (Q1–Q3) of climatic stability index (Fig. 7A; Table 1), i.e. species occur in regions that had more unstable climate. The top species rich- Paleoclimatic data analyses identified both climatic stable refugia areas (IT2, IT4, CB) showed medium–high values of climatic oscilla­ and highly variable regions along IT and surrounding territories tions (Fig. 7B). (Figs. 5–6). We detected local stable climatic conditions since middle- late Pliocene to present in the Mediterranean region, Anatolia, Meso­ potamia, Iranian Plateau, central-eastern Himalayas, and south-eastern 8 S. Herrando-Moraira and C. Roquet Molecular Phylogenetics and Evolution 189 (2023) 107928 Fig. 4. Estimation of ancestral states reconstructed with the fasAnc function of R package phytools for the most variable environmental traits across the Jurinea phylogeny. Loading values for each species are represented in vertical barplots besides the time-calibrated phylogeny. The asterisk indicates the clade containing species exclusively found in the Caucasus. 4. Discussion 4.1. Cold-arid climatic conditions as triggers of Jurinea radiation The present study offers new evidence on the tempo, mode, and The main evolutionary events of Jurinea were probably substantially drivers of evolution of the IT flora, which is one of the least studied influenced by climate-related turnovers, especially those that led to Eurasian floras despite its high diversity and rate of endemism. The colder and drier conditions: (1) the origin of the genus (10.7 Ma; Fig. 1) integrative approach applied here, combining data from multiple sour­ occurred during a major global climate change, turning from warm- ces (phylogeny, biogeography, geology, ecology, morphology, and tropical and subtropical to cool-temperate conditions at the onset of paleoclimate), revealed the importance of climate shifts and topo­ global cooling and marked seasonality during mid-late Miocene (Potter graphic heterogeneity as triggers or modulators of species radiations. and Szatmari, 2009; Pound et al., 2012; Herbert et al., 2016); (2) long- Importantly, the paleoclimatic analyses highlighted that either long- distance dispersals from IT to other floristic regions (~7–4 Ma; Fig. 1) term climatic stability or moderate climate oscillations could explain occurred around the Miocene-Pliocene boundary, a period characterized high species richness in two biogeographic regions depending on their by global biotic turnovers (LaJeunesse, 2005), such as the expansion of general topography. In the Iranian Plateau, the high climatic stability C4 grasslands (Shen et al., 2018) in response to a major cooling trend detected would explain its high richness of xerophytic elements and its between 7.0 and 5.5 Ma (Holbourn et al., 2018) and the aridification of role as a museum and a cradle for the IT flora. Contrarily, regional cli­ Central Asia (Miao et al., 2012); and (3) diversification dynamics of the matic fluctuations in the Circumboreal region and mountains of Central main lineages of Jurinea (i.e., the CB and CA clades) suffered a signifi­ Asia fostered allopatric diversification, especially in lineages adapted to cant diversification shift during the Pliocene-Pleistocene Transition (ca. semi-arid and cold areas. 3 Ma; Fig. 1), in which the net diversification rate almost doubled, from λ = 0.93 to λ = 1.77. During this later period, climate turned colder and drier at the onset of the Northern Hemisphere glaciations when large continental ice sheets appeared, resulting in important environmental 9 S. Herrando-Moraira and C. Roquet Molecular Phylogenetics and Evolution 189 (2023) 107928 Fig. 5. Mapped values of climatic stability index (CSI) for the study area from Pliocene (3.3 Ma) to present, at a resolution of 5 × 5 km. Colours range from blue for low standard deviation (SD) values, which repre­ sent areas with low climatic fluctuations from Plio­ cene to present, to red for high SD values, which show areas where high climatic fluctuations probably took place. Main geographic locations (e.g. mountains, deserts) are indicated in the zoomed maps. This index is based on data from six variables (bio1 = annual mean temperature; bio4 = temperature seasonality; bio9 = mean temperature of driest quarter; bio12 = annual precipitation; bio15 = precipitation season­ ality; and bio17 = precipitation of driest quarter). (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.) changes (Willeit et al., 2015). Lineages that could not adapt to novel 4.2. Insights into evolutionary and biogeographic hypotheses on the Irano- environments either dispersed to suitable habitats or faced extinction. Turanian flora Jurinea seems to have taken evolutionary advantage of these climate shifts towards dry-cold regimes. The diversification increase detected for In relation to the temporal origin of typical xerophytic species within the last 3 million years could be related to geographic range modifica­ the IT, we found that Jurinea originated in the middle-late Miocene, ca. tions or to adaptation to these new ecological opportunities and niches 10.7 Ma (9.5–12.0; Fig. 1), a more recent origin than that suggested for deserted by other plant groups unable to adapt to plummeting temper­ the IT elements (Manafzadeh et al., 2017). Spatially, Jurinea started to ature and precipitation. Indeed, during this period Jurinea lineages also diverge in Western Asia around Mesopotamia and the Iranian Plateau showed an increasing trend of form lability (Fig. 2) and ecological niche (region IT1–IT2 here; Fig. 1). Within IT, Western Asia is a centre of shifts (Figs. 3–4), probably as an adaptation to novel habitats and new origin for numerous plant radiations that also occurred during middle- geographic areas. Overall, Jurinea is an excellent example of the late Miocene (Barres et al., 2013; Karl and Koch, 2013; Lauterbach extensive impact of gradual global cooling trend over the past 15 Ma on et al., 2019; Peterson et al., 2019). The definitive closure of the Meso­ diversification history and expansion of the temperate biota of the potamian Seaway in the middle Miocene at 13.8 Ma (Bialik et al., 2019) Northern Hemisphere (Folk et al. 2019, 2020 and references therein). could have promoted this deep-lineage divergence, as new lands to be colonized were available (Potter and Szatmari, 2009; Herbert et al., 2016). The colonization of adjacent regions by IT elements could be the result of independent dispersion waves. In Jurinea, the first colonized 10 S. Herrando-Moraira and C. Roquet Molecular Phylogenetics and Evolution 189 (2023) 107928 Table 1 Summary of main results of five key questions in Jurinea radiation. Species formation and growth form transition processes were calculated according to retention (same inferred area or habit) or shift in internal nodes and tips respect their respective MRCA node. Abbreviations: BS = bootstrap; E-C = expansion–contraction model; IT = Irano-Turanian; LPP = local posterior probability; MRCA = most recent common ancestor; p% = probability in percentage; Phy-signal = phylogenetic signal; spp = species; SS = strong significant phylogenetic signal. Specific attributes or Jurinea Eastern Asia (EA clade) Circumboreal (CB clade) Central Asia (CA clade) measurable traits (whole tree) Number of species 187 19 65 89 Estimated sampling 77% 68% 74% 67% fraction Clade support (BS/ 100/1 100/1 100/1 100/1 LPP) Divergence time Late Miocene Late Miocene Pliocene Pliocene (crown age) (10.7 Ma) (5.9 Ma) (4.5 Ma) (3.8 Ma) Ancestral areas (p%) IT1–IT2 (64%) IT2–HIM (81%) IT2–CB (55%) IT2–IT2 (100%) Species formation in situ (60%) colonization (76%) in situ (70%) in situ (57%) process 6 6 6 6 Geographic scale (spp/ 5.9 × 10− 1.9 × 10− 3.3 × 10− 11.1 × 10− km2) Growth form conservation (81%) conservation (64%) conservation (90%) conservation (78%) conservation or transition Phy-signal growth 6.21 1.33 18.89 SS 4.15 form: δ Phy-signal niche PC1: 0.58 0.78 1.29 SS 0.80 K Phy-signal niche PC2: 0.52 0.96 1.19 SS 0.74 K Phy-signal altitude: K 0.53 0.62 1.22 SS 0.75 Climatic stability 16.05(14.70–18.68) 16.20(15.60–17.35) 17.8(14.63–20.70) 15.9(14.8–18.0) index:median (Q1- Q3) Radiation type Mixed model Adaptive Non-adaptive Adaptive proposed Diversification drivers - Climate changes - Niche partitioning on ecologically - Range and niche shift towards - Climate fluctuations and orogeny in Pamir- and history - Pre-adaptation to constrained adaptive landscapes- Low humid conditions Tian Shan mountains- Growth form and aridity- Geographic evolutionary potential in mesic habitats - Niche conservatism- Climatic niche lability fostered local adaptation speciation oscillations stimulated (longitudinal, allopatric speciation latitudinal, altitudinal) region out of IT was East Asia (late Miocene), the second the Circum­ favoured the persistence of Jurinea species as well as other xerophytic IT boreal region (Pliocene), and the last one the Mediterranean region elements. (Pleistocene), which was colonized repeatedly by different clades Takhtajan’s (1986) old hypothesis is once again verified here: the IT (Fig. 1). Among these regions, the genus notably diversified only in the flora originated along the Iranian Plateau and served as biodiversity Circumboreal region, whereas in Eastern Asia the net diversification rate source for the adjacent regions. For some millions of years during late was significantly lower than in other regions. Internal movements and Miocene (ca. 10.7–6 Ma) Jurinea probably remained within the Iranian range expansions between IT areas (IT1-IT5) occurred during the Plateau (Fig. 1). Its diversification was constant during the so-called Pliocene-Pleistocene boundary. period of the “Pikermian chronobiome”, when arid-seasonal climates and open environments dominated in Eurasia (Hurka et al., 2019). The 4.2.1. The Iranian Plateau as both a museum and a cradle for plant main transregional dispersal events took place around 6–4 Ma during diversity the Miocene-Pliocene transition, when the genus successfully colonized Our results suggest a double role of the Iranian Plateau as a three geographic areas: East Asia, the Circumboreal region, and Central “museum” (old-lineage persistence) and a “cradle” centre (in situ Asia. This period is coincident with an intense phase of uplift and speciation; cf. Moreau and Bell, 2013) for Jurinea, which has also been deformation along the Iranian mountain ranges: Zagros, Alborz, and documented for other IT radiations such as Acantholimon (Moharrek Kopeh Dagh (Mouthereau, 2011), which provided isolated geographic et al., 2019). For instance, most of the deepest nodes of Jurinea were areas that promoted allopatric speciation, as well as new habitats that distributed in this area, according to our biogeographic estimations; in created opportunities for ecological speciation. It is also coincident with contrast, some recent nodes gave rise to several endemics (Fig. 1). The the onset of the Central Asia aridification (7–5.3 Ma; Sun et al., 2015) dual role of this region as a museum and a cradle is likely linked to one of and a phase of intense cooling in Northern Hemisphere (the “Late our main findings: here we demonstrated for the first time that the Miocene Cooling”; 7–5.4 Ma; Zachos et al., 2001; Herbert et al., 2016). Iranian Plateau is a region characterized by long-term high climatic These climatic changes likely led to a decrease in native competitors and stability (Figs. 5 and 7B), one of the main causes determining high herbivores as shown by the regional extinction events of mammalian species richness and endemism of world’s biodiversity hotspots (Harri­ biochronological indicators (Casanovas-Vilar et al., 2010, 2011; Ge son and Noss, 2017). Specifically, our analysis of paleoclimatic data et al., 2012). Consequently, new lineages as Jurinea, other xerophytic revealed a low temporal variation in IT2 for annual precipitation and, plants and C4 grasses (Huang et al., 2007; Arakaki et al., 2011) could particularly, precipitation of driest quarter (bio 12 and bio17, respec­ have spread as a result of antagonists decrease (Yoder et al., 2010) and tively; Fig. 6; Supplementary Fig. 6); thus, the precipitation regime their pre-adapted niches to arid habitats. Here, we confirmed the pre- remained rather constant since the Pliocene, without abrupt regional adaptation hypothesis based on the phylogenetically conserved values changes even during the periods of strongest climatic changes such as obtained for aridity (Fig. 4), which may have resulted in higher rates of the LGM (time “5” in Fig. 6; Supplementary Fig. 5). This factor probably successful establishment in new available areas by Jurinea lineages 11 S. Herrando-Moraira and C. Roquet Molecular Phylogenetics and Evolution 189 (2023) 107928 Fig. 6. Representation of climatic variable oscillations through the 12 time intervals (see Supplementary Fig. S5) of PaleoClim dataset (Brown et al. 2018), as implemented in Herrando-Moraira et al. (2022). On the left, we show values for temperature variables (bio1, bio4, bio9) and precipitation variables (bio12, bio15, bio17) for each of the nine biogeographic areas defined. On the right, for each variable and biogeographic area, the standard deviation (SD) along the 12 time periods considered is mapped. The upper panel shows the evolution of global temperature since the Pliocene compared to the average temperature during the period 1960–1990. 12 S. Herrando-Moraira and C. Roquet Molecular Phylogenetics and Evolution 189 (2023) 107928 Fig. 7. Plots showing range of climatic stability index values for (A) the three main clades and (B) the nine biogeographic areas (see Fig. 1) based on mean index values for each species. The graphical representation type is a beanplot in (A), in which white lines represent species values, and a boxplot in (B) with a superposed line showing the species richness for each region. (Peterson, 2003), as documented in alpine radiations by cold-adapted stable area (Harrison and Noss, 2017). Reasons are probably two-fold: lowland clades (Uribe-Convers and Tank, 2015; Ye et al., 2019). on the one hand, species of Jurinea reached the Mediterranean only in the Pleistocene, a too short time for a radiation to occur. On the other 4.2.2. Successful dispersal pathways and geographic dead ends hand, Mediterranean vegetation would have been well-established at The first dispersal of Jurinea out of IT occurred at 5.9 Ma and led to the arrival of the Jurinea ancestral taxa and, thus, with a low availability the emergence of EA clade, extending the genus from the Iranian Plateau of empty niches. towards the Himalayas (probably across the mountain ranges connect­ The route towards Central Asia (CA clade) was established during the ing IT2 and HIM), Central Asia Mountains (IT4), and East Asia (sub­ Pliocene (at ca. 3.8 Ma). Members of this clade reached all IT regions, tropical China, Myanmar, Vietnam, and Taiwan). A similar dispersal and even the Mediterranean region at least once (J. cypria Boiss., a near- corridor has been postulated for other IT elements (e.g. Gagea Salisb. endemism of Cyprus, belongs to the CA clade). More than half of Jurinea [Liliaceae], Peterson et al., 2019) but, as in Jurinea, this dispersal event sampled species originated within this successful lineage. Unfortunately, originated a species-poor clade. Later, we discuss hypotheses that could a fine-scale discussion about the biogeographic history of this lineage is explain why some IT elements found a dead end along this East Asia hindered by low support values recovered within CA clade (Supple­ corridor. mentary Fig. 7–8). Further exploration of another type of high- During the Miocene-Pliocene transition, early diverging clades that throughput molecular markers may be necessary to unravel the CA ra­ coexisted in the Iranian Plateau gave rise to CB (at 4.9 Ma) and CA clades diation. COS target loci, which are conserved orthology loci identified (at 3.8 Ma), both of which colonized large parts of Eurasia (Fig. 1). The from expressed sequence tags (Mandel et al., 2014), do not seem the CB clade is comprised of two main subclades that are highly differen­ most adequate markers for disentangling rapid species radiations when tiated in terms of geographical expansion. One subclade contains species incomplete lineage sorting and contradictory signals among coding se­ exclusively found along the Caucasus mountains, which corresponds to a quences occur, as is probably the case for the CA clade. Such pattern of known radiation of ca. 20 endemic species (Szukala et al., 2019). The conflicting gene information has also been suggested for the radiation of other one is much more widespread geographically, as it includes spe­ Andean Lupinus L. by Nevado et al. (2016), who found that over 40% of cies distributed through vast CB territories (i.e. N Anatolia, Balkan genes had undergone adaptive evolution of coding sequences. Peninsula, the river basins of Dnieper, Don, and Volga rivers, the Cas­ pian Depression, and South-Central Russian uplands). Moreover, back- colonizations to IT2 and dispersal towards other IT regions (IT3, IT4 4.3. Parallel and independent diversification scenarios of three Jurinea and IT5) also occurred within this subclade. Two main causes may ac­ lineages count for the poor diversity of these IT lineages: (1) founder effects associated with long distance dispersal, which implies low potential Given that the three main clades (EA, CB, and CA) followed inde­ genetic plasticity to adapt to a new environmental space; or (2) high pendent evolutionary histories and were probably driven by distinct climatic instability of most of the newly colonized areas (especially in extrinsic forces, we separately discuss their diversification dynamics IT3 and IT5 regions; Figs. 5–6) that would have increased local extinc­ below. tion rates. These areas were later recolonized by recently-diverged The CB clade, dominated by circumboreal taxa, contains 65 species, Jurinea lineages, whose recent origin could explain the absence of representing 35% of the total explored richness. It originated with the large in situ radiations. expansion of Jurinea from IT2 to CB region ca. 4.5 Ma, coinciding with a Jurinea also colonized the Mediterranean Basin several times during climate reversal from cool to warm-humid conditions known as the the Pleistocene, coincident with the oldest xeric period of the region and “Pliocene climatic optimum” (Sniderman et al., 2016; Jiménez-Moreno the emergence of Mediterranean-type vegetation (ca. 2.3 Ma, Suc, et al., 2019) that probably facilitated northward expansions. In the first 1984). Such transregional dispersal is responsible for the so-called stages of diversification of the CB clade, ecological speciation may have “Kiermack disjunctions” between western Mediterranean and Central prevailed, given the wide range of ecological niche diversity detected at Asia (Ribera and Blasco-Zumeta, 1998). Contrary to other IT elements the tree backbone in the CB clade (Figs. 3–4). Later, such ecologically that actively diversified in the Mediterranean (Banasiak et al., 2013; divergent lineages further diversified, probably through allopatric Manafzadeh et al., 2014; Peterson et al., 2019), this region did not speciation given the high phylogenetic signal recovered for niche, become a secondary speciation centre for Jurinea (only ca. 6 endemic growth form, and altitude of CB clade (Table 1; Fig. 4). Indeed, the CB species are known), as it would be expected for a long-term climatically clade harbours a high number of species growing in regions that have been climatically unstable from the Pliocene to present (Figs. 5 and 7A). 13 S. Herrando-Moraira and C. Roquet Molecular Phylogenetics and Evolution 189 (2023) 107928 Drastic climatic oscillations in the circumboreal region during the 5. Conclusions Pleistocene probably lead to repeated cycles of range expansion and contraction that fostered further diversification, as shown for other We focused on the radiation of Jurinea as a proxy to understand the floras (e.g. Winkworth et al., 2005). evolutionary assembly of the Irano-Turanian flora. Our results depict a The diversification of the Central Asian clade, the most species-rich dual role of the Iranian Plateau for plant diversity both as a museum and lineage of Jurinea (89 species, ca. 50% of total richness), was probably cradle region, in which long-term climatic stability favoured old-lineage triggered by a combination of climatic and orogenic factors. During the persistence. This region has constituted an important source of xero­ cold-arid period in the Asian interior dated to 3.6 Ma (Rea et al., 1998; phytic lineages for other Eurasian hotspots such as the Mountains of An et al., 2001; Zheng et al., 2004), plant communities experienced Central Asia, where Jurinea highly diversified. Jurinea illustrates that marked turnovers to drought tolerant taxa (An et al., 2005; Koutso­ large plant radiations can be the result of both adaptive and non- dendris et al., 2019). Congruently, we observed niche changes towards adaptive processes (Donoghue and Sanderson, 2015). Indeed, Jurinea adaptation to higher aridity and open habitats in the origin of CA clade is constituted by three main clades whose diversity results from different (Fig. 4), and a growth form reversion to the most ancestral habit (syn­ macroevolutionary scenarios. The Circumboreal lineage is probably the drome “E”, foliose subshrubs; Fig. 2), likely as an adaptation to arid result of a non-adaptive radiation, in which allopatric processes pre­ climates in the Mesopotamian-Iranian region (IT1–IT2; Fig. 1). There­ dominated as a result of high climatic instability of the region during the fore, aridity could have acted as a speciation stimulus to CA radiation Pleistocene. In contrast, the Central Asian and the East Asian clades (Stebbins, 1952). Speciation was also probably fostered by orogenic often show divergent ecological niches and high growth form disparity processes in Central Asia resulting in new habitats. The Pamir-Tian Shan among closely related species. The evolutionary success of the Central convergence began at ca. 3–5 Ma (Fu et al. 2010; Thompson et al. 2015) Asian clade (the most diverse one) is probably both due to geographical when Pamir advanced about 300 km northwards (Burtman and Molnar, isolation as well as new ecological opportunities provided by orogenic 1993; Burtman, 2000), and parallel and subsequent pulses of deforma­ processes in the Pamir-Tian Shan mountains occurred since the Pliocene. tion and rapid uplifts of Pamir-Tian Shan mountain belts occurred Moreover, in this region, predominating cold and arid conditions during (Zhang et al., 2013). We hypothesize that CA ancestors diverged the Pleistocene likely favoured preadapted taxa as is the case of Jurinea. ecologically and occupied multiple habitats provided by mountain Contrarily, the East Asian clade is the most species-poor one, probably as ranges, partitioning their niche along altitudinal gradients differentiated a result of the low evolutionary potential of Jurinea in mesic habitats. In by topography, microclimate, predator pressure, pollinators, soils, or sum, the radiation of Jurinea resulted from the combined action of cli­ vegetation structure (Karl and Koch, 2013; Seehausen, 2015; Sundue matic, orogenic and ecological factors occurring at different temporal et al., 2015). This is supported by (1) the lack of phylogenetic signal scales, showing that large radiations are not the result of a single factor. detected for the global niche and the altitude variable, indicating niche divergence between sister species (Figs. 3 and 4, Table 1), and (2) CRediT authorship contribution statement growth form lability (Fig. 2; Table 1), which likely enabled rapid adaptation to new environments. Overall, the diversification of CA Sonia Herrando-Moraira: Investigation, Data Curation, Formal shows signs of adaptive radiation such as high morphological disparity analysis, Visualization, Writing – original draft. Cristina Roquet: and weak niche conservatism (Stroud and Losos, 2016). However, we Formal analysis, Writing – original draft, Writing – review & editing. suggest that diversification was also fostered by geographical isolation Juan-Antonio Calleja: Resources, Writing – review & editing. You- processes given the high climatic oscillations detected in Pamir/S Tian Sheng Chen:Resources. Kazumi Fujikawa: Resources. Mercè Gal­ Shan area (Figs. 5 and 6), which probably induced altitudinal and/or bany-Casals: Conceptualization, Investigation, Supervision, Writing – latitudinal migrations. review & editing. Núria Garcia-Jacas: Conceptualization, Supervision, The EA clade is the oldest and most species-poor one. In the Writing – review & editing. Jian-Quan Liu: Resources, Investigation. biogeographic history of Jurinea, the Himalayas acted as a bridge be­ Javier López-Alvarado: Resources. Jordi López-Pujol: Formal anal­ tween two climatic and floristic regions: the continental, summer-dry ysis, Funding acquisition, Writing – review & editing. Jennifer R Irano-Turanian and the summer-monsoon East Asian regimes (Takhta­ Mandel: Formal analysis, Writing – review & editing. Iraj Mehregan: jan, 1986; Djamali et al., 2012). However, Jurinea did not explosively Resources. Llorenç Sáez: Writing – review & editing. Alexander N radiate in either the Himalayan corridor or in East Asia (Fig. 1; Table 1), Sennikov: Resources, Investigation. Alfonso Susanna: Conceptualiza­ remaining there as a relatively depauperate lineage despite showing tion, Supervision, Funding acquisition, Writing – review & editing. wide ecological niche breadth (Figs. 3 and 4) and high morphological Roser Vilatersana: Resources. Lian-Sheng Xu: Resources, disparity (Fig. 2; Table 1). The low speciation rate detected (half of CB Investigation and CA clades) could be the result of the fact that the suitable adaptive landscape space was already almost filled when Jurinea arrived. Exam­ Declaration of Competing Interest ples of successful “native” competitors would be the closely-related, species-rich genus Saussurea DC., and other evolutionary successful The authors declare that they have no known competing financial groups such as Delphinium L. and Pedicularis L. (Wen et al., 2014). Given interests or personal relationships that could have appeared to influence its xerophytic ancestral background, Jurinea might not present enough the work reported in this paper. evolutionary potential to successfully diversify in humid environments. Analogous cases exist in two sister subtribes of Cardueae, Arctiinae and Data availability Onopordinae (Herrando-Moraira et al., 2019a). In Arctiinae, the explosive radiation of the mostly xerophytic Cousinia Cass. (ca. 600 spp.) Data will be made available on request. in the IT region contrasts with the depauperate genus Arctium L. (27 spp.), which is mainly adapted to mesophilous habitats (López-Vinyal­ Acknowledgements longa et al., 2009; Susanna and Garcia-Jacas, 2009). The subtribe Onopordinae shows a lineage with greater success in xeric habitats Authors thank Maria Luisa Gutiérrez, Fernando Castro, and Laia (Onopordum L., 60 spp.), in contrast to seven small genera in the mesic Barres for providing technical support during the laboratory process. We habitats of CA (Alfredia Cass., Ancathia DC., Lamyropappus Knorring & also thank the herbaria that provided material for the study: BC, DUSH, Tamamsch., Olgaea Iljin, Syreitschikovia Pavlov, Synurus Iljin, and Xan­ E, FRU, GDA, H, LE, MBK, MEM, MHA, MO, MU, OS, SANT, TAD, TARI, thopappus C.Winkl.). TI, TK, US, W, and WU. We gratefully acknowledge the University of Memphis computing facilities to allow us the utilization of its high- 14 S. Herrando-Moraira and C. Roquet Molecular Phylogenetics and Evolution 189 (2023) 107928 performance cluster (HPC). We also acknowledge funding support from Folk, R.A., Siniscalchi, C.M., Soltis, D.E., 2020. Angiosperms at the edge: Extremity, diversity, and phylogeny. Plant Cell Environ. 43, 2871–2893. the Ministerio de Ciencia e Innovación (Project CGL2015-66703-P Fragoso-Martínez, I., Salazar, G.A., Martínez-Gordillo, M., Magallón, S., Sánchez- MINECO and Ph.D. grant to Sonia Herrando-Moraira) and the Catalan Reyes, L., Lemmon, E.M., Lemmon, A.R., Sazatornil, F., Mendoza, C.G., 2017. A pilot government (“Ajuts a grups consolidats” 2017-SGR1116 and study applying the plant Anchored Hybrid Enrichment method to New World sages 2021SGR00315). (Salvia subgenus Calosphace; Lamiaceae). Mol. Phylogenet. Evol. 117, 124–134. Fu, B., Ninomiya, Y., Guo, J., 2010. Slip partitioning in the northeast Pamir-Tian Shan convergence zone. Tectonophysics 483, 344–364. Appendix A. Supplementary data Ge, D., Zhang, Z., Xia, L., Zhang, Q., Ma, Y., Yang, Q., 2012. Did the expansion of C4 plants drive extinction and massive range contraction of micromammals? Inferences from food preference and historical biogeography of pikas. Palaeogeogr. Supplementary data to this article can be found online at https://doi. Palaeoclimatol. Palaeoecol. 326, 160–171. org/10.1016/j.ympev.2023.107928. Harrison, S., Noss, R., 2017. Endemism hotspots are linked to stable climatic refugia. Ann. Bot. 119, 207–214. Herbert, T.D., Lawrence, K.T., Tzanova, A., Peterson, L.C., Caballero-Gill, R., Kelly, C.S., References 2016. Late Miocene global cooling and the rise of modern ecosystems. Nat. Geosci. 9, 843–847. Herrando-Moraira, S., The Cardueae Radiations Group (Calleja, J.A., Carnicero- An, Z., Kutzbach, J.E., Prell, W.L., Porter, S.C., 2001. Evolution of Asian monsoons and Campmany, P., Fujikawa, K., Galbany-Casals, M., Garcia-Jacas, N., Im, H.T., Kim, S. phased uplift of the Himalaya-Tibetan plateau since Late Miocene times. Nature 411, C., Liu, J.Q., López-Alvarado, J., López-Pujol, J., Mandel, J.R., Massó, S., Mehregan, 62–66. I., Montes-Moreno, N., Pyak, E., Roquet, C., Sáez, L., Sennikov, A., Susanna, A., An, Z., Huang, Y., Liu, W., Guo, Z., Clemens, S., Li, L., Prell, W., Ning, Y., Cai, Y., Vilatersana, R.)., 2018. Exploring data processing strategies in NGS target Zhou, W., Lin, B., Zhang, Q., Cao, Y., Qiang, X., Chang, H., Wu, Z., 2005. Multiple enrichment to disentangle radiations in the tribe Cardueae (Compositae). Mol. expansions of C4 plant biomass in East Asia since 7 Ma coupled with strengthened Phylogenet. Evol. 128, 69–87. monsoon circulation. Geology 33, 705. Herrando-Moraira, S., The Cardueae Radiations Group (Calleja, J.A., Galbany-Casals, M., Arakaki, M., Christin, P.A., Nyffeler, R., Lendel, A., Eggli, U., Ogburn, R.M., Spriggs, E., Garcia-Jacas, N., Liu, J.Q., López-Alvarado, J., López-Pujol, J., Mandel, J.R., Massó, Moore, M.J., Edwards, E.J., 2011. Contemporaneous and recent radiations of the S., Montes-Moreno, N., Roquet, C., Sáez, L., Sennikov, A., Susanna, A., Vilatersana, world’s major succulent plant lineages. Proc. Natl. Acad. Sci. USA 108, 8379–8384. R.)., 2019a. Nuclear and plastid DNA phylogeny of tribe Cardueae (Compositae) Banasiak, Ł., Piwczyński, M., Uliński, T., Downie, S.R., Watson, M.F., Shakya, B., with Hyb-Seq data: A new subtribal classification and a temporal diversification Spalik, K., 2013. Dispersal patterns in space and time: a case study of Apiaceae framework. Mol. Phylogenet. Evol. 137, 313–332. subfamily Apioideae. J. Biogeogr. 40, 1324–1335. Herrando-Moraira, S., The Cardueae Radiations Group (Calleja, J.A., Chen, Y.S., Barres, L., Sanmartín, I., Anderson, C.L., Susanna, A., Buerki, S., Galbany-Casals, M., Fujikawa, K., Galbany-Casals, M., Garcia-Jacas, N., Kim, S.C., Liu, J.Q., López- Vilatersana, R., 2013. Reconstructing the evolution and biogeographic history of Alvarado, J., López-Pujol, J., Mandel, J.R., Mehregan, I., Roquet, C., Sennikov, A., tribe Cardueae (Compositae). Am. J. Bot. 100, 867–882. Susanna, A., Vilatersana, R., Xu, L.S.)., 2020. Generic boundaries in subtribe Beck, H.E., Zimmermann, N.E., McVicar, T.R., Vergopolan, N., Berg, A., Wood, E.F., Saussureinae (Compositae: Cardueae): Insights from Hyb-Seq data. Taxon 69, 2018. Present and future Köppen-Geiger climate classification maps at 1-km 694–714. resolution. Sci. Data 5, 180214. Herrando-Moraira, S., Nualart, N., Herrando-Moraira, A., Chung, M.Y., Chung, M.G., Bialik, O.M., Frank, M., Betzler, C., Zammit, R., Waldmann, N.D., 2019. Two-step closure López-Pujol, J., 2019b. Climatic niche characteristics of native and invasive Lilium of the Miocene Indian Ocean Gateway to the Mediterranean. Sci. Rep. 9, 8842. lancifolium. Sci. Rep. 9, 1–16. Blomberg, S.P., Garland, Jr.T., Ives, A.R., 2003. Testing for phylogenetic signal in Herrando-Moraira, S., Nualart, N., Galbany-Casals, M., Garcia-Jacas, N., Ohashi, H., comparative data: behavioral traits are more labile. Evolution 57, 717–745. Matsui, T., Susanna, A., Tang, C.Q., López-Pujol, J., 2022. Climatic Stability Index Bolger, A.M., Lohse, M., Usadel, B., 2014. Trimmomatic: a flexible trimmer for Illumina maps, a global high resolution cartography of climate stability from Pliocene to sequence data. Bioinformatics 30, 2114–2120. 2100. Sci. Data 9, 48. Borges, R., Machado, J.P., Gomes, C., Rocha, A.P., Antunes, A., 2019. Measuring Hillis, D.M., Bull, J.J., 1993. An empirical test of bootstrapping as a method for assessing phylogenetic signal between categorical traits and phylogenies. Bioinformatics 35, confidence in phylogenetic analysis. Syst. Biol. 42, 182–192. 1862–1869. Holbourn, A.E., Kuhnt, W., Clemens, S.C., Kochhann, K.G., Jöhnck, J., Lübbers, J., Broennimann, O., Fitzpatrick, M.C., Pearman, P.B., Petitpierre, B., Pellissier, L., Andersen, N., 2018. Late Miocene climate cooling and intensification of southeast Yoccoz, N.G., Thuiller, W., Fortin, M.J., Randin, C., Zimmermann, N.E., Graham, C. Asian winter monsoon. Nat. Commun. 9, 1584. H., Guisan, A., 2012. Measuring ecological niche overlap from occurrence and Huang, Y., Clemens, S.C., Liu, W., Wang, Y., Prell, W.L., 2007. Large-scale hydrological spatial environmental data. Glob. Ecol. Biogeogr. 21, 481–497. change drove the late Miocene C4 plant expansion in the Himalayan foreland and Brown, J.L., Hill, D.J., Dolan, A.M., Carnaval, A.C., Haywood, A.M., 2018. PaleoClim, Arabian Peninsula. Geology 35, 531–534. high spatial resolution paleoclimate surfaces for global land areas. Sci. Data 5, Hughes, C.E., Atchison, G.W., 2015. The ubiquity of alpine plant radiations: from the 180254. Andes to the Hengduan Mountains. New Phytol. 207, 275–282. Burin, G., Alencar, L.R.V., Chang, J., Alfaro, M.E., Quental, T.B., 2019. How well can we Hughes, C.E., Nyffeler, R., Linder, P.H., 2015. Evolutionary plant radiations: where, estimate diversity dynamics for clades in diversity decline? Syst. Biol. 68, 47–62. when, why and how? New Phytol. 207, 249–253. Burtman, V.S., 2000. Cenozoic crustal shortening between the Pamir and Tien Shan and a Hurka, H., Friesen, N., Bernhardt, K.-G., Neuffer, B., Smirnov, S.V., Shmakov, A.I., reconstruction of the Pamir-Tien Shan transition zone for the Cretaceous and Blattner, F.R., 2019. The Eurasian steppe belt: Status quo, origin and evolutionary Palaeogene. Tectonophysics 319, 69–92. history. Turczaninowia 22, 5–71. Burtman, V.S., Molnar, P.H., 1993. Geological and geophysical evidence for deep Iljin, M.M., 1962. Jurinea Cass., in: Shiskin, B.K., Bobrov, E.G. (Eds.), Flora of the USSR, subduction of continental crust beneath the Pamir (Special Paper 281). The vol. 27. Izdatel’stvo Akademii Nauk SSSR, Moscow, pp. 538–704. Geological Society of America, Boulder. Jabbour, F., Renner, S.S., 2011. Consolida and Aconitella are an annual clade of Capella-Gutiérrez, S., Silla-Martínez, J.M., Gabaldón, T., 2009. trimAl: a tool for Delphinium (Ranunculaceae) that diversified in the Mediterranean basin and the automated alignment trimming in large-scale phylogenetic analyses. Bioinformatics Irano-Turanian region. Taxon 60, 1029–1040. 25, 1972–1973. Jiménez-Moreno, G., Pérez-Asensio, J.N., Larrasoaña, J.C., Sierro, F.J., Garcia- Casanovas-Vilar, I., García-Paredes, I., Alba, D.M., van den Hoek Ostende, L.W., Moyà- Castellanos, D., Salazar, A., Salvany, J.M., Ledesma, S., Mata, M.P., Mediavilla, C., Solà, S., 2010. The European Far West: Miocene mammal isolation, diversity and 2019. Early Pliocene climatic optimum, cooling and early glaciation deduced by turnover in the Iberian Peninsula. J. Biogeogr. 37, 1079–1093. terrestrial and marine environmental changes in SW Spain. Glob. Planet. Change Casanovas-Vilar, I., Alba, D.M., Garcés, M., Robles, J.M., Moyà-Solà, S., 2011. Updated 180, 89–99. chronology for the Miocene hominoid radiation in Western Eurasia. Proc. Natl. Acad. Johnson, M.G., Gardner, E.M., Liu, Y., Medina, R., Goffinet, B., Shaw, A.J., Zerega, N.J. Sci. USA 108, 5554–5559. C., Wickett, N.J., 2016. HybPiper: Extracting coding sequence and introns for Flora of Turkey and the East Aegean Islands vol. 5, 1975, 439–450. phylogenetics from high-throughput sequencing reads using target enrichment. App. Djamali, M., Brewer, S., Breckle, S.W., Jackson, S.T., 2012. Climatic determinism in Plant Sci. 4, 1600016. phytogeographic regionalization: a test from the Irano-Turanian region, SW and Jones, K.E., Fér, T., Schmickl, R.E., Dikow, R.B., Funk, V.A., Herrando-Moraira, S., Central Asia. Flora 207, 237–249. Johnston, P.R., Kilian, N., Siniscalchi, C.M., Susanna, A., Slovák, M., Thapa, R., Donoghue, M.J., Sanderson, M.J., 2015. Confluence, synnovation and depauperons in Watson, L.E., Mandel, J.R., 2019. An empirical assessment of a single family-wide plant diversification. New Phytol. 207, 260–274. hybrid capture locus set at multiple evolutionary timescales in Asteraceae. Appl. Drummond, A.J., Suchard, M.A., Xie, D., Rambaut, A., 2012. Bayesian phylogenetics Plant Sci. 7, e11295. with BEAUti and the BEAST 1.7. Mol. Biol. Evol. 29, 1969–1973. Karl, R., Koch, M.A., 2013. A world-wide perspective on crucifer speciation and Etienne, R.S., Haegeman, B., Stadler, T., Aze, T., Pearson, P.N., Purvis, A., Phillimore, A. evolution: phylogenetics, biogeography and trait evolution in tribe Arabideae. Ann. B., 2012. Diversity-dependence brings molecular phylogenies closer to agreement Bot. 112, 983–1001. with the fossil record. Proc. R. Soc. B Biol. Sci. 279, 1300–1309. Katoh, K., Standley, D.M., 2013. MAFFT multiple sequence alignment software version 7: FitzJohn, R.G., 2012. Diversitree: comparative phylogenetic analyses of diversification in improvements in performance and usability. Mol. Biol. Evol. 30, 772–780. R. Meth. Ecol. Evol. 3, 1084–11032. Köppen, W., Geiger, R., 1936. Handbuch der Klimatologie. Gebruder Borntraeger, Berlin. Folk, R.A., Stubbs, R.L., Mort, M.E., Cellinese, N., Allen, J.M., Soltis, P.S., Soltis, D.E., Koutsodendris, A., Allstädt, F.J., Kern, O.A., Kousis, I., Schwarz, F., Vannacci, M., Guralnick, R.P., 2019. Rates of niche and phenotype evolution lag behind Woutersen, A., Appeld, E., Berke, M.A., Fang, X., Friedrich, O., Hoorn, C., diversification in a temperate radiation. Proc. Natl. Acad. Sci. USA 116, Salzmann, U., Pross, J., 2019. Late Pliocene vegetation turnover on the NE Tibetan 10874–10882. 15 S. Herrando-Moraira and C. Roquet Molecular Phylogenetics and Evolution 189 (2023) 107928 Plateau (Central Asia) triggered by early Northern Hemisphere glaciation. Glob. Peterson, A.T., 2003. Predicting the geography of species’ invasions via ecological niche Planet. Change 180, 117–125. modeling. Q. Rev. Biol. 78, 419–433. Kozlov, A.M., Darriba, D., Flouri, T., Morel, B., Stamatakis, A., 2019. RAxML-NG: a fast, Peterson, A., Harpke, D., Peterson, J., Harpke, A., Peruzzi, L., 2019. A pre-Miocene Irano- scalable and user-friendly tool for maximum likelihood phylogenetic inference. Turanian cradle: Origin and diversification of the species-rich monocot genus Gagea Bioinformatics 35, 4453–4455. (Liliaceae). Ecol. Evol. 9, 5870–5890. Kück, P., Longo, G.C., 2014. FASconCAT-G: extensive functions for multiple sequence Potter, P.E., Szatmari, P., 2009. Global Miocene tectonics and the modern world. Earth- alignment preparations concerning phylogenetic studies. Front. Zool. 11, 81. Sci. Rev. 96, 279–295. LaJeunesse, T.C., 2005. “Species” radiations of symbiotic dinoflagellates in the Atlantic Pound, M.J., Haywood, A.M., Salzmann, U., Riding, J.B., 2012. Global vegetation and Indo-Pacific since the Miocene-Pliocene transition. Mol. Biol. Evol. 22, 570–581. dynamics and latitudinal temperature gradients during the Mid to Late Miocene Lauterbach, M., Veranso-Libalah, M.C., Sukhorukov, A.P., Kadereit, G., 2019. (15.97–5.33 Ma). Earth Sci. Rev. 112, 1–22. Biogeography of the xerophytic genus Anabasis L. (Chenopodiaceae). Ecol. Evol. 9, R Development Core Team., 2019. R: a Language and Environment for Statistical 3539–3552. Computing. R Foundation for Statistical Computing, Vienna, Austria http://www.R- Lisiecki, L.E., Raymo, M.E., 2005. A Pliocene-Pleistocene stack of 57 globally distributed project.org/. benthic δ18O records. Paleoceanography 20, PA1003. Rea, D.K., Snoeckx, H., Joseph, L.H., 1998. Late Cenozoic Eolian deposition in the North López-Vinyallonga, S., Mehregan, I., Garcia-Jacas, N., Tscherneva, O., Susanna, A., Pacific: Asian drying, Tibetan uplift, and cooling of the north hemisphere. Kadereit, J.W., 2009. Phylogeny and evolution of the Arctium-Cousinia complex Paleoceanography 13, 215–224. (Compositae, Cardueae-Carduinae). Taxon 58, 153–171. Flora Iranica vol. 139a, 1979, 180–209. Louca, S., Pennell, M.W., 2020. Extant timetrees are consistent with a myriad of Ree, R.H., Moore, B.R., Webb, C.O., Donoghue, M.J., 2005. A likelihood framework for diversification histories. Nature 580, 502–505. inferring the evolution of geographic range on phylogenetic trees. Evolution 59, Magyari, E.K., Chapman, J.C., Gaydarska, B., Marinova, E., Deli, T., Huntley, J.P., 2299–2311. Allen, J.R.M., Huntley, B., 2008. The ‘oriental’ component of the Balkan flora: Ree, R.H., Sanmartín, I., 2018. Conceptual and statistical problems with the DEC+ J evidence of presence on the Thracian Plain during the Weichselian late-glacial. model of founder-event speciation and its comparison with DEC via model selection. J. Biogeogr. 35, 865–883. J. Biogeogr. 45, 741–749. Mahmoudi Shamsabad, M., Moharrek, F., Assadi, M., Nieto Feliner, G., 2021. Revell, L.J., 2012. phytools: an R package for phylogenetic comparative biology (and Biogeographic history and diversification patterns in the Irano-Turanian genus other things). Methods Ecol. Evol. 3, 217–223. Acanthophyllum s.l. (Caryophyllaceae). Plant Biosyst. 155, 425–435. https://doi.org/ Ribera, I., Blasco-Zumeta, J., 1998. Biogeographical links between steppe insects in the 10.1080/11263504.2020.1756974. Monegros region (Aragón, NE Spain), the eastern Mediterranean, and central Asia. Manafzadeh, S., Salvo, G., Conti, E., 2014. A tale of migrations from east to west: the J. Biogeogr. 25, 969–986. Irano-Turanian floristic region as a source of Mediterranean xerophytes. J. Biogeogr. Rudov, A., Mashkour, M., Djamali, M., Akhani, H., 2020. A review of C4 plants in 41, 366–379. Southwest Asia: an ecological, geographical and taxonomical analysis of a region Manafzadeh, S., Staedler, Y.M., Conti, E., 2017. Visions of the past and dreams of the with high diversity of C4 eudicots. Front. Plant Sci. 11, 546518. future in the Orient: the Irano-Turanian region from classical botany to evolutionary Sales, F., Hedge, I.C., 2013. Generic endemism in South-West Asia: an overview. studies. Biol. Rev. 92, 1365–1388. Rostaniha 14, 22–35. Mandel, J.R., Barker, M.S., Bayer, R., Dikow, R.B., Gao, T.-G., Jones, K.E., Keeley, S., Salvo, G., Manafzadeh, S., Ghahremaninejad, F., Tojibaev, K., Zeltner, L., Conti, E., 2011. Kilian, N., Ma, H., Siniscalchi, C., Susanna, A., Thapa, R., Watson, L., Funk, V.A., Phylogeny, morphology, and biogeography of Haplophyllum (Rutaceae), a species- 2017. The Compositae Tree of Life in the age of phylogenomics. Journal of rich genus of the Irano-Turanian floristic region. Taxon 60, 513–527. Systematics and Evolution 55, 405–410. Sayyari, E., Mirarab, S., 2016. Fast coalescent-based computation of local branch support Mandel, J.R., Dikow, R.B., Funk, V.A., Masalia, R.R., Staton, S.E., Kozik, A., from quartet frequencies. Mol. Biol. Evol. 33, 1654–1668. Michelmore, R.W., Rieseberg, L.H., Burke, J.M., 2014. A target enrichment method Seehausen, O., 2015. Process and pattern in cichlid radiations - inferences for for gathering phylogenetic information from hundreds of loci: an example from the understanding unusually high rates of evolutionary diversification. New Phytologist Compositae. App. Plant Sci. 2, 1300085. 207, 304–312. Matzke, N.J., 2013. BioGeoBEARS: BioGeography with Bayesian (and likelihood) Shen, X., Wan, S., Colin, C., Tada, R., Shi, X., Pei, W., Tan, Y., Jiang, X., Li, A., 2018. evolutionary analysis in R Scripts: CRAN: The Comprehensive R Archive Network. Increased seasonality and aridity drove the C4 plant expansion in Central Asia since [WWW document] URL Available from: https://cran.r-project.org/. the Miocene-Pliocene boundary. Earth Planet. Sci. Lett. 502, 74–83. Mandel, J.R., Dikow, R.B., Siniscalchi, C., Funk, V.A., 2019. A fully resolved backbone Smith, S.A., O’Meara, B.C., 2012. treePL: divergence time estimation using penalized phylogeny reveals numerous dispersals and explosive diversifications throughout the likelihood for large phylogenies. Bioinformatics 28, 2689–2690. history of Asteraceae. Proc. Nat. Acad. Sci. USA 116, 14083–14088. Sniderman, J.K., Woodhead, J.D., Hellstrom, J., Jordan, G.J., Drysdale, R.N., Tyler, J.J., Miao, Y., Herrmann, M., Wu, F., Yan, X., Yang, S., 2012. What controlled Mid-Late Porch, N., 2016. Pliocene reversal of late Neogene aridification. Proc. Natl. Acad. Sci. Miocene long-term aridification in Central Asia?—Global cooling or Tibetan Plateau USA 113, 1999–2004. uplift: A review. Earth Sci. Rev. 112, 155–172. Stebbins, G.L., 1952. Aridity as a stimulus to plant evolution. Am. Nat. 86, 33–44. Middleton, N., Thomas, D., 1997. World atlas of desertification, second ed. United Stroud, J.T., Losos, J.B., 2016. Ecological opportunity and adaptive radiation. Annu. Rev. Nations Environment Programme/Edward Arnold, London. Ecol. Evol. Syst. 47, 507–532. Miller, M.A., Pfeiffer, W., Schwartz, T., 2010. Creating the CIPRES Science Gateway for Suc, J.P., 1984. Origin and evolution of the Mediterranean vegetation and climate in inference of large phylogenetic trees. Proceedings of the Gateway Computing Europe. Nature 307, 429–432. Environments Workshop (GCE). New Orleans, LA: IEEE, 1–8. Sun, J., Gong, Z., Tian, Z., Jia, Y., Windley, B., 2015. Late Miocene stepwise aridification Mittermeier, R.A., Turner, W.R., Larsen, F.W., Brooks, T.M., Gascon, C., 2011. Global in the Asian interior and the interplay between tectonics and climate. Palaeogeogr. biodiversity conservation: the critical role of hotspots. In: Zachos, F.E., Habel, J.C. Palaeoclimatol. Palaeoecol. 421, 48–59. (Eds.), Biodiversity hotspots distribution and protection of conservation priority Sundue, M.A., Testo, W.L., Ranker, T.A., 2015. Morphological innovation, ecological areas. Springer, Heidelberg, pp. 3–22. opportunity, and the radiation of a major vascular epiphyte lineage. Evolution 69, Moharrek, F., Sanmartín, I., Kazempour-Osaloo, S., Nieto Feliner, G., 2019. 2482–2495. Morphological innovations and vast extensions of mountain habitats triggered rapid Susanna, A., Garcia-Jacas, N., 2007. In: Kadereit, J.W., Jefferey, C. (Eds.), The families diversification within the species-rich Irano-Turanian genus Acantholimon and genera of vascular plants, 8. Springer, Berlin / Heidelberg, pp. 123–146. (Plumbaginaceae). Front. Genet. 9, 698. Susanna, A., Garcia-Jacas, N., 2009. Cardueae (Carduoideae). In: Funk, V.A., Moreau, C.S., Bell, C.D., 2013. Testing the museum versus cradle tropical biological Susanna, A., Stuessy, T.F., Bayer, R.J. (Eds.), Systematics, evolution, and diversity hypothesis: phylogeny, diversification, and ancestral biogeographic range biogeography of Compositae. IAPT, Vienna, pp. 293–313. evolution of the ants. Evolution 67, 2240–2257. Susanna, A., Garcia-Jacas, N., Hidalgo, O., Vilatersana, R., Garnatje, T., 2006. The Morlon, H., Lewitus, E., Condamine, F.L., Manceau, M., Clavel, J., Drury, J., 2016. Cardueae (Compositae) revisited: insights from ITS, trnL-trnF, and matK nuclear and RPANDA: an R package for macroevolutionary analyses on phylogenetic trees. Meth. chloroplast DNA analysis. Ann. Missouri Bot. Gard. 93, 150–171. Ecol. Evol. 7, 589–597. Szukala, A., Korotkova, N., Gruenstaeudl, M., Sennikov, A.N., Lazkov, G.A., Morlon, H., Robin, S., Hartig, F., 2022. Studying speciation and extinction dynamics Litvinskaya, S.A., Gabrielian, E., Borsch, T., von Raab-Straube, E., 2019. Phylogeny from phylogenies: addressing identifiability issues. Trends Ecol. Evol. 37, 497–506. of the Eurasian genus Jurinea (Asteraceae: Cardueae): Support for a monophyletic Mouthereau, F., 2011. Timing of uplift in the Zagros belt/Iranian plateau and genus concept and a first hypothesis on overall species relationships. Taxon 68, accommodation of late Cenozoic Arabia-Eurasia convergence. Geol. Mag. 148, 112–131. 726–738. Takhtajan, A., 1986. Floristic Regions of the World. University of California Press, Nee, S., Mooers, A.O., Harvey, P.H., 1992. Tempo and mode of evolution revealed from Berkeley. molecular phylogenies. Proc. Natl. Acad. Sci. USA 89, 8322–8326. Thompson, J.A., Burbank, D.W., Li, T., Chen, J., Bookhagen, B., 2015. Late Miocene Nevado, B., Atchison, G.W., Hughes, C.E., Filatov, D.A., 2016. Widespread adaptive northward propagation of the northeast Pamir thrust system, northwest China. evolution during repeated evolutionary radiations in New World lupins. Nat. Tectonics 34, 510–534. Commun. 7, 12384. Uribe-Convers, S., Tank, D.C., 2015. Shifts in diversification rates linked to Noroozi, J., Naqinezhad, A., Talebi, A., Doostmohammadi, M., Plutzar, C., Rumpf, S.B., biogeographic movement into new areas: An example of a recent radiation in the Asgarpour, Z., Schneeweiss, G.M., 2019. Hotspots of vascular plant endemism in a Andes. Am. J. Bot. 102, 1854–1869. global biodiversity hotspot in Southwest Asia suffer from significant conservation Wen, J., Zhang, J., Nie, Z.L., Zhong, Y., Sun, H., 2014. Evolutionary diversifications of gaps. Biol. Conserv. 237, 299–307. plants on the Qinghai-Tibetan Plateau. Front. Genet. 5, 4. Paradis, E., Schliep, K., 2019. ape 5.0: An environment for modern phylogenetics and White, F., Léonard, J., 1991. Phytogeographic links between Africa and Southwest Asia. evolutionary analyses in R. Bioinformatics 35, 526–528. Flora Veg. Mundi 9, 229–246. 16 S. Herrando-Moraira and C. Roquet Molecular Phylogenetics and Evolution 189 (2023) 107928 Wiens, J.J., Ackerly, D.D., Allen, A.P., Anacker, B.L., Buckley, L.B., Cornell, H.V., Yu, Y., Harris, A.J., Blair, C., He, X., 2015. RASP (Reconstruct Ancestral State in Damschen, E.I., Davies, T.J., Grytnes, J.A., Harrison, S.P., Hawkins, B.A., Holt, R.D., Phylogenies): a tool for historical biogeography. Mol. Phylogenet. Evol. 87, 46–49. McCain, C.M., Stephens, P.R., 2010. Niche conservatism as an emerging principle in Zachos, J.C., Pagani, M., Sloan, L., Thomas, E., Billups, K., 2001. Trends, rhythms, and ecology and conservation biology. Ecol. Lett. 13, 1310–1324. aberrations in global climate 65 Ma to present. Science 292, 686–693. Willeit, M., Ganopolski, A., Calov, R., Robinson, A., Maslin, M., 2015. The role of CO2 Zhang, Z., Han, W., Fang, X., Song, C., Li, X., 2013. Late Miocene-Pleistocene decline for the onset of Northern Hemisphere glaciation. Quat. Sci. Rev. 119, 22–34. aridification of Asian inland revealed by geochemical records of lacustrine-fan delta Winkworth, R.C., Wagstaff, S.J., Glenny, D., Lockhart, P.J., 2005. Evolution of the New sediments from the western Tarim Basin, NW China. Palaeogeogr. Palaeoclimatol. Zealand mountain flora: Origins, diversification and dispersal. Organisms Div. Evol. Palaeoecol. 377, 52–61. 5, 237–247. Zhang, C., Rabiee, M., Sayyari, E., Mirarab, S., 2018. ASTRAL-III: polynomial time Wu, S.D., Lin, L., Li, H.L., Yu, S.X., Zhang, L.J., Wang, W., 2015. Evolution of Asian species tree reconstruction from partially resolved gene trees. BMC Bioinformatics interior arid-zone biota: Evidence from the diversification of Asian Zygophyllum 19, 153. (Zygophyllaceae). PloS ONE 10 (9), e0138697. Zheng, H., Powell, C.M., Rea, D.K., Wang, J., Wang, P., 2004. Late Miocene and mid- Ye, X.Y., Ma, P.F., Yang, G.Q., Guo, C., Zhang, Y.X., Chen, Y.M., Guo, Z.H., Li, D.Z., 2019. Pliocene enhancement of the East Asian monsoon as viewed from the land and sea. Rapid diversification of alpine bamboos associated with the uplift of the Hengduan Glob. Planet. Change 41, 147–155. Mountains. J. Biogeogr. 46, 2678–2689. Zohary, M., 1981. On the flora and vegetation of the Middle East structure and evolution. Yoder, J.B., Clancey, E., Des Roches, S., Eastman, J.M., Gentry, L., Godsoe, W., In: Uerpmann, H.P., Frey, W. (Eds.), Contributions to the Environmental History of Harmon, L.J., 2010. Ecological opportunity and the origin of adaptive radiations. Southwest Asia, vol. 8. Reichert, Wiesbaden, pp. 1–17. J. Evol. Biol. 23, 1581–1596. Zohary, M., 1973. Geobotanical Foundations of the Middle East (Volume I and II). Amsterdam: Gustav Fischer Verlag, and Stuttgart : Swets & Zeitlinger. 17