Biogeography of Rhaponticoides, an Irano-Turanian element in the Mediterranean flora

Floristic relationships between the Irano-Turanian and Mediterranean regions have been known from old. However, only a few biogeographical analyses based on molecular data have evaluated the history of steppe plants within the Mediterranean basin. Our study aims to contribute to a better understanding of the migratory and diversification processes by reconstructing the biogeography of Rhaponticoides (Cardueae), distributed in the Mediterranean and Irano-Turanian regions. We generated nuclear and plastid sequences that were analyzed by Bayesian inference. We used the resulting phylogeny for dating the diversification of the genus and examining the dispersal pathways. Two clades were recovered, an Irano-Turanian clade and a Mediterranean clade. The origin of the genus was placed in the Anatolian plateau in the Middle Miocene. The genus experienced several diversifications and expansions correlated to the Messinian salinity crisis and the environmental changes in the Pliocene and the Quaternary. Rhaponticoides migrated following two routes reflecting the two souls of the genus: Irano-Turanian taxa colonized the steppes of Eurasia whilst Mediterranean species migrated via eastern and central Mediterranean and North Africa, leaving a trail of species; both pathways ended in the Iberian Peninsula. Our study also confirms that more work is needed to unravel phylogenetic relationships in Rhaponticoides.

. 50% majority-rule consensus tree obtained by Bayesian analysis of the combined nuclear dataset, indicating supported clades. Numbers occurring above branches are posterior probabilities (PP lower than 0.6 are not shown). Likelihood Bootstrap values figure below branches (BS lower than 60% are not shown). Capital letters following the names of species correspond to the countries of origin for species with more than one sample (see Table 1

Rhaponticoides aytachii
Rhaponticoides africana PT  Figure 3. 50% majority-rule consensus tree obtained by Bayesian analysis of the combined plastid dataset, indicating supported clades. Numbers occurring above branches are posterior probabilities (PP lower than 0.6 are not shown). Likelihood Bootstrap values figure below branches (BS lower than 60% are not shown). Capital letters following the names of species correspond to the countries of origin for species with more than one sample (see Table 1 www.nature.com/scientificreports/ Biogeography within a time-calibrated framework. The BayArea + J model performed best to reconstruct the ancestral area of Rhaponticoides (Table S2). The final time-calibrated biogeographic model revealed that Turkey is the most plausible ancestral area of Rhaponticoides in the middle of the Miocene (12 Ma). The genus diverged in two geographical lines, the Irano-Turanian and Mediterranean clades, in the Late Miocene (8 Ma). Both branches diversified again in the Miocene-Pliocene transition (5 Ma) and show a radiation process in the Pliocene-Quaternary transition and throughout the Quaternary period (Fig. 4, Table 3). The diversification processes entailed different expansion processes towards western and central Asia, on the one hand, and towards the Mediterranean basin on the other. This diversification occurred in different lineages and on different dates. The Irano-Turanian steppe line originated in the steppe environments of Turkey, i.e., the Anatolian plateau (red region in Fig. 4) and diversified (e.g. Rh. amasiensis) (nodes 6 and 10 in Fig. 4). This line also expanded towards Iran and the Caucasus (nodes 7 and 8, route 1a in Fig. 4, Table 3), and gave rise to the Iranian species complex of taxa (Rh. lachnopus) and the Caucasian group (Rh. hajastana). The same branch also originated the Rh. alpina/ruthenica complex, which integrates populations that occur from West India through the Alps and Mediterranean basin up to the Iberian Peninsula. These European populations (yellow and green areas in Fig. 4) represent a westward migration track (route 1b in Fig. 4) across the northern rim of the Mediterranean basin.
The Mediterranean line (node 11 in Fig. 4) also exhibits an Asian origin with ulterior diversification of species. This line shares with the steppe one the ancestral area in Turkey (red area in Fig. 4), and originated a first clade (nodes 11 and 12 in Fig. 4) Fig. 4).

Discussion
Origin of Rhaponticoides, tempo, and location. Steppe flora and vegetation are present in the Mediterranean basin (comprised Western Europe) since the Paleogene and gained relevance throughout the Neogene, especially during the Miocene 10,20 . The geographical origin of the main steppe lineages that colonized the Mediterranean basin is Central and Western Asia 15,30,51 and this is also the case for Rhaponticoides (Fig. 4) because Turkey appears as the most probable ancestral area. Turkey's landscape is dominated by the Anatolian plateau, which belongs to the Irano-Turanian region 15,52 , and a Mediterranean belt located in the south and west of the Anatolian Peninsula.
The Anatolian plateau, located at the intersection of Europe, Asia, and Africa, benefits from the diversity of three hotspots: East Mediterranean, Iran-Anatolia, and the Caucasus 29,50 . It is also a meeting place and dispersal corridor for different lineages that originated in Asia and colonized Europe and North Africa during the Cenozoic 15,53 , among them subtribe Centaureinae to which Rhaponticoides belongs 37 . Moreover, Anatolian plateau exhibits a high level of endemism and rich biodiversity probably driven by a complex paleogeographic history with dramatic topography and climate changes during the Cenozoic 54,55 .
The origin of Rhaponticoides should be placed in Middle Miocene as suggested for different taxa of vascular plants of Irano-Turanian origin 50,56 and many groups from the Compositae 33,34 . Most diversification events within Cardueae and especially in subtribe Centaureinae are related to recurrent connection and isolation episodes between Anatolian plateau and the Mediterranean basin throughout the Miocene 37 . In turn, these episodes are linked to environmental changes across the Irano-Turanian and Mediterranean regions such as climatic changes, normally tending to cooling and aridification, and collisions of tectonic plates and subsequent uplifts of plateaus and mountains comprising those close to Anatolia, e.g., Zagros 29,50 .
Diversification of Rhaponticoides. The diversification of Rhaponticoides took place in the late Cenozoic ( Fig. 4) as described for other vascular plant lineages 57 . It started at the end of the Miocene with subsequent radiation events in the Miocene-Pliocene and Pliocene-Quaternary transitions (Fig. 4). The first diversification events in the late Miocene and the Miocene-Pliocene transition especially implied dispersion. As registered in other steppe Irano-Turanian taxa 25,49 , Rhaponticoides steppe lineage expanded eastwards and originated species currently growing in Iran like Rh. lachnopus. In parallel, the genus shows a Mediterranean line that originated in the Mediterranean region of Turkey, which migrated to the west reaching southwestern Europe and northwestern Africa. These dispersal and diversification events coincide with those inferred in other steppe xerophytes that are present in the Mediterranean basin such as Anabasis 28 , Acantholimon 29 , and Haplophyllum 15 . Expansion of Irano-Turanian xerophytes had been favored by the extreme aridification and partial desiccation of the Mediterranean Sea during the Messinian crisis, which would have led to landmass connections, i.e., corridors towards the western end of the Mediterranean basin. Thus, the early diversification of Rhaponticoides fits the Messinian Model 14 .
The Messinian Model is compatible in turn with posterior dispersal and diversification processes that occurred during the Pliocene-Quaternary transition and Quaternary that may be related to Mediterranean climate onset 8,58 and also associated with the climatic and sea-level oscillations associated with glaciations 21,28,29 . The main lineages of Rhaponticoides which were previously originated in the Miocene-Pliocene transition experienced diversification and colonized the Mediterranean basin and western Asia. The Mediterranean line evolved probably from the steppe line into species that occurred in the Mediterranean belt of Anatolia (e.g. Rh. amasiensis, Rh. hierroi, Rh. mykalea) and in Armenia (Rh. hajastana), and then migrated to the eastern and central Mediterranean (Balkans and Italy, e.g., Rh. amplifolia, Rh. calabrica, Rh. centaurium, Rh. wagenitziana). In addition, within the steppe line, some species emerged and achieved wide and disjunct distributions either by fragmentation of a wide original range or by long-distance dispersal events 59 Figure 4. Molecular dating and biogeographic analyses. Maximum clade credibility tree from the relaxed molecular-clock analysis with exponential distribution and Birth and Death speciation process of ITS and ETS sequences in BEAST. Numbers refer to the supported nodes (Table 3). Pie charts (shown only for supported nodes, Table 3 Considering the environmental variations throughout the late Cenozoic, vicariance could explain, at least partially, the origin of lineages and species of Rhaponticoides. However, speciation by dispersal events seems to be a lot more frequent than vicariance in steppe taxa located in the Mediterranean and Irano-Turanian regions 15,29,30,51 . In our analyses, BayArea + J was selected as the best model, and this model embraces cladogenesis and sympatry as well as anagenesis and dispersion whereas vicariance has no relevance 60 . Rhaponticoides species of the steppe and Mediterranean lines would have originated in regions with biogeographical particularities and complex environmental history such as the Anatolian plateau, which harbors multiple endemic species 54,55 . Range expansions were favored by the corridors enabled by events such as the Messinian crisis 15,29 . The incongruence between plastid and nuclear data detected in Iberian Rh. fraylensis and Rh. africana (Figs. 2 and 3) is very illustrative on the origin of the Mediterranean line within the Irano-Turanian pool: both species appear in the plastid phylogeny nested in the Irano-Turanian clade. The usual explanation for these inconsistencies is ancient hybridization and subsequent plastid capture 61,62 . In our case, the most plausible hypothesis suggests that the ancestor of Rh. africana and Rh. fraylensis acquired the Irano-Turanian chloroplast by introgression in the contact zone between both clades in Anatolia. The other incongruences between plastid and nuclear phylogenies (i.e., Rh. amasiensis, Rh gokceoglui, Rh. iconiensis and Rh. phytiae) can be reduced to only one supported incongruence in Rh. gokceoglui since the position of Rh. amasiensis, Rh. iconiensis and Rh. phytiae is unsupported in the plastid tree (PP = 0.76, BS = 63%, Fig. 3). The conflicting position of Rh. gokceoglui should be also explained by an ancient event of hybridisation and plastid capture, which is a common fenomenon in tribe Cardueae 33,36 .
Migration routes toward the Mediterranean basin. To date, two routes have been usually reported for the Irano-Turanian steppe flora currently located within the Mediterranenan basin. Some taxa followed a "northern route" encompassing mountains and/or steppes of southern Europe and the northern rim of the Mediterranean Sea 11,27 . In contrast, other taxa tracked a "southern route" through North Africa, the Mediterranean islands and eventually through landmasses that emerged during the arid crises of the late Neogene and Quaternary 26,28 . Rhaponticoides comprises lineages that match both routes and reflect both biogeographic patterns (routes 1b and 2 in Fig. 4).
Regarding the north route (route 1b in Fig. 4), the steppe line that originated the Rh. alpina/ruthenica complex reached the western edge of the Mediterranean basin (Iberian Peninsula) through the Balkans and Alps, leaving in both massifs relictic populations of Rh. alpina as milestones 42 . This line also encompasses species from the western and central Asian steppes with a distribution centered on the Irano-Turanian region reaching the steppes of Iran in the south (represented in our analyses by Rh. lachnopus) (Fig. 4). As for the south route (route 2 in Fig. 4), the Mediterranean line integrates a series of species of narrow distribution, many of them from the Mediterranean part of Anatolia 41,46 extending westwards to North Africa and the Iberian Peninsula (Fig. 4).
The coexistence of both pathways and the arrival of the two lines-steppe and Mediterranean-to the Iberian Peninsula matches the case of Centaurea sect. Acrocentron (Cass.) DC. 25 . The pathway of Centaurea from Sicily to Spain is punctuated by some relict species that stand as stepping-stones, namely (from east to west) C. tauromenitana, C. carolipauana, and C. clementei. In contrast, stepping-stones in the south pathway of Rhaponticoides consist only in relict, isolated populations of a single species, namely Rhaponticoides africana, with small stands in Sicily, North Africa, and south and NW Iberia. The remarkable journey of Rh. africana ended in Galicia (north-west Spain), where some relict populations barely survive on the coastal scree 63 .
Taxonomic implications. The aims of this paper are not taxonomic, but our results show that the taxonomy of the genus is far from being well-known. The infrageneric taxonomy of Rhaponticoides has been studied for long but primarily based on morphological criteria. The most exhaustive and recent work 41 proposed a sectional classification of the genus in which virtually every species is a section or subsection, much in agreement with 39 . However, our results reveal discrepancies in the monophyly of some sections and species. Sections Ruthenicae and Aralocaspicae are not monophyletic and both of them include species recovered into different and unrelated clades. Such inconsistencies might arise from the species delimitation made by 39 . For instance, the morphological evidence for separating Rh. ruthenica from Rh. alpina are extremely weak and merely based on some vegetative characters 41 whereas our phylogenetic results confirm that this segregation lacks support from molecular evidence (Fig. 2).
In the same line, many taxonomic entities at species level of quite doubtful value in our opinion were described within the wide range Rh. ruthenica. Supporting this view, our results place Rh. razdorskyi, one of these segregate species, grouped with Rh. ruthenica. Within the same section Ruthenicae, the independence of Rh. linaresii, represented in our analyses by the Valencia population of Rh. alpina and sustained by 41 , has been rejected by all revisions of the group 42,64 , who merged it into Rh. alpina as confirmed by our results. Similarly, Rh. carrisoi is listed as separate species in 41 but it is actually a synonym of Rh. fraylensis 64 . Finally, Rh. centaurium and Rh. calabrica have been considered synonyms 65 which, according to our results, is highly probable. Instead, the species would be very narrowly connected to Rh. wagenitziana from the Balkans 45 .
In sum, the global richness and delimitation of some infrageneric entities, sections, and species in Rhaponticoides requires a deep, serious revision that should be carried out within the framework of the integrative taxonomy 66 . We suggest incorporating a wider representation of taxa and the use of the phylogenomic approach already used in the Cardueae 37 .

Concluding remarks
The biogeographical history of the genus set its origin in the Irano-Turanian part of Turkey (Anatolian plateau) in the Middle Miocene. The genus experienced different diversifications and westward and eastward expansions related not only to the Messinian salinity crisis but also to ulterior environmental changes in the www.nature.com/scientificreports/ Pliocene-Quaternary and Quaternary periods. Rhaponticoides, like other Irano-turanian steppe taxa, colonized the Mediterranean basin on different dates. However, in contrast to most of the previously studied taxa, the genus Rhaponticoides migrated and diversified across the Mediterranean basin following two different routes. The construction of the phylogenetic relationships within Rhaponticoides reveals the urgent necessity of a comprehensive integrative study on the genus to resolve the delimitation of some infrageneric taxa.

Material and methods
The methods comply with local and national guidelines.
Sampling. Sampling was designed to cover all the area and the nuclei of diversification of the genus: Caucasus, Iran, Balkans, Italian and Iberian peninsulas, Sicily, North Africa, and the steppes of Eurasia, reaching the extreme of the area in India (Fig. 1). Special focus was posed in Turkey, where the genus shows its peak of species according to 44 Table 1.
DNA extraction and amplification. Total genomic DNA was extracted from silica gel-dried tissue of one plant per population. The extraction of DNA followed the CTAB method of 67 with the modifications of 68 including three washing steps with sorbitol buffer. The ITS, ETS, rpl32-trnL UAG , and ycf3-trnS regions were amplified by polymerase chain reaction (PCR). The amplification primers for the nuclear regions were ITS1 and ITS4 69 for the ITS, and ETS1F 70 and 18SETS 71 for the ETS region. For the plastid rpl32-trnL UAG region, we used rpl32-F as the forward primer and trnL UAG as the reverse primer 72 . For the plastid ycf3-trnS region, we used SP43122F as the forward primer and SP44097R as the reverse primer 73 . The PCR reactions were performed under the conditions detailed in 74 . The amplified DNA segments were sequenced on an ABI 3730XL Analyser (Applied Biosystems, Foster City, CA, USA) following the manufacturer's protocol at Macrogen Inc., Korea.
Phylogenetic analyses. Nucleotide sequences were edited using Bioedit v7.0.5.3 75 and aligned visually by sequential pairwise comparison 76 . Basic data on ITS, ETS, rpl32-trnL UAG , and ycf3-trnS matrices are available in Table S3. The ITS plus ETS matrix was 1594 bp long (dataset 1) and the aligned plastid matrix was 1821 bp long (dataset 2). Likelihood analysis of both datasets was carried out by heuristic search using PAUP v4.0 77 using TBR branch swapping with character states specified as unordered and unweighted. The likelihood criterion was set to HKY85 (the default option in PAUP). Bootstrap analyses 78 were performed using 100 replicates of heuristic search with the default options. Internodes with a Bootstrap (BS) value > 75% were considered statistically significant. Bayesian inference of the two datasets was calculated using MrBayes v3.2.6 79 . The best available model of molecular evolution required for Bayesian estimations of phylogeny was selected using Akaike information criteria (AIC) for both datasets as implemented in the software MrModeltest v2.2 80 . The best fitting model was GTR + G + I for both datasets. Bayesian inference analyses were initiated with random starting trees and were run for 40 × 10 6 generations in two independent runs of four Metropolis-coupled chains. We saved one out of every 1000 generations, resulting in 40,000 sampled trees. Data from the first 1000 sampled trees were discarded as the ''burn-in'' period, after confirming that log-likelihood values had stabilized prior to the 1000th sampled tree. The stationarity of the runs and the convergence between the runs were checked with Tracer v.1.5.0 81 . Internodes with posterior probabilities (PP) > 0.95 were considered statistically significant.
Dating analyses. Divergence times were estimated using the nrDNA (ETS and ITS) sequences organized in a matrix with one-two accessions for each taxon of the ingroup (dataset 3, Table 1). The outgroups included Klasea coriacea, K. serratuloides (as in dataset 1) plus Plagiobasis centauroides Schrenk, Leuzea acaulis (L.) Holub., and L. conifera (L.) DC. Dating analyses were performed by using BEAST v.1.8.4. Four monophyletic groups were defined in BEAUti v.1.8.4 (included in BEAST package): (i) all species of dataset 2, (ii) Plagiobasis centauroides, (iii) Klasea plus Leuzea clade, (iv) and Rhaponticoides genus. These four groups were also implemented as secondary calibration points based on a previous phylogenomic study focused on the Cardueae tribe 37 , see Table 2.
To obtain the best time-calibrated phylogram, models with strict and uncorrelated log-normal relaxed clocks were run under two different speciation tree models, Yule and birth-death 82,83 . The options of lognormal and exponential distributions were also tested in the case of models with log-normal relaxed clocks (see Table S1). The six resulting BEAST models were run for four independent chains of 50 million generations each, sampling every 1000 generations. Their convergence was assessed by confirming that all parameters had reached stationarity and sufficient effective sample sizes (> 200) in all converged runs using Tracer v1.7 81 . All models and replicates were run in CIPRES Science Gateway 84 .
The six 6 different BEAST models were compared according to their respective log values of Marginal Likelihood Estimators (MLE) that were obtained with path-sampling (PS) and stepping-stone (SS) 85 implemented in BEAST v.1.8.4 (see Table S1). PS and SS log values were estimated with 100 path steps, a chain length of 106 generations and likelihoods saved every 1000 generations. The resulting log values of MLEs were averaged across  86 (i.e., twice the difference between the harmonic mean likelihoods of the two models). Values for 2BF those are greater than 2, 6, and 10 indicate positive, strong, and decisive support, respectively, for the generic hypothesis with minor marginal likelihood.
The best model with decisive support was that with a relaxed clock with exponential distribution and Birth and Death speciation process (Table S1). After discarding the burn-in steps (25%), tree files from the four independent runs of the selected model were combined using LogCombiner 1. 8 Ancestral area estimation. We inferred the origin of the genus Rhaponticoides and its possible routes of expansion and speciation from the interspecific phylogenetic relationships supported by a time-calibrated tree that was subjected to a biogeographical analysis, BioGeoBEARS 87 . We defined 9 geographic regions based, mainly, on the richness and endemicity of species in Rhaponticoides genus (Fig. 4, above). The time-calibrated tree resulting from Bayesian inference of BEAST was used as an input file to estimate the probabilities of ancestral ranges (dataset 3). BioGeoBEARS calculates maximum-likelihood estimates of the ancestral states at internal nodes by modeling transitions between geographical ranges along phylogenetic branches as a function of time. BioGeoBEARS encompasses six different biogeographical models (DEC, DEC + J, DIVA, DIVA + J, BAYAREA-LIKE, BAYAREA-LIKE + J) as implemented in the R package BioGeoBEARS 87 . All models entail dispersal (d) and extinction (e) as free parameters. Three models also comprise an additional parameter "j" (+ J) to embrace Table 2. Summary of the calibration points and the prior distribution applied in dating BEAST analyses represented in Fig. 4.  Table 3. Age estimates and reconstructed ancestral ranges for each of the nodes in the chronogram represented in Fig. 4. Bayesian posterior probabilities (BPP > 0.9), 95% highest posterior density (HPD) intervals (millions of years = Ma) based on a relaxed molecular-clock analysis of ITS-ETS sequences in BEAST. Letters correspond to the following ancestral areas or combination of areas with a relative probability (BayArea + J model) equal or greater 0.05. A Alps, B Balkans, C Transcaucasus, E Eurasia, I Iberian Peninsula, N North Africa, R Iran, T Turkey, Y south Italy.