# General destabilizing effects of eutrophication on grassland productivity at multiple spatial scales

## Abstract

Eutrophication is a widespread environmental change that usually reduces the stabilizing effect of plant diversity on productivity in local communities. Whether this effect is scale dependent remains to be elucidated. Here, we determine the relationship between plant diversity and temporal stability of productivity for 243 plant communities from 42 grasslands across the globe and quantify the effect of chronic fertilization on these relationships. Unfertilized local communities with more plant species exhibit greater asynchronous dynamics among species in response to natural environmental fluctuations, resulting in greater local stability (alpha stability). Moreover, neighborhood communities that have greater spatial variation in plant species composition within sites (higher beta diversity) have greater spatial asynchrony of productivity among communities, resulting in greater stability at the larger scale (gamma stability). Importantly, fertilization consistently weakens the contribution of plant diversity to both of these stabilizing mechanisms, thus diminishing the positive effect of biodiversity on stability at differing spatial scales. Our findings suggest that preserving grassland functional stability requires conservation of plant diversity within and among ecological communities.

## Introduction

Humans are altering global nutrient cycles via combustion of fossil fuels and fertilizer application1. We have more than doubled preindustrial rates of nitrogen (N) and phosphorus (P) supply to terrestrial ecosystems2. Terrestrial N and P inputs are predicted to reach levels that are three to four times preindustrial rates by 2050 (ref. 3). This pervasive global eutrophication will have dramatic consequences on the structure and functioning of terrestrial and aquatic ecosystems3. In grasslands, nutrient enrichment usually increases primary productivity, but reduces plant diversity, and alters the ability of ecosystems to reliably provide functions and services for humanity4,5,6,7.

Concerns that eutrophication compromises both the diversity and stability of ecosystems have led to a growing number of theoretical and empirical studies investigating how these ecosystem responses may be mechanistically linked4,6,8,9,10,11. These studies have repeatedly shown that the positive effect of plant species richness on the temporal stability of community productivity in ambient (unfertilized) conditions is usually reduced with fertilization4,5,6. However, these studies have primarily focused on plant responses at relatively small scales (i.e., within single local communities). Whether fertilization reduces the positive effect of diversity on temporal stability at larger scales (i.e., among neighboring local communities) remains unclear. Filling this knowledge gap is important because the stable provision of ecosystem services is critical for society12. This is especially true, given an increasing concern for large variability of environmental conditions due to multiple anthropogenic influences, including eutrophication and climate change13.

A recent theoretical framework allows the quantification of the processes that determine the stability of ecosystem functioning at scales beyond the single local community (Fig. 1)14,15,16. Stability at any given scale is defined as the temporal mean of primary productivity divided by its standard deviation17. Higher local scale community stability (alpha stability) can result from two main processes. First, a higher average temporal stability of all species in the community (species stability) can stabilize community productivity due to lower variation in individual species abundances from year to year (Fig. 1b). Second, more asynchronous temporal dynamics among species in response to environmental fluctuations (species asynchrony) can stabilize community productivity because declines in the abundance of some species through time are compensated for by increases in other species (Fig. 1c). Higher stability at the larger scale (gamma stability) can result from higher alpha stability and more asynchronous dynamics across local communities (spatial asynchrony; Fig. 1d). Thus, the stabilizing effect of spatial asynchrony on productivity at the larger scale (spatial insurance hypothesis)14,18 mirrors the stabilizing effect of species asynchrony on productivity at the local scale (species or local insurance hypothesis)8,16,19,20. Higher species asynchrony and species stability can result from higher local species diversity through higher species richness9,21,22, higher species evenness8, or both (e.g., higher values of diversity indices—such as the Shannon index—that combines the two23; Fig. 1e). Higher spatial asynchrony can result from greater local species diversity or higher variation in species composition among communities (beta diversity)16.

