Changes in plant species richness due to land use and nitrogen deposition across the globe

The effects of land use and atmospheric nitrogen (N) deposition on plant species richness are typically studied in isolation. Here, we quantified the combined effects of these two pressures on terrestrial plant species richness at a 0.25º spatial resolution across the globe.


| INTRODUC TI ON
Terrestrial plant communities are key to global biogeochemical cycles, the distribution of terrestrial animal biodiversity and the provision of benefits to people, including carbon sequestration and water regulation (Franklin et al., 2016;Vellend et al., 2017). Hence, it is crucial to understand how plant communities will change in response to human pressures on the environment. Land use and atmospheric nitrogen (N) deposition are two important pressures that affect plant biodiversity worldwide (Bobbink et al., 2010;Franklin et al., 2016).
The human use of land, for example for residential, industrial, or agricultural purposes, typically comes with the conversion or disturbance of the original natural habitat. This in turn may decrease the number of plant species compared to undisturbed habitats Newbold et al., 2015;Vellend et al., 2017). Some lowintensity types of land use (e.g. occasional logging, cities with a low human population density) may be characterized by plant species richness similar to or even higher than plant species richness in undisturbed vegetation. However, conversion of natural vegetation to agricultural land, which is the dominant human land use by far, typically results in species richness declines (Newbold et al., 2015;Vellend et al., 2017).
The effects of atmospheric N deposition (i.e. the input of reactive N from the atmosphere to the biosphere) on plant species richness are more equivocal, as both increases and decreases in species richness have been observed. In general, increases in plant species richness tend to occur in response to relatively low amounts of N deposition in regions where N is scarce, for example in remote areas (Bobbink et al., 2010). The number of species may also increase in well-buffered soils with high pH and cation exchange capacity, where the decline of sensitive species is slower and new species might colonize (Dise et al., 2011;Stevens et al., 2011;Van den Berg et al., 2011). In contrast, higher values of N input may lead to decreases in plant species richness, with soil eutrophication and acidification as important underlying mechanisms (Stevens et al., 2010;Vellend et al., 2017). Specifically, eutrophication may disproportionally promote the growth of nitrophilous plant species, which then outcompete others for light and resources, leading to an overall decrease in plant species richness (Dise et al., 2011;Lawrence, 2003;Midolo et al., 2019;Stevens et al., 2018). Soil acidification may lead to reductions in plant species richness particularly in soils with already low pH, as the pool of species adapted to low pH is relatively small (Bobbink et al., 2010;Dise et al., 2011;Stevens et al., 2010Stevens et al., , 2018. The emission and deposition of N are typically enhanced by land use-related activities such as fertilizer application (and volatilized livestock excreta), fossil fuel combustion (e.g. vehicle emissions), electricity production, or seasonal fires (Bobbink et al., 2010;Dise et al., 2011;Fowler et al., 2013;McClean et al., 2011). Hence, land use and associated activities act as a source of N emissions, while (semi-)natural areas are the sinks (Deng et al., 2019;Lv et al., 2019), resulting in spatial correlations between both pressures (Bowler et al., 2020;McClean et al., 2011). Thus, landscapes characterized by high-intensity human land use might be characterized by disproportional cumulative declines of plant species richness, because the direct on-site land-use impacts are complemented by impacts of increased N deposition in remaining (semi)-natural habitats. Yet, the landscape-level effects of both pressures combined are rarely studied and have, to our knowledge, never been assessed at a continental to global extent.
Here, we aimed to quantify the relative and combined effects of land use and N deposition on landscape-level plant species richness across the globe. To that end, we developed an approach based on the countryside species-area relationship (cSAR) model, which has been used to study impacts of land use on biodiversity (Gerstner et al., 2017;Pereira et al., 2014). The cSAR model quantifies changes in species richness due to land use, accounting for the affinity of species to different human-modified habitats (Pereira et al., 2014).
We extended the cSAR to account also for the effect of nitrogen deposition, arriving at a new multi-pressure species-area relationship (mp-SAR) model. We accounted for the effect of nitrogen deposition by quantifying the proportion of plant species able to persist under different levels of N addition, similar to how land-use impacts are quantified in the cSAR. We determined the proportional changes in plant species richness for different land-use types and N deposition values through meta-analyses of local monitoring data obtained from the literature. We then combined the site-level responses with global land use and N deposition maps in the new mp-SAR framework to provide estimates of changes in plant species richness at a resolution of 0.25º (about 25 km at the equator) worldwide.

| General approach
Our multi-pressure SAR (mp-SAR) model calculates the change in species richness within a spatial unit (grid cell) based on the areas within that spatial unit affected by different pressures combined with the affinity of the species group to these pressures, as follows (where j and k denote two different pressures): where RSR new is the relative species richness compared to a reference situation (natural habitat), A tot is the total area of the spatial unit, A j,k is the area within that unit that is subject to pressures j and k, h j,k represents the affinity of the species group to pressures j and k, and z is the slope of the species-area relationship. Following the countryside SAR model, the affinity h j,k is derived from the ratio between species richness in a human-modified habitat and an original or natural situation (Pereira et al., 2014): where S j,k represents the average number of species in the habitat modified by the pressures j and k, and S n represents the average number of species in the natural habitat. The affinity value is higher than 1 if the number of species is higher in the modified than the natural habitat. Ideally, the affinity of species groups to the combination of pressures j and k is determined based on observational or experimental data covering both pressures simultaneously (e.g. data from multifactorial experiments). As these data are scarce, we propose to retrieve the affinity to pressure combinations from single-pressure studies based on the assumption that the pressures act independently, implying that we can apply response addition (Plackett & Hewlett, 1952;Vinebrooke et al., 2004) as: Note that this is a generic equation that also applies if one of the two pressures is absent, as then the corresponding species richness ratio equals 1.
To quantify the impacts of land use and N deposition on plant diversity across the globe, we first assembled local monitoring and experimental data representing relative plant species richness in relation to different land-use types (S j ∕S n in Eq. 3) and to different levels of N addition (S k ∕S n in Eq. 3). We then applied the mp-SAR (Eq. 1) by quantifying the areas affected by land use, N deposition and both pressures combined within each 0.25º grid cell (A j,k in Eq. 1) and the combined affinity to land use and N deposition in that grid cell (h j,k in Eq. 1, 3), using slope values (z) specific to the biome that the grid cell belongs to. We further detail each of these steps below.

| Quantifying plant species richness responses to land use
To quantify the response of plant species richness to land use (S j ∕S n in Eq. 3), we first merged four existing databases with comparisons of plant species richness between specific land-use types and natural habitats (de Baan et al., 2013;Elshout et al., 2014;Gerstner, Dormann, Václavík, et al., 2014;Hudson et al., 2014). As these datasets together cover publications up to 2015, we searched for additional, more recent publications (from 2015 to 2020) using the following search key within the ISI Web of Knowledge database: (TS = (("land use") AND (plant OR plants) AND ("species richness" OR "species composition" OR "species abundance"))) AND  S1). A list of the data sources is found in Appendix S1. We categorized the different land use types and intensities following Hudson et al. (2014), whereby we merged some categories (Table 1).
To estimate the plant species richness response to different landuse types, we performed a meta-analysis of the natural logarithm of the response ratios (RR) as extracted from the studies (Hedges et al., 1999): where S j,i and S n,i represent the mean species richness in a land-use type j and the corresponding natural habitat n for a study i, respectively. Thus, the log response ratio quantifies the proportionate response of plant species richness to land use, based on an effect size measure standardized across studies (Hedges et al., 1999). We analysed the ln RR in relation to land-use type using linear mixedeffects models fitted with the package "metafor" (Viechtbauer, 2010). Following Konstantopoulos (2011), we included a random intercept with each observation nested within study and publication to account for the possibility that the underlying true effects within grouping level are not homogeneous. Moreover, we accounted for non-independence of observations sharing a common natural reference through the variance-covariance matrix (Lajeunesse, 2011).
To account for uncertainties in the observations, we used the inverse of the variance of each observation per study i to weight the ln RR, with the variance calculated as (Hedges et al., 1999): where SD 2 h,i and SD 2 n,i represent the standard deviations of S h,i and S n,i , respectively, and N h,i and N n,i are the sample sizes (Hedges et al., 1999).
For studies that did not report estimates of variation (34% of the studies), we imputed the SD with the "Bracken1992" method using the coefficient of variation from all complete cases (Bracken & Sinclair, 1992).
Finally, we investigated whether our results are influenced by publication bias, that is a bias in peer-refereed journals to publishing statistically significant results (Nakagawa & Santos, 2012). To assess possible publication bias, we used a funnel plot and an Egger's test to assess the funnel plot's asymmetry using the meta-analytic residuals with our multi-level random effect structure and the precision (1/SE) as covariate (Fernández-Castilla et al., 2021;Nakagawa & Santos, 2012). We did not find signals of publication bias (Egger's test p-value = 0.26, Fig. S2) suggesting that our model outcomes are representative (Fig. S2). For prediction purposes, we used the average estimates from the model ( Figure 1, Table S1) to calculate the land-use affinities (Eq. 2).

| Quantifying plant species richness responses to nitrogen deposition
To quantify the response of plant species richness to N deposition

| Application of the multi-pressure SAR
We applied the mp-SAR to quantify the combined effects of land use and N deposition on plant species richness at a 0.25º x 0.25º resolution (~25 km) globally, using biome-specific slopes (z-values,  Figure 1). We assigned the use intensity of croplands using the amount of N application (kg ha −1 yr −1 ) as retrieved from the LUH2 data for the year 2015 (Schipper et al., 2020). To that end, we classified cropland use intensity as minimal (<100 kg ha −1 yr −1 ), light (100-250 kg ha −1 yr −1 ), or high (>250 kg ha −1 yr −1 ; Temme and Verburg, (2011)).
TA B L E 1 Classification of land use types and intensities (modified from Hudson et al. (2014))

Land use Description
Primary vegetation Natural (original) vegetation with no evidence of prior destruction and minor (if any) disturbance by humans. It can be forest or non-forest.

Secondary vegetation
Regenerating vegetation upon the destruction of the primary vegetation with currently minor disturbance by humans.

Plantations
Plantations of trees or shrubs for subsistence or commercial use. Including extensively managed or mixed timber, fruit/coffee, oil-palm, or rubber plantations in which native understorey and/or other native tree species are tolerated, which are not treated with pesticide or fertilizer, and which have not been recently (< 20 years) clear-felled; and monoculture fruit/coffee/rubber plantations with pesticide input, or mixed species plantations with pesticide inputs. Present signs of clear-felling (trees of different or same ages).
To quantify the amount of N deposition, we used the recently  1987-1993, 1997-2003, 2007-2013). We then resampled the new long-term average N deposition map to a 0.25º resolution using a bilinear interpolation (Fig. S4). To include the modifying effects of CEC and mean annual temperature on plant species richness responses to N deposition, we obtained CEC values at a 250 m resolution as averages across soil depths of 0-5, 5-15 and 15-30 cm from SoilGrids (Fig. S4, Hengl et al., 2017) and resampled the map to a 0.25º resolution by averaging the values.
Similarly, we averaged monthly temperature values for each year in the period 1984-2015 from the global Climate Research Unit database, then averaged over the years and resampled the result to a 0.25º resolution using a bilinear interpolation (Fig. S4, Harris et al., 2020).

