Nitrogen and phosphorus fertilization consistently favor pathogenic over mutualistic fungi in grassland soils

1Ecosystems across the globe are receiving elevated inputs of nitrogen (N) and phosphorus (P) from agriculture and urban activities, and atmospheric N deposition has increased threefold from pre-industrial levels and will likely increase in the foreseeable future. Many plant communities are N limited, and additional N can promote plant productivity if P is non-limiting,. However, growth responses to nutrient addition are not always positive and can decline over time. This is often attributed to altered plant community composition, involving loss of dominant species and an increase in ruderal plants, but there is growing recognition that nutrient addition may also alter the composition and function of soil microbial communities, which may feedback to limit plant community productivity,.

2As mutualists, pathogens, and saprotrophs, soil fungi regulate key ecosystem processes, including plant primary productivity and carbon mineralization and sequestration,. Global patterns in fungal guild abundance are related to climate, but guilds are also affected by soil fertility. For example, most herbaceous plants, shrubs, and trees in temperate and tropical habitats form a root symbiosis with the putative mutualists, arbuscular mycorrhizal fungi (AMF). Inside roots, AMF exchange nutrients (especially P) acquired in soil for plant carbon (C). Where P is plentiful, AMF abundance is often lower due to reduced C allocation from plants. The addition of N can benefit AMF if it exacerbates plant P limitation, but may be suppressive if nitrophilic, ruderal plants replace plants that allocate more C to AMF,. Thus, responses by AMF likely depend on the extent to which nutrient addition alleviates plant deficiencies and alter plant communities. Like AMF, interactions between plants and fungal pathogens are complicated and depend on both host and pathogen responses to environmental conditions. Most research to date has been conducted on cultivated plants, and less is known about the role of pathogens in natural communities, although recent work highlights their importance for maintaining plant diversity. One emerging pattern from work in natural systems is that pathogens often thrive in resource-rich environments,, and fertilizer addition has been linked to increased disease incidence in plants,. This effect is likely greater with N than P, because P addition can enhance plant vigor, which tends to decrease susceptibility to pathogens,. Fungal saprotrophs have received more attention recently because of their key role in soil C flux and storage– where relatively small changes in abundance and activity could have outsized consequences for C budgets. While specific responses to P addition are poorly understood, N addition can affect saprotrophic activity and decomposition rates. However, such effects are hard to predict, can change over time, and depend on the saprotrophic community composition, litter quality, soil fertility, and N supply rate,–. For example, N addition may accelerate initial decomposition rates but retard turnover of more recalcitrant litter due to shifts in oxidative enzymes,. As such, it has been suggested that anthropogenic N deposition can promote C sequestration,.

3Most previous studies have focused on individual guild responses at single locations or have spanned multiple vegetation types, which complicates direct comparisons among locations within a single vegetation type. To better predict current and future functioning of soil microbial communities, we need to understand how altered nutrient availability influences fungal guilds, if responses are consistent across locations, and the nature of the underlying drivers.

4In this work, we combined a previously published dataset that documented phylum-level responses (Ascomycota, Basidiomycota, Glomeromycota, and Zygomycota) to N and P addition with novel data on AMF colonization in roots collected from the same soil cores. Because only Glomeromycota can be assigned into a functional guild at the phylum level of resolution, these higher taxonomic groups offer limited information about possible functional consequences of nutrient addition. We therefore used taxonomic information more closely related to species, and assigned potential function to taxa using FUNGuild. This allowed us to quantify responses by AMF, putative pathogens, and saprotrophs after 1–4 years of experimental N and P addition at 25 Nutrient Network (NutNet) grassland sites across four continents (Supplementary Fig. 1). We predicted that AMF would be suppressed by P and respond less consistently to N, whereas pathogens would be promoted by N and possibly be suppressed by P. Given the context-dependent nature of saprotroph responses in the literature, we expected less consistent responses by saprotrophic fungi. To identify the underlying mechanisms, we used structural equation modeling (SEM) to determine if guild responses were directly driven by N and P (while accounting for initial soil conditions), or indirectly mediated by plant community responses. We predicted that responses would primarily be driven by shifts in plant communities given the dependency on plants by all three guilds. The novel data on AMF colonization in roots allowed us to determine if N or P addition shifted AMF biomass allocation between roots and soil. Finally, we used co-occurrence networks to assess how nutrient addition may change potential interactions within and among guilds. To explore if differences in guild abundance across sites were related to inherent differences in soil properties, plant communities, and climate, we used a mediation test restricted to unfertilized control plots.

