Assessing the vulnerability of habitats through plant rarity patterns in the Pyrenean range

Biodiversity is being severely affected by global change worldwide, but consequences differ across individual species and environments. Since rarity is a condition associated with vulnerability and extinction risk, it is a concerning factor itself that can also amplify the detrimental impact of global change drivers. Here, we assess differences in vulnerability derived from species rarity across environments in a South European mountain range, the Pyrenees, a complex landscape sheltering more than 3400 vascular plants. We first analyzed patterns of alpha taxonomic (richness) and phylogenetic diversity (standardized phylogenetic diversity, sesPD) from more than 18,000 plant surveys across 14 different habitats to find out which ones shelter the highest diversity. Then, we classified plant species according to four criteria of rarity (Pyrenean endemic, regional geographic range, ecological specialization and local abundance) and estimated their frequency among habitats and their contribution to PD. Generalized linear models were used to compare richness and sesPD among habitats, and the relationship between each diversity metric and rare species proportion was assessed. Grasslands harbor the greatest number of species but are poor in sesPD, while inland surface waters contain fewer species but are the most phylogenetically diverse. Deciduous, evergreen and Mediterranean forests show both high species richness and sesPD. Rare species account for half of the regional species pool, and are more abundant in habitats with fewer species and higher sesPD. Diversity metrics showed opposite tendencies between them and with rare species proportion, emphasizing that they are not interchangeable. Our results highlight the vulnerability of rocky and water‐related habitats due to the rarity of the species they shelter, but only the latter match other assessments of habitat vulnerability to global change drivers. The analysis of rarity patterns can guide conservation efforts by identifying priority targets when information about direct threats to habitats is scarce or incomplete.


| INTRODUCTION
Mountains around the globe are considered major plant diversity hotspots due to the large amount of habitats and species they harbor, along with a high level of endemism (e.g., Mittermeier et al., 2004). At the same time, these regions, and particularly those in or around the Mediterranean Basin, are among the most threatened areas by global change (Hock, Rasul, Adler, C aceres, & Gruber, 2019). Mountain biodiversity is widely acknowledged as highly vulnerable due to two main factors: temperatures increasing faster than the global average (Bravo, Araújo, Lasanta, & Moreno, 2008;Cramer et al., 2018;Pachauri et al., 2014), and the dramatic habitat transformation caused by land use changes during the last century (García-Vega & Newbold, 2020;Mottet, Ladet, Coqué, & Gibon, 2006;Newbold et al., 2015).
The vulnerability of biodiversity within a region, however, is not defined solely by its exposure to external hazards like global change drivers. Some features relative to species distribution, specialization, or population abundance may result in rarity, prompting an intrinsic sensitivity to factors such as demographic stochasticity (Davies, Margules, & Lawrence, 2004;Matthies, Bräuer, Maibom, & Tscharntke, 2004), or deterministic perturbations such as slow habitat succession (Weißhuhn, Müller, & Wiggering, 2018). Rarity may arise from a variety of factors, which are biological, ecological and historical in nature (Stebbins, 1980), resulting in each type of rarity having different risks associated with them. For example, species with limited geographic range have been linked to higher extinction risk, likely due to a reduction in the buffering effect of range size against abiotic and biotic stressors (Harnik, Simpson, & Payne, 2012;Payne & Finnegan, 2007). Habitat specialists are constrained by certain environmental conditions and in turn are more vulnerable to habitat fragmentation and loss (Fontaine et al., 2007). Finally, species with smaller populations are more vulnerable to stochastic (Matthies et al., 2004) and denso-dependent phenomena like the Allee effect (Kuparinen, Keith, & Hutchings, 2014). Since the type of rarity affects species vulnerability to particular factors or processes, a comprehensive analysis of the frequency of different kinds of rarity within a geographic region can give us a broader knowledge about the possible effects of global change drivers and other disturbances on biodiversity patterns of such areas (Caro & Girling, 2010).
The loss of rare species might, in turn, have strong implications for ecosystems beyond a reduction in taxonomic diversity. On the one hand, rare species can be considered important assets for ecosystem functioning as they sometimes perform unique or novel functions despite their relative low abundances (Jain et al., 2014;Lyons, Brigham, Traut, & Schwartz, 2005;Mouillot et al., 2013). On the other hand, biodiversity goes beyond the number and abundance of species present in a community and thus the consequences of their loss also extend to the evolutionary and functional aspects of diversity (Craven et al., 2018;Naeem et al., 2016). In consequence, exploring and comparing the frequency of various types of rare species across different communities, habitats or ecosystems, along with their contribution to taxonomic and phylogenetic diversity (PD), may help us assess not only their vulnerability to intrinsic factors, but also the consequences of their loss for those habitats.
The aim of this study is to explore how vulnerability of biodiversity varies across different habitats based on the frequency and type of rare species they shelter throughout a large and environmentally heterogeneous mountain area in Southern Europe: The Pyrenees. We used the largest available dataset of plant communities in the region (~18,000) and classified each plant species according to four important features associated with rarity based on the proposal of Rabinowitz (1981): amplitude of their geographic distribution, regional abundance, ecological specificity, and local abundance. In particular, we strive to answer the following questions: (1) which are the kinds of rarity more prevalent in the Pyrenean range, and how is the frequency of rare plants distributed among habitats? We expect certain kinds of rare species, for instance habitat specialists, to accumulate in locations with highly distinct environmental conditions such as rocky habitats and alpine environments or wetlands (Boulangeat, Lavergne, Van Es, Garraud, & Thuiller, 2012;Enquist et al., 2019;Harrison & Noss, 2017). (2) What is the relationship between the frequency of rare plants, and taxonomic and phylogenetic diversity across habitats? We expect the proportion of rare species to increase in species-rich communities, as described by Heegaard, Gjerde, and Saetersdal (2013), and in turn, with sesPD. (3) Do rare species contribute more to the phylogenetic diversity than other species, and does this contribution differ across habitats and rarity types? Here, we hypothesize that the contribution of rare species to PD will be higher than their less rare counterparts as observed by Mi et al. (2012) across a global network of forest communities.

