Phylogenetic Patterns of Extinction Risk in the Endemic Flora of a Mediterranean Hotspot as a Guiding Tool for Preemptive Conservation Actions

Climate change is a major driver of biodiversity decline with pervasive effects in biodiversity hotspots, where many endemic and threatened species thrive. However, the biological drivers of extinction susceptibility remain largely elusive, which hampers the implementation of effective conservation policies. Here, we advocate for the use of phylogenies as a complementary tool to inform policy makers. If we assume that the traits that determine extinction susceptibility are somewhat evolutionarily conserved, identifying the clades that accumulate a disproportionate amount of threatened species may help to mitigate potential increases in extinction risk among currently unthreatened species in these clades, even if the underlying biological drivers are unknown. We focused on the complete endemic angiosperm flora of a Mediterranean hotpot (Iberian Peninsula) to examine phylogenetic patterns in extinction risk expressed as IUCN categories (Least Concern “LC”, Near Threatened “NT”, Vulnerable “VU”, Endangered “EN” and Critically Endangered “CR”) using alpha and beta diversity metrics, comparative methods and a “hot node” approach. Phylogenetic diversity was significantly low for EN species and marginally significant for NT and CR, while LC and VU categories showed random pattern. Phylogenetic beta diversity (PBD) between IUCN categories was intermediate (0.40 – 0.61) and predominantly due to the “true” turnover component of PBD. Phylogenetic turnover was significantly low between NT – VU and VU – EN, suggesting that closely related species tend to show different threat status. In contrast, the comparisons involving the CR category sit toward the higher tail of the distribution, indicating a somewhat higher degree of clade specificity for CR species. In line with these patterns, phylogenetic signal in extinction risk was rather low (lambda = 0.23). Several of the “hot” clades that accumulated a significantly high number of species with the same threat status were specific to certain IUCN categories, yet few of them were observed across the categories. Most notably, the Caryophyllales stood out as the main threat-accumulating lineage, particularly within the Plumbaginaceae. All in all, our results indicate that few phylogenetic clades concentrate a great fraction of the extinction-risk gradient in the endemic flora of the western Mediterranean, and monitoring programs should pay particular attention to these extinction-prone lineages.