5We find that nutrient addition promotes pathogens, suppresses AMF, and reduces co-occurrences within and among fungal guilds in ways that appear primarily mediated by shifts in plant communities. We also document global distribution patterns of fungal guilds, with soil properties influencing pathogens, plant community characteristics influencing mutualists, and climate affecting saprotrophs. Our results demonstrate that fungal guilds respond predictably to nutrient addition despite substantial differences in plant and fungal communities, climatic conditions, and edaphic properties across grasslands. The generality of these patterns contributes to a growing body of knowledge that may ultimately help us build better models to identify areas where disease prevalence or severity could be particularly high, and where soils are more likely to be carbon sources than sinks under anthropogenically altered environmental conditions.

6Independent and combined N and P addition promoted fungal pathogens as we observed an average increase of 140% in N + P plots relative to control plots in relative pathogen abundances. The effect of N aligns with previous work in mostly cultivated systems where N addition has increased disease severity,. The independent promotion of pathogens by P was surprising given that P has more often had no effect or has been associated with improved plant health and reduced disease. Future work should assess if this is largely driven by a subset of taxa, as pathogens differ in their response to fertilizer. Unfortunately, this could not be done here due to an insufficient number of sequences for robust community analyses.

7We predicted that P would suppress AMF, but this only occurred when P was added together with N (−33% relative to control plots), whereas independent addition of N or P had no effect (Fig. 1). Because AMF rely on plants for all their C, this general suppression of AMF with N and P indicates that plants can somewhat adjust C flow to AMF and reduce the risk of parasitism when nutrient availability is high. Direct links between fungal guilds and plant–soil feedbacks have been shown previously. The promotion of pathogens and—to a lesser extent—the suppression of AMF with nutrient addition observed here may help explain why plant biomass responses to fertilizer decline with time,. It also supports predictions that the effects of soil biota on plant growth are more negative in resource-rich environments,,.

8To acquire C, AMF must colonize roots, but relative biomass allocation between roots and soil differs among fungal taxa and environmental conditions,. AMF occupancy of both roots and soil differed across sites here, but unlike AMF soil colonization, root colonization was not suppressed by combined N and P addition. On the contrary, root colonization appeared to have been slightly promoted by N (P = 0.105, Fig. 1), which suggests a shift in fungal allocation with N addition where a greater proportion of total fungal biomass was inside, not outside, the root. Nitrogen addition can alter AMF communities from taxa that grow extensively in soils to those that preferentially occupy roots–, but see ref.. Due to the generally low AMF sequence counts, we could not assess if shifts in composition accompanied the slight shift in allocation, here. Regardless, AMF that allocate biomass preferentially inside rather than outside roots are sometimes considered less beneficial with a lower capacity to acquire P,, but they can provide protection against pathogens,. Thus, we speculate that an intriguing alternative hypothesis to nutrient-induced parasitism is that N promotes pathogens, which initiates a switch in AMF function from resource acquisition to pathogen protection. Testing the pathogen protective ability of AMF that are promoted by N in soils where P is non-limiting would be a fruitful area for future research, because, if true, it would alter how we view the mutualism-to-parasitism continuum of this symbiosis.

9Compared with pathogens and mutualists, nutrient addition did not alter the relative abundance of fungal saprotrophs (Fig. 1 and Supplementary Fig. 1), concordant with an earlier study conducted along a soil fertility gradient. This lack of response was evident despite an increase in productivity and may be because litter biomass, and thus substrate inputs, did not increase with N or P addition (Supplementary Table 1). This lack of response is also in agreement with the observation that soil C concentration has not changed in response to nutrient addition in many of the same grasslands. It is possible that the context-dependency of responses to N addition outlined earlier may limit directional responses due to differences in climate, or that measurable responses may simply be slow to manifest. However, responses across all three guilds did not depend on the number of years elapsed since the start of the nutrient addition (P > 0.05), suggesting that factors other than time are stronger determinants of fungal responses to nutrient addition. Nonetheless, some of sites in the Nutrient Network have been maintained for >10 years, and repeated sampling could assess whether differences in the first 1–4 years decrease over time as communities acclimate to the novel nutrient conditions, stabilize and remain the same, or become more pronounced.

