Threats of land use to the global diversity of vascular plants

Land use is a main driver of biodiversity loss worldwide. However, quantifying its effects on global plant diversity remains a challenge due to the limited availability of data on the distributions of vascular plant species and their responses to land use. Here, we estimated the global extinction threat of land use to vascular plant species based on a novel integration of an ecoregion‐level species‐area model and the relative endemism richness of the ecoregions.


IUCN Red
List, are mostly available for mammals, birds, and amphibians (Cazalis et al., 2022;Davidson et al., 2017;Santini et al., 2019;Tracewski et al., 2016). For the majority of the vascular plant species, however, there is not enough data for assessing extinction risk at the species level. For instance, only 14% of the flowering plant species, a subgroup of vascular plants, have had their conservation status assessed by the IUCN Red List (IUCN, 2021).
In the absence of species-specific data on distribution and population size, extinction threats from land use can be estimated with species-area relationship (SAR) models (Chaudhary et al., 2015;Chaudhary & Brooks, 2018;Pereira & Daily, 2006). These models can be used to assess potential species losses due to land use in specific biogeographical regions (Van Vuuren et al., 2006). However, extinction within a given biogeographical region does not necessarily equate global extinction, as a species may occur in multiple regions.
The latter also implies that simply summing up regional extinctions may overestimate global extinctions because of double counting (e.g., Van Vuuren et al., 2006). Estimating global extinction threat based only on species endemic to the biogeographical regions of concern comes with the reverse problem: the estimated global extinction might be underestimated, as it excludes non-endemic species that could go extinct globally if they would become extirpated across all regions of occurrence (e.g., Brooks et al., 2002;Chaudhary & Brooks, 2019;Jantz et al., 2015).
Recent studies have attempted to tackle this issue by combining SAR-based estimates of potential total species loss per region with a weighting factor representing the aggregated global extinction risk of the species within the region based on their range size and IUCN threat status (Chaudhary et al., 2015;Chaudhary & Brooks, 2019;Kuipers et al., 2019). However, this approach requires species-specific data on extinction vulnerability and occurrence range, which is not available for the majority of the vascular plant species. Another approach proposed recently is to incorporate compositional similarity (beta diversity) in the estimation of extinction threats (Di Marco et al., 2019). This permits local or regional extinction threats to be added up to global extinction threats because it accounts for species shared among different locations. Di Here, we present a novel approach to estimate the threats of land use to the global diversity of vascular plants. To that end, we integrate countryside SAR-models with endemism richness. We use the countryside SAR models to quantify extinction threat due to land use per biogeographical unit and then scale the regional extinction threats to global extinction threat using endemism richness. Endemism richness is a metric that combines the species richness and endemism of an area or biogeographical region (Kier & Barthlott, 2001). A region has a high endemism richness if there is a large number of species occurring exclusively in that region. Endemism richness is especially useful because the sum of the endemism richness across all regions equals the total global number of species. Therefore, we can interpret endemism richness as the contribution of an area or region to the global biodiversity (Kier & Barthlott, 2001). We based our analysis on ecoregions as biogeographical units because they are relatively homogeneous regions characterized by specific flora (Olson et al., 2001). In addition, the use of ecoregions as a biogeographical unit simplifies calculations and thus reduces the model's computational requirements. We calculated the endemism richness of 818 ecoregions and combined it with the ecoregion-specific extinction threat estimates to estimate the global species extinction threat due to land use. Next, we allocated the ecoregion-specific contribution to global extinction threat to land use types using a high-resolution land-use map. Our approach is unique in quantifying the contribution of each land use type and intensity in each ecoregion to the global extinction threat of vascular plants. Knowing these contributions and how they vary across ecoregions is a prerequisite for designing effective policy measures to achieve national and international objectives regarding biodiversity conservation.