Climate change is a major driver of biodiversity decline with pervasive effects in biodiversity hotspots, where many endemic and threatened species thrive. However, the biological drivers of extinction susceptibility remain largely elusive, which hampers the implementation of effective conservation policies. Here, we advocate for the use of phylogenies as a complementary tool to inform policy makers. If we assume that the traits that determine extinction susceptibility are somewhat evolutionarily conserved, identifying the clades that accumulate a disproportionate amount of threatened species may help to mitigate potential increases in extinction risk among currently unthreatened species in these clades, even if the underlying biological drivers are unknown. We focused on the complete endemic angiosperm flora of a Mediterranean hotpot (Iberian Peninsula) to examine phylogenetic patterns in extinction risk expressed as IUCN categories (Least Concern "LC", Near Threatened "NT", Vulnerable "VU", Endangered "EN" and Critically Endangered "CR") using alpha and beta diversity metrics, comparative methods and a "hot node" approach. Phylogenetic diversity was significantly low for EN species and marginally significant for NT and CR, while LC and VU categories showed random pattern. Phylogenetic beta diversity (PBD) between IUCN categories was intermediate (0.40 -0.61) and predominantly due to the "true" turnover component of PBD. Phylogenetic turnover was significantly low between NT -VU and VU -EN, suggesting that closely related species tend to show different threat status. In contrast, the comparisons involving the CR category sit toward the higher tail of the distribution, indicating a somewhat higher degree of clade specificity for CR species. In line with these patterns, phylogenetic signal in extinction risk was rather low (lambda = 0.23). Several of the "hot" clades that accumulated a significantly high number of species with the same threat status were specific to certain IUCN categories, yet few of them were observed across the categories. Most notably, the Caryophyllales stood out as the main threat-accumulating lineage, particularly within the INTRODUCTION Climate change is a major driver of biodiversity decline (Sala et al., 2000) and particularly in biodiversity hotspots where many endemic and threatened species thrive (Bellard et al., 2014). With nearly 22500 plant species and 11700 of them as endemics, the Mediterranean Basin is among the most outstanding plant biodiversity hotspots of the world (Mittermeier et al., 2004), covering over 60% of the world Mediterranean climate extent (Klausmeyer and Shaw, 2009). This extraordinary concentration of plant diversity has urged conservation biologists to find solutions to minimize the projected loss of biodiversity in the region (Galli et al., 2012).
Understanding the association between species biological attributes and extinction susceptibility could help to make conservation planning more efficient (Purvis et al., 2000;Cardillo and Meijaard, 2012). However, and despite few traits have been identified as potential biological drivers of extinction risk in both plants (Lughadha et al., 2017;Mankga and Yessoufou, 2017) and animals (Cardillo et al., 2008), there is a noticeable lack of consistency among studies (Davies, 2019), and the relationship between the vast majority of traits analyzed and extinction risk remains largely elusive (Chichorro et al., 2019). Critically, this lack of empirical consistency hampers the identification of the species that might be most in need of conservation action.
Many traits that are linked to plant extinction risk such as phenology and biotic pollination (Lughadha et al., 2017) are evolutionarily conserved Cirtwill et al., 2020), suggesting that climate-driven extinctions may be nonrandom with respect to phylogeny. For example, Willis et al. (2008) found that flowering-time responses to climate change are shared among closely related species and strongly correlated with species abundances, which explained phylogenetically clumped patterns of species loss. On the other hand, the association between extinction risk and species phenotypes appears to vary geographically (Sodhi et al., 2008;Godefroid et al., 2014), suggesting that the clades that are threatened by climatic stressors may be region-specific if we assume the biological drivers of extinction susceptibility are evolutionarily conserved. For example, if pseudanthia (a strongly conserved floral trait) were associated with high extinction risk in humid regions, Asteraceae representatives could be threatened in the tropics but not in drylands. Therefore, identifying the clades that accumulate a disproportionate amount of threatened species in specific regions may help to guide preemptive conservation actions to mitigate potential increases in extinction risk among currently unthreatened (or near threatened) species in these clades (Yessoufou et al., 2012), even if the underlying biological drivers are unknown. Moreover, considering extinction within a phylogenetic framework allows quantifying its impacts on the tree-of-life as the loss of evolutionary history (Davies and Yessoufou, 2013;Faith, 2018), an increasingly appreciated perspective within conservation goals (Laity et al., 2015;Gumbs et al., 2018;Carta et al., 2019).
Here, we focused on the complete endemic angiosperm flora of the Iberian Peninsula in the western Mediterranean to explore phylogenetic patterns of extinction risk using IUCN assessments and time-calibrated molecular phylogenies. The Iberian Peninsula is one of the two major centers of narrow endemics within the Mediterranean Basin biodiversity hotspot (Medail and Quezel, 1997), and thus represents an ideal system to evaluate the usefulness of phylogenetic tools for the conservation of endemic plants in the Mediterranean. Specifically, we used the phylogenetic diversity metric (Faith, 1992) to assess whether extinction risk is phylogenetic clumped in this flora (phylogenetic diversity of IUCN categories), and further quantified the degree of clade-specificity in threat categories using the "true" turnover component of phylogenetic beta diversity (Leprieur et al., 2012). Complementarily, we used a "hot node" approach (e.g., Saslis-Lagoudakis et al., 2011;Molina-Venegas et al., 2020b) to elucidate the identity of the clades with a significant overabundance of species in each category. Finally, we evaluated the degree of phylogenetic signal in extinction risk using comparative methods.