10In accordance with our predictions and previous work,, the SEM supported plant community mediation of fungal guild responses in most cases. AMF were suppressed by an increased plant biomass, saprotrophs were promoted by greater root biomass, and both AMF and pathogens responded to shifts in community dissimilarity. Here AMF abundance was suppressed as the phylogenetic distance increased within plant communities, whereas pathogen abundance was promoted when plant communities became more dissimilar from control communities (Fig. 2). How plant communities may drive the observed shifts is unknown. But if nutrient addition changes plant communities toward taxa with more ruderal traits, this can promote pathogens and suppress mutualists, because ruderal plants allocate less C to AMF and are more susceptible to disease,–. Two previous studies that focused on many of the same sites show that N and P addition reduces root:shoot ratios, promotes nutrient concentrations in leaves, and increases specific leaf areas,; traits that are all associated with faster life histories. Given that plant richness and functional groups are often highlighted as potential drivers of plant-soil feedbacks, it is noteworthy that plant functional groups played only a minor role compared to overall community descriptors here, and that the effect of plant richness was mostly mediated by changes in phylogenetic distance within communities (MNTD). Whether shifts in fungal guilds precede—and possibly cause—shifts in plant communities, or if fungal guilds respond to altered plant species and traits will require detailed analyses of plant and fungal community changes over time. We also cannot rule out a direct effect of N and P on pathogens, because a model replacing plant community similarity with the direct effect of treatment (while keeping soil properties) performed similarly (Supplementary Table 2). Careful experimentation could disentangle the potential bidirectional interactions of plant community and soil chemical properties on fungal guild abundances, especially pathogens.

11One unexpected finding was the relatively strong effect of pre-treatment soil nutrient concentrations (Fig. 2 and Supplementary Table 3). For instance, iron was associated with a higher saprotroph abundance and root biomass but a lower AMF soil abundance and root colonization (Supplementary Fig. 2). Even though we found directional responses to N and P addition by putative pathogens and AMF across sites, among-site variability and pre-treatment soil conditions explained more variability than the N and P effects or the differences among plant communities at the time of sampling. Several soil nutrients were correlated (Supplementary Figs. 3 and 4), so the specific effect of individual nutrients is difficult to determine. Further, pre-treatment soil conditions are likely correlated with initial plant communities, so we cannot isolate plant effects from soil effects on the fungal guilds with the available information. It is safe to say, however, that soil chemistry is an important mediator of fungal guild abundances. Soil pH is a known driver of fungal communities,, and pre-treatment pH influenced pathogens but had little effect on saprotrophs and AMF (Fig. 2 and Supplementary Fig. 2). Why pH—or other soil properties that correlate with soil pH—affects one guild more than another, and in different directions, is unknown but could have consequences for identifying sites where disease prevalence or severity could be particularly high.

12We found more complex networks in control plots than plots receiving nutrients (Fig. 3 and Supplementary Table 4 for identities). Network complexity can be reduced by perturbation, and fungal communities may respond to fertilization as a disturbance (possibly short-term), where the enriched soil nutrient environment is completely altered and homogenized in ways that could alter the available niches. Our results could be a product of habitat filtering where taxa respond similarly to environmental shifts, but do not necessarily interact. As such, more research is needed to disentangle the cause and consequences of these changing patterns. If nutrient addition reduces interspecific interactions, it could have functional consequences, as more connected soil networks are more efficient at C uptake and nutrient cycling.

13We observed clear biogeographical patterns of fungal taxa within all guilds (Fig. 4), suggesting that responses to N and P addition occurred despite substantial differences in taxonomic composition. Also, sites differed in their relative guild abundances (Supplementary Table 5). Some sites contained many saprotrophs, but very few AMF, and vice versa. When restricting analyses only to control plots (those receiving no fertilizer inputs) to assess potential underlying drivers of these distribution patterns, we found that each guild followed a distinctive pattern (Table 1). Similar to the SEM, plant communities were the most important predictor of AMF (root and soil colonization), whereas pathogens were more affected by soil properties (Supplementary Table 6 and Supplementary Fig. 2). In contrast, climate was the most important driver of saprotroph abundance. It is also possible that dissimilarities in plant community composition could be important for predicting differences in pathogen abundances, similar to our SEM and earlier findings targeting the whole fungal community. Unfortunately, we could not include a measure of plant community dissimilarity in these analyses as, per definition, it relies on pairwise comparisons and a “reference” community, analogous to the control plots used to calculate the treatment effect in the SEM. Regardless, our results build on previous work across biomes showing differential responses among fungal guilds to plant community characteristics, soil edaphic properties, and climate,. De-aggregating bulk fungal responses and examining differences among fungal guilds will help us better predict soil function and responses to future perturbations.

