Local nutrient addition drives plant diversity losses but not biotic homogenization in global grasslands

Experimental setup

The experimental sites used in this study are part of the Nutrient Network (NutNet, Fig. S1 and Table S1). The experimental design includes a factorial manipulation of nutrients (N, P, and K) plus two fences to exclude herbivores, see ref. ³⁹ for more details. For the analyses here, we used plots under two treatments: Ambient (Control) and fertilization by nitrogen, phosphate, and potassium together (i.e., NPK). Treatments were randomly assigned to 5 m × 5 m plots and were replicated in three or more blocks. A micronutrient mix consists of Fe (15%), S (14%), Mg (1.5%), Mn (2.5%), Cu (1%), Zn (1%), B (0.2%), and Mo (0.05%) was added once only at the start of the experiment (i.e., year 1) for the nutrient addition plots, but not in subsequent years to avoid toxicity. Nitrogen, phosphate, potassium were added annually before the growing season of each treatment year at most sites. Nitrogen was added as 10 g m⁻² yr⁻¹ time-release urea [(NH₂)₂CO], phosphate was added as 10 g m⁻² yr⁻¹ triple-super phosphate [Ca(H₂PO₄)₂], while potassium was added as 10 g m⁻² yr⁻¹ sulfate [K₂SO₄].

Data were retrieved from the NutNet database in November 2023. We analyzed data from 72 sites where 1) nutrients were applied for at least four years; and 2) each site had at least three blocks. These sites are distributed across six continents and include a wide range of grassland types. See Fig. S1 and Table S1 for details of geolocation, grassland types, and experimental years used.

Sampling protocol

Scientists at NutNet sites followed standard sampling protocols³⁹. Specifically, a 1 m × 1 m subplot within each plot was permanently marked for annual recording of plant species composition. Species cover (%) was estimated visually for individual species in the subplots; thus the total cover of living plants may sometimes exceed 100% for multilayer canopies. At most sites, cover was recorded once per year at peak biomass. At some sites with strong seasonality, cover was recorded twice per year to include a complete list of species. For those sites, the maximum cover for each species and total biomass were used in the analyses. When taxa could not be identified to the species level, they were aggregated at the genus level but referred to as “species” for simplicity.

Quantifying changes in α, γ, and β diversity

We measured α and γ diversity using species richness (i.e., number of species) because it is the most commonly examined diversity metric⁴³. At each site, α diversity was estimated as the number of species in each permanent subplot (1 m × 1 m), and γ diversity as the total number of species occurring in three permanent subplots (for each treatment separately). To standardize sampling effort, for sites with more than three blocks, we selected the first three blocks according to the block number recorded by site PIs. The framework relies on Whittaker’s multiplicative β diversity partition, and it quantifies β diversity using the effective number of communities¹². As such, if all subplots share the same species, then β diversity would equal to one. In contrast, if each subplot has unique species, then β diversity would equal to three. We calculated ∆α as the richness difference in local communities (subplots) and ∆γ as the difference in the sum of the subplots under nutrient addition relative to that of control treatment on the log scale. That is, ∆α = log(α_NPK/α_Control) and ∆γ = log(γ_NPK/γ_Control). We calculated ∆β as ∆γ minus \(\overline{\Delta \alpha }\), where\(\,\overline{\Delta \alpha }\) is the average of ∆α over three blocks. A decrease in ∆β indicates nutrient addition causes species composition to be more similar among three subplots than that among control subplots. Because sites are not evenly distributed around the world, many sites are aggregated in North America, we checked spatial autocorrelation of diversity change under nutrient addition using Moran’s I⁴⁴. We found that\(\,\overline{\Delta \alpha }\), ∆γ, and ∆β did not appear to be more similar for sites that are closer to each other (Table S2).

We fitted multilevel (also referred as mixed effects or hierarchical) models for ∆α, ∆γ, and ∆β (as the response variable; all on the log scale) separately. We included random intercept for each site, model was coded as: richness change ~ 1 + (1 |sites) to estimate site-level variation. We used Bayesian analysis because it yields full posterior distributions of parameters rather than point estimates and p-values, which provides a deeper understanding of the uncertainty and variability in the results⁴⁵. Models described above were fitted using the Hamiltonian Monte Carlo (HMC) sampler in Stan and coded using the package ‘brms’ (version 2.21.0) in R (version 4.4.1)^46,47. Models were fitted without explicitly specifying priors, allowing brms to assign its default priors. Models were fitted with 6 chains and 3000 iterations (1000 iterations for warm up). Visual inspection of the HMC chains and Rhat summaries showed model convergence (all Rhats <1.03; Tables S3, S5 and S6). We visually checked posterior predictive plots to determine how well models can reproduce the data (Fig. S2).