According to this framework, fertilization can affect the links between diversity, asynchrony, and stability across spatial scales (Fig. 1e and Table 1). At the local scale, fertilization can decrease niche dimensionality, and favor a few dominant plant species by affecting the competitive balance among species, potentially reducing the insurance effects of local diversity7,22. At the larger scale, fertilization can reduce spatial heterogeneity in community composition, and decrease variations among local plant community structure, potentially reducing the spatial insurance effect of beta diversity16. Moreover, fertilization often reduces plant diversity, which could in turn reduce asynchrony and stability at multiple scales4,9,17,24. However, the role of fertilization in mediating the functional consequences of biodiversity changes (variations in the number, abundance, and identities of species) and compensatory mechanisms (variation and compensation in species responses) that can affect the stable provisioning of ecosystem functions at larger spatial scales remains to be elucidated25.

To our knowledge, only one recent study has assessed the effect of nutrient enrichment on stability within and among interconnected communities in a temperate grassland26. By adding different nitrogen treatments to communities in ten blocks spread out within a single site, that study found that 5 years of chronic nitrogen addition reduced alpha stability through a decline in species asynchrony, but had no effect on spatial asynchrony. However, these conclusions were based on a single grassland site manipulating a single nutrient, with the implicit assumption that the relationship between diversity and stability was unaffected by eutrophication. This argues for multisite comparative studies assessing the generality of the mechanistic links between these ecosystem responses to eutrophication.

Here, we use a coordinated, multisite and multiyear nutrient enrichment experiment (±chronic nitrogen, phosphorus, and potassium addition, Nutrient Network (NutNet)27) to assess the scale dependence of fertilization impacts on plant diversity and stability. Treatments were randomly assigned to 25 m2 plots and were replicated in three blocks at most sites (Supplementary Data 1). Samples were collected in 1 m2 subplots across 243 communities from 42 grassland sites on six continents and followed a standardized protocol at all sites27. We selected these sites as they contained between 4 and 9 years of experimental duration (hereafter “period of experimental duration”), and three blocks per site, excluding additional blocks from sites that had more than three (Supplementary Data 1). Sites spanned a broad range of seasonal variation in precipitation and temperature (Supplementary Fig. 1), and a wide range of grassland types (Supplementary Data 1). In our analysis, we treated each 1 m2 subplot as a “community” and the replicated subplots within a site as the “larger scale” sensu Whittaker28. We computed diversity, asynchrony, and stability within a community (local “alpha” scale) and across the three replicated communities within a site (larger “gamma” scale) (see “Methods”). We then used bivariate analysis and structural equation modeling (SEM)29 to assess fertilizer impacts, and disentangle the relative contributions of diversity and asynchrony to stability (Fig. 1e).

## Results and discussion

### Fertilization effects on diversity, asynchrony, and stability

Analyses of variance revealed the negative effects of nutrient inputs on biodiversity and stability at the two scales investigated, consistent with recent findings from a single site26. Fertilization consistently reduced species richness, alpha, and gamma stability, but had no effect on beta diversity (Supplementary Fig. 2). Bivariate analyses further revealed the negative effects of nutrient inputs on biodiversity–stability relationships at the two scales investigated (Fig. 2). Relationships were generally consistent across the different periods of experimental duration considered (Supplementary Table 1). Under ambient (unfertilized) conditions, species richness was positively associated with alpha and gamma stability (Fig. 2a, b), but fertilization weakened the positive effect of species richness on stability at the two scales (Fig. 2c, d). Fertilization reduced local stability of grassland functioning by increasing temporal variability in species-rich communities (Supplementary Fig. 3). Similarly, high beta diversity (variation in species composition among communities) was positively associated with spatial asynchrony and gamma stability under ambient conditions (Fig. 2e, f), but again fertilization weakened the positive effect of beta diversity on spatial asynchrony and gamma stability (Fig. 2g, h). These results remained when accounting for variation in climate using residual regression (Supplementary Fig. 4), when using local diversity indices accounting for species abundance (Supplementary Fig. 5), and when data were divided into overlapping intervals of 4 years (Supplementary Fig. 6). Our results extend previous evidence of the negative impact of fertilization on the diversity–stability relationship obtained within local plots and over shorter experimental periods4,6,26. Importantly, they show that these negative effects propagate from within to among communities. To our knowledge, our study is the first to report the negative impacts of fertilization on the relationships of beta diversity, with spatial asynchrony and gamma stability.

### Mechanisms linking diversity and stability