14While we show clear directional effects of N and P addition on fungal guilds, these findings are based on relatively short-term responses and rely on sequence counts and their link to potential function. This approach assumes that sequence counts accurately reflect abundance or biomass and that potential function informs actual function, which may not always be the case,. For example, AMF can sometimes be parasitic when costs outweigh benefits, and pathogens can occur in host tissues without causing disease. We argue that these potential discrepancies, including potential PCR biases, apply equally to all treatments and in comparisons across sites, and that our method can be used to assess relative shifts across treatments and global distributions. The use of relative sequence counts is another potential issue where responses by one guild may constrain responses by another given the proportional nature of these data. While the SEM indicated that some guilds are correlated (Fig. 2), in accordance with previous surveys, ~40% of the fungal sequence data could not be functionally assigned, so responses among the guilds we tested should therefore be relatively decoupled. Interactions between the producer and decomposer food webs also are expected, and correlations among fungal guilds have been documented previously, suggesting that the correlations observed here represent actual ecological responses. Most importantly, the shifts in pathogen to mutualist ratios with nutrient addition we observed is consistent with previously documented changes in plant responses to soil biota,. For example, plants inoculated with soil from fertilized fields were 50% smaller than plants grown in soil collected from unfertilized fields (estimated from Revillini et al., Fig. 4). Two of the four sites in Revillini et al. (Konza LTER and Short Grass Steppe LTER) showed a more than tenfold increase in pathogen to mutualist ratios between control and N + P addition plots in our analyses. Thus, while we cannot unequivocally support a direct link between pathogen abundance and disease development or growth suppression here, this has been shown previously,, supporting the likely ecological relevance of the results presented here.

15We used a published fungal dataset generated from 25 NutNet grassland sites that are distributed world-wide and where researchers follow the same treatment and sampling protocols. In brief, each grassland site is situated in a relatively homogeneous ~1000 m2 area divided into three blocks. Each block is made up of 5 ×5 m plots that are surveyed annually for plant community composition and productivity, and a subset of plots receive applications of 10 g of either N [(NH2)2CO] or P [Ca(H2PO4)2] or both m−2 per year. Soil samples for fungal community analyses were collected from control plots and plots with N and P added alone and in combination in 2011 or 2012, 1–4 years after the initial nutrient addition. DNA was extracted and amplified using fungal-specific primers (ITS1F/ITS2) targeting the internal transcribed spacer (ITS1) region. For more details on experimental design as well as soil sampling, DNA extraction, amplification, and bioinformatics see Borer et al., Leff et al., Prober et al., and supplemental material.

16We matched taxonomic identities of fungal sequences randomly rarefied to 485 sequences per plot with ecological guilds using the expert-curated database FUNGuild. This rarefaction level represents a trade-off between keeping as many sites as possible in our analyses while characterizing most taxa within each treatment at each site. This rarefaction was also chosen as it was the same as in Leff et al., thus allowing direct comparisons with that publication. While this rarefaction level did not result in an exhaustive characterization of all taxa at all sites, it did capture abundant taxa in all sites (Supplementary Fig. 6). We then assessed if sequence numbers belonging to different guilds varied across sites and among control plots and those receiving N, P, or N + P. Of the 164,900 total sequence reads, 60% belonged to taxa annotated with a guild assignment. These assignments were further subset into three guilds: highly probable and probable arbuscular mycorrhizal (AMF, putative mutualists), plant pathogenic fungi, and saprotrophic fungi. For the latter, both the soil saprotroph and undefined saprotroph guilds were included. The specific substrate or habitat of many saprotrophic fungi is undefined, and by only including soil and undefined saprotrophs, we sought to exclude those taxa known to associate with non-representative substrates such as wood or dung. Total sequence counts were 5659 for AMF, 13,263 for plant pathogens, and 24,333 for saprotrophs. To complement sequence data from soils, we also quantified the abundance of AMF in roots using trypan blue and the gridline intersect method. For this additional analysis, we included four additional sites that contained roots but not sequence data (n = 29).

