General destabilizing effects of eutrophication on grassland productivity at multiple spatial scales
Study sites and experimental design
The study sites are part of the NutNet experiment (Supplementary Data 1; http://nutnet.org/)27. Plots at each site are 5 × 5 m separated by at least 1 m. All sites included in the analyses presented here included unmanipulated plots and fertilized plots with nitrogen (N), phosphorus (P), and potassium and micronutrients (K) added in combination (NPK+). N, P, and K were applied annually before the beginning of the growing season at rates of 10 gm−2 y−1. N was supplied as time-release urea ((NH2)2CO) or ammonium nitrate (NH4NO3). P was supplied as triple super phosphate (Ca(H2PO4)2), and K as potassium sulfate (K2SO4). In addition, a micronutrient mix (Fe, S, Mg, Mn, Cu, Zn, B, and Mo) was applied at 100 gm−2 y−1 to the K-addition plots, once at the start of the experiment, but not in subsequent years to avoid toxicity. Treatments were randomly assigned to the 25 m2 plots and were replicated in three blocks at most sites (some sites had fewer/more blocks or were fully randomized). Sampling was done in 1 m2 subplots and followed a standardized protocol at all sites27.
Site selection
Data were retrieved on 1 May 2020. To keep a constant number of communities per site and treatment, we used three blocks per site, excluding additional blocks from sites that had more than three (Supplementary Data 1). Sites spanned a broad envelope of seasonal variation in precipitation and temperature (Supplementary Fig. 1), and represent a wide range of grassland types, including alpine, desert and semiarid grasslands, prairies, old fields, pastures, savanna, tundra, and shrub-steppe (Supplementary Data 1).
Stability and asynchrony measurements are sensitive to taxonomic inconsistencies. We adjusted the taxonomy to ensure consistent naming over time within sites. This was usually done by aggregating taxa at the genus level when individuals were not identified to species in all years. Taxa are however referred to as “species”.
We selected sites that had a minimum of 4 years, and up to 9 years of posttreatment data. Treatment application started at most sites in 2008, but some sites started later resulting in a lower number of sites with increasing duration of the study, from 42 sites with 4 years of posttreatment duration to 15 sites with 9 years of duration (Supplementary Data 1). Longer time series currently exist, but for a limited number of sites within our selection criteria.
Primary productivity and cover
We used aboveground live biomass as a measure of primary productivity, which is an effective estimator of aboveground net primary production in herbaceous vegetation36. Primary productivity was estimated annually by clipping at ground level all aboveground live biomass from two 0.1 m2 (10 × 100 cm) quadrats per subplot. For shrubs and subshrubs, leaves and current year’s woody growth were collected. Biomass was dried to constant mass at 60 °C and weighed to the nearest 0.01 g. Areal percent cover of each species was measured concurrently with primary productivity in one 1 × 1 m subplot, in which no destructive sampling occurred. Cover was visually estimated annually to the nearest percent independently for each species, so that total summed cover can exceed 100% for multilayer canopies. Cover and primary productivity were estimated twice during the year at some sites with strongly seasonal communities. This allowed to assemble a complete list of species and to follow management procedures typical of those sites. For those sites, the maximum cover of each species and total biomass were used in the analyses.
Diversity, asynchrony, and stability across spatial scales
We quantified local scale and larger-scale diversity indices across the three replicated 1-m2 subplots for each site, treatment and duration period using cover data37,38. In our analysis, we treated each subplot as a “community” and the collective subplots as the “larger scale” sensu Whittaker28. Local scale diversity indices (species richness, species evenness, Shannon, and Simpson) were measured for each community, and averaged across the three communities for each treatment at each site resulting in one single value per treatment and site. Species richness is the average number of plant species. Shannon is the average of Shannon–Weaver indices39. Species evenness is the average of the ratio of the Shannon–Weaver index and the natural logarithm of average species richness (i.e., Pielou’s evenness40). Simpson is the average of inverse Simpson indices41. Due to strong correlation between species richness and other common local diversity indices (Shannon: r = 0.90 (95% confidence intervals (CIs) = 0.87–0.92), Simpson: r = 0.88 (0.86–0.91), Pielou’s evenness: r = 0.62 (0.55–0.68), with d.f. = 324 for each), we used species richness as a single, general proxy for those variables in our models. Results using these diversity indices did not differ quantitatively from those presented in the main text using species richness (Supplementary Fig. 5), suggesting that fertilization modulate diversity effects largely through species richness. Following theoretical models15,16, we quantified abundance-based gamma diversity as the inverse Simpson index over the three subplots for each treatment at each site and abundance-based beta diversity, as the multiplicative partitioning of abundance-based gamma diversity: abundance-based beta equals the abundance-based gamma over Simpson28,42, resulting in one single beta diversity value per treatment and site. We used abundance-based beta diversity index because it is directly linked to ecosystem stability in theoretical models15,16, and thus directly comparable to theories. We used the R functions “diversity”, “specnumber”, and “vegdist” from the vegan package43 to calculate Shannon–Weaver, Simpson, and species richness indices within and across replicated plots.
Stability at multiple scales was determined both without detrending and after detrending data. For each species within communities, we detrended by using species-level linear models of percent cover over years. We used the residuals from each regression as detrended standard deviations to calculate detrended stability17. Results using detrended stability did not differ quantitatively from those presented in the main text without detrending. Stability was defined by the temporal invariability of biomass (for alpha and gamma stability) or cover (for species stability and species asynchrony), calculated as the ratio of temporal mean to standard deviation14,17. Gamma stability represents the temporal invariability of the total biomass of three plots with the same treatment, alpha stability represents the temporal invariability of community biomass averaged across three plots per treatment and per site, and species stability represents the temporal invariability of species cover averaged across all species and the three plots per treatment14. The mathematical formula are:
$${mathrm{Species}},{mathrm{stability}} = frac{{sum _{i,k}m_{i,k}}}{{sum _{i,k}sqrt {w_{ii,kk}} }},$$
(1)
$${mathrm{Alpha}},{mathrm{stability}} = frac{{sum _kmu _k}}{{sum _ksqrt {v_{kk}} }},$$
(2)
$${mathrm{Gamma}},{mathrm{stability}} = frac{{sum _kmu _k}}{{sqrt {sum _{k,l}nu _{kl}} }},$$
(3)
where mi,k and wii,kk denote the temporal mean and variance of the cover of species i in subplot k; μk and vkk denote the temporal mean and variance of community biomass in subplot k, and vkl denotes the covariance in community biomass between subplot k and l. We then define species asynchrony as the variance-weighted correlation across species, and spatial asynchrony as the variance-weighted correlation across plots:
$${mathrm{Species}},{mathrm{asynchrony}} = frac{{sum _{i,k}sqrt {w_{ii,kk}} }}{{sum _ksqrt {sum _{ij,kl}w_{ij,kl}} }},$$
(4)
$${mathrm{Spatial}},{mathrm{asynchrony}} = frac{{sum _ksqrt {v_{kk}} }}{{sqrt {sum _{k,l}nu _{kl}} }},$$
(5)
where wij,kl denotes the covariance in species cover between species i in subplot k and species j in subplot l.
These two asynchrony indices quantify the incoherence in the temporal dynamics of species cover and community biomass, respectively, which serve as scaling factors to link stability metrics across scales14 (Fig. 1). To improve normality, stability, and asynchrony measures were logarithm transformed before analyses. We used the R function “var.partition” to calculate asynchrony and stability across spatial scales14.
Climate data
Precipitation and temperature seasonality were estimated for each site, using the long-term coefficient of variation of precipitation (MAP_VAR) and temperature (MAT_VAR), respectively, derived from the WorldClim Global Climate database (version 1.4; http://www.worldclim.org/)44.
Analyses
All analyses were conducted in R 4.0.2 (ref. 45) with N = 42 for each analysis unless specified. First, we used analysis of variance to determine the effect of fertilization, and period of experimental duration on biodiversity and stability at the two scales investigated. Models including an autocorrelation structure with a first-order autoregressive model (AR(1)), where observations are expected to be correlated from 1 year to the next, gave substantial improvement in model fit when compared with models lacking autocorrelation structure. Second, we used bivariate analyses and linear models to test the effect of fertilization and period of experimental duration on biodiversity–stability relationships at the two scales investigated. Again, models including an autocorrelation structure gave substantial improvement in model fit (Supplementary Table 1)46,47,48. We ran similar models based on nutrient-induced changes in diversity, stability, and asynchrony. For each site, relative changes in biodiversity, stability, and asynchrony at the two scales considered were calculated, as the natural logarithm of the ratio between the variable in the fertilized and unmanipulated plots (Supplementary Fig. 9). Because plant diversity, asynchronous dynamics, and temporal stability may be jointly controlled by interannual climate variability22, we ran similar analyses on the residuals of models that included the coefficient of variation among years for each of temperature and precipitation. Results of our analyses controlling for interannual climate variability did not differ qualitatively from the results presented in the text (Supplementary Fig. 4). In addition, to test for temporal trends in stability and diversity responses to fertilization, we used data on overlapping intervals of four consecutive years. Results of our analyses using temporal trends did not differ qualitatively from the results presented in the text (Supplementary Fig. 6). Inference was based on 95% CIs.
Second, we used SEM29 with linear models, to evaluate multiple hypothesis related to key predictions from theories (Table 1). The path model shown in Fig. 1e was evaluated for each treatment (control and fertilized), and we ran separate SEMs for each period of experimental duration (from 4 to 9 years of duration). We generated a summary SEM by performing a meta-analysis of the standardized coefficients across all durations for each treatment. We then tested whether the path coefficients for each model differed by treatment by testing for a model-wide interaction with the “treatment” factor. A positive interaction for a given path implied that effects of one variable on the other are significantly different between fertilized and unfertilized treatments. We used the R functions “psem” to fit separate piecewise SEMs49 for each duration and combined the path coefficients from those models, using the “metagen” function50.
Reporting summary
Further information on research design is available in the Nature Research Reporting Summary linked to this article. More