Endemic Flora of the Iberian Peninsula
We compiled all angiosperm species and subspecies included in the published volumes of Flora iberica (Castroviejo, 1986(Castroviejo, -2019 specifically registered as endemics for said work, which covers continental Portugal, continental Spain (including the Balearic Islands, which are an extension of the Baetic range in southern Spain; see Figure 1A) and Andorra (a small country in the eastern Pyrenees). To make up for the absence of volumes dedicated to the Poaceae family (still unpublished), we took advantage of the comprehensive monography by Romero Zarco (2015) and different generic treatments published in recent years (e.g., Devesa et al., 2013;Ortega-Oliviencia and Devesa, 2018). As the Flora iberica project has been under development for forty years, several taxa have been described after their first volumes were published, many of them being endemics to the territory. Thus, an exhaustive search for these taxa (8.4% of the species in the dataset) was undertaken using the option "new taxa" on the Flora iberica website (floraiberica.es/miscelania/nuevos_taxones.php), revising published generic monographs and descriptions of new species and subspecies during these years (e.g., Salvà, 2006;Zozomová-Lihová et al., 2014;Marques et al., 2017).
The Pyrenees, a mountain range between Spain and France, comprises the northern limit of the Iberian Peninsula ( Figure 1A). Therefore, species whose distribution is restricted to the French part of the Pyrenees (a region representing 1.6% of the area of the Iberian Peninsula) can be considered as endemics of the Iberian Peninsula. Such species are very low in number (3-6 species; see Villar and García, 1989), and they were not included in our dataset because we lack detailed floristic information of this region as to rigorously confirm their Iberianendemic status. The resulting list of Iberian endemics consisted of 2,138 angiosperm species and subspecies. Endemic gymnosperms and ferns were not considered in the analyses because (1) they are comparatively rare in the study area (only 9 species), and (2) the much older most recent common ancestor (MRCA) of vascular plants compared to that of the angiosperms can obscure phylogenetic patterns among the latter (Letcher, 2010;Qian et al., 2017).

IUCN Assessments
The allocation of risk categories according to the IUCN scheme had recently been addressed for all the Spanish flora in a previous work (Muñoz-Rodríguez et al., 2016). These authors combined the last published Spanish Red List (Moreno-Saiz, 2008) with additional information appeared in regional Red Books, and the threat status for the species that were not included in these sources was assigned based on their area of occupation, number of localities, and future distribution trends using environmental niche models (IUCN Standards Petitions Committee, 2019). Some updates on the threat status of the species have been carried out after the publication of the new volume of the Spanish Red Book (Moreno-Saiz et al., 2019). For Iberian endemics whose distribution is restricted to Portugal, the categories adjudged on the Lista Vermelha project website were used (LVFVPC, 2020). Finally, several new taxa described in the last years included a proposed IUCN category in their descriptions, which were used in our analysis (see Supplementary Table 1 for a list of Iberian endemics with IUCN categories).

Phylogenetic Information
We obtained a time-calibrated phylogeny of DNA for the endemic angiosperm flora of the Iberian Peninsula using the R package V.PhyloMaker (Jin and Qian, 2019). This software uses the largest time-calibrated mega-phylogeny of vascular plants available to date (GBOTB.extended; see Jin and Qian, 2019) to generate a subtree from a given species list (standardized to The Plant List, 2013 v 1.1 nomenclatural and spelling criteria) following a three-steps procedure. First, V.PhyloMaker finds for each genus in the list the MRCA of all the tips in the largest monophyletic cluster within GBOTB.extended, and defines it as the crown node of the genus. Second, the species in the list that are missing in GBOTB.extended are bound to the tree using a sticking scenario that is specified by the user. Here, we used scenario 2, which binds the new tips to a randomly selected node at and below the corresponding genus crown node (or family node if the genus is missing in GBOTB.extended). Finally, the resultant tree is pruned to include only the species in the input list. The GBOTB.extended phylogeny does not include infraspecific taxa, which are particularly abundant in biodiversity hotspots. Thus, phylogenetic relationships among accepted subspecies (after The Plant List, 2013 v. 1.1 nomenclatural criteria) were randomly resolved in the corresponding terminal tips (i.e., species) using the R package phytools (Revell, 2012). No infra-subspecific taxa were considered in the study. In order to account for phylogenetic uncertainty (i.e., random binding of species and subspecies), we repeated this procedure iteratively to obtain 1000 different phylogenetic hypotheses, and all subsequent analyses were replicated and results averaged over the 1000 trees (Rangel et al., 2015).

Phylogenetic Patterns of Extinction Risk
First, we explored whether extinction risk is phylogenetic clumped in the endemic angiosperm flora of the Iberian Peninsula. To do so, we grouped species in the following categories: Least Concern (LC, 885 species), Near Threatened (NT, 297 species), Vulnerable (VU, 357 species), Endangered (EN, 180 species), Critically Endangered (CR, 132 species), Data Deficient (DD, 86 species) and Not Evaluated (NE,196 species). Species categorized as Extinct in the Wild (EW, 2 species) or Extinct (EX, 3 species) were not considered in the analyses. For each category, we quantified whether the phylogenetic diversity (PD sensu Faith, 1992) encapsulated by their constituent species was lower (phylogenetic clustering), random or higher (phylogenetic overdispersion) than expected using a null-model approach. The observed PD values were compared to a null distribution of values generated by shuffling taxa labels in the phylogeny 999 times using the function ses.pd implemented in Picante R package v.1.8 (Kembel et al., 2010) with default settings. For a nominal alpha of 5%, ses.pd values (Z-scores) below and above −1.96 and +1.96 can be considered as significant phylogenetic clustering and overdispersion, respectively (Kembel, 2009). This standardization corrects for differences in species richness between the categories. It is important to note that DD and NE categories do not represent actual threat categories but incomplete sampling or taxonomic knowledge. Thus, in order to assess the impact of these categories in the computation of Z-scores, we conducted the analyses using two different species pools. Pool 1 included the species categorized as LC, NT, VU, EN, and CR (n = 1851) and Pool 2 included also DD and NE species (n = 2133). We also explored the phylogenetic structure of threatened species (i.e., CR + EN + VU as a unique category) following the procedure described above.
Complementarily to phylogenetic diversity analyses, we used a "hot node" approach (e.g., Saslis-Lagoudakis et al., 2011;Molina-Venegas et al., 2020b) to elucidate the identity of the clades with a significant overabundance of species in each category. For each IUCN category i, the number of species assigned to i that descended from each node j of the phylogeny k was counted. Then, this number was compared to a null distribution of values generated by shuffling trait values (i.e., 1 if the species is included in category i, 0 otherwise) across the tips of phylogeny k 999 times to compute a Z-score for each node j of k. Nodes were considered as "hot" if their associated Z-scores were higher than +1.96 (5% nominal alpha). We only evaluated those clades that included 10 species or more, since previous studies have documented unacceptable rates of statistical errors for smaller lineages (Parra et al., 2010). Note that the identity of the hot nodes may change across phylogenies due to (1) different randomization schemes for the missing species (specially below genus-level, where most missing species were inserted) and (2) Z-scores that sit very close to the significance threshold. To account for such uncertainty, we only considered as "hot" those nodes that were observed and showed statistical significance in at least 50% of the trees analyzed.
We quantified the degree of clade-specificity in threat categories (i.e., whether phylogenetic clades include only species with the same threat status) using the PhyloSor similarity index (Bryant et al., 2008) as a metric of phylogenetic beta diversity (PBD). The PhyloSor index computes the fraction of evolutionary units (typically branch-length) that is shared between two samples, and it ranges between 0 (no branch-length is shared) and 1 (all branch-length is shared). Thus, PBD can be defined as 1 -PhyloSor index. PBD can be decomposed into two additive components, namely "true" turnover (pβsim) and nestedness, which represent different aspects of phylogenetic beta diversity (Leprieur et al., 2012). While the nestedness component represents the fraction of PBD that is simply due to differences in PD, the turnover component implies the replacement of an exact amount of branch-length, the branch-length that is replaced being exclusive to each sample. The observed pβsim values between IUCN categories were standardized following the same procedure described earlier for PD. Lower and higher than expected pβsim values indicate that phylogenetic "true" turnover (hereafter "turnover") between the samples (here IUCN categories) tends to occur toward the tips and the root of the phylogeny, respectively (Molina-Venegas et al., 2015, 2020a. Thus, significantly low turnover between IUCN categories indicates that closely related species tend to show either conservation status (i.e., low clade-specificity in threat categories), while higher than expected turnover indicates that close relatives tend to show the same status (i.e., high clade-specificity).
In addition, we assessed the degree of phylogenetic signal in extinction risk using models of trait evolution. To do so, we assumed that transitions between threat categories (codified as LC = 1, NT = 2, VU = 3, EN = 4, and CR = 5) occur in a stepwise fashion (e.g., from LC to NT, CR to EN, and VU to EN, etc.) without skipping intermediate steps (meristic transitions between neighboring states), and tested three different models: (1) Brownian motion, (2) Pagel's lambda, and (3) white noise (evolution independent of phylogeny). Pagel's lambda equates Brownian motion when the estimated parameter lambda is equal to 1 (strong phylogenetic signal), and white noise if lambda = 0 (lack of phylogenetic signal). Model accuracy was estimated using likelihood ratio tests. Although more complex models are possible (e.g., delta time-dependent model, early burst), we are not interested in the specific mode of evolution of the IUCN categories but in providing empirical evidence for phylogenetic signal in extinction risk (i.e., whether closely related species are similarly threatened), and thus the evolutionary models used here suffice to fulfill the objectives of the study. These comparative analyses were conducted using only the species in Pool 1 (LC, NT, VU, EN, and CR) because meristic transitions between DD or NE and any other category would be meaningless.
In order to examine potential variation of the patterns across phylogenetic scales (see Graham et al., 2018), we also conducted the analyses within the eudicot (1887 species) and monocot (244 species) clades, respectively. All the analyses were conducted in R (R Core Team, 2020) using the packages picante (Kembel et al., 2010), geiger (Pennell et al., 2014), betapart (Baselga et al., 2018) and the R code provided in Molina-Venegas et al. (2020b) for the "hot nodes" analysis.

Phylogenetic Diversity of IUCN Categories
Phylogenetic diversity (PD) of LC and VU categories (i.e., averaged Z-score) did not depart from random expectation, while it was significantly low (i.e., phylogenetic clustering) for EN ( Figure 1B). The averaged Z-score for NT and CR sit very close to the threshold of significance toward the lower tail of the distribution (i.e., marginally significant clustering, Figure 1B). DD and NE categories showed a strong phylogenetic clustering, yet the averaged Z-scores for the actual threat categories (i.e., LC, NT, VU, EN, and CR) were only very slightly affected when DD and NE species were included in the analysis, the general pattern remaining virtually unchanged (Supplementary Figure 1). Phylogenetic diversity patterns within the eudicots were very similar, the only difference being that NT and CR became even more clumped and EN slightly less (Supplementary Figure 2A,B). In contrast, FIGURE 2 | Phylogenetic beta diversity analyses (PBD) among IUCN categories (LC, NT, VU, EN, and CR) for the endemic angiosperm flora of the Iberian Peninsula.
The barplot (A) shows the fraction of PBD due to phylogenetic "true" turnover (in red) and nestedness (in blue) for each pairwise comparison. The box-and-whisker plot (B) shows standardized values (Z-scores) for the turnover component. The horizontal black lines within the boxes represent the mean value of the distribution of Z-scores in each pairwise comparison, and the horizontal gray dashed lines are visual references at y = ± 1.96 (significance thresholds for a 5% nominal alpha).
Frontiers in Ecology and Evolution | www.frontiersin.org phylogenetic diversity was randomly distributed within the monocot clade with the only exception of the NE category, which remained clustered (Supplementary Figure 2C,D). The subset of threatened species (i.e., CR + EN + VU) was phylogenetically more clumped than unthreatened species regardless of the phylogenetic scale, although differences were not significant (Supplementary Figure 3). Phylogenetic uncertainty had a low impact in the analyses, as the distribution of Z-scores showed strong kurtosis and very narrow 95% confidence intervals (Supplementary Figure 4).

Phylogenetic Beta Diversity Between IUCN Categories
Phylogenetic beta diversity (PBD, 1 -Phylosor index) was overall intermediate (ranged between 0.40 -0.61) and predominantly due to the turnover component except in LC -EN and LC -CR comparisons, where PBD was largely due to differences in PD between the categories (i.e., nestedness) (Figure 2A). Phylogenetic turnover was significantly low between NT -VU and VU -EN (Figure 2B), while the averaged Z-score of the comparisons involving the CR category sit toward the higher tail of the distribution. Specifically, phylogenetic turnover was significantly high and marginally significant between CR -NT and CR -VU, respectively. This general pattern of significance persisted when DD and NE species were included in the analyses (Supplementary Figure 5). PBD patterns were very similar for the eudicots, although the CR -NT and CR -VU comparisons shifted away from the significance threshold (Supplementary Figure 6). The only significant turnover within the monocot clade occurred between NT -EN (lower than expected, Supplementary Figure 7). The general pattern of significance persisted for both the eudicots and the monocots when DD and NE species were included in the analyses (Supplementary Figures 8,9), and phylogenetic uncertainty had also a low impact in the PBD analyses (Supplementary Figures 10-12).

Phylogenetic Hot Nodes
While several of the hot nodes detected in the analysis were specific to certain IUCN categories, few of them were observed across the categories (full details are shown in Table 1). Most notably, the Caryophyllales stood out as the main threat-accumulating clade for CR, EN, and VU categories, particularly within the Plumbaginaceae (Figure 3). The Alchemilla + Potentilla (Rosaceae) clade was also detected as a shared hot node between EN and NT categories (although no Potentilla species were observed in the former category), and the Primulaceae hot node was observed in both VU and NT categories. The eudicots and eudicots excluding Ranunculales clades were identified as deep hot nodes in the NT and EN categories, respectively. LC species were significantly numerous within few deep clades such as monocots, asterids and Saxifragales (Figure 4). No new hot nodes emerged when CR, EN, and VU species were analyzed as a unique thread category (Supplementary Figure 13). We found a strong overdominance of Rosaceae, monocots, Malpighiales, Asterales, and Orobanchaceae representatives in the DD and NE categories (Supplementary Figure 14).

Phylogenetic Signal in Extinction Risk
We found no higher support for the Brownian meristic model of evolution than the non-phylogenetic one (p-value of the likelihood ratio test = 1 in all cases). However, the Pagel's lambda model fit well to our data compared to both the non-phylogenetic (p-value < 0.001 for 99.9% of the trees) and Brownian (pvalue < 0.001 in all cases) models, although the signal was rather weak (Pagel's lambda = 0.23 ± 0.001; mean lambda and 95% confidence interval). Results were very similar for the eudicot clade (lambda = 0.26 ± 0.001), however neither the Brownian, As a general rule, only the most inclusive hot clades in each category are shown. Hot clades that were nested within the Caryophyllales (ID 1), Fabales + Rosales + Fagales (ID 7) and Lamiales (ID 12) hot clades are also shown (highlighted with dots). Most inclusive hot clades that did not contribute with any species to the total number of species at risk observed within their respective child nodes are not shown. For example, the Euphorbia + Linum clade was detected as a most inclusive hot node in the CR category, but no Linum species are categorized as CR. Thus, the Euphorbia hot node (ID 5) is shown instead. Note that the eudicots (ID E1) and the eudicots excluding Ranunculales (ID E2) clades were hot clades for NT and EN, respectively, and all hot nodes detected in the analysis are nested within these clades. Clade ID links table entries to phylogenic nodes in Figure 3. NPAAA: non-protein amino acid-accumulating clade.
FIGURE 3 | One randomly selected tree topology for the endemic angiosperm flora of the Iberian Peninsula that included all the hot clades detected in the analysis (CR, EN, VU, and NT categories). Hot nodes are marked with pie charts, representing level of support for the nodes, this is, the number of times they were observed with statistical significance (5% nominal alpha) across all the trees analyzed. The circle symbols on the phylogenetic tips represent the species of each category. The numbers correspond to entries in Table 1, where the identity of the nodes is shown.
not the Pagel's lambda models were supported for the monocots (lambda = 0 in all cases).

DISCUSSION
Minimizing expected plant diversity loss in biodiversity hotspots is a major challenge for conservation biologists, and phylogenetic tools may provide a useful handle to assess current and future threats in these endemic species-rich regions (Purvis, 2008).
Here, we show that extinction risk in the endemic angiosperm flora of the Iberian Peninsula is largely clumped in the phylogeny (Figure 1), and that specificity in the relationship between IUCN categories and phylogenetic clades is overall low (Figure 2). This lack of association between clades and threat categories is consistent with the low phylogenetic signal detected in our comparative analysis, suggesting that closely related species tend to differ in their threat status (although CR species showed a somewhat higher degree of clade specificity, Figure 2). All in all, our results indicate that few phylogenetic clades concentrate FIGURE 4 | One randomly selected tree topology for the endemic angiosperm flora of the Iberian Peninsula that included all the hot clades detected in the analysis (LC category). Hot nodes are marked with pie charts, representing level of support for the nodes, this is, the number of times they were observed with statistical significance (5% nominal alpha) across all the trees analyzed. The circle symbols on the phylogenetic tips represent LC species. 1: monocots, 2: Saxifragales, 3: Genisteae (Fabaceae), 4: asterids.
a great fraction of the extinction-risk gradient in the endemic flora of the western Mediterranean, and monitoring programs should pay particular attention to these extinction-prone lineages ( Table 1 and Figure 3). The Caryophyllales clade, and particularly the subclade Plumbaginaceae, stood out as the main threat-accumulating lineage for CR, EN, and VU categories. This pattern is supported by previous non-phylogenetic studies that showed a significant association between several Caryophyllales representatives such as Limonium (Plumbaginaceae), Armeria (Plumbaginaceae), Silene (Caryophyllaceae), or Petrocoptis (Caryophyllaceae) and high extinction risk (Buira et al., 2020). Most species in these genera are associated to coastal, rocky habitats, salt marshes and edaphic islands that are either naturally or anthropogenically fragmented, which may have largely contributed to their current threat status (Lienert, 2004). The genus Limonium (Plumbaginaceae) is well-known for its ability to produce apomictic off-spring (Róis et al., 2016), which may represent a particular adaptation to maintain small populations in fragmented landscapes. Many Alchemilla representatives (Rosaceae) show a similar apomictic reproduction mode and low-sized populations in restricted rocky habitats, and the clade was also identified as an endangered group (Table 1 and Figure 3). Delimiting taxonomic species within these speciesrich genera is often difficult due to the overwhelming number of locally restricted microspecies they contain, which has resulted in many taxa being left out of current species-oriented conservation legislation (Ennos et al., 2005;Domínguez Lozano et al., 2007). The phylogenetic perspective presented here may help to bridge this gap between conservation law and threatened biodiversity, as many taxa that are taxonomically difficult to distinguish would be preemptively protected under the phylogenetic "umbrella".
While few clades concentrated a great fraction of the extinction-risk gradient in the study area (Table 1 and Figure 3), others stood out as accumulators of species categorized as "least concern" instead, particularly the monocots and the asterids (Figure 4). Both clades include species-rich genera (e.g., Teucrium, Thymus, Linaria, Centaurea, Narcissus, and Festuca) that unlike the ones discussed earlier (e.g., Limonium, Armeria, and Potentilla) they do not appear to be particularly threatened (with the exception of Taraxacum, Table 1). A considerable number of these unthreatened endemic species occur in mountainous habitats of the study area (Buira et al., 2020), which might be acting as refugia against extinction. However, increasing global warming cast reasonable doubts on the persistence of these putative refugia (Ohlemüller et al., 2008), which makes long-term monitoring programs in mountainous areas such as the GLORIA, 2020 project 1 an indispensable reference for policy makers.
Although threatened species were overall more clumped in the phylogeny than unthreatened ones (Figure 1), we found that IUCN categories showed contrasting phylogenetic diversity values. As such, EN and CR species were phylogenetically clumped, but VU species were not. This suggests that the trend toward phylogenetic clustering in extinction risk reported here (Figure 1 and Supplementary Figure 3) is driven by the most threatened species. In a previous study conducted in the Eastern Arc Mountains biodiversity hotspot in Africa, Yessoufou et al. (2012) found just the opposite pattern, this is, phylogenetic clustering in VU and random pattern in EN and CR categories. These authors conducted their analyses on the native flora of Tanzania regardless of the endemic status of the species, while we focused on the endemic flora of the Iberian Peninsula, which is largely concentrated in species-rich lineages that have experienced recent evolutionary diversification during the late Pliocene and Quaternary (Rundel et al., 2016). Some authors have reported a tight relationship between young and fast-evolving plant lineages (neo-endemics) and high extinction risk (Davies et al., 2011;Tanentzap et al., 2020), which may explain the discrepancies between our findings and the study by Yessoufou et al. (2012). Further, our results support the idea that evolutionary history itself might be an important predictor of extinction risk, to the extent that the latter could be simply determined by differential diversification rates between lineages (Davies et al., 2011;Tanentzap et al., 2020). Nonetheless, such hypothesis is yet to be confirmed for the endemic flora of the western Mediterranean.
It is remarkable that although the gaps of information in the conservation status of Iberian endemics did not affect the phylogenetic patterns reported here, species categorized as "Data Deficient" (DD) and "Not Evaluated" (NE) were strongly clumped in the phylogeny (Supplementary Figure 1). The significant relatedness among DD species may respond to the taxonomic challenge posed by the clades where these species are largely clumped, namely Rosaceae (mostly Alchemilla and Rubus) and Taraxacum (Supplementary Figure 14). As such, these species are often referred to as aggregates in the botanical literature (i.e., group of similar species) due to extreme lack 1 www.gloria.ac.at of consistency in diagnostic traits (Graham and Woodhead, 2009) or simply because they have remained understudied. On the other hand, many Poaceae and Asteraceae species that are endemic to the study area have been described very recently without being assigned a threat status, which may explain the overwhelming accumulation of NE species in the monocot and Asterales clades, respectively. Finally, lack of taxonomic consensus between specialists may have also contributed to the accumulation of NE species in certain groups that are currently evolving, such as Narcissus (monocots). Therefore, future conservation efforts should be directed to fill in these gaps of data if we are to take well-informed conservation actions to preserve the endemic flora of this emblematic hotspot in the western Mediterranean.

DATA AVAILABILITY STATEMENT
The list of Iberian endemics and their IUCN status is provided as part of the Supplementary Material for this article.

AUTHOR CONTRIBUTIONS
RM-V and JCM-S conceived the ideas. RM-V conducted the analyses and led the writing. JCM-S and IR-G compiled the dataset and contributed to the writing. IR-G contributed to Figure 1 (map). All authors approved the submitted version.

FUNDING
RM-V was supported by the TALENTO program (2018-T2/AMB-10332), and IR-G by the Garantía Juvenil program (PEJ-2018-AI/AMB-9865), both from the Regional Government of the Community of Madrid, Spain. This paper contributes to the project REMEDINAL TE (P2018/EMT-4338) from the Regional Government of the Community of Madrid.