The concurrent increase in drought and atmospheric nitrogen deposition has profoundly impacted multitrophic biodiversity and ecosystem functioning in grasslands. Despite the well-documented individual effects of reduced precipitation and nitrogen enrichment, their interactive effects, especially on multitrophic cascading responses (e.g., plant, nematode, and arthropod communities), remain largely unknown.

Using a four-year field experiment in a typical steppe of Inner Mongolia, we explored the effects of three drought scenarios (intense drought, excluding 100% of rainfall in June; reduced precipitation frequency, reducing rainfall events by 50% without changing total rainfall from June to August; and chronic drought, excluding 50% of each rainfall event from June to August) and nitrogen addition (+10 g N m–2 yr–1) on species diversity, functional group abundance, and functional group associations within and between trophic levels, including plant, ground-dwelling arthropod and soil nematode communities, as well as their relationships with grassland productivity.

We found that: (1) Drought and nitrogen addition had contrasting effects on multitrophic species diversity, functional group abundance, and grassland productivity. Chronic drought significantly reduced productivity independent of nitrogen, while nitrogen addition enhanced it. Intense drought increased the abundance of bacterivorous and fungivorous nematodes, but this trend was absent with nitrogen addition. Reduced precipitation frequency had no significant effect on multitrophic communities or productivity under any nitrogen condition. (2) Drought, particularly chronic drought, enhanced positive associations between ground-dwelling arthropod and soil nematode functional groups, whereas nitrogen addition was accompanied by a weakening of these functional group interconnections. (3) Increased productivity with nitrogen addition was associated with reduced positive associations within plant functional groups and between arthropod and nematode functional groups, along with increased soil nitrogen availability. Drought was related to lower productivity overall, although it also coincided with reduced associations within plant functional groups, which were related to higher productivity.

Synthesis. Our results indicate that nitrogen deposition under drought scenarios is linked to adverse effects on the abundance and interconnections of multitrophic functional components, highlighting the importance of species interactions across trophic levels for understanding grassland responses to environmental change. 

Species richness, Hill–Shannon, Hill–Simpson (Chao et al., 2014), and Pielou evenness indices were used to assess the species diversity of plants, soil nematodes, and ground-dwelling arthropods. Hill–Shannon diversity was calculated as eH′ = exp (−Σpi lnpi), Hill–Simpson diversity as 1/ Σpi2, and Pielou evenness as J′ = H′/ ln(S), where pi is the relative abundance of species i, H′ is the Shannon index, and S is species richness. Data from repeated measurements of soil moisture and available nitrogen were averaged prior to statistical analysis to facilitate pairing with other data.

First, linear mixed-effects models ("nlme" package in R) were used to evaluate the effects of drought scenarios, nitrogen addition and their interaction on ANPP, soil available nitrogen, soil moisture, and species diversity and the abundance of functional groups of plants, soil nematodes, and ground-dwelling arthropods, with drought and nitrogen addition as fixed effects and block as a random effect. Data on the abundance of fungivorous nematodes were square root-transformed to meet the assumption of normality and homogeneity of variance. Duncan’s multiple range tests were used for post hoc comparisons among drought scenarios within each nitrogen level, while paired t-tests were used for assessing nitrogen addition effects within each drought scenario.

Second, correlation networks were constructed to examine the associations among species from multiple functional groups across different trophic levels (including perennial tall bunch grasses, perennial short bunch grasses, perennial rhizome grasses, perennial forbs, and annual and biennial herbs; plant-parasitic, bacterivorous, fungivorous, and omnivorous-predatory nematodes; and phytophagous, predatory, omnivorous, saprophagous, and parasitic arthropods) using the R packages "ggplot2", "igraph", and "visNetwork". Within each treatment, pairwise associations were calculated using Spearman's rank correlations among species from each pair of groups. The number of significant positive correlations (P< 0.05) between species in each pair of functional groups was divided by the total number of all possible correlations to obtain the average association strength. Positive and negative associations within and between trophic levels were also quantified based on the correlation networks. Network topological properties, including number of edges, connectance, mean degree, modularity, and average clustering coefficient, were calculated to characterise network structure (Ma et al., 2016). Sensitivity analyses were performed using the false discovery rate (FDR) correction to control for multiple comparisons. Pearson correlations were used to analyse the relationships between positive associations among multitrophic functional groups and ANPP.

Finally, a piecewise structural equation model (SEM) was used to explore the direct and indirect effects of drought and nitrogen addition on ANPP. A priori model was based on the bivariate relationships between ANPP and the diversity and functional group abundance of different biological communities, as well as on correlations among multitrophic functional groups (see Figure S2 for the explanation and illustration of each pathway). The initial model was simplified by removing insignificant pathways and variables according to regression weight estimates (Lefcheck, 2016), retaining key pathways while maintaining a reasonable sample-to-parameter ratio. Model fit was evaluated using Fisher's C statistic (P> 0.05 indicates good fit) and AIC values (Akaike Information Criteria, lower AIC indicates better fit) (Shipley, 2013). The model was constructed using the R package "piecewiseSEM". To test robustness, SEM was also rerun using FDR-corrected correlation matrices. All analyses were conducted using R version 4.4.1 (R Core Team, 2022).