| RE SULTS
Worldwide, land use and N deposition combined resulted in a species richness decline of 26 ± 12% (mean ±standard deviation) across the analysed grid cells (Figure 2; Table S3). Europe was the most impacted continent with an average species decline of 34% (± 8%) due to both pressures combined, while South America was the least impacted continent (16% ±11%) (Figure 2). In North America, Africa and Asia, the average impact of the two pressures combined was 25% (± 11%), 22% (± 12%) and 26% (± 12%), respectively ( Figure 2; Table S3). Land use was the main driver of species decline globally and per continent, with a global average impact of 19% (± 11%)  Table S4).

| DISCUSS ION
In this study, we developed and applied a novel SAR-based model to quantify the effects of land use and N deposition on plant species richness in 0.25° grid cells across the globe. As for any SAR-based model, our results reflect projected changes in plant species richness rather than actual (instantaneous) losses or gains (Pereira et al., 2014), because losses or gains of species in response to environmental change typically occur with a lag time (leading to "extinction debt" or "colonization credit," respectively; Hanski and Ovaskainen, (2002)). This delay is currently not accounted for by SAR models (Lewis, 2006), which implies that our estimates are not directly comparable to empirical observations. However, our model is based on comprehensive meta-analyses of site-level effects of both N deposition (Midolo et al., 2019) and land use (de Baan et al., 2013;Newbold et al., 2015), providing a solid basis for the mp-SAR and our results.
We estimated that on average across the globe, landscape-level plant species richness might eventually be reduced by 26%  However, the increase in species richness due to N deposition is typically caused by the colonization of generalist nitrophilous species (Dise et al., 2011;Lawrence, 2003;Perring et al., 2018). These colonization events and the overall increase in species richness may mask the decrease of specialist species that are outcompeted, as well as the resulting ecosystem homogenization (Perring et al., 2018;Stevens et al., 2010). Moreover, possible future agricultural intensification as well as an increasing prevalence of seasonal fires to create and prepare agricultural land may increase N deposition levels in tropical regions (Bauters et al., 2018;Chen et al., 2010), which may turn plant species richness increases to declines. To better understand the effects of nitrogen addition on plant biodiversity, future research should include additional characteristics of the plant communities (e.g. abundance of individual species, or species turnover), because species richness may mask diverging responses of individual species (Fleishman et al., 2006;Schipper et al., 2016).
Our approach is not without uncertainties. First, we assumed a time frame of 32 years of cumulative N deposition based on the N deposition data available (Ackerman et al., 2018). However, impacts may accumulate over longer times, resulting in larger declines (Fig. S3;Midolo et al. (2019)). Further, we assumed a constant level of N deposition across the 32 years, thus not accounting for recent reductions in N deposition as observed for example in Europe. The potential recovery of vegetation communities upon reductions in N inputs is, however, typically slow or may not take place at all as long as the remaining inputs are above critical levels (Dise et al., 2011;Stevens, 2016). Second, we made several assumptions with regard to the combined effects of land use and nitrogen deposition. For  secondary vegetation, we assumed that land use and N deposition act independently and that species' sensitivities to both pressures are uncorrelated, which justifies calculating the combined effect based on the assumption of response addition (Plackett & Hewlett, 1952;Vinebrooke et al., 2004). However, the assumption of uncorrelated sensitivity might be too simplistic. For example, N-tolerant generalist plant species might be tolerant also to land cover change, while certain specialist species might be sensitive to both pressures. In case of positively correlated sensitivities to both pressures, the response addition approach leads to an overestimation of their combined effect (Vinebrooke et al., 2004). If we assume no further decline of species richness due to N deposition in secondary vegetation (i.e. the land cover type assumed to be affected by both pressures), the average decline changes from 6% for the combined effect of both pressures to 5% for land use only (Fig. S6). For cropland, pastures, plantations and urban areas, we assumed that the direct effects of land use override the effects of atmospheric nitrogen deposition. We made this assumption because experimental data were not available to quantify nitrogen impacts in these land-use types (Midolo et al., 2019).
However, N deposition may still have an additional impact on plant species richness in these land-use types, particularly if nitrogen input from alternative sources (e.g. manure or artificial fertilizers) is low (Kleijn et al., 2009;Wüst-Galley et al., 2021). To improve our model, it would be highly beneficial to collect and analyse experimental or monitoring data that allow for quantifying and disentangling the separate and combined effects of N enrichment (from both deposition and fertilization) and land use (excluding fertilization) on plant species richness in comparison to natural habitat. This need is most prominent for croplands or pastures, which make up a large share of the anthropogenic land. Finally, we note that we applied our model only to grid cells within its applicability domain (i.e. 33% of all the 0.25º cells in the world), which is not necessarily a representative sample of the globe as a whole. Applied to the entire globe, our model predicts lower species declines on average, due to the inclusion of areas with low N deposition levels and little land use (global mean species decline of 18% instead of 26%; Fig. S7 and Fig. S8). However, the global mean impact of nitrogen deposition remains similar (Fig. S8) and land use remains the main driver of plant species richness decline.
Our study provides a new methodology to quantify impacts of multiple pressures on landscape-level plant species richness. Our new mp-SAR approach can be expanded to predict the combined effect of other human pressures, such as climate change or invasive species, which are currently lacking in the SAR framework (Hanski et al., 2013;Pereira et al., 2014). The required affinity values can be quantified based on experimental data, based on field observations obtained along a pressure gradient, or a combination of experimental and monitoring data, provided that the data include also records in natural habitat (Koricheva & Gurevitch, 2014). Given the paucity of studies quantifying impacts of multiple pressures (Bonebrake et al., 2019;Bowler et al., 2020;Franklin et al., 2016), our current methodology may fill a gap and, ultimately, help to underpin conservation policies. Tackling multiple pressures would certainly increase the effectiveness of conservation policies and can even be imperative for success, as removing a dominant pressure could simply reveal the impacts of remaining pressures without a net biodiversity gain (Bonebrake et al., 2019;Bowler et al., 2020). For example, our results suggest that restoring agricultural land to natural vegetation may not be fully effective in regions with high levels of N deposition, such as eastern Asia, the United States and Europe. The consistent negative impact of land use together with the long-term effect of N deposition in those areas point at a need to tackle both pressures simultaneously if nature is to be conserved or restored, for example through a combination of alternative farming practices (e.g. organic farming), dietary changes (notably a reduction of meat consumption) and targeted N emission reduction measures (Dise et al., 2011;Mace et al., 2018;Pe'er et al., 2020). These efforts may need to be supplemented by active restoration measures to speed up or even enable the recovery of damaged habitats (Stevens, 2016). Mitigating or reducing the combined impacts of N deposition and land use on plant communities is not only needed to halt the ongoing decline of biodiversity, but also the associated services to human society.

ACK N OWLED G EM ENTS
This study was conducted as part of the GLOBIO project (www.globio.info). We are thankful to Martijn Schaap for pointing us at global N deposition maps and to Gabriele Midolo for helping us to model the impact of N deposition on plants. We thank Melinda de Jonge for helping us to double-check our classification of observations into different land-use classes.

DATA AVA I L A B I L I T Y S TAT E M E N T
All data and R code necessary to reproduce the analysis are available from the Dryad Digital Repository: https://doi.org/10.5061/dryad.  1984-1986, 1994-1996, 2004-2006, and 2014-2016. https