17To understand overall trends in fungal communities in response to nutrient addition, we compared fungal guild sequences across nutrient addition treatments using a linear mixed-effects model and incorporating site and block within sites as random effects using lmerTest package. Fungal guild sequences were log +1 transformed to meet assumptions of normality. We calculated the partial effect of the treatment and the partial residuals using the visreg package (Fig. 1). To assess the potential mechanisms driving the treatment effect, we included three sets of covariates to determine if any of them could drive the observed patterns: post-treatment vegetation, pH and root biomass, and pre-treatment soil conditions (Supplementary Table 5). For each guild, we eliminated one of each pair of highly correlated variables identified using a variance inflation factor (VIF) and, to obtain a more parsimonious model, we retained only significant variables. Post-treatment vegetation was described using functional group descriptors or whole community descriptors, and the best model in each case was compared between them. For every regression, we required at least two control plots and eight total samples per site for inclusion. Not enough samples had all the information needed so we dropped the random block factor and excluded root colonization in the model. Also, because some predictor or response variables were missing for some sites, the SEM analyses were restricted to 15 sites with complete data matrices.

18The significant predictors (p < 0.05) from each of the three groups of covariates were pulled together to build the final regression model for each fungal guild. We combined the three regressions obtained for AMF, pathogens and saprotrophs into a structural equation model (SEM), including the effect of treatment and each significant predictor (p < 0.05) using the piecewiseSEM package. We modified these final models and obtained the marginal and conditional R2 to test if adding treatment effects or replacing the plant community information with treatments improved the models. Finally, to assess the overall explanatory power of each of those variables, we estimated the marginal and conditional R2 of models including plant community and pre-treatment soil conditions only.

19Because the relative abundance within each guild differed among sites (Supplementary Table 5), we ran a mediation test to assess whether differences in soil (including total N range: 0.05–1.6% and available P range: 7–248 mg kg−1soil), climate, or plant communities better predict the observed abundance trends. We used control plots only and averaged observed values across years (years 0–4 after establishment where available) to reduce the effects of plot-to-plot and year-to-year variation. Thus, each site was represented by a single value for each predictor and response. We identified the best predictors for each guild from each set of covariates (plant communities, soil, climate, location), selecting variables using a stepwise backward selection model, using BIC to identify the most parsimonious model. Then, for each combination of fungal guild and predictor set, we tested if the best predictors of another set could improve the model (provide new information) or decrease the significance of some of the original predictors (which could suggest that one of the variables added mediates the effect of the original predictor). See Supplementary Information for detailed methods.

20To examine co-occurrence patterns, first, we grouped all taxa that were identified to guild by genus and looked for significant correlations between different genera’s sequence abundances using the SparCC method, which accounts for compositional data (‘sparccboot’ function from the SpiecEasi package in R). Only those genera with significant correlations (both positive and negative, with p < 0.05) by bootstrapping (100 times) were included in the networks. We created a separate network for each of the four treatments (N, P, N + P, and control), where a node represents a genus, and an edge between two genera represents a significant correlation in abundances between the pair. To calculate p-values for the network statistics in Fig. 3, we performed 10,000 permutations of the networks. To calculate p-values for number of links, mean degree, and connectance, we took the edge lists for all four of the networks, which consist of each possible pairwise combination of genera (those present in any of the four networks) and 0 if no link, and 1 if linked. We then randomly shuffled the edge weights (0 or 1) across pairwise treatments, keeping the nodes constant, and reformed the two focal networks from these shuffled weights. We calculated the metrics for each of the networks for each permutation. P-values were calculated as the proportion of permutations at least as extreme as the observed differences between the network metrics. Superscripts in Fig. 4 indicate treatments that were significant at p < 0.05.

21Using the ‘cmdscale’ function in the vegan package, we performed principal coordinate analyses (PCoAs) on Bray–Curtis distances of Hellinger transformed relative abundance data for each fungal guild and plot. PCoAs were plotted using ggplot2. To detect significant differences among treatments and continents within each guild, we performed permutational analyses of variance (Permanova), using the ‘adonis’ function. Permanova analyses were run on both Bray-Curtis distances as well as on Raup-Crick transformed presence/absence data.

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