| Study area
The Pyrenean range is located in the southwest of Europe, and separates the Iberian Peninsula from the rest of the continent (Figure 1). Due to its location in the transition between the Alpine, Mediterranean, and Atlantic biogeographic regions, together with a broad altitudinal gradient (highest summit at 3404 m a.s.l.), the Pyrenees show a wide variety of climatic conditions (Ninot, Carrillo, & Ferré, 2017). This, combined with intensive land use (García-Ruiz et al., 1996), and geological heterogeneity (García-Ruiz et al., 2015), results in a wide diversity of natural and semi-natural habitats. The Pyrenean flora is composed of more than 3600 native vascular plants with a broad range of environmental requirements, from Mediterranean species adapted to dry, hot conditions; to boreoalpine species inhabiting the cold environments of summits (G omez, García, Font Castell, & Aizpuru Oiarbide, 2017).

| Plant data collection
We downloaded and validated 18,608 plant inventories or relevés from SIVIM, a database of Ibero-macaronesian vegetation (Font, Balcells, & Mora, 2017). They were located between 400 and 3300 m a.s.l. and included over 400,000 records of a total of 2550 species, which represent around 74.2% of the native mountain species in the Pyrenean flora (those occurring over 400 m a.s.l.; G omez, . Each inventory was categorized according to 14 European EUNIS habitats, a hierarchical classification system for the terrestrial and marine habitats of Europe and its surrounding waters (Moss, 2008). We used different EUNIS levels in order to better accommodate habitat F I G U R E 1 Map of the study region. The grid marks the UTM cells the limits established by the Working Community of the Pyrenees (CTP). Numbers within cells indicate the number of plant inventories located in that cell frequency in the Pyrenees and obtain a more balanced sample of communities. Habitats with very high anthropic influence such as ruderal communities and irrigated meadows were excluded, and only inventories with more than five species were used. The final list of habitats was: surface waters (C), mires and fens (D), dry grasslands (E1), seasonally wet and wet grasslands (E3), alpine and subalpine grasslands (E4), woodland fringes and clearings and tall forb stands (E5), arctic, alpine and subalpine scrub (F2), temperate and Mediterranean-montane scrub and heathlands (F3-F4), garrigues (F6), broadleaved deciduous woodland (G1), broadleaved evergreen woodland (G2), coniferous woodland (G3), screes (H2), and inland cliffs, rock pavements and outcrops (H3). F3 and F4 habitats were joined together due to their low abundance in our dataset as well as ecological and floristic similarities. Some of these habitats are clearly associated to a specific altitudinal interval (the high altitude alpine grasslands, and the low and mid altitude broadleaved evergreen woodlands for example) while others extend over a wider altitudinal gradient (alpine and subalpine scrubs) or are relatively independent of altitude because they are intimately associated to specific abiotic elements (inland cliffs, rock pavements and outcrops, screes, inland surface waters and mires, bogs and fens).
All species names of this dataset were validated using the Atlas Flora Pyrenaea (http://www.atlasflorapyrenaea. eu/src/home/index.php?idma=0) with the exception of those genera for which species identification was not explicit in the inventory (i.e., the apomictic Alchemilla or Hieracium). Most inventories followed the classical Braun-Blanquet format, but the scale used to record plant coverabundance varied between inventories. Consequently, we standardized all the values to the extended Braun-Blanquet scale, which ranges from 1 to 9 (van der Maarel, 1979).

| Phylogenetic inference
We built a genus-level phylogeny of the Pyrenean flora ( Figure 2) using the workflow proposed by Roquet, F I G U R E 2 Example of the genus-level phylogeny used in our study after randomly resolving species-level branches. Dots indicate the location of rare species in the tree. Endemic: species exclusive to the Pyrenees; HS: habitat specialists; LA: species with low local abundance; RGR: species with limited regional geographic range , and based on the species present in our data and in the atlas of the Pyrenean flora (FLORAPYR). We downloaded from Genbank three conserved chloroplastic regions (rbcL, matK, and ndhF) plus the ITS region for a subset of families, which we aligned separately by taxonomic clustering. We aligned all coding sequence clusters with MACSE (Ranwez, Harispe, Delsuc, & Douzery, 2011) and noncoding ones with MAFFT (Katoh & Standley, 2013), and trimmed all alignments with TrimAl (Capella-Gutierrez, Silla-Martinez, & Gabaldon, 2009). We concatenated all alignments to obtain a supermatrix. We then conducted maximumlikelihood (ML) phylogenetic inference analyses with RAxML (Stamatakis, 2014), applying the most appropriate partitioning scheme and substitution model obtained with PartitionFinder (Lanfear, Calcott, Kainer, Mayer, & Stamatakis, 2014) and a supertree constraint at the family-level obtained with the online software Phylomatic v.3 (tree R20120829). Specifically, we performed 100 independent tree searches. The best ML tree (the one with the highest probability) was dated applying the penalized likelihood method in treePL (Smith & O'Meara, 2012) and the following node calibrations: we fixed the node corresponding to the ancestor of eudicots at 125 Ma based on the earliest eudicot fossil (Hughes & McDougall, 1990), and applied minimum age constraints to 15 nodes based on fossil information extracted from Smith, Beaulieu, and Donoghue (2010) and Bell, Soltis, and Soltis (2010). To deal with unknown within-genera relationships, we simulated 10 scenarios of within-genera random branchings using a Yule process as implemented in the R package apTreeshape (Bortolussi, Blum, Durand, Francois, & Maliet, 2020). These 10 trees represent a distribution of possible hypotheses about evolutionary relationships in our dataset sensu Rangel et al. (2015). The only species lacking in our tree that were present in the inventories were those belonging to monotypic genera Cytinus, Ptychotis and Xatardia as the data available in GenBank was insufficient for proper phylogenetic inference. In order to explore and avoid potential biases derived from the limitations of our phylogenetic tree, we conducted our analysis with both our phylogeny and the ALLOTB tree published by Smith and Brown (2018), a publicly available and highly inclusive phylogeny of worldwide seed plants (See Supplementary files).
As previous studies have shown, inferences based on phylogenies with deep initial bifurcations such as the one found between the major groups of terrestrial plants may be inconsistent depending on how species are distributed among the branches of the phylogeny (Molina-Venegas, Aparicio, Lavergne, & Arroyo, 2015). In order to avoid this, all the analyses were conducted using only angiosperms, excluding gymnosperms and ferns from our communities and phylogeny.

| Taxonomic and phylogenetic alpha diversity
Taxonomic diversity of the studied communities was measured as species richness. Phylogenetic diversity (PD) was estimated as the sum of the length of the branches from the phylogeny associated with the species in each inventory (Faith, 1992). As the latter measure is highly correlated with species richness and biases comparisons between communities (Tucker et al., 2017), we standardized the observed values by subtracting from them the expected mean phylogenetic diversity value for communities of a certain species richness and dividing the result by the standard deviation of said value (hereafter referred to as standardized phylogenetic diversity or sesPD). The calculation and standardization of phylogenetic diversity values was done using the PhyloMeasures R package (Tsirogiannis & Sandel, 2017). In order to consider all possible evolutionary hypotheses, we computed sesPD using all 10 phylogenetic trees and averaged the results.

| Species rarity: Classification and contribution to habitat vulnerability
We classified each species present in our dataset based on the three criteria proposed by Rabinowitz (1981): geographic range, habitat specificity (HS) and local population size (LA). Given the difficulty of estimating total geographic range, we split it into endemic status (hereafter End, taxa restricted to the Pyrenean range), and regional geographic range (hereafter RGR, number of 10 km 2 UTM cells in which each species was found based on the location of the above mentioned relevés). Then, different thresholds were tentatively used to avoid overrepresentation of rare plants due to their regional geographic range. Finally, a species was considered to have a small regional geographic range (narrowly distributed) in the Pyrenees if its RGR was <5% of the maximum number of UTM cells occupied by any species in our study area. Habitat specificity was assessed from the number of habitats in which a species was found based on the habitat types associated with each relevé, with species having at least 90% of their occurrences in only one habitat being classified as habitat specific. In order to estimate each species' local abundance, we first calculated the specific average cover-abundance across all the inventories, along with its 95% confidence interval. We classified species as having low local abundance (locally scarce) if the upper limit of their cover-abundance confidence interval was lower than 2 in the extended Braun-Blanquet scale. Finally, we calculated the proportion of rare species per inventory from the number of species that fell within each of the rarity categories separately and any combination of them (hereafter referred to as rare species), and the total species richness of that inventory.
The contribution of rare species to the phylogenetic diversity of plant communities was calculated as the difference between the observed PD of an inventory and the PD resulting from removing any rare species from it, following a similar approach to Pool, Grenouillet, and Villéger (2014). We tested if this contribution differed from random expectation by comparing it with a null model. The null distribution of expected values for the difference in PD was obtained by randomly reassigning the rarity categories among the species in each inventory, removing any rare species from them and recomputing the difference in PD, and repeating this process 999 times. We then calculated the standard size effect by subtracting the expected difference in PD to the observed value and dividing the result by the expected standard deviation. A two-tailed test was applied, with values outside the range (À1.96, +1.96) considered as statistically higher or lower than expected by chance with a 95% confidence (Mazel et al., 2016). For this analysis we only used inventories containing at least one rare species in order to avoid the influence of communities without rare species, which comprise the majority of our inventories.

| Statistical analysis
Differences between habitats regarding species richness, standardized phylogenetic diversity, proportion of each rarity type and overall rare species as well as their contribution to phylogenetic diversity were analyzed with generalized linear models using a Poisson distribution for species richness, a Gaussian distribution for community sesPD and the contribution of rare species to it, and a binomial distribution for the proportion of rare species. All analyses were conducted using R version 4.0.4 (R Core Team, 2021). The statistical significance of general differences between habitats were assessed using analysis of variance with the car package (Fox & Weisberg, 2019), and pairwise comparisons were done by estimating marginal means and their confidence intervals via the emmeans package (Lenth, 2021). In addition, we tested the relation between both species richness and sesPD and the proportion of each type of rare species using Pearson's correlation coefficient. For comparative purposes, we tested the correlation between the results obtained with our phylogeny and those from the ALLOTB phylogenetic tree.

| Patterns of taxonomic and phylogenetic diversity across habitats
Species richness varied greatly between habitats, ranging from 483 species in inland surface waters to 1230 in woodland fringes, clearings and tall forb stands (Table 1). Mean species richness per inventory was 19.7 (SD 9.9), with the richest 1% of inventories found in alpine and subalpine grasslands, broadleaved deciduous forests and dry grasslands. Differences in species richness per inventory were statistically significant between habitats (χ 2 = 6843.1, df = 13, p < 0.0001; Figure 3a), with dry grasslands being the richest (26.34, 95% CI: 25.86 to 26.82) and inland surface waters the poorest (9.27, 95% CI: 8.82 to 9.75).

| Rare species: Differences in abundance across habitats, and relationship with diversity
Around 52% of the observed species showed some kind of rarity, and about one third of those were represented by nonendemic, nonspecialized and locally scarce species with small regional geographic ranges. Only seven species (0.3% of the species total) were included in the rarest category composed by endemic, specialized and locally scarce taxa with small regional geographic ranges.
Although rare species accounted, on average, for 22.64% (SD 6) of species found in each habitat, the number and type of rare plants varied widely between them, with dry grasslands having the highest proportion of rare species and arctic and alpine scrubs the lowest. Endemic species were more prevalent in arctic, alpine and subalpine grasslands, screes and inland cliffs, rocky outcrops and pavements. Species with small regional geographic ranges were more abundant in garrigues, dry grasslands, and inland surface waters. Habitat specialists concentrated in inland surface waters; inland cliffs, rocky outcrops and pavements and woodland fringes and clearings and tall forb stands. Locally scarce species were most prominent in dry grasslands, garrigues and inland surface waters (Table 1).
The relationship between diversity measures and the proportion of rare species was always significant, but coefficients were low and varied considerably between richness, sesPD and rarity types ( Figure 5). The proportion of rare species related negatively with species richness. The same applied for each type of rarity except for locally scarce species. The opposite trend was observed for sesPD and rare species proportion, with both variables relating positively except for locally scarce species.
The patterns of phylogenetic diversity we calculated with the ALLOTB phylogeny of Smith and Brown (2018) were highly correlated (r = 0.72) and consistent with our phylogeny ( Figure S1). Although the distribution of T A B L E 1 Number of inventories (N), rare species by rarity type, and total and rare species richness (pooled rarity types) per habitat; in parenthesis the percentage of rare species considering total richness sesPD in each habitat differed slightly between phylogenies ( Figure S2A), the main pattern of species-poor habitats like inland surface waters, screes and inland cliffs, rock pavements and outcrops holding great phylogenetic diversity, and species-rich habitats like dry grasslands having low sesPD was still present. Regarding the contribution of rare species to the phylogenetic diversity of their communities, the results varied between phylogenies (Table S1). However, this contribution did not differ significantly from randomness, supporting the results from our phylogeny ( Figure S2B).

| DISCUSSION
Here, we presented the first comprehensive description and comparison of taxonomic and phylogenetic plant diversity patterns, as well as the abundance of rare species and their contribution to habitat species richness and PD in a temperate mountain range, the Pyrenees. Our analysis revealed that rocky and aquatic habitats shelter the highest proportion of rare species. We also found contrasting patterns in the relationship between diversity and rarity, with the richest habitats harboring the least amount of rare species but the most phylogenetically diverse habitats showing the highest proportion of rare species. Interestingly, the loss of rare species would not translate into a stronger reduction of phylogenetic diversity compared with the loss of more common plants. The phylogenetic analyses were consistent across the two different plant phylogenies, which support the robustness of the patterns we observed.

| Diversity and rarity patterns across habitats
The observed patterns of species richness are partially consistent with the description of the Pyrenean flora made by G omez, , who found grasslands to harbor the largest number of species, followed by wetlands. We, on the other hand, found the highest number of species in grasslands, broadleaved deciduous forests and shrublands. These differences between our study and theirs might arise from the use of regional herbarium dataset instead of a systematic sampling, as well as slight differences in habitat classification (G omez,  used broader habitat categories, especially for grasslands). Concerning sesPD, the general patterns we found are in line with the phylogenetic patterns observed by Pardo et al. (2017)  Although more than half of the species present in our study were rare according to at least one criterion, these F I G U R E 4 Estimated marginal mean and its 95% confidence interval per habitat of rare species proportion per inventory (a) and standardized contribution of rare species to PD in said inventories (b). In both plots black lines and symbols represent data coming from species with one or more rarity types and the ones in color represent data from each type individually. HS: habitat specialists; LA: low local abundance; RGR: small regional geographic range F I G U R E 5 Relationship between proportion of rare species, richness and phylogenetic diversity based on Pearson's correlation coefficient. HS, habitat specialists; LA, Low local abundance; RGR, small regional geographic range. Brackets show 95% confidence intervals for each correlation index and asterisks indicate statistically significant values species were unevenly distributed among habitats. We observed a higher incidence of endemic species in rocky habitats and alpine grasslands, consistent with previous studies of the distribution of endemic plants in the Pyrenees Tejero, García, & G omez, 2017). Regarding the regional abundance of plant species in our study area, G omez, Lorda, Garmendia, and  observed that these species were located mostly in grasslands and wetlands, followed by rocky habitats. We found a similar pattern for wetlands and rocky habitats but not for grasslands, probably because the latter are much more abundant in the Pyrenees and thus species associated with them tend to have a broader distribution in the range. Another possibility is that we analyzed different types of grasslands separately while G omez, Lorda, et al. (2017) pooled all of them together. Regardless of this, the abundance of species with narrow areas of distribution could be attributed to multiple factors (Schemske et al., 1994). Recent speciation promoted by isolation in particular habitats, for example, can lead to species accumulating in specific areas if their dispersal capabilities are limited (Lesica, Yurkewycz, & Crone, 2006). Another opposite case scenario could be the reduction of a species' geographic range due to past environmental changes, which culminates in that species being secluded in refugial areas that are still suitable for them or have more stable conditions (Postigo Mijarra, Barr on, G omez Manzaneque, & Morla, 2009).
Habitat specialists in the Pyrenees tend to concentrate in distinct habitats with unique abiotic conditions compared with their surroundings like rocky cliffs, pavements and outcrops; screes; mires, fens and bogs, and tall forb stands. G omez, Lorda, et al. (2017) showed that specialist plants in the Pyrenees were mostly restricted to grasslands and rocky or humid habitats, and Boulangeat et al. (2012) found that specialist taxa in the Alps were also prone to accumulate in distinctive environments such as wetlands. With the exception of forb stands, these rocky and aquatic habitats tend to be poor in species, which is usually a consequence of harsh environmental conditions such as cold temperatures in alpine areas, damp or poorly developed soils in mires and fens, and low nutrient availability in rocky cliffs and pavements (Adamidis, Kazakou, Baker, Reeves, & Dimitrakopoulos, 2014;Kleidon, Adams, Pavlick, & Reu, 2009;Niedrist, Cantonati, & Füreder, 2018). The filters that these environments impose promote the presence of species exclusively adapted to them, leading to high proportions of specialized taxa (De ak et al., 2018;Pandit, Kolasa, & Cottenie, 2009). The phylogenetic patterns derived from these filters would depend on the distribution along the phylogenetic tree of the traits that allow plants to survive in such conditions. On the one hand, if these traits had evolved in just a few, phylogenetically close lineages, they would be located in specific areas of the phylogenetic tree, and thus the standardized phylogenetic diversity of habitats harboring these species would be low (Cavender-Bares, Kozak, Fine, & Kembel, 2009). This seems to be the case for dry grasslands, which are the only habitat with significantly low sesPD, although they have high species richness. On the other hand, if these traits had evolved several times in separate lineages through a convergent evolution process, these communities would be phylogenetically overdispersed, as the species able to withstand those environmental filters would be scattered throughout the phylogenetic tree (Cavender-Bares, Ackerly, Baum, & Bazzaz, 2004). None of the habitats in the Pyrenees shows significant phylogenetic overdispersion, however, those with the highest values of sesPD, namely inland surface waters and rocky habitats, tend to have a limited number of species due to strong environmental filters, suggesting that the adaptations to these environments may have evolved several times across the phylogenetic tree, leading to higher sesPD.

| Connecting species rarity and habitat vulnerability
Since species vary in their vulnerability to different factors depending on their rarity, the communities they inhabit will be more or less susceptible to diversity loss according to the type and abundance of rare species they harbor, although the consequences can differ between aspects of diversity. One of the most interesting results of our study is the contrasting contribution of rare species in habitats where they are frequent: their loss would translate into an important reduction of richness but would not translate into a significantly higher reduction of phylogenetic diversity than the loss of co-occurring and more common species. Nonetheless, their loss can lead to dramatic changes in ecosystem functioning and stability depending on its magnitude.
Studies on the relation between community stability and biodiversity highlight the importance of the latter in ensuring the functioning of ecosystems: those with less species tend to show reduced asynchrony in species responses to abiotic factors and less functional redundancy, which lead to decreased stability (Schäfer et al., 2019;Xu et al., 2021). While more diverse and redundant communities can still function after the loss of part of its diversity thanks to the functional redundancy of the remaining species (Thibaut & Connolly, 2013;Yachi & Loreau, 1999), less diverse communities have more difficulties for compensating the possible loss of keystone species and in consequence are more vulnerable to disturbances. If, in addition, less diverse communities have a higher proportion of rare species, which are considered to be more vulnerable to extinction (Kempel et al., 2020), and tend to play important roles in ecosystem functioning (Jain et al., 2014;Mouillot et al., 2013), disturbances in those communities could have severe consequences for their stability and function.

| Intrinsic vulnerability versus external hazards
It is noteworthy to mention that habitat vulnerability inferred from species rarity does not necessarily match habitat's exposure to global change drivers (Wilson et al., 2005). Our approach gives insight on a particular aspect of the intrinsic vulnerability of habitats due to the solely factor of species rarity. However, their vulnerability is still dependent on the combination of multiple internal and external factors like land use or climate change (Weißhuhn et al., 2018). According to the European Red List of Habitats (ERLH; Janssen et al., 2016), about half of the freshwater and grassland habitats in the European Union, along with more than 80% of mires and bogs are threatened by global change. Our rarity-based approach to vulnerability adds a new layer of concern to this assessment, given the high incidence of rare species in inland surface waters, mires, fens and bogs, and tall forb stands (which are included within the grassland category in the ERLH). Such coincidence translates into a double risk for these kinds of habitat, as they are highly vulnerable to both external and internal factors. Conversely, screes and inland cliffs, rocky pavements and outcrops, which our analyses pointed out as vulnerable, were deemed of least concern by the ERLH regardless of the sensitivity of the species in them. This highlights the importance of following a multifaceted course of action for conservation practices, integrating both intrinsic and extrinsic aspects of habitats and species (Dawson, Jackson, House, Prentice, & Mace, 2011).
In conclusion, our analyses provided insight on the patterns of taxonomic and phylogenetic plant diversity across habitats in a South European mountain range. They highlighted the role of distinct environments such as tall forb stands, aquatic, and rocky habitats as hotspots of both phylogenetic diversity and rare species. When taking into account the high proportion of rare species as an approximation to the intrinsic vulnerability of habitats, we found notable differences in habitat sensitivity. On the one hand, inland cliffs, rocky pavements and outcrops along with screes are particularly sensitive to demographic stochasticity due to small populations and low local abundances (Matthies et al., 2004). On the other hand, aquatic habitats and tall forb stands are sensitive to risks associated with low regional abundance and habitat specialization in addition to being among the most sensitive habitats to current global change drivers, as stated by the European Red List of Habitats. Through the study of both internal and external factors we can better identify the most vulnerable and priority habitats under the current climatic and land use change scenario. Our study highlights the importance of taking an integrative approach towards habitat vulnerability assessment, one that considers both internal and external drivers of vulnerability. It also shows how, in the absence of information about direct threats to species and habitats, accounting for rarity patterns could be a useful tool to guide conservation managers and policy makers.