To examine whether diversity changes were sensitive to species relative covers, we redid the above analyses (i.e., based on species richness) using Shannon diversity and Simpson diversity (both converted to effective numbers)⁴⁸ (Fig. S4). Species richness is most sensitive to rare species, followed by Shannon diversity, and Simpson diversity is more sensitive to the numbers of relatively abundant species. We calculated the exponential of Shannon diversity and the inverse form of Simpson diversity using the R package vegan (version 2.6-6.1)⁴⁹. These three diversity metrics equal to diversity with order q = {0, 1, 2}, where increasing q decreases the influence of rare species, and D_q = \({\left({\sum }_{i=1}^{s}{p}_{i}^{q}\right)}^{1/\left(1-q\right)}\), where p is the relative cover of species i, s is the total number of species. These diversity metrics are also referred to as Hill numbers^48,50.

Site covariates

We investigated whether the effects of nutrient addition on \(\overline{\Delta \alpha }\), γ, and β diversity based on species richness were mediated by site characteristics. We included site characteristics that have been shown in previous literature to influence ∆α, ∆γ, and ∆β in grasslands: site species pool, site productivity, drought intensity, and grazing intensity^24,25,34,40. We quantified drought intensity as the sum of annual evapotranspiration/precipitation, and averaged it from year 0 to 4 at each site. Precipitation and potential evapotranspiration used were downloaded from https://crudata.uea.ac.uk/cru/data/hrg/cru_ts_4.07/. We quantified the site species pool as the total number of species and site productivity as the average aboveground biomass from year 0 to 4 under the control treatment in the three blocks at each site. Aboveground biomass was harvested within two 1 × 0.1 m strips (in total 0.2 m²), strips were moved from year to year to avoid resampling the same location. For subshrubs and shrubs occurring within strips, we collected all leaves and current year’s woody growth. All biomass was dried at 60 °C (to constant mass) before weighing to the nearest 0.01 g. We used published methods to quantify an integrated grazing intensity metric from vertebrate herbivores at each site. Specifically, herbivore species (>2 kg) that consume grassland biomass were documented at each site by site PIs, and each species was assigned an importance value from 1 (present but low impact and frequency) to 5 (high impact and frequency). An index value was calculated for each site as the sum of herbivore importance values for all herbivores following refs. ^51,52. We also investigated relationships between change in diversity and distance among blocks, because species composition may become less similar as the distance between sampled communities increases. The average pairwise distance among the three blocks within sites ranged from 23.04 to 12538.09 m, with a mean of 513.01 m and a median of 118.7 m across 54 sites that have geolocation data for each block. We first calculated three Euclidean distances between pairs of blocks, we then used the mean of these pairwise distances as the average distance among blocks. We used the average distance among blocks instead of area, because blocks are arranged in parallel at some sites. We fitted linear regression models with\(\,\overline{\Delta \alpha }\), ∆γ, and ∆β as the response variable separately, and each of the site characteristics was used as a predictor variable.

Species groups

We then investigated the effects of nutrient addition on α, γ, and β diversity within groups of species with similar characteristics following the method for changes in α, γ, and β diversity in the entire communities. We eliminated sites where no species occurred in control, nutrient addition, or both plots for a particular group because the value of the log (0) is undefined. We ran the analyses separately for native and non-native species. Native and non-native species were classified by site PIs. Then, we investigated effects of nutrient addition on species richness for different life forms including forb, graminoid, legume, and woody species because previous studies have shown that different life forms may show distinct responses to nutrient addition^6,11,53.

Sensitivity test

We tested whether effects of nutrient addition on species richness across spatial scales depend on experimental duration because a few single-site experiments have shown that the effects of nutrient additions on changes in diversity, especially β diversity, may take several years to emerge^29,31. To that end, we used a subset of 14 sites that had data 14 years after nutrient additions began. Also, because three blocks may be limited in spatial extent, we tested whether combining more blocks to create the γ scale would alter our results. We redid the analyses using data from 11 sites that had five spatial blocks.

Reporting summary

Further information on research design is available in the Nature Portfolio Reporting Summary linked to this article.