| Global extinction threat
Inspired by Pereira and Daily (2006) and Kier and Barthlott (2001), we estimated global extinction threats to vascular plants due to land use by summing the ecoregion-specific extinction threats corrected for the relative endemism richness, over all ecoregions: where GET is the global extinction threat due to land use; RE j is the relative endemism richness in ecoregion j ( Figure S1), and the part in parentheses is the extinction threat for ecoregion j according to the countryside SAR model, where A new,j is the original habitat area remaining in the ecoregion j after conversion to land use; h i,j is the affinity of vascular plants to the land use type i in ecoregion j; A i,j is the area of land use type i in ecoregion j; z j is the slope of the species-area curve in the biome that ecoregion j is located, and A org,j is the original habitat area in ecoregion j. GET is a dimensionless indicator ranging between 0 and 1, where 0 means that there are no species threatened with global extinction and 1 means all species are threatened with global extinction.
The relative endemism richness is the sum of the endemic weights of all species occurring in ecoregion j relative to the global species richness: where S org,g is the global species richness and w s,j is the species endemic weight of a species s occurring in ecoregion j, which is defined as the inverse of the total number of ecoregions in which species s occurs.
Because we do not have data on the distributions of species within ecoregions, we prefer calculating the endemism weights based on the presence and absence of species in an ecoregion rather than the areas of the ecoregions that the species occurs in.
The affinity (h i,j ), introduced in Equation (1), ranges between 0 and 1, where 0 means that the taxonomic group has no affinity to the land use type i and will not occur there, whereas an affinity value of 1 means that the land use type is equivalent to natural habitat (which has an affinity value of 1 by definition). The affinity value is calculated as (Martins & Pereira, 2017;Pereira et al., 2014): where RR i is the local relative species richness in land use type i; S i is the local species richness in land use type i, and S org is the local species richness in a natural habitat.

