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. 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).

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 well-supported 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 Muehlenbeckia-type 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).

We estimated ancestral ranges for Limonium on “Supermatrix-ITS-like” and “Supermatrix-cpDNA-like” MCC trees using the R package BioGeoBEARS 1.1.2 (Matzke, 2013; R Core Team, 2018). Considering current species distributions, patterns of endemism, and the floristic regions occupied by the species, we identify nine major biogeographic areas: South Africa, Euro-Mediterranean, North Africa, Irano-Turanian, Macaronesia, East Asia-Australia, Circumboreal, America, and Arabia-NE Africa. 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-within-area” 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, DIVA-like, 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.

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 birth-death). 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 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.

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 BiSSE-like 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 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 and 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.

Phylogenies inferred with ML and BI for each dataset recovered very similar topologies. Only a small number of conflicts were found between cpDNA and ITS trees, all located within Limonium (i.e., “rogue clades”: subclades within “Mediterranean lineage,” Circinaria clade and Jovibarba-Ctenostachys clade, and six “rogue taxa”; see Supplementary Figures S1, S2 and Malekmohammadi et al., 2017; Koutroumpa et al., 2018). The 50% majority-rule BI trees for all three Supermatrices showed better resolution (“Supermatrix”: 181/300, “Supermatrix-ITS-like”: 177/293, and “Supermatrix-cpDNA-like”: 176/300 nodes resolved, i.e., have posterior probabilities ≥0.5) than ITS and cpDNA trees (154/259 and 142/298 nodes resolved, respectively). Our results strongly support monophyly for Plumbaginaceae subfamilies (Plumbaginoideae and Limonioideae), tribes (Aegialitideae and Limonieae) and genera (except for Plumbago and Acantholimon; both non-monophyletic), and for major clades of Limonium corresponding to infrageneric subdivisions (except for section Schizhymenium). Hereafter, results for “Supermatrix-ITS-like” tree and “Supermatrix-cpDNA-like” tree will be reported first and second, respectively, separated by a vertical bar, unless otherwise specified.

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) within-area 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).

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 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). 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, range-dependent 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).

State-dependent diversification analyses in Pteroclados s.l. clade showed no effect of habit (herbaceous habit vs. woodiness) 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).

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; 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-ITS-like” 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., 2008, 2009; Salvo 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 the increase in diversification rates of the Mediterranean lineage, giving rise to its current high diversity of more than 400 species (see below).

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, 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., 2005, 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).

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 (Warren et al., 2018). 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 inter-archipelago 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 within-archipelago 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 post-colonization loss of dispersal ability in colonists of oceanic islands followed by extensive in situ diversification (Carlquist, 1966; MacArthur and Wilson, 1967).

Woodiness is very common in oceanic islands (Darwin, 1859; Wallace, 1878; Carlquist, 1965, 1974; Lens et al., 2013; Burns, 2019). Two competing hypotheses about the origin of insular woodiness have been proposed, which differ in whether woodiness evolved before or after island colonization. The “islands-as-museums” hypothesis views woodiness as a relictual trait already present in continental taxa that colonized islands and later went extinct from the continent (Cronk, 1992, 1997). Alternatively, woodiness has been interpreted as a derived trait that evolved in situ from herbaceous colonists subsequent to island colonization (e.g., Nürk et al., 2019). Our results in Limonium support the latter hypothesis by showing that woodiness is a derived trait in Pteroclados s.l. that emerged after colonization of the Canaries (woody Nobiles clade; Figure 6). The derived nature of woodiness is also supported by most studies on woody insular plants (e.g., Echium: Böhle et al., 1996; Lavatera: Fuertes-Aguilar et al., 2002; Sideritis: Barber et al., 2002; Sanctambrosia: Kool and Thulin, 2017), while fewer studies support its relictual nature (e.g., Tolpis: Moore et al., 2002; Descurainia: Goodson et al., 2006).

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.

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.