To understand the relative role of local vs. larger scale community properties in determining asynchrony and stability at different spatial scales, we conducted SEM analyses, including all measures in a single causal model (Fig. 3, Supplementary Fig. 7 and Supplementary Table 2). Under ambient conditions, SEM revealed that higher plant species richness contributed to greater alpha and gamma stability largely through higher asynchronous dynamics among species (species asynchrony, standardized path coefficient = 0.39), and not necessarily through greater species stability (standardized path coefficient = 0.01; Fig. 3a and Supplementary Fig. 8a, b). The positive association between species richness and alpha stability is consistent with existing experimental17,24 and shorter-term observational evidence4,30,31. Our results confirm that the stabilizing effects of species richness in naturally assembled grassland communities is largely driven by species asynchrony, but not species stability4,6,22,26. In addition, they show that the positive impact of species richness on the stability of community productivity via species asynchrony in turn leads to greater stability of productivity at the larger spatial scale.

While correlated with species richness, higher beta diversity also contributed to greater gamma stability through an independent pathway, namely via higher asynchronous dynamics among local communities (spatial asynchrony, standardized path coefficient = 0.20, Fig. 3a). While theoretical studies have suggested a role for beta diversity in driving spatial asynchrony15,16, previous empirical studies conducted along a nitrogen gradient at a single site26 or across 62 sites with non-standardized protocols21 did not find an association between these two variables. Here, we show that the presence of different species among local communities is linked to higher variation in dynamics among them, demonstrating the stabilizing role of beta diversity at larger spatial scales through spatial asynchrony. This also indicates the need for multisite replication with standardized treatments and protocols to detect such effects.

Importantly, fertilization acted to destabilize productivity at the local and larger spatial scale through several mechanisms (Fig. 3 and Table 2). At the local scale, fertilization weakened the positive effects of plant species richness on alpha and gamma stability (Fig. 2a–d) via a combination of two processes (Fig. 3b and Supplementary Fig. 8c, d). First, the positive relationship between species richness and species asynchrony in the control communities (standardized path coefficient = 0.39, Fig. 3a), was weaker in the fertilized communities (standardized path coefficient = 0.20, Fig. 3b). Moreover, this general positive effect of richness on asynchrony was counteracted by a second stronger negative relationship of richness with species stability (standardized path coefficient = −0.37). Such negative effect of fertilization on species stability was not observed under ambient conditions, and could be due to shifts in functional composition in species-rich communities from more stable conservative species to less stable exploitative species in a temporally variable environment32,33. Together, these two effects explain the overall weaker alpha stability at higher richness with fertilization. We did not find evidence that the loss of diversity caused by fertilization (an average of −1.8 ± 0.5 species m−2, Supplementary Fig. 2a and Supplementary Fig. 9a) was related to the decline of alpha stability, confirming results from other studies5,6 and earlier NutNet results4 obtained over shorter time periods. This could be because the negative feedback of the loss of richness caused by fertilization on stability requires a longer experimental duration, or greater loss of plant diversity, to manifest9,34. Another possible explanation is that fertilization may have a direct positive effect on stability, by increasing community biomass (t = 2.41, d.f. = 326, P = 0.016) and enhancing stability via overyielding effects35, a formal test that would require monocultures.

At the larger scale, fertilization reduced the strength of the relationship between beta diversity and gamma stability by reducing the strength of the relationship between beta diversity and spatial asynchrony (standardized path coefficient = 0.20 in Fig. 3a vs. standardized path coefficient = 0.03 in Fig. 3b). This result provides evidence that fertilization can reduce the stabilizing role of spatial asynchrony among initially dissimilar communities. We did not find evidence that this was due to a negative feedback of changes in beta diversity caused by fertilization on gamma stability (Supplementary Fig. 2b and Supplementary Fig. 9b). The positive relationship between beta diversity and spatial asynchrony, and the negative impact of fertilization on that relationship, suggests that the spatial insurance effect caused by variation in species composition among local communities may be disrupted in a eutrophic world.

### Implications