| Global extinction threat per land use type
The ecoregion-specific contribution to global extinction threat (term in Equation 1  where GET i is the global extinction threat imposed by land use type i, and GET j is the contribution of ecoregion j to the global extinction threat (i.e., Equation 1 without the summation over the ecoregions).
The reasoning for preferring allocating impacts to land use types based on the local relative richness (RR i ) rather than affinity (h i,j ), which is also sometimes used (Chaudhary & Kastner, 2016), can be found in the Appendix S1.

| Model parameterization
To determine the vascular plant species richness per ecoregion (S org,j ) and globally (S org,g ), we extracted vascular plant species occurrence data from the Global Biodiversity Information Facility (GBIF, 2021b).
We removed records without geographic coordinates, records of cultivated and invasive species, and records of fossil material and records with spatial issues, such as a mismatch between coordinates and verbatim location (country name). Next, we substituted synonyms by the corresponding accepted scientific name using the GBIF Backbone Taxonomy (GBIF, 2021a), in order to avoid doublecounting species. We overlaid the geographical coordinates of the records with a map (shapefile) of the ecoregions (Olson et al., 2001) and with a land use map (see below). Finally, we calculated the number of unique vascular plant species for all ecoregions with land use data available and for the entire globe. The full description of how we cleaned and processed the GBIF data can be found in the Appendix S1.  (Table S3).

| Sensitivity analysis
We analysed the sensitivity of the estimated global species extinction threat to variability in the z values, the local relative richness in relation to land use, the relative endemism richness, and the total area of primary vegetation by a data randomization procedure. We randomized each of the four variables 100 times by randomly reordering the values of each variable (without replacement of elements) and calculated corresponding global extinction threats. We then calculated the Pearson correlation (r) between species extinction threats estimated by the model with the input variable of concern randomized and the default model (i.e., non-randomized data) and calculated the importance of the variable as 1−r (Thuiller et al., 2021), averaged across the 100 repetitions. Thus, the smaller the correlation between the results from a model with the variable randomized and the default results (resulting in a high value for 1−r), the more important the variable. We randomized each parameter separately as follows: we randomly swapped z-values (z j ) among biomes; we randomly swapped local relative richness (RR i ) among land use types, with exception for primary vegetation because it is the reference to which the other land uses are compared; we randomly swapped the relative endemism richness (RE j ) across ecoregions; and we randomly swapped the percentage of natural (% A org ) and nonnatural (% A i ) land among ecoregions while keeping the original total area of each ecoregion unchanged.
To test the sensitivity of our results to the land use input map, we performed the calculations also based on land use data from the Land use Harmonization (LUH2) data set (year 2015), with a spatial resolution of 0.25° (Hurtt et al., 2020). We harmonized the land use categories from the GLOBIO and LUH2 maps, resulting in the following categories: primary vegetation, secondary vegetation, urban, forestry, cropland of minimal intensity, cropland of high intensity, and pasture of minimal and high intensity (Table S1). We note that secondary vegetation is not present in the GLOBIO map and forestry is not present in LUH2. As the LUH2 land use map does not distinguish land use intensities, we classified pasture and cropland intensities following the same procedure as used to compile the GLOBIO map (Schipper et al., 2020). Specifically, we identified pasture and cropland of high-intensity use based on the inorganic nitrogen fertilization rate, assuming that a rate above 100 kg/ha characterizes high-intensity use (Temme & Verburg, 2011). Pasture and cropland with a nitrogen fertilization under 100 kg/ha were classified as minimal intensity.

| Global extinction threats
Overall

| Sensitivity analysis
According to our sensitivity analysis, the difference in global extinction threat between ecoregions depends more on the differences in relative endemism richness and percentage of natural versus F I G U R E 1 Percentage of vascular plant species threatened with global extinction due to land use per ecoregion. The sum of the percentages across the ecoregions equals the global total percentage of species threatened with extinction (11%). Note that intervals differ among classes for visualization purposes. The map is displayed based on the Mollweide equal area projection.
non-natural habitat area than on the differences in the local relative richness per land use type. The estimate of relative endemism richness clearly has the greatest influence on the estimated extinction threat. The model results were least sensitive to differences in z-values ( Figure 3a).
Overall, the results based on the LUH2 instead of the GLOBIO land use map were similar, with 13% instead of 11% of global vascular plant species extinction threat, respectively. Further, the estimated global extinction threats per ecoregion were highly correlated between the models using GLOBIO and LUH2 maps (r = .84; Figure 3b). However, with LUH2 as land use map input, the most threatened regions differ partially: Madagascar, South Africa, and the Western European broadleaf forests stand out due to estimated extinctions threats ( Figure S2). Based on the LUH2 land use map, secondary vegetation and pasture of high intensity pose the highest threat ( Figure S3).

| Global extinction threats to vascular plants due to land use
We estimate that 11%  that has identified these areas, among others, as regions with a high estimated extinction risk for plants and vertebrates (Chaudhary & Brooks, 2018;Davidson et al., 2017;Di Marco et al., 2019;Tracewski et al., 2016). In order to reduce biodiversity loss, it is important to take conservation measures particularly in these areas. In the Global South, where cropland of minimal intensity contributes the most to extinction threats, cropland abandonment is a widespread phenomenon (Laue & Arima, 2015;Wang et al., 2015;Yin et al., 2020). If well-managed, abandoned cropland can offer an opportunity to promote biodiversity conservation (Chazdon et al., 2020). Abandoned croplands can be restored to native landscapes or be re-used as pastures (Knoke et al., 2014), which could avoid the deforestation of new areas. Other measures to avoid land conversion include protecting the remaining wilderness areas, improving the management of existing protected areas, and targeting demands and supply chains (Maxwell et al., 2020;Winkler et al., 2021).

| Challenges
Validating global extinction threat models is challenging as the results cannot be checked against field observations. Some of the estimated extinctions might have already happened, while other species extinctions have not yet occurred because of the lag time between impact and extinction (relaxation time) (Cronk, 2016). A tentative validation can be obtained by comparing the estimated proportion of extinctions with the proportion of species categorized as threatened or extinct by the IUCN (Chaudhary & Brooks, 2018, 2019De Baan et al., 2013;Koh & Ghazoul, 2010;Kuipers, Hilbers, et al., 2021). According to the IUCN (2021), 6% of vascular plant species are categorized as "Extinct," F I G U R E 3 (a) Sensitivity of global extinction threat per ecoregion to the variability in the relative endemism richness, percentages of natural versus non-natural habitat area per ecoregion, local relative richness (RR) per land use type, and z-values. The variable importance is the average of 1−r across all simulations, where r is the Pearson correlation between estimated species extinction threats from a model using randomized input values for the variable of concern and the model outputs based on the default (i.e., non-randomized) data. (b) Comparison of estimated extinction threats (%) per ecoregion between models using LUH2 and GLOBIO as land use input maps. The sum of the percentages across all the ecoregion equals the global total estimated extinction threat (11% for GLOBIO and 13% for LUH2). The solid black line represents the line of equality.
"Extinct in the wild," "Critically endangered," "Endangered," and "Vulnerable." Of these, two-thirds (4% in total) are threatened or extinct due to threats related to land use. This suggests that our result of 11% is an overestimate. However, the IUCN estimate is based on only 15% of the vascular plant species that underwent an extinction risk assessment, which is not necessarily a representative sample. The IUCN Sampled Red List Index, which gives an extinction risk indicator for a whole group based on the sampling of assessed species (Baillie et al., 2008), estimates that 22% of vascular plant species are threatened with extinction (Brummitt et al., 2015). Assuming that two-thirds of these are threatened due to land use, this would imply that 15% of vascular plant species are threatened with extinction due to land use, which is close to our estimate of 11%.
We encountered several challenges related to species data availability. We used GBIF data as input to determine the endemism richness and ecoregion-and global-level vascular plant species richness, yet GBIF does not include all known species. For instance, in our dataset, only 6% of the plant species are endemic to the Cerrado region in Brazil, whereas, according to floristic inventories, the percentage of endemic vascular plant species is around 17% (Forzza et al., 2012). The fact that our data set is incomplete means that we likely underestimate the global extinction threats due to land use per ecoregion. In addition, spatial accuracy issues and the inclusion of non-native species in the GBIF data may lead to overestimating species richness in certain ecoregions (Maldonado et al., 2015). Although we excluded records with spatial issues and of cultivated and non-native species, we acknowledge that there is a certain level of uncertainty inherent to opportunistic biodiversity data, which may have affected our species richness and endemism richness estimates. This, in turn, may contribute to explaining why our estimates of relative endemism richness differ from those of Kier et al. (2009), who used a combination of species lists and species distribution data to calculate endemism richness, as well as a different delineation of regions (90 biogeographical regions instead of our 818 ecoregions). Finally, our approach is based on the average sensitivity of plant species assemblages to land use (RR i ), ignoring the variability of responses among biomes (Newbold et al., 2020) and species (Dalle Fratte et al., 2019). Our results likely underestimate the global extinction threat in Tropical and Mediterranean biomes, as these biomes are more sensitive to land use (Newbold et al., 2020). In order to reflect the variability of plant responses to land use, future research could derive local relative species richness for different land use types at the biome or ecoregional level.

| Future directions
The method we developed to estimate global extinctions due to land use intensities can also be applied to estimate global extinctions threats due to other pressures that affect habitat quality and for which affinity values are available, such as nitrogen deposition (Gallego-Zamorano et al., 2022) and climate change (Newbold et al., 2020). In addition, our method can be used to evaluate future scenarios, like the Shared Socioeconomic Pathways (SSPs), or to update the results of the present study with improved species richness estimates, as many new plant species are described each year (Cheek et al., 2020) and new data are constantly added to GBIF (Cornwell et al., 2019). Finally, the method can be applied to other taxonomic groups that are not comprehensively assessed by IUCN but are sufficiently represented in GBIF and in databases needed for calculating the affinity of the species assemblage to the pressure of interest. Applying our method to multiple taxonomic groups, pressures, and scenarios would allow for a systematic prediction and comparison of global extinction risks in accordance with a consistent framework grounded in ecological theory. This would be an asset for global biodiversity assessments, in particular for species groups where data limitations do not allow for performing specieslevel extinction risk assessments.

ACK N O WLE D G E M ENTS
We thank C.J.M. Musters for his contribution to the overall idea and the methodology development and for reviewing multiple versions of the manuscript. We are grateful to J. Gallego-Zamorano for sharing his code and results regarding land-use impacts on plant species assemblages before his work was published. AMS acknowledges funding from PBL Netherlands Environmental Assessment Agency as part of the GLOBIO project. We gratefully thank LP, MCZ, MH, and MAJH for acquiring funds for the PhD project of HM, of which this research is part.

CO N FLI C T O F I NTE R E S T S TATE M E NT
None.

DATA AVA I L A B I L I T Y S TAT E M E N T
The raw data that support the findings of this study are openly avail-