We compiled all angiosperm species and subspecies included in the published volumes of Flora iberica (Castroviejo, 1986–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 Iberian-endemic 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).

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

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

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

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

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 (p-value < 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, not the Pagel’s lambda models were supported for the monocots (lambda = 0 in all cases).