Our results support the idea that asynchronous dynamics among species in species-rich communities play a stabilizing role, and show that this effect propagates to larger spatial scales21,26. Furthermore, to our knowledge, our study is the first to report the positive association between beta diversity and gamma stability through spatial asynchrony in real-world grasslands. Importantly, fertilization reduced the contribution of biodiversity to these stabilizing mechanisms at both scales, diminishing the local and spatial insurance of biodiversity on stability. Such diminished insurance effects lead to a reduced ecosystem stability at larger scales. Future climate will be characterized by more variability, including more frequent extreme events13. Our results indicate that preserving ecosystem stability across spatial scales in a changing world requires conserving biodiversity within and among local communities. Moreover, policies and management procedures that prevent and mitigate eutrophication are needed to safeguard the positive effects of biodiversity on stability at multiple scales.

## Methods

### 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 _k\mu _k}}{{\sum _k\sqrt {v_{kk}} }},$$
(2)
$${\mathrm{Gamma}}\,{\mathrm{stability}} = \frac{{\sum _k\mu _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 _k\sqrt {\sum _{ij,kl}w_{ij,kl}} }},$$
(4)
$${\mathrm{Spatial}}\,{\mathrm{asynchrony}} = \frac{{\sum _k\sqrt {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.

## Data availability

The data that support the findings of this study are available via GitHub (https://github.com/YannHautier/NutNetStabilityScaleUp). Data sources are provided with this paper. WorldClim global climate database is freely available through the World Data Center for Climate (WDCC; cera-www.dkrz.de), as well as through the CCAFS-Climate data portal (http://ccafs-climate.org).

## Code availability

R code of all analyses are available via https://github.com/YannHautier/NutNetStabilityScaleUp.

## References

1. 1.

Erisman, J. W. et al. Consequences of human modification of the global nitrogen cycle. Philos. Trans. R. Soc. B Biol. Sci. 368, 20130116 (2013).

2. 2.

Galloway, J. N. The global nitrogen cycle: past, present and future. Sci. China Ser. C. Life Sci. 48, 669–677 (2005).

3. 3.

Tilman, D. et al. Forecasting agriculturally driven global environmental change. Science 292, 281–284 (2001).

4. 4.

Hautier, Y. et al. Eutrophication weakens stabilizing effects of diversity in natural grasslands. Nature 508, 521–525 (2014).

5. 5.

Xu, Z. W. et al. Environmental changes drive the temporal stability of semi-arid natural grasslands through altering species asynchrony. J. Ecol. 103, 1308–1316 (2015).

6. 6.

Zhang, Y. H. et al. Nitrogen enrichment weakens ecosystem stability through decreased species asynchrony and population stability in a temperate grassland. Glob. Change Biol. 22, 1445–1455 (2016).

7. 7.

Harpole, W. S. et al. Addition of multiple limiting resources reduces grassland diversity. Nature 537, 93–96 (2016).

8. 8.

Thibaut, L. M. & Connolly, S. R. Understanding diversity-stability relationships: towards a unified model of portfolio effects. Ecol. Lett. 16, 140–150 (2013).

9. 9.

Hautier, Y. et al. Anthropogenic environmental changes affect ecosystem stability via biodiversity. Science 348, 336–340 (2015).

10. 10.

Koerner, S. E. et al. Nutrient additions cause divergence of tallgrass prairie plant communities resulting in loss of ecosystem stability. J. Ecol. 104, 1478–1487 (2016).

11. 11.

Yang, H. J. et al. Diversity-dependent stability under mowing and nutrient addition: evidence from a 7-year grassland experiment. Ecol. Lett. 15, 619–626 (2012).

12. 12.

Millennium Ecosystem Assessment. Ecosystems and Human Well-being: Synthesis (Island Press, 2005).

13. 13.

Shukla, P. R. et al. in Climate Change and Land: an IPCC Special Report on Climate Change, Desertification, Land Degradation, Sustainable Land Management, Food Security, and Greenhouse Gas Fluxes in Terrestrial Ecosystems (IPCC, 2019).

14. 14.

Wang, S. P., Lamy, T., Hallett, L. M. & Loreau, M. Stability and synchrony across ecological hierarchies in heterogeneous metacommunities: linking theory to data. Ecography 42, 1200–1211 (2019).

15. 15.

Wang, S. P. & Loreau, M. Ecosystem stability in space: alpha, beta and gamma variability. Ecol. Lett. 17, 891–901 (2014).

16. 16.

Wang, S. P. & Loreau, M. Biodiversity and ecosystem stability across scales in metacommunities. Ecol. Lett. 19, 510–518 (2016).

17. 17.

Tilman, D., Reich, P. B. & Knops, J. M. H. Biodiversity and ecosystem stability in a decade-long grassland experiment. Nature 441, 629–632 (2006).

18. 18.

Loreau, M., Mouquet, N. & Gonzalez, A. Biodiversity as spatial insurance in heterogeneous landscapes. Proc. Natl Acad. Sci. USA 100, 12765–12770 (2003).

19. 19.

Lamy, T. et al. Species insurance trumps spatial insurance in stabilizing biomass of a marine macroalgal metacommunity. Ecology 100, e02719 (2019).

20. 20.

Loreau, M. & de Mazancourt, C. Biodiversity and ecosystem stability: a synthesis of underlying mechanisms. Ecol. Lett. 16, 106–115 (2013).

21. 21.

Wilcox, K. R. et al. Asynchrony among local communities stabilises ecosystem function of metacommunities. Ecol. Lett. 20, 1534–1545 (2017).

22. 22.

Gilbert, B. et al. Climate and local environment structure asynchrony and the stability of primary production in grasslands. Glob. Ecol. Biogeogr. 29, 1177–1188 (2020).

23. 23.

Zhang, Y., Loreau, M., He, N., Zhang, G. & Han, X. Mowing exacerbates the loss of ecosystem stability under nitrogen enrichment in a temperate grassland. Funct. Ecol. 31, 1637–1646 (2017).

24. 24.

Hector, A. et al. General stabilizing effects of plant diversity on grassland productivity through population asynchrony and overyielding. Ecology 91, 2213–2220 (2010).

25. 25.

Mori, A. S., Isbell, F. & Seidl, R. beta-Diversity, community assembly, and ecosystem functioning. Trends Ecol. Evolution 33, 549–564 (2018).

26. 26.

Zhang, Y. H. et al. Nitrogen addition does not reduce the role of spatial asynchrony in stabilising grassland communities. Ecol. Lett. 22, 563–571 (2019).

27. 27.

Borer, E. T. et al. Finding generality in ecology: a model for globally distributed experiments. Methods Ecol. Evolution 5, 63–73 (2013).

28. 28.

Whittaker, R. H. Evolution and measurement of species diversity. Taxon 21, 213–225 (1972).

29. 29.

Grace, J. B. et al. Guidelines for a graph-theoretic implementation of structural equation modeling. Ecosphere 3, art73 (2012).

30. 30.

Bai, Y., Han, X., Wu, J., Chen, Z. & Li, L. Ecosystem stability and compenatory effects in the inner Mongolia grassland. Nature 431, 181–184 (2004).

31. 31.

Tilman, D. Biodiversity: population versus ecosystem stability. Ecology 77, 350–353 (1996).

32. 32.

Polley, H. W., Isbell, F. I. & Wilsey, B. J. Plant functional traits improve diversity-based predictions of temporal stability of grassland productivity. Oikos 122, 1275–1282 (2013).

33. 33.

Majekova, M., de Bello, F., Dolezal, J. & Leps, J. Plant functional traits as determinants of population stability. Ecology 95, 2369–2374 (2014).

34. 34.

Isbell, F. et al. Nutrient enrichment, biodiversity loss, and consequent declines in ecosystem productivity. Proc. Natl Acad. Sci. USA 110, 11911–11916 (2013).

35. 35.

de Mazancourt, C. et al. Predicting ecosystem stability from community composition and biodiversity. Ecol. Lett. 16, 617–625 (2013).

36. 36.

Oesterheld, M. & McNaughton, S. J. Herbivory in terrestrial ecosystems, in Methods in ecosystem science (eds Sala, O. E., Jackson, R. B., Mooney, H. A. & Howarth, R. W.) 151–157 (Springer, New York, 2000).

37. 37.

Tuomisto, H. An updated consumer’s guide to evenness and related indices. Oikos 121, 1203–1218 (2012).

38. 38.

Jost, L. et al. Partitioning diversity for conservation analyses. Diversity Distrib. 16, 65–76 (2010).

39. 39.

Shannon, C. E. A mathematical theory of communication. Bell Syst. Tech. J. 27, 379–423 (1948).

40. 40.

Pielou, E. C. Measurement of Diversity in different types of biological collections. J. Theor. Biol. 13, 131-& (1966).

41. 41.

Simpson, E. H. Measurement of diversity. Nature 163, 688–688 (1949).

42. 42.

Olszewski, T. D. A unified mathematical framework for the measurement of richness and evenness within and among multiple communities. Oikos 104, 377–387 (2004).

43. 43.

Dixon, P. VEGAN, a package of R functions for community ecology. J. Veg. Sci. 14, 927–930 (2003).

44. 44.

Hijmans, R. J., Cameron, S. E., Parra, J. L., Jones, P. G. & Jarvis, A. Very high resolution interpolated climate surfaces for global land areas. Int. J. Climatol. 25, 1965–1978 (2005).

45. 45.

R Core Team. R: A language and environment for statistical computing. R Foundation for Statistical Computing, Vienna, Austria http://www.R-project.org/ (2013).

46. 46.

Pinheiro, J. C. & Bates, D. M. in Mixed-Effects Models in S and S-Plus (Spinger, New York, 2000).

47. 47.

Trikalinos, T. A. & Olkin, I. Meta-analysis of effect sizes reported at multiple time points: a multivariate approach. Clin. Trials 9, 610–620 (2012).

48. 48.

Hedges, L. V. & Olkin, I. Statistical Methods For Meta-analysis (Academic, 1985).

49. 49.

Lefcheck, J. S. PIECEWISESEM: piecewise structural equation modelling in R for ecology, evolution, and systematics. Methods Ecol. Evol. 7, 573–579 (2016).

50. 50.

Viechtbauer, W. Conducting meta-analyses in R with the metafor package. J. Stat. Softw. 36, 1–48 (2010).

51. 51.

Loreau, M. & de Mazancourt, C. Species synchrony and its drivers: neutral and nonneutral community dynamics in fluctuating environments. Am. Nat. 172, E48–E66 (2008).

## Acknowledgements

The research leading to these results has received funding from the European Union. Seventh Framework Programme (FP7/2007-2013) under grant agreement no. 298935 to Y.H. (with A.H. and E.W.S.). This work was generated using data from the Nutrient Network collaborative experiment, funded at the site scale by individual researchers and coordinated through Research Coordination Network funding from NSF to E.B. and E.W.S. (grant #DEB-0741952). Nitrogen fertilizer was donated to the Nutrient Network by Crop Production Services, Loveland, CO. We acknowledge support from the LTER Network Communications Office and DEB-1545288. M.L. was supported by the TULIP Laboratory of Excellence (ANR-10-LABX-41), and by the BIOSTASES Advanced Grant funded by the European Research Council under the European Union’s Horizon 2020 research and innovation programme (grant agreement no. 666971). S.W. was supported by the National Natural Science Foundation of China (31988102). We thank Rita S. L. Veiga and George A. Kowalchuk for suggestions that improved the manuscript.

## Author information

Authors

### Contributions

Y.H., P.Z., K.R.W., M.L., and S.W. developed and framed research questions. Y.H. and S.W. analyzed the data with help from P.Z., K.R.W., E.W.S., J.E.K.B., S.E.K., K.J.K., and J.S.L. Y.H. wrote the paper with contributions and input from all authors. E.W.S. and E.T.B. are Nutrient Network coordinators. The author contribution matrix is provided as Supplementary Data 2.

### Corresponding authors

Correspondence to Yann Hautier or Shaopeng Wang.

## Ethics declarations

### Competing interests

The authors declare no competing interests.

Peer review information Nature Communications thanks the anonymous reviewer(s) for their contribution to the peer review of this work. Peer reviewer reports are available.

Publisher’s note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

## Rights and permissions

Reprints and Permissions

Hautier, Y., Zhang, P., Loreau, M. et al. General destabilizing effects of eutrophication on grassland productivity at multiple spatial scales. Nat Commun 11, 5375 (2020). https://doi.org/10.1038/s41467-020-19252-4

• Accepted:

• Published: