Multidimensional responses of grassland stability to eutrophication

Global eutrophication remains widespread, with most nations failing to meet targets to limit excess nutrients, threatening biodiversity, ecosystem functioning, and nature’s contributions to people. Concurrently, climate extremes such as droughts and floods are increasing in frequency and intensity. How eutrophication affects ecosystem stability under increasing climatic variability is unclear. Ecological stability is multifaceted, encompassing temporal invariability (mean-to-variance through time), resistance to perturbations during extremes, and recovery after extremes. Prior work often assessed single facets (commonly invariability or drought resistance) and mainly focused on biomass, despite stability also being definable for composition and species richness. Correlations among stability facets and among community aspects can indicate low (strongly correlated) or high (weakly correlated) dimensionality. Understanding whether eutrophication alters the magnitudes of distinct stability facets and the correlations among them is essential to predict and manage grassland responses to global change. This study tests how nutrient addition (NPK) influences five stability facets across three community aspects and whether eutrophication changes the dimensionality (correlation structure) of stability.

Previous studies show that different stability facets can be uncorrelated and that their relationships may shift under global change, including climate extremes, eutrophication, consumer removal, light addition, and heatwaves. While biomass stability has been the primary focus, compositional and richness stability can strongly influence ecosystem functions. High compositional invariability can elevate biomass invariability, but compensatory dynamics may allow high biomass stability despite compositional variability. Greater species richness can enhance biomass invariability via asynchronous species responses, though eutrophication and aridity can modulate this relationship. Evidence on eutrophication effects on biomass invariability, resistance, and recovery is mixed, and effects during and after wet extremes are rarely examined. Eutrophication typically reduces species richness and shifts composition toward fast-growing or invasive species, but impacts on the stability of richness and composition remain poorly resolved. Theory suggests eutrophication could either strengthen or weaken correlations among stability facets by altering the balance between deterministic and stochastic assembly processes. A recent grassland study found eutrophication did not alter relationships between biomass and compositional stability, highlighting the need for a multi-site, multi-facet assessment.

Study system and design: The work leveraged 55 grassland sites from the Nutrient Network (NutNet) spanning five continents. Plots (5 × 5 m) in randomized blocks (typically three per site) received either Control or NPK fertilization. Fertilized plots received N, P, and K annually at 10 g m−2 yr−1 (N as time-release urea, P as triple superphosphate, K as potassium sulfate). A micronutrient mix (Fe, S, Mg, Mn, Cu, Zn, B, Mo) was applied once at 100 g m−2 at the start of the experiment. Data were retrieved November 2022. Inclusion criteria for sites: three blocks; at least four years of post-treatment measurements; and at least one dry or wet growing season during experimental years. Sampling: Within each plot, a permanently marked 1 × 1 m subplot was used for annual plant cover surveys (visual percent cover; multilayer canopies could exceed 100% cover). Aboveground biomass was harvested annually adjacent to the subplot from two 1 × 0.1 m strips (0.2 m2 total), dried at 60 °C to constant mass, and weighed (g m−2). For shrubs/subshrubs, leaves and current-year woody growth were collected. At seasonal sites, cover was recorded twice per year; maximum cover values and total biomass were used. Taxonomy was standardized within sites (aggregating to genus where necessary, treated as “species” for analyses). Climate extremes: The standardized precipitation-evapotranspiration index (SPEI) was computed for each site as the growing-season standardized water balance (precipitation minus potential evapotranspiration) using CRU TS 4.06 data (1901–2021). PET followed the Penman–Monteith method. Growing seasons were classified as dry (≤25th percentile), normal (25–75th), or wet (≥75th) using ±0.67 SD thresholds; robustness checks used ±1.28 SD (dry/wet ≈10th/90th percentiles). To minimize confounding from consecutive extremes, rules excluded certain years when different extremes occurred consecutively or when multiple same-type extremes occurred; recovery was calculated only when followed by a normal season (details in Supplementary Table 2 and Figs. 4–5). Stability facets and community aspects: Five stability facets were quantified for three community aspects (biomass, composition, richness): temporal invariability, resistance during dry and wet growing seasons, and recovery after dry and wet growing seasons. For biomass and richness, temporal invariability was mean/SD over experimental years; SD used residuals after detrending by linear regression on year to remove directional trends. Following Isbell et al. (2015), resistance was defined as the inverse proportional deviation relative to normal levels during an extreme, operationalized as (Ye − Yn)/Yn; recovery as (Ye+1 − Ye)/Ye (with higher absolute values transformed so that higher stability corresponds to higher metric; biomass and richness metrics were log-transformed for variance homogeneity). For composition, Bray–Curtis dissimilarity based on cover was converted to similarity (1 − dissimilarity) so higher values indicate higher stability. Temporal invariability in composition was the multi-year similarity (betapart::beta.multi.abund). Resistance for composition was the similarity between the community during an extreme and the average community of normal years (reference) within treatment and block (betapart::beta.pair.abund). Recovery for composition was the ratio of post-extreme similarity to during-extreme similarity. Resistance and recovery were averaged over years per site-treatment to match invariability scale. For all facets, higher values represent greater stability. Statistical analyses: Linear mixed-effects models (nlme::lme) tested nutrient addition effects on each stability facet within each community aspect, with treatment as fixed effect and site and block nested within site as random effects. To evaluate sensitivity of rare vs. common species in richness metrics, effective species diversity (Hill numbers Q = 0, 1, 2) weighted by cover was calculated. Dimensionality analyses computed within-site Pearson correlations among the five facets (10 pairs) within each community aspect and tested treatment effects on correlation coefficients using mixed models with site as random effect. Similarly, correlations among stability across the three community aspects were calculated for each facet and tested for treatment effects. Robustness checks used stricter SPEI thresholds (±1.28 SD), detrended SPEI time series, and a subset of 22 sites with ≥10 years of data.

• Scope: 55 grassland sites across five continents with ≥4 years of standardized NPK addition; across sites, 150 dry, 247 normal, and 131 wet growing seasons were identified (0.67 SD SPEI thresholds).
• Effects on stability facets:
  • Species richness and community composition: Nutrient addition significantly reduced temporal invariability and resistance during both dry and wet growing seasons. For composition, NPK altered normal-season abundance distributions and increased compositional deviations during extremes, lowering resistance. For richness, NPK reduced resistance during dry extremes via lower normal levels and greater deviations, and during wet extremes primarily via reduced normal levels.
  • Aboveground biomass: No significant NPK effect on any stability facet (invariability, resistance during dry/wet, recovery after dry/wet). Although NPK increased biomass in normal seasons and increased absolute deviations during extremes, proportional deviations were similar between treatments, yielding no net effect on resistance or invariability.
  • Recovery: Nutrient addition did not significantly affect recovery after dry or wet extremes for biomass, composition, or richness. However, absolute deviations during and one year after extremes were larger under NPK, especially for composition and biomass.
• Dimensionality of stability (correlation structure):
  • Within community aspects: Under control conditions, only a minority of facet pairs were significantly correlated (biomass: 4/10; composition: 4/10; richness: 2/10), indicating high dimensionality. This pattern persisted under NPK, with some pairwise correlations weakened or strengthened. Temporal invariability was positively correlated with resistance but not with recovery across all community aspects and treatments. Resistance and recovery tended to be negatively correlated (trade-off) for biomass and composition under both treatments; this trade-off was not evident for richness.
  • Among community aspects: For any given facet, correlations of stability among biomass, richness, and composition were generally weak under both treatments, indicating that stability of different community aspects represents separate dimensions.
• Dominance effects: During wet seasons, resistance weighted more heavily toward abundant species (Hill number Q increases) showed weaker NPK effects, suggesting dominant species were more resistant than rarer species under eutrophication in wet extremes.
• Robustness: Results were qualitatively consistent using stricter SPEI thresholds (±1.28 SD), after removing SPEI trends, and when restricting to 22 long-term sites (≥10 years). Specific significant correlation pairs varied among robustness analyses due to site inclusion differences.

The study shows that eutrophication primarily destabilizes biodiversity-related facets (species richness and community composition) by reducing both their long-term invariability and short-term resistance to dry and wet climate extremes, while biomass stability remains largely unaffected. This decoupling indicates that ecosystem function (biomass) can appear stable even as biodiversity stability erodes under nutrient enrichment, likely due to compensatory dynamics and dominance by resistant species. The positive association between invariability and resistance across aspects, coupled with a general lack of association with recovery, implies that long-term temporal stability depends more on withstanding extremes than on post-extreme rebound. The pervasive trade-off between resistance and recovery for biomass and composition suggests systems more impacted during extremes tend to recover faster, indicating reversible changes, whereas species richness lacks this trade-off, implying slower or more uncertain recovery. Weak correlations among stability of biomass, richness, and composition under both ambient and eutrophic conditions reveal a consistently high dimensionality of stability, so mechanisms regulating one aspect or facet provide limited predictability for others. Management should prioritize strategies that enhance resistance—such as incorporating drought- and flood-tolerant species or genotypes—especially under eutrophication, and recognize that preserving biodiversity stability may require different, context-specific actions beyond maintaining biomass production.

By synthesizing data from 55 globally distributed grasslands, this study demonstrates that nutrient enrichment reduces the temporal invariability and resistance of species richness and community composition during both dry and wet extremes, while leaving biomass stability largely unchanged. Stability is highly multidimensional: most facets are weakly correlated within aspects, and stability among biomass, richness, and composition is generally decoupled. These findings highlight that biodiversity stability is particularly vulnerable to eutrophication and that enhancing resistance is key to maintaining long-term invariability. Future work should develop theory and methods to quantify resistance and recovery under consecutive and repeated extremes (to account for legacy effects), identify trait- and dominance-based mechanisms underpinning facet-specific responses, and design management interventions tailored to distinct community aspects and stability facets.

• Quantifying resistance and recovery under consecutive extreme seasons (e.g., dry–wet or repeated extremes) is challenging; existing metrics assume single pulse perturbations, and legacy effects may confound estimates.
• While qualitative conclusions were robust across alternative SPEI thresholds, detrending, and longer site durations, the identity of significantly correlated facet pairs varied due to differing site subsets.
• Site durations varied (4–15 years), potentially affecting the power to detect recovery dynamics, especially for species richness.