Geo-Climatic Changes and Apomixis as Major Drivers of Diversification in the Mediterranean Sea Lavenders (Limonium Mill.)

The Mediterranean realm, comprising the Mediterranean and Macaronesian regions, has long been recognized as one of the world’s biodiversity hotspots, owing to its remarkable species richness and endemism. Several hypotheses on biotic and abiotic drivers of species diversification in the region have been often proposed but rarely tested in an explicit phylogenetic framework. Here, we investigate the impact of both species-intrinsic and -extrinsic factors on diversification in the species-rich, cosmopolitan Limonium, an angiosperm genus with center of diversity in the Mediterranean. First, we infer and time-calibrate the largest Limonium phylogeny to date. We then estimate ancestral ranges and diversification dynamics at both global and regional scales. At the global scale, we test whether the identified shifts in diversification rates are linked to specific geological and/or climatic events in the Mediterranean area and/or asexual reproduction (apomixis). Our results support a late Paleogene origin in the proto-Mediterranean area for Limonium, followed by extensive in situ diversification in the Mediterranean region during the late Miocene, Pliocene, and Pleistocene. We found significant increases of diversification rates in the “Mediterranean lineage” associated with the Messinian Salinity Crisis, onset of Mediterranean climate, Plio-Pleistocene sea-level fluctuations, and apomixis. Additionally, the Euro-Mediterranean area acted as the major source of species dispersals to the surrounding areas. At the regional scale, we infer the biogeographic origins of insular endemics in the oceanic archipelagos of Macaronesia, and test whether woodiness in the Canarian Nobiles clade is a derived trait linked to insular life and a biotic driver of diversification. We find that Limonium species diversity on the Canary Islands and Cape Verde archipelagos is the product of multiple colonization events followed by in situ diversification, and that woodiness of the Canarian endemics is indeed a derived trait but is not associated with a significant shift to higher diversification rates. Our study expands knowledge on how the interaction between abiotic and biotic drivers shape the uneven distribution of species diversity across taxonomic and geographical scales.


INTRODUCTION
Biodiversity on Earth is unevenly distributed. Species richness varies across habitats, geographic regions, and taxonomic groups, raising long-standing questions about the ecological and evolutionary mechanisms underpinning this variation. Factors that drive speciation and extinction (i.e., species diversification) have altered biodiversity in space and time. Drivers of diversification include abiotic factors, for example, tectonic processes and climatic events (Barnosky, 2001), and biotic factors, for example, reproductive traits, ploidy levels, hybridization, and habit (Soltis et al., 2019). Numerous studies have focused on identifying a single key trait and linking it to shifts in diversification rates (e.g., Mayhew et al., 2008;de Vos et al., 2014;Howard et al., 2020). However, both biotic and abiotic factors can act synergistically toward changes in diversification (e.g., Bouchenak-Khelladi et al., 2015;Donoghue and Sanderson, 2015;Condamine et al., 2018). Analyses of spatio-temporal evolution and drivers of diversification in species-rich lineages are crucial to clarify the role of past biotic and abiotic changes in the origins of species diversity and predict how lineages will be affected by ongoing environmental changes.
Flowering plant diversity is partitioned taxonomically, geographically, and environmentally. Angiosperms comprise more than 13,000 genera (Christenhusz and Byng, 2016) ranging in size from a single to thousands of species, yet only about 70 genera are characterized as species rich (with ≥500 species ;Frodin, 2004;Mabberley, 2017). Furthermore, plant diversity is concentrated in 36 global biodiversity hotspots that cover only 16% of Earth's surface but harbor more than 50% of endemic vascular plants and are undergoing remarkable loss of habitat (Myers et al., 2000;Mittermeier et al., 2011). The Mediterranean has been identified as one of the world's biodiversity hotspots by Mittermeier et al. (2011) andCritical Ecosystem Partnership Fund (2019). These authors defined the Mediterranean hotspot as comprising the Mediterranean region (i.e., mainland areas surrounding the Mediterranean Sea and the islands in it) and Macaronesia (i.e., Canaries, Cape Verde, Azores, Madeira, and Selvagen Islands). The Mediterranean region and Macaronesia combined are also often referred to as the Mediterranean realm (e.g., Comes, 2004;Mansion et al., 2009). This hotspot covers only 1.6% of the Earth's surface, yet accommodates 10% of its total plant species richness, representing the third richest hotspot with approximately 25,000 species, more than half of which are endemic (Medail and Quezel, 1997;Blondel et al., 2010;Critical Ecosystem Partnership Fund, 2019). Its geographic location at the crossroads of three continents (Europe, Africa, and Asia) makes the Mediterranean a large contact zone for taxa of different biogeographic origins (e.g., Eurasian Circumboreal, Irano-Turanian, and Saharo-Arabian), which, together with taxa that originated and diversified in situ, form its remarkably diverse flora (Blondel and Aronson, 1999).
The Mediterranean region has undergone multiple geological and climatic upheavals (Thompson, 2005). Geologically, the region originated from two ancient, independent ocean basins: the Alpine Tethys Ocean (opened during the Middle to Late Jurassic and related to the opening of the Central Atlantic) in the Northwest and the Neotethys Ocean (opened from the Triassic to the Jurassic between Laurasia and Gondwana) in the Southeast (van Hinsbergen et al., 2020 and references therein). From the Cretaceous to the early Miocene, a continuing convergence of tectonic plates brought Europe and Africa progressively closer (Rosenbaum et al., 2002). In the late Miocene, uplift at the continental margins of Iberia and Africa triggered extensive basin desiccation (Duggen et al., 2003;Gargani and Rigollet, 2007). This period, known as the Messinian Salinity Crisis (MSC, from ca. 5.96 to 5.33 Ma; Krijgsman et al., 1999), has been described as "one of the most dramatic events on Earth during the Cenozoic era" (Hsü et al., 1973;Duggen et al., 2003).
The MSC and the onset of the Mediterranean climate (3.2-2.8 Ma; Suc, 1984) were landmark events in the evolution of diversity in the Mediterranean region (Fiz- Palacios and Valcárcel, 2013). The creation of saline deserts during the MSC (Hsü et al., 1973) produced land bridges between islands and continental areas that potentially facilitated migrations of plants with the necessary dispersal properties and salt-tolerance (e.g., halophytes). The MSC is considered to have facilitated speciation in arid-adapted lineages and extinction in sub-tropical Tertiary lineages (Rodríguez-Sánchez et al., 2008;Jiménez-Moreno et al., 2010;Crowl et al., 2015). The refilling of the basin at the end of the MSC disrupted previously formed land bridges, thus promoting vicariance, and mitigated aridity (García-Castellanos et al., 2009), thus possibly causing extinction of arid-adapted lineages (Fiz-Palacios and Valcárcel, 2013). While the effects of MSC on the Mediterranean flora are still debated, the positive effects of the emergence of the Mediterranean climatic regime on diversification are corroborated by multiple studies (e.g., Valente et al., 2010;Fiz-Palacios and Valcárcel, 2011). Furthermore, several plant lineages show a temporal period of reduced diversification rate from the Messinian event to the onset of the Mediterranean climate that has been variably attributed to either mass extinction, rate stasis, or a combination of the two (Fiz- Palacios and Valcárcel, 2013).
During the late Pliocene-early Pleistocene, cooler and dryer conditions were implicated in several extinctions (Bessedik et al., 1984), while Pleistocene glacial cycles and eustatic sea-level changes (2.58-0.01 Ma; Lisiecki and Raymo, 2007) further impacted Mediterranean plant diversification and distributions. Pleistocene geoclimatic oscillations caused species range contractions and expansions as plant populations fragmented and merged in response to the appearance and disappearance of dispersal barriers through time, with contrasting effects on diversification of different Mediterranean plant lineages (Nieto Feliner, 2014). Furthermore, different types of islands (oceanic and continental), substrates, and microclimates provided opportunities for adaptation and speciation driven by both ecological and geographical barriers. Overall, geomorphological and climatic processes, coupled with a long history of human activities in the Mediterranean, created a mosaic of heterogeneous habitats, where a diversity of abiotic factors had profound impacts on diversification.
Frontiers in Plant Science | www.frontiersin.org 3 January 2021 | Volume 11 | Article 612258 In addition to extrinsic environmental factors, inherent biological features also play a key role in diversification. For example, hybridization, polyploidy, plant habit, and reproductive strategies have all been invoked to explain species divergence and eco-geographical differentiation in plants (Rieseberg et al., 2003;Goldberg et al., 2010;Goldberg and Igić, 2012;Soltis et al., 2013). In particular, asexual reproduction via apomixis (i.e., cloning through seeds; Asker and Jerling, 1992) can promote rapid diversification by enabling reproductively isolated genotypes to form rapidly, providing reproductive assurance in the absence of pollinators and/or mates, and offering an escape from sterility in newly formed polyploid hybrid individuals (Baker, 1955;Darlington, 1958;Majeský et al., 2017). Apomixis can enable the establishment of a new population from a single individual, an analogous effect to that of selfing as proposed by Baker (1955). Furthermore, in oceanic island systems, such as Macaronesia, the evolution of woodiness from herbaceous ancestors following island colonization has been proposed as a key driver of insular radiations (Nürk et al., 2019). Thus, elucidating the timing, space, and rates of diversification of Mediterranean and Macaronesian lineages and attempting to correlate them with likely biotic and abiotic drivers is essential to explain the evolution of biodiversity in the Mediterranean hotspot (Comes, 2004).
To achieve a deeper understanding of diversification dynamics in the Mediterranean realm, it is necessary to focus on widely distributed plant groups with high species diversity in this region and biotic traits that could act as triggers of diversification given the unique geo-climatic history of the area. Limonium Mill. (sea lavender; Plumbaginaceae) qualifies as such a group. It is a species-rich genus (ca. 600 species; Koutroumpa et al., 2018) in the top 0.005% of angiosperm genera for size (according to Frodin's, 2004 criteria) and has a cosmopolitan distribution with center of diversity in the Mediterranean (ca. 70% of the species, mostly endemics, occur in the region). Limonium species, characterized as facultative halophytes (i.e., salt tolerant), grow predominantly on saline and metal-rich soils of mainland and coastal areas (Erben, 1978). They represent an important component of coastal vegetation in Mediterranean ecosystems and contribute the dominant species in several vegetation types (e.g., Crithmo-Limonietea class; Brullo et al., 2017). Limonium displays marked variation of chromosome numbers, ranging from 12 to 18 in diploids and from 24 to 72 in polyploids (Erben, 1979;Brullo and Erben, 2016). Sea lavenders can reproduce both sexually and asexually via apomixis, with most sexual species characterized by pollen-stigma dimorphism and sporophytic self-incompatibility (Baker, 1966). Apomixis occurs exclusively in polyploid taxa, some of presumed hybrid origin (Erben, 1979). The combined effects of polyploidy, hybridization, and apomixis have been suggested as having shaped Limonium diversity in the Mediterranean region (e.g., Ingrouille, 1984;Palacios et al., 2000;Lledó et al., 2005). Previous phylogenetic and taxonomic analyses concluded that all described apomictic species of Limonium, together with some sexual species, belong to a single large clade formed by the vast majority of Mediterranean endemics, named the "Mediterranean lineage" by Koutroumpa et al. (2018). The same study placed the Macaronesian endemics in four clades of the Limonium tree, with the majority of them included in the Nobiles and Jovibarba-Ctenostachys clades. The Nobiles clade consists entirely of Canarian endemics with a woody (suffruticose) habit. The Jovibarba-Ctenostachys clade comprises endemics of the Canaries and Cape Verde archipelagos together with endemics in Morocco and Hispaniola (Malekmohammadi et al., 2017;Koutroumpa et al., 2018). A previous phylogenetic study, limited to ca. 8% of Limonium species and based on a single biogeographic calibration, inferred a late Miocene (ca. 6-7 Ma) origin for the Mediterranean clade of Limonium with most diversification during the Messinian and Pliocene (Lledó et al., 2005). The same study inferred the diversification of Macaronesian clades between the late Miocene and Pliocene. Owing to its large size, worldwide distribution, and uneven diversity among regions, Limonium is an ideal group to elucidate the factors shaping the partitioning of biodiversity through space and time, warranting a new study with enhanced sampling and methodology.
Here, we generate and time calibrate the largest phylogeny to date of Limonium and Plumbaginaceae. By reconstructing ancestral areas of distribution, inferring the tempo of diversification, and estimating diversification rate dynamics, we address the following questions. At the global scale we ask: (1) When and where did the genus originate and diversify?
(2) Can we detect significant shifts of diversification rates and link them to specific abiotic (major geological and/or climatic events) and/or biotic factors (apomixis)? More specifically, we test whether diversification rates are constant or heterogeneous in Limonium, and whether potential rate changes are concomitant with changes in extrinsic and intrinsic variables. At the regional scale, we focus on island biogeography, trait evolution, and diversification of Macaronesian Limonium by asking: (1) What are the biogeographic origins of island endemics in Macaronesia?
(2) Did the transition from herbaceousness to woodiness precede or follow island colonization in the Canarian Nobiles clade? Specifically, we hypothesize that woodiness is a derived state linked to insular life, as suggested by several studies (e.g., Carine et al., 2010;Lens et al., 2013), rather than an ancestral state preserved in islands, as proposed by the "islands-as-museums" hypothesis (Cronk, 1997). (3) Did the transition to woodiness trigger an increase in diversification rate in the Canarian Nobiles clade, as found in other insular clades (Nürk et al., 2019)? Our findings highlight the importance of both extrinsic and intrinsic factors in generating the remarkable plant diversity of a global biodiversity hotspot, the Mediterranean realm.

Taxon and Molecular Sampling
We sampled more than one-third of all Limonium species (i.e., 216 taxa representing 214 species: 212 species were identified at the species level, one species was identified at the subspecies level, and one species was represented by three varieties) together with 66 species of 20 other Plumbaginaceae genera.
Frontiers in Plant Science | www.frontiersin.org Additionally, 20 species from the sister family of Polygonaceae (both subfamilies sampled) were used as outgroup taxa for phylogenetic inference (Supplementary Table S1). Our sampling of Limonium includes representatives from all taxonomic subdivisions and phylogenetic clades (Malekmohammadi et al., 2017;Koutroumpa et al., 2018) and spans its cosmopolitan distribution. Sampling of Plumbaginaceae and Polygonaceae provided appropriate nodes for fossil and secondary calibrations. Overall, our molecular dataset comprised 302 taxa and a total of 1,039 sequences from one nuclear (261 ITS sequences) and three chloroplast (291 trnL-F, 243 rbcL, and 244 matK sequences) markers. Here, we increase taxon sampling by 15 taxa compared to the previously published phylogeny of Limonium and Plumbaginaceae of Koutroumpa et al. (2018). Molecular methods for 70 newly generated sequences and alignments followed Koutroumpa et al. (2018).

Phylogenetic and Dating Analyses
Phylogenies were estimated using maximum likelihood (ML) and Bayesian inference (BI) as implemented in RAxML 8.2.10 (Stamatakis, 2014) and MrBayes 3.2.6 (Ronquist et al., 2012), respectively, following conditions described in Koutroumpa et al. (2018; except for BI of supermatrices for which chains ran for 30 million generations). Trees were inferred for each locus separately and in concatenation (concatenated datasets: 3-loci cpDNA matrix and 4-loci "Supermatrix"). We inspected the resulting cpDNA and ITS trees to detect incongruences between well-supported clades (as defined in Koutroumpa et al., 2018) and identify "rogue taxa" (i.e., taxa placed in different wellsupported clades in the cpDNA vs. the ITS tree). Given the small number of conflicts thus identified, we additionally inferred phylogenies for a Supermatrix from which we excluded ITS sequences of "rogue clades and taxa" (hereafter "Supermatrix-cpDNA-like") and a Supermatrix from which we excluded cpDNA sequences of "rogue clades and taxa" (hereafter "Supermatrix-ITS-like"). Trees from these two trimmed supermatrices have higher resolution than those from the three-loci cpDNA matrix and ITS matrix.
Divergence time estimates were performed in BEAST 2.5.1 (Bouckaert et al., 2019) for ITS, cpDNA, "Supermatrix, " "Supermatrix-ITS-like, " and "Supermatrix-cpDNA-like" datasets, using five calibration points. Following a comprehensive review of literature for the Plumbaginaceae-Polygonaceae fossil record and dated angiosperm phylogenies, and following guidelines for nodal assignment of fossils (Magallon and Sanderson, 2001;Rutschmann et al., 2007), we calibrated four nodes (two within Plumbaginaceae and two within Polygonaceae), and the stem node of Plumbaginaceae, forming a total of five calibration points (see Supplementary Figures S1, S2); the details of the calibration process are described below. Although fossils for Plumbaginaceae are sparse, two internal nodes could be calibrated: an Armeria-type pollen fossil from upper Miocene (Van Campo, 1976;Muller, 1981) was used as a minimum age constraint (5.333 Ma) for the Limonieae crown-node (the origin of Armeria-type pollen according to Costa et al. 's, 2019 ancestral state estimates), and a pollen fossil of Ceratolimon cf. feei (former Limoniastrum; Beucher, 1975;Muller, 1981) from the Pliocene was used as minimum age constraint (2.58 Ma) for the Ceratolimon stem-node. We also used two fossils from Polygonaceae: a Coccoloba pollen fossil from the upper Miocene (Graham, 1976;Muller, 1981) was used to assign a minimum age (5.333 Ma) to the Coccoloba stem-node and Muehlenbeckiatype pollen ( †Rhoipites muehlenbeckiaformis; Macphail, 1999) from the upper Eocene was used as minimum age constraint (33.9 Ma) for the Muehlenbeckia stem-node. Finally, we used a secondary calibration for the stem of Plumbaginaceae based on age estimates from the dated angiosperm phylogeny of Magallón et al. (2015) [mean: 67.9 Ma and 95% Highest Posterior Density (HPD): 65.63-78.21 Ma]. We employed uniform priors for fossil calibrations with a hard maximum bound of 78.21 Ma (Magallón et al. 's, 2015 HPD max age for Plumbaginaceae-Polygonaceae split) and minimum bounds mentioned above. Employing the conservative approach of uniform priors, we consider that fossils can provide only minimum ages and that paucity of fossil records hinders the use of more informative priors (similar to other angiosperm studies, e.g., Bouchenak-Khelladi et al., 2014;Boucher et al., 2016). For the secondary calibration, we ran analyses using either a normal (mean age of 67.9 Ma and SD of 5.25 Myr) or a uniform prior (65.63-78.21 Ma), which produced very similar age estimates (Supplementary Table S2). Thus, we present results from analyses with uniform priors following Schenk (2016), who demonstrated reduced error in uniform vs. normal priors for secondary calibrations.
Divergence times were inferred using a relaxed uncorrelated lognormal clock and a nucleotide substitution model-averaging method (bModelTest tool; Bouckaert and Drummond, 2017) for each of the two partitions (ITS vs. cpDNA). bModelTest integrates over 203 time-reversible models while simultaneously estimating other parameters (estimates weighted by the support of each model). Independent runs of 200 million generations were combined (after removing up to 22.5% as burn-in) using LogCombiner 2.5.1, and maximum clade credibility (MCC) trees were constructed in TreeAnnotator 2.5.1 (Drummond et al., 2012). Chain convergence was verified in Tracer 1.7.1 [Effective Sample Size (ESS) >200 for all parameters; Rambaut et al., 2018]. The 95% HPD intervals of inferred ages for the five datasets overlapped with each other (Supplementary Table S2). Thus, we performed all subsequent analyses on the Limonium clade pruned from the "Supermatrix-ITS-like" and "Supermatrix-cpDNA-like" MCC trees (see also justification above), unless otherwise specified. BEAST and MrBayes analyses were performed using the CIPRES Science Gateway (Miller et al., 2010).
Taxa were coded as present or absent in these areas (Figure 1 and Supplementary Table S3). We tested DEC (Dispersal-Extinction-Cladogenesis, Ree and Smith, 2008) and DIVA-like models (Dispersal-Vicariance Analysis, Ronquist, 1997). We chose not to use the BayArea-like model (Landis et al., 2013) since it does not allow for vicariance or "subset-withinarea" speciation (allowed in DIVA-like and/or DEC models), which are both plausible considering species distributions and the broad areas defined in our study.
Time-stratified analyses with dispersal multiplier matrices were also employed to account for differences in dispersal probabilities due to changes in distances and connectivities between areas over time (Supplementary Table S4). In all biogeographic analyses, Macaronesia was allowed as an ancestral range only after 27 Ma, corresponding to the age of the oldest extant volcanic islands, i.e., the Selvagen Islands (Borges et al., 2008;Carine et al., 2010;Triantis et al., 2010;Fernández-Palacios et al., 2011). Considering the uncertainty surrounding the reconstruction of past sea levels (Miller et al., 2005;Rohling et al., 2014), we chose the age of 27 Ma as the most widely accepted, hence conservative estimate for the oldest emerged island in Macaronesia. Unlikely disjunct ranges (i.e., non-neighboring, geographically unconnected areas over the inferred existence of Limonium, starting ca. 40 Ma) were removed from the list of possible states. All models were also run with j parameter, which accounts for founder-event speciation. However, due to criticism by Ree and Sanmartín (2018) concerning "degenerate" inferences of +j models and inappropriate statistical comparisons of models with and without j, and considering that our inferences with and without j parameter closely match each other, we present only results of models without j and compare their likelihoods using AIC and AICc. We additionally estimated the type and number of biogeographical events using Biogeographical Stochastic Mapping (BSM; Matzke, 2016;Dupin et al., 2017) in BioGeoBEARS 1.1.2, providing the best-fitting model for each MCC tree and summarizing event frequencies over 200 discrete stochastic maps.
At a finer biogeographic scale, we investigated the origins of island endemics in Macaronesia and analyzed a clade (Jovibarba-Ctenostachys clade) comprising endemics in the Canaries, Cape Verde, Morocco, and Hispaniola in the Caribbean (four coded areas). Since the Jovibarba-Ctenostachys clade shows topological differences between datasets, we pruned it from cpDNA, ITS and "Supermatrix" MCC trees, ran DEC, DIVAlike, and BayArea-like models and compared results for all three trees. In contrast to the genus-wide biogeographic analyses, we included also the BayArea-like model for the biogeographic analyses of the Jovibarba-Ctenostachys clade because the absence of vicariance and "subset-within-area" speciation of this model was not considered a limitation for a lineage consisting exclusively of endemics, most of which are insular. Their biogeographic history could only be explained by dispersal, extinction, and sympatric speciation, all of which are accounted for in the BayArea-like model.

Diversification Analyses
We estimated speciation and extinction rate dynamics of Limonium applying recently developed methods to both "Supermatrix-ITS-like" and "Supermatrix-cpDNA-like" MCC trees. First, we tested for shifts in diversification rates along tree branches and through time, using Bayesian Analysis of Macroevolutionary Mixture (BAMM 2.5; Rabosky et al., 2014). Sampling fractions per Limonium clade/section (following the classification of Koutroumpa et al., 2018: Supplementary Table S3) were provided in order to account for missing species (Supplementary Table S5). The R package BAMMtools 2.1.6 (Rabosky et al., 2014) was used to choose appropriate priors for BAMM analyses, analyze, and visualize results. BAMM was run with four independent chains of 100 million generations each and convergence was confirmed with ESS >200 for all parameters. The absence of recent, comprehensive monographic treatments for Limonium prevents us from performing analyses that account for potential discrepancies among different taxonomic treatments (e.g., Faurby et al., 2016). Moreover, taxonomic uncertainties could especially affect the apomictic species of the "Mediterranean lineage" (Haveman, 2013;Majeský et al., 2017). Acknowledging that taxonomic uncertainties could affect estimates of diversification rate dynamics (Faurby et al., 2016;Fernández-Mazuecos et al., 2019), we ran preliminary sensitivity analyses in BAMM, where we arbitrarily assumed that the number of species for the "Mediterranean lineage" is only half of the currently accepted number of "Mediterranean lineage" species (i.e., ca. 200 vs. 400 species), and thus assigned a higher sampling fraction to this lineage to account for missing species in analyses of both MCC trees (i.e., ½ instead of ¼; the former being the fraction used in the preliminary analyses). The signal in our data was so strong that even in these preliminary analyses that underestimated the number of species in the "Mediterranean lineage, " a shift in rates was recovered for this lineage in both MCC trees, despite the extreme and rather unrealistic taxonomic scenario.
Second, we tested for episodic tree-wide rate shifts caused, for example, by a mass extinction event, such as the one hypothesized for the temporal gap of diversification from the MSC to the onset of the Mediterranean climate in other Mediterranean plant lineages. We used the CoMET algorithm implemented in the R package TESS 2.1 (Höhna et al., 2015) to test for signals of mass extinction events and for shifts in speciation and extinction rates over time. Changes in speciation and extinction are modeled to occur at discrete time points and affect all clades in a tree simultaneously (episodic birthdeath). We ran models assuming one or two rate shifts and allowing mass extinction events. A global sampling fraction of Limonium was provided to account for sampling probability of extant taxa ("Supermatrix-ITS-like": 0.35 and "Supermatrix-cpDNA-like": 0.36). We also used the automatic stopping functionality and specified a minimum ESS of 1,000 to allow MCMC simulations to reach convergence.
Third, we fitted different diversification models with constant or varying speciation and/or extinction rates through time and in relation to paleoenvironmental variables (Morlon et al., 2010;Condamine et al., 2013) using the R package RPANDA 1.6 (Morlon et al., 2016). Analyses were repeated using different initial parameters and choosing the solution with the highest likelihood for each model. Models were compared using AICc Frontiers in Plant Science | www.frontiersin.org 6 January 2021 | Volume 11 | Article 612258 and a multimodel inference approach was employed for parameter estimates of time-dependent analyses (Burnham and Anderson, 2002). According to Morlon et al. (2011), rate heterogeneity can affect diversification rate estimates if such heterogeneity is not accounted for. Thus, apart from fitting models to the entire tree, RPANDA allowed us to test and account for rate heterogeneity by fitting models to subclades for which shifts were inferred in BAMM analyses (see section Results) and to the remaining phylogeny after removing these subclades. Models with or without rate shifts were then compared using the AICc. Global or clade-specific sampling fractions were assigned in analyses of either the entire phylogeny or the subclades defined above (in heterogeneity analyses), respectively, to account for incomplete taxon sampling. Paleoenvironmental variables used in RPANDA analyses as potential correlates of diversification rates were (i) global temperature throughout the Cenozoic (Zachos et al., 2008), fitted to the entire tree or all parts of the tree in heterogeneity analyses and (ii) past sea-level estimates (Miller et al., 2005;Rohling et al., 2014), fitted to Mediterranean subclade(s) that comprise coastal species (heterogeneity analyses; see environmental curves in Supplementary Figure S3). Specifically, the recent sea-level reconstructions of Rohling et al. (2014) in the Mediterranean are limited to the past 5.3 Myr and thus were fitted to the temporally corresponding part of the Mediterranean subclade, requiring that a few early-diverging taxa were removed from the subclade.

State-Dependent Diversification
We tested for the effect of reproductive system on diversification rates using several state-dependent speciation and extinction (SSE) models. We assembled available information on reproductive strategies, coding 175 Limonium taxa (out of the 216 taxa included in the trees) as sexual or apomictic based on the literature (see Supplementary Table S6). Tips with missing reproductive information were dropped from "Supermatrix-ITS-like" and "Supermatrix-cpDNA-like" MCC trees before analyses. First, we used BiSSE (Maddison et al., 2007) implemented in the R package diversitree 0.9-11 (FitzJohn et al., 2009;FitzJohn, 2012) to evaluate 10 models in which speciation, extinction, and transition rates varied or remained equal between states; two of these models had zero extinction (pure-birth models). Since we lack information on reproductive system for many Limonium species, we conducted three analyses with three different sampling fractions for each state (sexual vs. apomictic) to account for missing species in the phylogeny and tested the robustness of the results. In the first analysis, we assumed that the same proportion of sexuals vs. apomicts sampled in the entire Limonium phylogeny also applies to unsampled taxa (i.e., global sampling fraction). In the other two analyses, we incorporated prior knowledge that all described apomicts are included in the "Mediterranean lineage" (while species in all other lineages are sexual) and assumed either 50-50 or 40-60% of sexuals vs. apomicts in this lineage; we also accounted for uneven sampling in major lineages of Limonium by using the "make.bisse.uneven" function. We consider the latter sampling scenario as more realistic, since the predominance of apomicts over sexuals in the Mediterranean has been documented by Brullo and Erben (2016). Finally, we accounted for confidence intervals of parameter estimates by running a Bayesian MCMC with 50,000 steps using an exponential prior (following FitzJohn, 2012) and the best-fitting model for each analysis. Chain convergence was checked using the R package coda (Plummer et al., 2006) and results were plotted after removing 5,000 steps as burn-in.
Second, due to criticism of BiSSE concerning incorrect assignments of rate differences for neutral observed states in phylogenies, where rate heterogeneity is caused by other factors (Maddison and FitzJohn, 2015;Rabosky and Goldberg, 2015), we applied the Hidden SSE (HiSSE; Beaulieu and O'Meara, 2016) model to specifically account for other unmeasured factors impacting diversification rates. We tested seven models (model details in Beaulieu and O'Meara, 2016): a full BiSSElike model, two HiSSE models with equal or varying transition rates (one rate for transitions among hidden states and two rates for transitions between sexual reproduction and apomixis), and four character-independent diversification (CID) models assuming that diversification is trait-independent but not constant across the tree (CID-2 and CID-4 models match complexity of BiSSE and HiSSE models, respectively). We also used the same three sampling scenarios for sexuals vs. apomictic states explained above, but without assigning uneven lineage-specific sampling probabilities (not implemented in HiSSE). Furthermore, we inferred the model-averaged rate estimates for all tips in both MCC trees and for all analyses to examine whether diversification rates varied between sexuals and apomicts. In BiSSE and HiSSE analyses, the tested models were compared with AIC.
Third, we used FiSSE (Rabosky and Goldberg, 2017) as a nonparametric test for reproductive system-dependent diversification in Limonium. This test is used for additional validation of results obtained from the other two state-dependent models (Rabosky and Goldberg, 2017). FiSSE shows the lowest "false positive" rates compared to BiSSE and HiSSE, but its power in detecting state-dependent diversification is much lower than BiSSE for trees with <300 tips (Rabosky and Goldberg, 2017). Although our trees are relatively small for the power limits of FiSSE, we expect that a significant state-dependent FIGURE 1 | Time calibrated phylogeny and biogeography of Limonium based on "Supermatrix-ITS-like" dataset. The nine areas used for biogeographic inference are color coded as in the world map and corresponding list in the upper left inset. Colored squares on nodes indicate the most likely ancestral areas; due to space limitations, some nodes have no colored squares: in these cases, ranges are the same as in their respective ancestral nodes. Triangles and X's along branches indicate range expansions and contractions, respectively, that occurred when geological events significantly impacted Limonium's distribution and diversification (i.e., in the early Oligocene, early to middle Miocene and during the Messinian). These periods are marked by vertical color bands. In two of these periods, namely in the early Oligocene and the early to middle Miocene, paleomaps show the Mediterranean geographic configuration in each time frame (paleomaps are from the article of Rögl, 1999 published in Geologica Carpathica journal under the CC BY-NC-ND 4.0 license). Uncertainties in ancestral range estimates are in Supplementary Figure S4. diversification, if detected, would show that the signal of our data is strong enough to overcome FiSSE's power limitations. Conversely, finding non-significant state-dependent diversification does not allow the possibility that this result is due to FiSSE's power limitations to be excluded, especially when the other two methods support state-dependent diversification. Here, we did not account for incomplete taxon sampling because this is not implemented in FiSSE.
Since most Limonium species occur in the Euro-Mediterranean area, we also tested the effect of this range on diversification rates using the fast Hidden Geographical SSE model (fGeoHiSSE; Caetano et al., 2018). Limonium species were coded as present only in the Euro-Mediterranean area or only in other areas or widespread in both the Euro-Mediterranean and other areas. Missing species were accommodated by providing sampling fractions for each of these three range categories. We fitted four models: two range-independent and two range-dependent diversification models, each either with or without a hidden area. Models were compared via AIC.
At shallower phylogenetic levels, we pruned the Limonium subgen. Pteroclados s.l. lineage (clade A in Figure 1; see also Koutroumpa et al., 2018) from "Supermatrix-ITS-like" and "Supermatrix-cpDNA-like" MCC trees, and tested whether woodiness (i.e., suffruticose habit) of Canarian Limonium endemics in Nobiles clade Bramwell, 1974, 1990;Mesa et al., 2001;Marrero and Almeida, 2003;Kunkel, 2012;Jiménez et al., 2017) had an impact on diversification rates. In this fully sampled Pteroclados s.l. lineage, all other clades apart from Nobiles consist of herbaceous species. We used BiSSE and HiSSE analyses, with models described above, to test for habit-dependent diversification in the Pteroclados s.l. lineage. Simulation studies have suggested that while it is feasible to use SSE methods for analyses of small-size clades such as Pteroclados s.l. (e.g., Gamisch, 2016), their power in detecting state-dependent diversification depends on the strength of speciation rate asymmetry in the tree. If state-dependent diversification were detected, this would mean that speciation rate asymmetry between woody and herbaceous taxa is sufficiently high (i.e., >2.5-fold) to be detected despite small clade size, as demonstrated by simulations (Gamisch, 2016;Kodandaramaiah and Murali, 2018). Conversely, if state-independent diversification were supported, this would not automatically exclude the existence of moderate levels of speciation rate asymmetry (i.e., <2.5-fold), but might simply mean that such levels are insufficient to be recovered by the SSE methods in small clades. We also tested whether woodiness is a derived state for the Pteroclados s.l. lineage employing model-averaged ancestral state estimations implemented in HiSSE package.
Divergence time estimates were almost identical between "Supermatrix-ITS-like" and "Supermatrix-cpDNA-like" MCC trees (detailed age estimates and HPDs in Table 1;  Supplementary Table S2 and Supplementary Figures S1, S2). The origin of Plumbaginaceae is placed between the Upper Cretaceous and the early Eocene (95% HPD: 52-77 Ma), with a median age at early Paleocene (ca. 65.7 Ma). The median crown ages for Plumbaginoideae and Limonioideae are ca. 29.5 Ma (95% HPD: ca. 18-43 Ma) and ca. 57 Ma (95% HPD: ca. 43-71 Ma), respectively. Combining divergence time and biogeographic estimates for Limonium (DEC without dispersal-multiplier matrices selected as best-fitting model; Figure 1; Supplementary Table S7 and Supplementary Figures S4, S5), we infer a late Paleogene origin (ca. 33 Ma median crown age; 95% HPD: ca. 24-44 Ma) in the proto-Mediterranean region (i.e., Euro-Mediterranean and North-Africa|North Africa: most likely states, yet with low probability) for the genus. The two subgenera Pteroclados s.l. and Limonium originated during Miocene, ca. 12 Ma (95% HPD: ca. 6-20 Ma) and ca. 21 Ma (95% HPD: ca. 15-29 Ma), respectively, within widespread regions. A widespread range was also inferred (with low to moderate probabilities) for the most recent common ancestor (MRCA) of clades B2 (comprising mostly non-Mediterranean species) and B3 ("Mediterranean lineage") in the middle Miocene (ca. 15.5 Ma; 95% HPD: ca. 11-21 Ma). This widespread ancestor gave rise to an endemic Euro-Mediterranean MRCA of the "Mediterranean lineage" (B3 crown-node; ca. 12 Ma) and either an endemic Irano-Turanian MRCA for the mostly non-Mediterranean B2 clade according to "Supermatrix-ITS-like" tree (ca. 14.4 Ma; Figure 1; Supplementary Figure S4) or a widespread MRCA for the B2 clade in the same ancestral widespread range according to "Supermatrix-cpDNA-like" tree (ca. 14 Ma; Supplementary Figure S5). In clade B2, dispersals, vicariance events and/or range contractions from middle Miocene through Plio-Pleistocene produced ancestors with endemic ranges in the Irano-Turanian, Arabia-NE Africa, East Asia-Australia, South Africa, America, and Macaronesia. Early divergence events of the "Mediterranean lineage" (clade B3) in late Miocene-early Pliocene were accompanied by dispersal from Euro-Mediterranean to North Africa, while the origin of the larger subclade of B3 (ca. 6 Ma median crown age) coincides with the MSC. Extensive speciation within the Euro-Mediterranean area during the Pleistocene was inferred for clade B3, with only little dispersal from Euro-Mediterranean to Macaronesia, North-Africa, Circumboreal and South Africa. Some of these biogeographic results should be interpreted with caution given the limited phylogenetic resolution within the "Mediterranean lineage. " Detailed nodal biogeographic reconstructions and uncertainties in ancestral range estimates for both trees are in Supplementary Figures S4, S5.
BSM analyses estimated 191|196 (SD: ±8| ± 7.79) withinarea speciation events, 67|63 (±2.71| ± 2.06) range expansion events and 16|18 (±2.18| ± 1.71) vicariance events. The large number of sympatric speciation events (ca. 70% of all biogeographic events) was expected, given the large size of the defined areas. Most dispersal events occurred from Euro-Mediterranean to North Africa (ca. 14) and from Euro-Mediterranean to Circumboreal (ca. 5), while all other combinations of areas had <3.5 events (Figure 2;  Supplementary Table S8). Overall, the major source area for dispersals was the Euro-Mediterranean (ca. 40% of estimated dispersal events), followed by North Africa (ca. 15%) and Irano-Turanian (ca. 14%), while the most common sink was North Africa (ca. 28%), followed by Circumboreal (16%). The longest-distance dispersals occurred from Macaronesia to the Americas, and from Arabia-NE Africa and North Africa (or Euro-Mediterranean region) to South Africa (Supplementary Table S8). Movements between areas were largely asymmetrical for most pairs of areas. For example, dispersals from Euro-Mediterranean to North Africa were ca. eight times more common than those in the opposite direction.
Focusing on island biogeography in Macaronesia, origins of island endemics included both neighboring continents and archipelagos (Supplementary Figures S4, S5; Figure 3). Macaronesian endemics were placed in five distinct clades, indicating at least five independent long distance dispersal (LDD) events (Supplementary Figures S4, S5). The Nobiles clade, comprising 16 Canarian endemics, originated ca. 2.3 Ma, while colonization of the Macaronesian area predated the origin of this clade and occurred from a widespread ancestor (Euro-Mediterranean and North Africa|Euro-Mediterranean, North and South Africa, Circumboreal and Arabia-NE Africa) between late Oligocene and early Pliocene. The Canarian endemic L. dendroides (clade B1) represents an early diverged lineage (ca. 21 Ma) that originated via dispersal from neighboring areas (i.e., North Africa or Euro-Mediterranean). The same dispersal pattern holds true for the closely related L. bollei and L. lowei (nested within "Mediterranean lineage") endemic to Canaries and Madeira, respectively, which diversified in the late Pleistocene. The Azorean endemic (L. eduardi-diasii) originated in the late Pleistocene in a clade comprising American endemics, either via a recent colonization from the Americas ("Supermatrix-cpDNA-like") or an earlier colonization (late Miocene-Pliocene) following range expansion from Irano-Turanian toward Euro-Mediterranean, Circumboreal, Macaronesia, and the Americas ("Supermatrix-ITS-like"). Fine-scale biogeographic estimations for the Jovibarba-Ctenostachys clade (DIVA-like selected as best-fitting model; Figure 3; Supplementary Table S9) that diversified during Plio-Pleistocene revealed two independent colonizations of Cape Verde (one probably from North Africa and the other from Canaries) and a LDD event from Cape Verde to Hispaniola. In this clade, the biogeographic origin of Canarian endemics is either from North Africa or Cape Verde and for Moroccan endemics either from Canaries or Cape Verde (inferences varied due to unresolved relationships). BSM reconstructs an overall total of seven dispersal events to Macaronesia (giving rise to native Macaronesian species as well as endemics) including two from North Africa and two from the Euro-Mediterranean. Furthermore, Macaronesia is inferred to be the source pool for four colonizations of continental areas, including one to the Americas and one to North Africa (Figure 2; Supplementary Table S8).

Diversification Dynamics
BAMM analyses strongly reject constant diversification along clades and through time for Limonium (Supplementary Figure S6). Instead, they reveal significant shift(s) toward higher net diversification rates within the "Mediterranean lineage" (Figure 4;  Supplementary Figure S7). In both trees, there is a shift close to the origin of the larger subclade of this lineage (crown age of subclade: ca. 6 Ma), while an additional shift within the 1 | Divergence time estimates for major clades in the "Supermatrix-ITSlike" and "Supermatrix-cpDNA-like" maximum clade credibility (MCC) trees.

Median ages and 95% Highest Posterior Density (HPD) intervals (in parentheses) for the crown-nodes of clades are in million years (Ma). Dating results are from analyses with uniform priors assigned to Plumbaginaceae stem-node (secondary calibration). * Coding of clades within Limonium follow
Frontiers in Plant Science | www.frontiersin.org FIGURE 2 | Summary of dispersal events from Biogeographical Stochastic Mapping (BSM) analyses. Arrows and their thickness represent the direction and frequency of dispersals between areas, respectively. Only event counts with a mean of ≥0.9 in both trees ("Supermatrix-ITS-like" and "Supermatrix-cpDNAlike") are depicted as arrows (detailed results of BSM analyses in Supplementary Table S8). Gray arrows represent dispersal events with a mean ≥0.9 in only one of the two trees. The nine areas used in BSM analyses are color coded as in the world map in the upper left inset (same color coding of areas as in Figure 1). AM, America; AR, Arabia-NE Africa; CB, Circumboreal; EA, East Asia-Australia; EM, Euro-Mediterranean; IT, Irano-Turanian; MA, Macaronesia; NA, North Africa; SA, South Africa. smaller subclade for "Supermatrix-cpDNA-like" tree is explained by the conflicting topologies obtained from plastid and nuclear genomes for some Aegean species. Further evidence supports an 18-20-fold increase in net diversification rate for Mediterranean subclades compared to background rate (1.07 vs. 0.06|1.18 vs. 0.06). The rate-through-time plots show an increase in net diversification rates starting ca. 6 Ma that intensifies during the past 2-3 Myr (Figure 4).
TESS analyses found no significant mass extinction event or episodic change of speciation or extinction rates that could concurrently affect all clades of Limonium (Bayes Factors: 2lnBF<6; Supplementary Figure S8). RPANDA analyses of rate heterogeneity further corroborated results from BAMM. Models with rate shift(s), allowing distinct patterns of rate variation for the "Mediterranean lineage" subclade(s), were strongly supported over models assuming no rate shift for Limonium (ΔAICc>30; Table 2). FIGURE 4 | Diversification dynamics of Limonium across clades and through time. Results of BAMM analyses for "Supermatrix-ITS-like" and "Supermatrix-cpDNAlike" trees are presented in the top and bottom row, respectively. The first figure at the left of each row shows rate heterogeneity across clades for Limonium (phylorate plot) and the other three figures show the net diversification rates through time (rate-through-time plots). In the phylorate plots, branches are colored from blue to red to indicate low to high speciation rates, respectively, and the best rate-shift configuration is denoted by a red circle (the 95% credible set of rate shift configurations are in Supplementary Figure S7). The rate-through-time plots show the net diversification rates for the entire Limonium phylogeny ("Entire tree"), the Mediterranean clade(s) for which a shift was detected by BAMM ("Med clade(s)"), and the remaining Limonium phylogeny after removing the Mediterranean clade(s) ("Non-Med tree"). Significant shift(s) in diversification rates were found within the "Mediterranean lineage" of Limonium, while the steep increase in net diversification rates over the past 2-3 Myr is concomitant with the Plio-Pleistocene climatic oscillations and the establishment of the Mediterranean climate. MSC, Messinian Salinity Crisis.
We inferred a ca. 30-fold higher net diversification rate for the larger "Mediterranean lineage" subclade compared to the background rate (1.21 vs. 0.04|1.11 vs. 0.04; Table 2). For this subclade, a model of constant speciation and exponentially decreasing extinction over time ("Supermatrix-ITS-like") or a pure-birth model with exponentially increasing speciation over time ("Supermatrix-cpDNA-like") had the lowest AICc values (Supplementary Table S10). In the slightly reduced version of this subclade (with Rohling et al. 's, 2014 sea-level data fitted), the highest support was assigned to a model with constant speciation and extinction positively correlated to past sea-level for "Supermatrix-ITS-like, " while for "Supermatrix-cpDNA-like" that the model was the second best-fitting after a pure-birth model (Supplementary Table S10). For the backbone tree [i.e., remaining tree after removing "Mediterranean lineage" subclade(s)], a constant birth-death model received the lowest AICc values, while for the entire tree a model with exponentially varying speciation and extinction, both positively correlated to paleo-temperature, was best-fitting. In all analyses, we observed model uncertainty, with several alternative models explaining diversification dynamics (i.e., ΔAICc<2; Supplementary Table S10). Paleoenvironment-dependent models always received some support, with past temperature and sea-level having similar impacts on speciation and extinction rates.
BiSSE, HiSSE, and FiSSE analyses consistently supported the reproductive strategy-dependent diversification for Limonium, with apomixis leading to higher speciation and net diversification rates. In BiSSE, a model with varying speciation and transition rates and equal extinction rates between sexual reproduction and apomixis received the highest support in all analyses (i.e., for both MCC trees and different sampling scenarios), except for one analysis that supported a full BiSSE model (Supplementary Table S11).
In most BiSSE analyses, up to two additional models received some statistical support (i.e., ΔAIC<2), yet all these alternative models supported the state-dependent diversification. Bayesian parameter estimations revealed higher speciation and net diversification rates for apomixis, and higher transition rates from apomixis to sexual reproduction than in the opposite direction ( Figure 5). Additionally, higher extinction rate for apomixis was observed in the analysis supporting a full BiSSE model (Figure 5; Supplementary Table S11). In HiSSE, either a BiSSE-like model or a HiSSE model with three varying transition rates (see Section Materials and Methods) was supported in all analyses (Supplementary Table S12). Model-averaged parameter estimates for species in our trees recovered higher mean speciation, extinction, and net diversification rates for apomicts across all trees and analyses ( Table 3). Given that a precise estimation of extinction from phylogenies of extant taxa is challenging (Pyron and Burbrink, 2013; but see Morlon et al., 2011 showing inference of realistic extinction rates), we focus our discussion on the general patterns rather than the estimates of extinction rates per se that need to be taken with caution. Complementary FiSSE analyses reported a ca. 2-fold higher speciation rate for apomixis (2.2|1.9 vs. 1.2|1.3), with the difference in rates being either significant (p = 0.016; "Supermatrix-ITS-like" tree) or only marginally non-significant (p = 0.0599; "Supermatrix-ITS-like" tree). Conversely, rangedependent diversification analysis (fGeoHiSSE) found no effect of presence in the Euro-Mediterranean area on diversification rates. Instead, analyses supported a model with diversification rate shifts across the tree, yet independent of the ranges (Supplementary Table S13).

(A) Likelihoods and AIC comparisons between a model assuming no shifts across the tree ("No rate shift") and a model assuming one and two rate shifts for "Supermatrix-ITS-like" and "Supermatrix-cpDNA-like" tree, respectively ("Heterogeneity analysis"); (B) Multimodel averaged speciation, extinction and net diversification rate estimates of time-dependent diversification analyses on the entire tree, the Mediterranean subclade(s) ("Med subclade" and "Med1 subclade": the larger Mediterranean subclade in "Supermatrix-ITS-like" and "Supermatrix-cpDNA-like" trees, respectively, for which a shift was detected by BAMM; "Med2 subclade": the small Mediterranean subclade for which a shift was detected by BAMM in "Supermatrix-cpDNA-like" tree) and the tree after removing the Mediterranean subclade(s) ("Non-Med tree").
on diversification rates (Supplementary Tables S14, S15). BiSSE analyses reported the lowest AIC for the constant rates model, although some state-dependent models with varying rates between woodiness and herbaceousness received partial support (i.e., ΔAIC<2). This result could be attributed to state-independent rate heterogeneity within Pteroclados s.l. clade, which BiSSE cannot test/account for, as revealed by HiSSE (i.e., the highest support for the CID-2 model with equal transition rates; Figure 6). However, caution should be taken when interpreting these results due to known limitations of SSE models in analyses of small-size clades such as the Pteroclados s.l. (see sections Materials and Methods and Discussion). Model-averaged ancestral state estimations combined with biogeographic results inferred woodiness as a derived state linked to insular life (Figure 6).

DISCUSSION
Identifying biotic and abiotic factors that drive the compartmentalization of biodiversity is central to biology (Benton, 2009(Benton, , 2015. With its large size (ca. 600 species), wide geographic range (six continents), and species richness concentrated in the Mediterranean region (ca. 70% of all species), Limonium is a paramount example of uneven taxonomic and geographic plant diversity. Our macroevolutionary study on sea lavenders is one of the few to test the impact of both species-extrinsic and -intrinsic factors on diversification by using carefully selected methods of biogeographic and diversification-rate analyses (Lagomarsino et al., 2016;Condamine et al., 2018;Letsch et al., 2018), thus expanding the knowledge of the tempo and mode of plant diversification. We found significant increases of diversification rates associated with both geo-climatic events (namely, the MSC, onset of Mediterranean climate, and Pleistocene sea-level fluctuations) and intrinsic species traits, such as apomixis. The significant role of apomixis in shaping the Mediterranean radiation of Limonium had been previously proposed and is here explicitly tested for the first time. Our study provides new insights into the origins of Mediterranean biodiversity and highlights a significant role and interplay of both biotic and abiotic factors in promoting species diversification.

Geological and Climatic Processes Shaping Species Distributions and Diversification
Our biogeographic and dating analyses of the largest Limonium phylogeny to date, combined with paleo-geological, paleo-climatic, and paleo-vegetational evidence (e.g., Rögl, 1999;Böhme, 2003; FIGURE 5 | Diversification rate analyses in Limonium dependent upon apomictic vs. sexual reproduction. Results for "Supermatrix-ITS-like" and "Supermatrix-cpDNA-like" trees are presented in the top and bottom rows, respectively. (A) Phylogenetic distribution of sexual reproduction and apomixis. (B) Bayesian parameter inferences from BiSSE analyses employing the best fitting model under the sampling scenario that assumes a 40-60% of sexuals vs. apomicts in the "Mediterranean lineage." For "Supermatrix-ITS-like" tree, a model with equal extinction rates was the best (see the single parameter estimate for both reproductive modes in the Extinction Rate graph in the top row) and for "Supermatrix-cpDNA-like" tree, a model with all rates varying was the best (Supplementary Table S11). The results show that speciation and net diversification rates are higher for apomictic than sexual taxa in both trees, extinction rates are also higher for apomicts but only in "Supermatrix-cpDNA-like" tree, while transitions from apomixis to sexual reproduction were more common that in the opposite direction. Thompson, 2005;Pound and Salzmann, 2017), allow us to propose likely processes that shaped the spatio-temporal evolution of the genus. Our analyses support the origin of Limonium in the late Paleogene (ca. 33 Ma) most probably in the proto-Mediterranean region ( Figure  1; Table  1; Supplementary Figures S4, S5). During this period, major tectonic activities changed the configuration of Eurasia. Europe, formerly an archipelago during the Eocene, became more continental, and the Mediterranean Sea separated from the newly formed intercontinental Paratethys Sea (Rögl, 1999 ; Figure 1). The newly created land corridors may have facilitated the initial range expansion of Limonium from the Euro-Mediterranean into the Irano-Turanian region (Figure 1). Furthermore, the origin of Limonium in the early Oligocene is marked by a climatic transition from sub-tropical to more seasonal conditions along the Euro-Mediterranean coasts. At that time, vegetation representatives of seasonal climates (i.e., warm temperate conifer forests and sclerophyll woodland) spread from the Iberian Peninsula to the coasts of France, and xerophytic shrubland occurred for the first time in the central Mediterranean coasts partially replacing subtropical paleo-biomes formerly dominating the Mediterranean area (Pound and Salzmann, 2017).
The divergence of the "Mediterranean lineage" (clade B3) of Limonium from its mostly non-Mediterranean sister lineage (clade B2) during the Miocene (Table 1; Figure 1) follows a temporal pattern of divergence found in many other Mediterranean angiosperm lineages (Vargas et al., 2018) and is concomitant with paleo-geological and -climatic changes. These two largest lineages of Limonium, clades B2 and B3, diverged during the middle Miocene (ca. 15.5 Ma) from an ancestor that was widespread in the Euro-Mediterranean, Irano-Turanian, and other regions (Figure 1; Supplementary Figure S5). Prior to this split, in the early Miocene and specifically late Burdigalian, the Arabian plate that was connected to Africa collided with the Anatolian plate for the first time, cutting off the seaway between the Mediterranean Sea and the Indian Ocean (Figure 1). In the early Langhian and concomitant with the divergence between B2 and B3 clades, a marine transgression flooded the entire Mediterranean and Paratethys causing the re-opening of the seaway between South Anatolia and Arabia and the formation of other seaways in Eastern Anatolia and probably also in the Aegean or Western Anatolian regions (Rögl, 1999 ; Figure 1). These short-lived seaways may have enabled vicariant speciation followed by range contraction giving origin to a Euro-Mediterranean ancestor for clade B3 and an Irano-Turanian ancestor for clade B2 ("Supermatrix-ITSlike" tree; Figure 1), or peripatric speciation (i.e., temporary peripheral isolate in Euro-Mediterranean) giving origin to a Euro-Mediterranean ancestor for clade B3 and a widespread ancestor in the same ancestral range for clade B2 ("Supermatrix-cpDNA-like" tree; Supplementary Figure S5). Moreover, the early Langhian split between B2 and B3 clades occurred during the Mid-Miocene Climatic Optimum (i.e., a global warming event at ca. 17-15 Ma; Böhme, 2003) and specifically in a period of increased seasonality in Central Europe and the Mediterranean area (up to 6 dry months; Böhme, 2003;Thompson, 2005). At that time, the Mediterranean flora began to resemble the current vegetation due to the extinction of several tropical and subtropical elements in the region (Thompson, 2005).
The diverse Mediterranean flora is composed of a combination of elements that immigrated to the region from other areas, and those that arose in situ (Thompson, 2005). Elucidating the relative roles of repeated colonization vs. in situ speciation is necessary to understand the processes that shaped Mediterranean diversity. Our results show that the great majority of Euro-Mediterranean Limonium species have an autochthonous origin (ca. 100 sympatric speciation vs. ca. 7 dispersal events; see BSM results above), supporting Thompson's (2005) hypothesis that in situ diversification contributed the main component of the heterogeneous Mediterranean flora. Additionally, the Euro-Mediterranean area received a few Limonium species, mostly from the neighboring Irano-Turanian, North African, and Circumboreal regions, consistent with previous studies on other taxa (e.g., Mansion et al., 2008Mansion et al., , 2009Salvo et al., 2010;Manafzadeh et al., 2014). The Euro-Mediterranean area also served as a major source of dispersals to other areas (see also tribe Antirrhineae: Gorospe et al., 2020), especially North Africa (Figure 2). Accordingly, diversity in this latter region was shaped mostly by dispersals (ca. 18 events, of which 14 from the Euro-Mediterranean area) and less by in situ speciation (ca. 3 events). Strong dispersal directionalities between areas (Figure 2) were also documented in Solanaceae by Dupin et al. (2017), who found a close relationship between species richness of an area and the number of dispersals from this area, a pattern congruent with our results. The inferred higher number of species dispersals from the Euro-Mediterranean "source" to the North African "sink" than to any other area could be explained by geoclimatic events that connected these two areas in the past, and by the availability of suitable habitats on their coastlines. Overall, the directional asymmetry of dispersals between the Euro-Mediterranean and other areas, consisting in more emigration than immigration to the region, suggests that extensive in situ speciation is the pump that promotes species movement from one main area (here the Euro-Mediterranean) to many others. Proposed abiotic landmarks in the evolution of Mediterranean biodiversity include the Messinian Salinity Crisis (ca. 6-5.33 Ma), the onset of the Mediterranean climate (3.2-2.8 Ma), and Pleistocene geo-climatic fluctuations (Fiz-Palacios and Valcárcel, 2013). Indeed, we found that these events played a key role in the diversification of the "Mediterranean lineage" of Limonium. Although Mediterranean sea lavenders originated already in the middle Miocene (ca. 12 Ma; Figure 1; Table 1), all major extant lineages started to diversify only about 6 Ma concomitantly with the onset of the MSC, as a product of a significant, clade-specific shift toward higher diversification rates in the largest subclade of the "Mediterranean lineage" (phylorate plots; Figure 4). Thus, the salinity crisis that massively increased the availability of saline habitats across the Mediterranean represented an ecological opportunity for a salt tolerant genus like Limonium to experience a major episode of increased diversification rates. Conversely, the MSC apparently caused the extinction of other genera in the region, including the subtropical mangrove Avicennia (Thompson, 2005). Furthermore, the land corridors available across the Mediterranean at that time (see Mediterranean Paleogeographic Reconstruction During the MSC in Anzidei et al., 2014) might have provided venues for the colonization of North Africa from the Euro-Mediterranean (Figures 1, 2), or vice versa, as documented for Borago (Mansion et al., 2009). The refilling of the Mediterranean Sea at the end of the MSC may have promoted further vicariant speciation of Limonium in the Euro-Mediterranean and North Africa (Figure 1; see also Romeiras et al., 2016).
Diversification continued to increase exponentially also in the past 2-3 Myr, after the establishment of the Mediterranean climate and during the Pleistocene climatic oscillations (rate-through-time plots; Figure 4), corroborating results from other Mediterranean plant taxa (e.g., Valente et al., 2010;Fiz-Palacios and Valcárcel, 2011). Here, we explicitly tested the impact of paleo-temperature and past sea-level on diversification rates for the larger Mediterranean subclade. A pattern supported in all analyses was that extinction rates were higher when sea-level and temperature were higher (Supplementary Table S10). Additionally, our results supported heterogeneous, yet range-independent diversification of Limonium (Supplementary Table S13), demonstrating that occurrence in the Euro-Mediterranean region per se cannot explain the burst of diversification rates. Rather, we conclude that it was the combination of environmental changes experienced by Limonium in this region and intrinsic biotic traits (i.e., apomixis) that triggered FIGURE 6 | Evolution of habit in relation to island colonization for Pteroclados s.l. clade of Limonium from "Supermatrix-ITS-like" and "Supermatrix-cpDNAlike" MCC trees (model-averaged HiSSE analyses). Yellow at the tips and nodes denotes herbaceous habit and blue denotes woody habit. The pies at the nodes represent the relative probabilities of the inferred ancestral states. Net-diversification rates increase from black to red along the branches. Shift to woodiness followed island colonization but did not drive a significant increase of diversification rates. Photos of the herbaceous L. lobatum (top photo) and the woody-suffruticose L. arboreum (bottom photo) from Ares Jiménez.
the increase in diversification rates of the Mediterranean lineage, giving rise to its current high diversity of more than 400 species (see below).

Apomixis as a Biotic Driver of Diversification and the Interplay Between Biotic and Abiotic Factors in Bolstering Diversification Rates
In explaining the origins of the Mediterranean flora, questions have been raised about the relative role of geological and climatic history, and intrinsic biological processes such as polyploidization and hybridization (Thompson, 2005). In Limonium, as in other angiosperms (Hörandl and Hojsgaard, 2012), polyploidy and hybridization are usually associated with apomixis (Erben, 1978(Erben, , 1979, and the combined effects of these three biotic factors were proposed to have shaped the Mediterranean radiation species (e.g., Palacios et al., 2000;Lledó et al., 2005Lledó et al., , 2011. While several apomictic Limonium species are characterized as obligate apomicts (e.g., Palacios and González-Candelas, 1997;Arrigoni and Diana, 1999;Palacios et al., 1999;Khan et al., 2012;Róis et al., 2016), facultative apomixis (i.e., occasional sexual reproduction in apomicts) has also been documented for some species (e.g., D' Amato, 1949;Hjelmqvist and Grazi, 1964;Artelari, 1989;Artelari and Georgiou, 2002;Georgakopoulou et al., 2006). At the macroevolutionary level, our analyses show for the first time that apomixis is associated with an acceleration of diversification rates. Speciation rates are consistently higher in apomictic than sexual taxa, while extinction rates are either equal or higher in apomicts ( Figure 5; Table 3; Supplementary Table S11).
The higher extinction rate inferred for apomicts than sexuals is congruent with the labile nature of apomicts proposed by Darlington (1958) and Holsinger (2000), who regarded apomixis as a "blind alley of evolution" due to reduced recombination, hence lower genetic variation in apomicts, eventually driving them to extinction. This view was based on the assumption that apomixis is an irreversible, derived trait. However, our results suggest that transitions from apomixis to sexuality are very common in Limonium ( Figure 5) and corroborate recent studies that support apomixis-to-sex reversals enabled by the usually facultative nature of apomixis in angiosperms (Hörandl et al., 2007;Hörandl and Hojsgaard, 2012;Hand and Koltunow, 2014;Hojsgaard et al., 2014;Hojsgaard and Hörandl, 2015;Carman et al., 2019;Sharma and Bhat, 2020). While the result that transitions from apomixis to sexuality are more common than vice versa is strongly supported for both MCC trees in Limonium, it should be interpreted with caution because the limited phylogenetic resolution in the "Mediterranean lineage" could affect the inferred transition rates. The molecular genetic basis of apomixis is unknown in Limonium. However, in some plant genera, the molecular pathway of apomixis seems to be superimposed onto the pathway of sexual reproduction, rather than being completely independent, thus apomixis can revert to sexuality relatively easily (e.g., Hand and Koltunow, 2014 and references therein). This explanation for the molecular basis of apomixis is consistent with the occurrence of both apomictic and sexually developed seeds in facultative apomicts, suggesting that, if an ovule fails to initiate the apomictic pathway, sexual reproduction is activated, since the sexual pathway remains available. Importantly, in addition to enabling apomixis-to-sex reversals, the occurrence of "a little bit of sex" in apomicts can prevent the genomic decay caused by the absence of meiotic recombination, which can lead to extinction (Hörandl and Hojsgaard, 2012;Hojsgaard and Hörandl, 2015;Hodač et al., 2019).
Our results point toward a synergistic relationship between apomixis and sea-level fluctuations in driving the diversification of Mediterranean sea lavenders. In the Euro-Mediterranean region, the majority of Limonium species occur in coastal areas (which are vast, including many islands). Geo-climatic changes caused intense sea-level oscillations that had the highest frequency during the Plio-Pleistocene (Supplementary Figure S3), directly impacting coastal habitats. High sea levels triggered inundation and loss of available coastal habitats for Limonium, thus increasing extinction (see positive correlation of extinction rates with past sea levels for the "Mediterranean lineage" in Supplementary Table S10). Additionally, when sea levels are high, some populations may split and diversify allopatrically. Conversely, low sea levels create new areas that allow populations to expand, occasionally coming in contact and hybridizing. The newly created polyploid hybrids can establish new populations by single individuals through apomixis, and further diversify (see Negative Correlation of Speciation Rates With Past Sea Levels for the "Mediterranean lineage" in Supplementary Table S10). Apomixis provides both the required escape from sterility for hybrid polyploids and reproductive assurance in the absence of pollinators or mates in newly colonized habitats (Baker, 1955;Darlington, 1958). Thus, sea-level fluctuations cause repeated events of species-range contraction and fragmentation when sea level is high, and expansion and reconnection when sea level is low, promoting the origin of incipient species that can become established through apomixis. Long-term survival of apomictic species can be enhanced through occasional sex enabling the filtering of deleterious mutations via purging selection (Hojsgaard and Hörandl, 2015;Hodač et al., 2019).

Island Biogeography of Macaronesian Endemics and Insular Woodiness
Oceanic islands emerge from the sea empty of life. Their biodiversity is attributed to a combination of dispersals across the sea of species from neighboring areas and local diversification, producing high levels of endemicity (Cowie and Holland, 2006;Warren et al., 2015). Limonium endemics in Macaronesian archipelagos are of recent Plio-Pleistocene origin, except for L. dendroides, which diverged much earlier during the early Miocene (Supplementary Figures S4, S5). Limonium dendroides is an endangered Canarian endemic that is morphologically and taxonomically unique. It is characterized by at least two striking autapomorphies, i.e., arborescent habit and salt-glands in spikelet, and is the sole species of L. sect. Limoniodendron (Figure 1). The phylogenetic distinctiveness of L. dendroides, indicated by its isolated long branch sister to the large clade comprising all other species of subgenus Limonium, suggests that L. dendroides is either an old species or it is the only surviving species of a once larger clade that has undergone extensive extinction . Many plant taxa in Macaronesia comprise multiple endemic species that stemmed from a single long-distance dispersal followed by in situ diversification (e.g., Böhle et al., 1996;Francisco-Ortega et al., 1996;Barber et al., 2007;Kim et al., 2008;Caujapé-Castells, 2011). In Limonium, however, Macaronesian diversity was shaped by multiple LDD events (ca. 7), mostly from the neighboring Euro-Mediterranean and North African areas, followed by in situ speciation (ca. 22 events; BSM results). The Canarian and Cape Verde archipelagos were colonized repeatedly (at least four times and twice, respectively), as found in several studies of other Macaronesian taxa (Carine et al., 2004;Sanmartín et al., 2008 and references therein;Navarro-Pérez et al., 2015;Jaén-Molina et al., 2020). The availability of different ecological niches in oceanic islands may facilitate independent colonization by different species, possibly via reducing competition. For example, the two Limonium endemic lineages stemming from independent dispersals to Cape Verde occupy distinct habitats, namely wet mountainous cliffs and coasts (Romeiras et al., 2015). Thus, habitat heterogeneity of oceanic island systems favors colonization by a diversity of species pre-adapted to contrasting ecological niches.
While oceanic islands have been traditionally regarded as major sinks of biodiversity (Carlquist, 1974), their role as sources of biodiversity for neighboring continents or archipelagos has been rarely documented (e.g., Mort et al., 2002;Carine et al., 2004). Our results support four LDD events for Limonium from Macaronesia including at least one to North Africa, but also Trans-Atlantic LDD to the Americas (Figure 2; Supplementary Table S8). Fine-scale analysis of Jovibarba-Ctenostachys clade reveals a dispersal event from Macaronesia to Hispaniola and another to Morocco, and one or two interarchipelago dispersals between Canaries and Cape Verde (Figure 3). These results highlight the significance of oceanic archipelagos in Macaronesia as a source flora for LDD to archipelagos within and outside the region, and to neighboring continents (see also Gruenstaeudl et al., 2017;Jaén-Molina et al., 2020). Furthermore, the noteworthy long-distance dispersal from Cape Verde across the Atlantic Ocean to Hispaniola (Caribbean) was probably facilitated by both the light diaspores of these species and sea currents, such as the North Equatorial Current, as suggested for other angiosperms (Renner, 2004). Our results further suggest that different dispersal properties in Macaronesian lineages may explain the geographically broad range of the Jovibarba-Ctenostachys clade vs. the withinarchipelago range of the Canarian Nobiles clade. Indeed, reduced dispersal abilities have been documented for the latter clade (Jiménez et al., 2017), supporting the hypothesis of postcolonization loss of dispersal ability in colonists of oceanic islands followed by extensive in situ diversification (Carlquist, 1966;MacArthur and Wilson, 1967).
Woodiness on oceanic islands has been regarded as a key innovation linked to higher diversification rates when it is associated with disparity in growth forms (e.g., arborescent shrubs, subshrubs, trees, cushion forms, woody lianas, and giant rosette plants), because it enables the exploration of a broader niche space (Nürk et al., 2019). However, our results suggest that insular woodiness in the Canarian Nobiles clade of Limonium is either not linked to accelerated diversification (Figure 6; Supplementary Tables S14, S5) or a shift in rates for woody vs. herbaceous taxa is too moderate for the SSE methods to detect it, considering their limited power in analyses of small-size clades (Gamisch, 2016;Kodandaramaiah and Murali, 2018). A possible explanation for the result that the evolution of woodiness on the Canary Islands did not trigger a significant shift to much higher diversification rates may be that, in Canarian Limonium, woodiness is not associated with a diversity of woody growth forms (as found by Nürk et al., 2019), since all species in Nobiles clade have similar subshrub suffruticose habit. The lack of diversity in woody forms observed in species of the Nobiles clade may have limited their ability to radiate in a way that can be detected as a significant shift of diversification rates.

CONCLUSION
At the global scale, our study shows that the evolution of Limonium in space and time was shaped by major geologic and climatic processes that resulted in its cosmopolitan distribution and extensive diversification in the Mediterranean area. The increased diversification of Mediterranean sea lavenders was enabled by the ability of these taxa to reproduce asexually, via apomixis, which allowed them to survive and diversify during climatic oscillations. Our results show that the joint effect of biotic and abiotic factors is responsible for the current diversity of Limonium. At the regional scale, focusing on the insular endemics of Macaronesia, diversity on oceanic islands was shaped by multiple colonization events from neighboring continents and archipelagos, and in situ speciation. The Macaronesian islands have also served as sources for dispersals to other archipelagos (Caribbean) and continents (North Africa and America). In addition, woodiness in the Canarian sea lavenders (Nobiles clade) is a derived trait linked to insularity but not too much higher diversification rates. Our study highlights the importance of analyzing multiple abiotic and biotic factors, and their interactions, to achieve an in-depth understanding of evolution in hotspots of biodiversity.

DATA AVAILABILITY STATEMENT
The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found in the article/Supplementary Material.

AUTHOR CONTRIBUTIONS
EC and KK contributed to conception and design of the study. AJ, BW, KK, and MR provided material and data for analyses. KK performed all analyses with help from BW, MC, and ST and wrote the first draft of the manuscript. All authors contributed to manuscript revision, read, and approved the final version.

FUNDING
Financial support was given by the University of Zurich (Department of Systematic and Evolutionary Botany) and Georges-und-Antoine-Claraz-Schenkung.