1 Introduction
Fertilization is a widely used and an essential practice to increase crop yields and maintain soil fertility (Steiner et al., 2007), and is also one of the governing factors influencing soil microbial communities and functions in agroecosystems (Geisseler and Scow, 2014). As heterotrophs, soil fungi not only mediate organic matter turnover and nutrient transformation (Aguilar-Trigueros et al., 2015), but also respond to addition of nutrients by changes in their biomass and community composition (Sun et al., 2016). Compared to bacterial communities, fungal communities have less functional redundancy due to their lower resistance in face of environmental shifts (Wohl et al., 2004), which leads to greater response of fungal communities to fertilization disturbances (Banerjee et al., 2016). Hence, better understanding and effective management of fungal resources are critical for promotion of crop nutrient uptake and optimization of crop productivity (Ellouze et al., 2014). To date, the majority of researches in this area have focused on the influences of organic and chemical fertilization on fungi (Ullah et al., 2019; Xiang et al., 2020). Applied organic fertilization generally has positive effect on fungal communities, such as increase of soil fungal diversity and decrease of the abundance of soil-borne fungal pathogens (Sun et al., 2016; Wen et al., 2020). Conversely, it was reported that the richness and abundance of soil fungi would be reduced, and the growth of potentially pathogenic fungi would be promoted following long-term addition of chemical fertilizers (Wang et al., 2018; Ye et al., 2020). Nevertheless, no study has specifically addressed the negative effects of application of chemical N, P or K alone or their combinations on fungal communities within an agroecosystem.
As the essential plant nutrients, the availability of soil N, P and K regulates soil microbial communities and functions in various ecosystems (Su et al., 2014; Turner and Wright, 2014). N addition has previously been reported to promote fungal growth in a few cases (Zhao et al., 2020; Li et al., 2021), but negative influences were more frequently documented (Treseder, 2008; Morrison et al., 2016; Spohn et al., 2016), and the main reasons were ascribed to the aggravation of soil N saturation, decrease of soil pH and mobilization of aluminum (Lu et al., 2014; Tang et al., 2018). Comparatively, weaker effects following P addition on fungal communities were reported (Camenzind et al., 2018). For example, regardless of P addition, N addition significantly changed soil fungal community composition in an alpine grassland of the Qinghai-Tibetan Plateau in China (Chen et al., 2020). Nevertheless, the significant decreases in diversity and richness of arbuscular mycorrhizal fungi (AMF) induced by long-term P fertilization are well documented (Lin et al., 2012; Liu et al., 2012c), and this negative effect of P addition on mycorrhizal colonization may consequently favor the growth of fungal pathogens (Mujica et al., 2020). Compared with the studies on N and P fertilizers, there is scanty understanding on the effects of K addition on microbial biomass and community. Kaspari et al. (2008) reported that the K-limitation of microbial decomposition in a tropical forest, and Luizão et al. (2007) presented evidence of positive K effects on soil microbial respiration in central Amazonian forest soils. In this context, due to the limited and inconsistent data, the response of fungal communities to addition of N, P and K fertilizers and the comparison of their effects need to be explored by more field experiments. Furthermore, in addition to N fertilizers, most studies of P addition were conducted in forest ecosystems, in which P rather than N limitation is widely recognized (Camenzind et al., 2018). Thus, it is urgent to call for a systematic investigation on the impacts of routine chemical fertilizations (N, P, or K alone and their combination) on the structure and function of soil fungal community in agroecosystems, which is an essential prerequisite for fertilization strategy optimization.
The black soil region in North-east China is one of the four famous black soil areas in the world (Liu et al., 2012b). With inherently fertile and productive soils, this region is well known as an important commodity grain bases and plays a pivotal role in maintaining food security in China. Chemical fertilizers, i.e., N, P, K and their combinations are widely applied fertilization practices in this region. In our previous studies (Hu et al., 2019; Yu et al., 2019), we found that addition of N fertilizers rather than P and K fertilizers significantly decreased bacterial and diazotrophic community diversity, and N fertilizations exerted stronger effects on bacterial and diazotrophic community structure than those of P and K fertilizers in black soils. However, it is still unclear whether similar shifts are also found in fungal communities and how fungal taxa respond to long-term addition of N, P and K fertilizers. To fill this gap, using the same soil samples as reported previously (Hu et al., 2019), the changes of soil fungal communities were investigated in this study. We hypothesized that (1) addition of N fertilizers would exert greater effects on fungal communities and co-occurrence network pattern relative to those of P, K fertilizations, (2) abundance of responding potential plant pathogens and their interactions would be stimulated by addition of N fertilizer that was ultimately detrimental to crop yields.
2 Materials and methods
2.1 Field experiment, soil sampling and soil properties analysis
Soil samples were collected from a long-term fertilization experimental station established since 1980 in Minzhuxiang of Harbin City, Heilongjiang Province, China (45°50.53′N, 126°51.13′E), as described in our previous study (Hu et al., 2019). Briefly, the site is subjected to wheat-maize-soybean annual crop rotation and soil samples were collected under soybean cultivation in 2014. Eight fertilization treatments, i.e., non-fertilization (NoF); addition of N, P, and K alone; addition of combined fertilizers NK, NP, PK and NPK, were randomly arranged in plots with three replicates, and each replicate contained 8 rows with 6 m long and 70 cm wide. The N, P and K fertilizers were applied in the form of urea, calcium superphosphate and potassium sulfate, respectively. Five individual soil cores were randomly collected at a soil depth of 0–20 cm and pooled into one sample to minimize variation within a plot. A total of 24 soil samples (8 treatments × 3 reps.) were obtained. The soil sample processing and soil properties determination were reported in our previous study (Hu et al., 2019). The soil properties and soybean yields were stated in the supplementary results (Text S1).
2.2 Soil DNA extraction and quantitative PCR
For each sample, soil total DNA was extracted from 0.5 g fresh soil using a Fast DNA®SPIN Kit for Soil (MP Biomedicals, USA) and stored at −20°C until use (Hu et al., 2019). The fungal abundance of all samples was examined in a LightCycler®480 (Roche Applied Science) with primers ITS1F (5′-CTT GGT CAT TTA GAG GAA GTA A-3′) and ITS2 (5′-GCT GCG TTC TTC ATC GAT GC-3′) (White et al., 1990). The quantitative PCR (qPCR) reaction and amplification conditions were previously described in Hu et al. (2017). The fungal ITS copy numbers were calculated using a standard curve that was generated from a clone with an ITS gene insert.
2.3 MiSeq sequencing and bioinformatics analyses
Amplicon libraries were created by using a primer pair ITS1F/ITS2 targeting the fungal ITS1 region. A 6-bp barcode was fused to each forward primer for distinguishing the amplified products. The amount of purified PCR products was examined using a Qubit fluorometer (Invitrogen, USA), and mixed in equimolar amounts to generate one library. The qualified library was sequenced using the Illumina MiSeq sequencing platform (Majorbio, Shanghai, China). The raw sequence data were submitted into NCBI with the accession number SRP 110636.
The raw sequences were processed in QIIME pipeline (http://qiime.sourceforge.net/). The barcode, primers and the low-quality reads with mean base quality score<20 were discarded. High-quality sequences were assigned to 1820 OTUs at 97% similarity level using the UPARSE (Edgar, 2013). The representative sequence of each OTU was blasted against the UNITE database (ver. 8.0). All samples were rarefied to 25604 sequences per sample to meet the even sequencing depth for comparison of community diversity. Alpha diversity of the fungal community was displayed by Shannon diversity with calculation using the alpha_diversity.py function in QIIME.
Fungal OTUs was classified into functional guilds with the FUNGuild database ver. 10 (Nguyen et al., 2016) with 795 OTUs allocated in this study. Only results with confidence levels of “probable” or “highly probable” were accepted for further analyses. Three main fungal trophic modes and eight guilds with some overlaps were concentrated in this study: saprotroph included dung saprotroph, wood saprotroph and plant saprotroph; pathotroph included fungal parasite, animal pathogen and plant pathogen; symbiotroph included ectomycorrhiza and endophyte.
2.4 Microbial network construction
Three fungal networks were analyzed for non-fertilization (NoF), for combination of P, K fertilizations (P, K, PK treatments) and for combination of N fertilizations (N, NK, NP, NPK treatments) by processing fungal OTUs within all samples. The similarity of OTUs was measured by calculating Spearman correlations, and the adjacency matrix was calculated based on the random matrix theory (RMT) method through the MENA pipeline (Zhou et al., 2010; Deng et al., 2012). The similarity threshold was automatically determined, and network topology characteristics were processed in the pipeline. The network graphs were visualized by Gephi software (Ver. 0.9.2).
2.5 Statistical analyses
The effects of different chemical fertilization treatments on soil properties, fungal abundance, relative abundances of fungal taxa, fungal diversity and community resistance were examined by analysis of variance (ANOVA, Turkey’s test) in SPSS software (ver. 21.0). The principal coordinates analysis (PCoA) based on the Bray-Curtis distance was conducted to reveal the dissimilarity among different fertilization treatments via the “vegan” package in R (ver. 4.0.2). Volcano plots, generated via the “ggplot2” package in R, was used to identify the significantly enriched or depleted taxa after addition of fertilizers. A comparison of the fungal taxa was conducted between fertilization treatments and the reference (NoF) by the Wilcox-test via the “stats” package in R, and then corrected by a false discovery rate (FDR). The response ratio (RR) analysis was conducted to quantify significant responses of potential plant pathogens in fertilization treatments compared to those of NoF, which was calculated at the confidence interval of 95% and P<0.05 was regarded as level of significance (Lajeunesse, 2011). The Spearman correlation coefficient was used to test the effects of soil properties on fungal taxa and community diversity, which was conducted in SPSS. The partial Mantel test was applied to analyze the relationships between soil properties and fungal community and network structures via the “vegan” package in R. A heatmap graph was generated to display the primary soil factors dominating the variations of fungal taxa via the “pheatmap” package in R. Structural equation modeling (SEM) (Grace, 2006) was used to evaluate the direct or indirect relationships among responding taxa (RT, the changeable taxa induced by fertilizers addition compared to NoF), responding plant pathogens (RP, the changeable taxa induced by fertilizers addition compared to NoF), responding plant pathogens in network (RPN), fungal alpha and beta diversity, community resistance, network connectivity and soybean yields, and was performed in SPSS AMOS (ver. 21.0). All the variables used in SEM were identified with the average standardized abundance of these data.
2.6 Calculation of resistance
The resistance indices of fungal communities were calculated according to the methodology of Orwin and Wardle (2004) as:
where C0 is fungal community of NoF indicated by PCoA1 value, and D0 is the difference of fungal community between NoF and fertilizers addition treatments indicated by PCoA1 value. The values of 1 represented fertilization with no effect or a full recovery of fungal communities, whereas values lower than 1 represented less resistance response to fertilization or a slower recovery rate of fungal communities (Cabrerizo et al., 2019).
3 Results
3.1 Fungal abundances and communities under long-term chemical fertilizations
After MiSeq sequencing, a total of 852 252 high-quality sequences with 46 649 to 25 604 per sample were obtained from 24 soil samples. The high-quality sequences, except unclassified sequences, were assigned to 5 fungal phyla. Ascomycota was the overwhelmingly dominated phylum with relative abundance ranging from 44.5%–74.7% across all soil samples, followed by Zygomycota (9.13%–23.9%) and Basidiomycota (4.82%–44.5%) (Fig. S3a). Glomeromycota and Chytridiomycota were occasionally detected at low frequencies. No significant difference was observed in fungal phyla among fertilization treatments. However, at the genus level, the relative abundances of several genera were distinctly shifted by different chemical fertilizations (Fig. S3b). For example, addition of P, K fertilizers significantly increased the Mortierella abundances while decreased the Guehomyces abundances compared to NoF. In contrast, addition of N fertilizers, especially NPK greatly enhanced the Guehomyces abundances.
Furthermore, compared to NoF, all the chemical fertilizations significantly increased fungal abundances, especially N, NK, NP and NPK addition (Fig. S4a). However, addition of P, K, PK, NP and NK fertilizers had no significant effect on alpha diversity of fungal community, while NPK clearly reduced this index (Fig. S4b). For the beta diversity, it is interesting to find that the fungal communities of P, K and PK treatments were clustered with that of NoF, which were separated from N fertilization treatments along with the PCoA1 axis (Fig. S4c). Based on this finding, the fungal abundances and communities were reanalyzed for NoF, for P, K fertilizations and for N fertilizations, which showed that N fertilizations largely increased fungal abundances, but decreased fungal diversity index when compared with NoF and P, K fertilizations (Fig. 1A, B). A similar phenomenon was detected in fungal community structure and the variation induced by addition of N fertilizers was clearly larger than that of P, K fertilizations when compared to NoF (Fig. 1C). It should be noted that application of P, K fertilizers led to the most resistant fungal community with the resistance values over 0.8, which is significantly greater than those of N fertilization treatments (Fig. 1D).
Fig.1 Fungal abundance for NoF, for P, K fertilizations and for N fertilizations assessed by qPCR (A). Shannon diversity index for NoF, for P, K fertilizations and for N fertilizations (B). Dissimilarity of fungal community structure between and within fertilization regimes based on their Bray-Curtis distance (C). Index values of fungal community resistance after different fertilizers addition (D). P, K fertilizations include P, K and PK fertilization treatments, and N fertilizations include N, NK, NP and NPK fertilization treatments. Bars with different letters indicate significant differences (P<0.05, ANOVA). |
3.2 Fungal responding taxa to long-term chemical fertilizations
To clarify the variation of fungal taxa in response to different chemical fertilizations, the significantly enriched and depleted taxa after addition of fertilizers was calculated and displayed by volcano plots (Fig. 2A). The proportions of responding taxa in N, NK, NP and NPK were significantly greater than those in P, K and PK. The count of responding taxa (CRT) relative to NoF, the variation of their relative abundances (VRA), and the variation amplitude of per responding OTU (VRA/CRTOTU) were further calculated based on the pairwise comparison between fertilization treatments and NoF (Table 1). Most fungal taxa were silent microbes in all fertilization treatments, but there were 44 genera and 330 OTUs with significant variation in N fertilizations vs NoF, which were significantly greater than those of P, K fertilizations vs NoF (15 genera and 88 OTUs). The total relative abundances of varied genera and OTUs accounted for 10.9% and 17.1% in N fertilizations, and 0.58% and 0.81% in P, K fertilizations. The responding amplitude per OTU in N fertilizations vs NoF presented with a substantial variation in a mass of fungi, while relatively lower responding amplitude in P, K fertilizations vs NoF presented with a moderate variation in a small number of fungi (Fig. 2B).
Fig.2 Volcano plots showing the significantly enriched and depleted fungal OTUs after addition of chemical fertilizers (A). Percentage proportions of the significantly varied OTUs. Variation patterns of responding fungal OTUs after P, K fertilizers and N fertilizers addition (B). The x-axis indicates the number of fungal OTUs, and the y-axis indicates the variation of the relative abundance (VRA). P, K fertilizations include P, K and PK fertilization treatments, and N fertilizations include N, NK, NP and NPK fertilization treatments. |
Tab.1 Variation of fungal communities based on responding taxa after P, K fertilizations and N fertilizations addition relative to NoF. |
| P, K fertilizations vs NoF | N fertilizations vs NoF | |||
|---|---|---|---|---|
| CRTa | Phylum | 0 | 0 | |
| Genus | 15 | 44 | ||
| OTU | 88 | 330 | ||
| VRAb | Phylum | 0 | 0 | |
| Genus | 0.58% | 10.9% | ||
| OTU | 0.81% | 17.1% | ||
| VRA/CRTOTU | 0.01% | 0.05% | ||
a CRT, count of responding taxa. b VRA, variation of the relative abundance. P, K fertilizations include P, K and PK fertilization treatments, and N fertilizations include N, NK, NP and NPK fertilization treatments. |
The top five responding fungi with most enrichment in P, K fertilizations vs NoF were OTU1189 (Mortierella), OTU1443 (Mortierella), OTU1407 (Lasiosphaeriaceae), OTU1570(Tetracladium) and OTU1646 (Purpureocillium), and the most depleted fungi were OTU1703 (Cryptococcus), OTU513 (Stachybotrys), OTU685 (Pseudogymnoascus), OTU493 (Schizothecium) and OTU1803 (Guehomyces). The most enriched fungi in N fertilizations vs NoF were OTU849 (Gibberella), OTU205 (Setophoma), OTU1630 (Chaetomiaceae), OTU78 (Helotiales) and OTU1038 (Microascaceae), and the most depleted fungi were OTU1703 (Cryptococcus), OTU1307 (Mortierella), OTU1057 (Chaetomiaceae), OTU1734 (Mortierella) and OTU64 (Schizothecium) (Fig. S5).
3.3 Prediction of fungal functional profiles
Based on the FUNGuild database, saprotroph, pathotroph and symbiotroph with eight guilds were categorized in NoF, P, K fertilizations and N fertilizations (Fig. 3A). N fertilizations induced the highest relative abundances of main functional guilds compared to NoF and P, K fertilizations. Specifically, the relative abundance of endophytic fungi increased from 5.34% in NoF to 40.19% in N fertilizations, wood saprophytic fungi increased from 1.69% to 16.9%, and plant pathogens and animal pathogens sharply increased from 12.4% to 32.1% and from 1.71% to 19.6%, respectively. In contrast to N fertilizations, addition of P, K fertilizers resulted in relatively lower abundance of endophytic fungi than NoF, and similar relative abundance of animal pathogens were detected between P, K fertilizations and NoF. Although slightly higher relative abundance of plant pathogens was observed in P, K fertilizations than that of NoF, the value was much smaller when compared to those of N fertilizations. Compared to NoF, the responding functional guilds were further calculated for P, K fertilizations and N fertilizations, and all the responding function guilds in N fertilizations presented higher relative abundances than those in P, K fertilizations, especially in wood saprotroph, animal pathogens and plant pathogens (Fig. 3B).
Fig.3 Relative abundance of fungal functional modes and guilds assigned using FUNGuild database for NoF, for P, K fertilizations and N fertilizations (A). Relative abundance of responding functional guilds under P, K fertilizations and N fertilizations (B). P, K fertilizations include P, K and PK fertilization treatments, and N fertilizations include N, NK, NP and NPK fertilization treatments. |
The OTUs of responding potential plant pathogens (RP) in P, K fertilizations and N fertilizations were further selected for comparison with those of NoF (Fig. S6). Obviously, the number of RP was much more in N fertilizations (21 taxa) than those in P, K fertilizations (5 taxa), and similar phenomena were observed in the relative abundances of these taxa. Curvularia and Leptosphaeria were concurrently detected in both P, K fertilization and N fertilization treatments, but all the fertilizations decreased the relative abundance of Leptosphaeria when compared to NoF. Additionally, Gibberella with high abundance in NoF (2.82%) was significantly enriched in N fertilizations (6.04%), and Microascaceae, Chloridium and Cyphellophora also presented a similar trend in N fertilizations.
3.4 Fungal networks under long-term chemical fertilizations
Three fungal networks for NoF, for P, K fertilizations and for N fertilizations were displayed (Fig. S7a, b, c), and the fundamental topological properties were listed in Fig. S7d. Similar number of nodes and edges were detected in networks of P, K fertilizations and N fertilizations, which were both higher than those of the NoF network. Although positive interactions dominated across all networks, the proportion of positive interaction was observed to be lower in N fertilizations than those in either NoF or P, K fertilizations. Similarly, addition of N fertilizers also resulted in lowest values of average degree (number of edges of a node to other nodes), and connectedness (including average clustering coefficient, network density and network centralization) when compared to NoF and P, K fertilizations.
Moreover, the subnetwork of responding potential plant pathogens (RPN) was generated for P, K fertilizations and N fertilizations (Fig. 4). N fertilizations harbored a much more complex pathogenic subnetwork, with 71 nodes and 319 edges, than those of P, K fertilizations, with 12 nodes and 27 edges. Positive interactions accounted for 100% in the subnetwork of P, K fertilizations, while 77.1% were observed in the subnetwork of N fertilizations. N fertilizations induced higher average degree relative to P, K fertilizations. Especially, many potential plant pathogens, such as Rhodotorula, Periconia, and Chloridium, presented intensive interactions with other fungal members in N fertilizations. The detailed interaction information of responding potential plant pathogens in subnetworks was listed in Table S1. It should be noted that some pathogens in N fertilizations, such as Microascaceae, Periconia and Cyphellophora presented positive interactions with Rhodotorula (OTU890), Gibberella (OTU1749) and Aspergillus (OTU1310), but some of them had negative interactions with Mortierella (OTU1443 and OTU1734).
Fig.4 Subnetwork of responding plant pathogen under P, K fertilizations (A) and N fertilizations (B), and relevant subnetwork parameters shown in table (C). P, K fertilizations include P, K and PK fertilization treatments, and N fertilizations include N, NK, NP and NPK fertilization treatments. Blue and red lines indicate positive and negative interactions between two individual nodes, respectively. The responding nodes assigned to potential plant pathogen are colored in red and identified at the genus level or higher. The size of each node is proportional to its degree (i.e., the number of edges of a node to other nodes). |
3.5 Soil properties and fungal communities associated with soybean yields
The variations of soil properties induced by different chemical fertilization treatments significantly affected fungal communities and network patterns. Spearman correlations revealed that only soil pH had significantly positive and negative correlations with the relative abundances of Ascomycota and Basidiomycota, respectively (Fig. S8b). At genus level, soil pH, TC, TN and NO3–-N were the most influencing factors in variation of relative abundances of several fungal genera (Fig. S8b). Soil pH, TC, TN and NO3–-N explained the most important variations in responding taxa (RT), responding plant pathogens (RP) and the responding pathogens detected in network (RPN), which were clearly displayed in heatmap plots (Fig. S9). These four soil properties also presented significant influence on fungal diversity and community resistance examined by Spearman correlations, and fungal community structure and network connectivity revealed by partial Mantel test (Table 2). Specifically, soil pH was the most influencing factor on these fungal communities and network structure.
Tab.2 Effects of soil properties on fungal Shannon diversity, community structure, community resistance and network connectivity. |
| Variables | Shannon diversity | Community structure | Community resistance | Network connectivity | |||||||
|---|---|---|---|---|---|---|---|---|---|---|---|
| ra | P | ra | P | rb | P | rb | P | ||||
| pH | 0.625** | 0.001 | 0.938** | <0.001 | 0.918** | <0.001 | 0.368** | 0.001 | |||
| TC | –0.404* | 0.048 | 0.704** | <0.001 | 0.547** | 0.009 | 0.097* | 0.010 | |||
| TN | –0.469* | 0.021 | 0.845** | <0.001 | –0.742** | <0.001 | 0.294** | 0.001 | |||
| TP | 0.130 | 0.543 | –0.242 | 0.254 | 0.439* | 0.046 | –0.176 | 1.000 | |||
| TK | –0.183 | 0.393 | 0.365 | 0.080 | 0.040 | 0.864 | –0.043 | 0.867 | |||
| NH4+-N | –0.238 | 0.262 | 0.319 | 0.129 | –0.349 | 0.121 | –0.095 | 1.000 | |||
| NO3–-N | –0.514** | 0.009 | 0.8107** | <0.001 | –0.737** | <0.001 | 0.258** | 0.001 | |||
| AP | 0.090 | 0.677 | –0.189 | 0.376 | 0.421 | 0.057 | –0.183 | 1.000 | |||
| AK | –0.183 | 0.392 | 0.100 | 0.642 | 0.025 | 0.913 | –0.128 | 1.000 | |||
a Correlations between fungal Shannon diversity, community resistance and soil properties examined by Spearman correlation coefficient. b Correlations between fungal community structure, network connectivity and soil properties examined by partial Mantel test. **, P<0.01; *, P<0.05. |
All the responding taxa, fungal community and network properties were further linked to crop yields, which was evaluated by SEM analysis (Fig. 5). RT had significantly positive and negative correlations with fungal diversity and community structure, respectively. RP was positively correlated with fungal diversity, community structure and network connectivity, but negatively correlated with community resistance. Negative correlations were observed between RPN and fungal diversity, community resistance and network connectivity. Nevertheless, all the responding taxa had no significant and direct effect on soybean yields, and the similar phenomena were observed in the relationships with fungal diversity and network connectivity associated with soybean yields. It should be noted that fungal community resistance presented significant and positive correlation with soybean yields (r = 0.424; P<0.01).
Fig.5 Structural equation modeling (SEM) relating responding taxa, fungal community indices and network structure affect soybean yields. Solid and dashed lines indicate significant and insignificant relationships, respectively. The widths of arrows represented the effect strength, with blue and red indicating positive and negative relationships, respectively. RT, responding taxa; RP, responding plant pathogen; RPN, responding plant pathogen detected in network. Alpha diversity represented by Shannon diversity, beta diversity represented by PCoA1 index, network connectivity as an indicator of network structure which is equal to average node degree. R2 values indicate proportion of the variance explained for the endogenous variable. , Chi-square value; df, degree of freedom; GFI, goodness of fit index; RMSEA, root mean square error of approximation. ***, p<0.001; **, p<0.01. |
4 Discussion
4.1 Stronger responses of fungal abundance, community and network structures to N fertilizations
Although soil pH decreased by N addition generally had negative effects on bacterial abundances (Shen et al., 2010), it exerted little inhibition on fungal abundance due to the stronger cell walls and pervasive mycelia of fungi (Ma et al., 2017). Contrarily, soil fungi generally benefited from addition of mineral N fertilizers as reviewed by Geisseler and Scow (2014), and similar phenomenon had been found in this study. The effect might be due to the favorable soil C/N ratio for fungal growth.
Similar to the changes of bacterial and diazotrophic communities examined in this location previously (Hu et al., 2019; Yu et al., 2019), the results of this study supported our first hypothesis that N fertilizations would exert stronger influences on fungal community than P, K fertilizations. Effects of N on microbes would become more pronounced over time (Zhang et al., 2018), and 35-year continuous addition of chemical fertilizers had been applied in this study. The significant decrease of fungal diversity and shifts in fungal community structure under N fertilizations resulted in lower community resistance and slower recovery from environmental disturbance (Fig. 1D, Girvan et al., 2005). N fertilizations significantly affecting fungal communities was likely due to the less N requirement of fungi per unit biomass C accumulation (Zhou et al., 2017) and the accumulation of soil organic matter caused by long-term N addition disrupt the habitats of fungal taxa (Morrison et al., 2016). In contrast, other past researches have indicated that P limitation for fungal growth predominated in forest ecosystems (Liu et al., 2012a; Turner and Wright, 2014), and P rather than N addition significantly changed fungal community structure through relieving fungal P limitation (Ma et al., 2021). The inconsistent influences of different nutrients on soil fungal processes largely depended on soil types and land-uses (Hayden et al., 2010). Although soil P and K contents were elevated by addition of long-term P, K fertilizers (Fig. S1), this accumulation might lead to physiologic adaptations of soil fungal communities in cropland with high community resistance to nutrient variations (Fig. 1D).
In this study, addition of N fertilizers significantly enriched the relative abundances of saprophytic fungi (Fig. 3A), which has been reported in various ecosystems as the increase of resources (i.e., N and C) contents favor these copiotrophic fungi (Song et al., 2015; Morrison et al., 2016; Yao et al., 2021). Abundant endophytes were detected in N fertilization treatments (Fig. 3A), which have previously been found to confer benefits to their crop host through insect resistance/repellence due to the secretion of alkaloids (Malinowski and Belesky, 2000). However, the degree of mutual benefit between host and endophytes is dependent on environmental conditions (Saikkonen et al., 2006), and negative effects of endophytes on their hosts have been observed under certain circumstances (Malinowski and Belesky, 2006). In support of our second hypothesis, the proportion of potential pathogenic fungi was obviously increased following N addition rather than P, K fertilizations (Fig. 3A), which probably exerted negative impacts on crop health as plant productivity respond strongly to the richness of pathogenic fungi (Peay et al., 2013). This phenomenon suggested that reasonable nutrient types and amounts are crucial to increase soil health and reduce the morbidity of crops.
Long-term chemical fertilization complicated fungal network structure (Fig S6), indicating the great potential for fungal interactions or niche-sharing stimulated by addition of fertilizers (Shi et al., 2016). In this work, positive interactions predominated across three fungal networks, however, relatively lower cooperative relationships among fungi were detected in network of N fertilizations (Fig. S7). This phenomenon suggested a decrease in fungal niche breadths under N fertilizations, which possibly aggravated taxa competition and led to production of toxins (Becker et al., 2012; Banerjee et al., 2016). Moreover, compared to NoF and P, K fertilizations, lower average degree and connectedness were detected in the network of N fertilizations, which suggested that the less connected network structure in N fertilizations might contribute to inefficient carbon utilization and resource transfer (Morriën et al., 2017). Compared to N fertilizations, higher connectedness in network of P, K fertilizations suggested more harmonious and stable fungal community that was moderately affected by P, K fertilizers (Yao et al., 2021). Additionally, the subnetwork of responding plant pathogens was reconstructed, which showed that N fertilizations harbored a more complex pathogenic subnetwork than that of P, K fertilizations (Fig. 4). Average degree represented the number of neighbors of a node, and higher value of this index in pathogenic subnetwork of N fertilizations probably promoted the propagation of fungal pathogenicity across the whole network (Wang et al., 2020). Nevertheless, it should be noted that different fungal co-occurrence networks under fertilizer addition depend on the statistically determined interactions among the relative abundance of OTUs, and whether these positive or negative interactions would occur under field conditions is debatable and requires further targeted confirmation.
4.2 Large alternations of fungal taxa under N fertilizations
Compared to P, K fertilizations, substantial variations of fungal taxa and their abundances were induced by addition of N fertilizers (Fig. 1A, Fig. 2). Consistent with the results of Spohn et al. (2016), high N availability not only suppressed microbial respiration but had a negative effect on utilization efficiency of soil microbial carbon. P, K fertilizations, on the other hand, presented no significant effect on microbial nutrient cycling. This suggested that long-term N fertilization strengthened the niche selection for microbes by increasing their competition for resources (Zhou et al., 2014), which was also confirmed by the high proportion of competitive interactions among fungi in network of N fertilizations (Fig. S7). Despite the inference of agroecological function from OTU data should be interpreted cautiously, the responding taxa with potentially important known functions were inspected for P, K fertilizations and N fertilizations (Fig. 3B). Clearly, Gibberella was found to be much responsive to addition of N fertilizers with a high relative abundance of 6.04% (Fig S4). This genus is a notably serious and prevailing pathogen of agricultural crops, causing crop root rot and tissue necrosis with its mycotoxins (Pioli et al., 2004), and its growth stimulated by N sources had been reported in cultivation-based study previously (Bonn and Cappellini, 1970). Differently, significant decrease of Mortierella abundance was observed in N fertilizations, and members of this genus have benefits for crop growth and health by degrading hemicellulose and chitin, meanwhile producing a range of antibiotic compounds to against some plant pathogens (Tagawa et al., 2010; Xiong et al., 2017). The negative effects of addition of N fertilizers on Mortierella growth were possibly due to its sensitivity to the amount and balance of C and N (Koike et al., 2001). In contrast, previous studies have identified the effectiveness of Mortierella in soil phosphate utility from various soil types (Osorio Vega et al., 2015), which might be the reason that a high proportion of Mortierella was observed in P, K fertilizations with possible promotion of the P uptake by crops (Fig. S5, Tamayo-Velez and Osorio, 2017).
Strong interactions among microbial taxa reflected their potential ecological synergistic or antagonistic relationships in diverse ecosystems (Fuhrman, 2009; Barberán et al., 2012). In this study, abundant potential plant pathogens belonging to the guilds of responding fungal taxa (RPN) were found in subnetwork of N fertilizations, while just solitary members of RPN were found in subnetwork of P, K fertilizations (Fig. 4). Species of Rhodotorula are common environmental yeast and recognized as emerging opportunistic pathogens in soil, water and air (Wirth and Goldani, 2012). Although no direct evidence discovered in plant disease is associated with this genus, growing concern on food contamination in humans infection causing by Rhodotorula fungemia (Lunardi et al., 2006) requires more attentions on its broad interactions with soil fungal taxa in agroecosystems (Table S1). Pathogenic Periconia is one of the major threats to the crop cultivation due to potential production of peritoxins, which could cause rotting and blacken of the roots and stem base of crops (Churchill et al., 2001). More seriously, in this study, this plant pathogen was positively and significantly correlated with Rhodotorula in N fertilizations that probably strengthened their pathogenicity with detriment to crop production and quality. Similar synergetic relationship was found between species of Gibberella and Microascaceae (Fig. 4), and the latter have been recognized as plant pathogenic fungi that might cause brown spots on plant leaves with progressive leaf lesions (Mirzaee et al., 2010). Additionally, some potential plant pathogens in N fertilizations presented competitions with plant-beneficial fungi, such as Mortierella (Table S1), which probably reduced the Mortierella abundance (Fig. S5) and inhibited its benefit on crop nutrient uptake. Nevertheless, the interactions and functions of responding taxa need to be further examined by more methods, such as metagenome or metatranscriptome sequencing or functional assays with isolated strains to promote microbes that benefit soil fertility and crop health.
4.3 Indirect effects of N fertilizations on soybean yields via fungal community resistance
Sustainable agricultures need nutrient managements and improvements of microbial ecological services that finally lead to profitable crop yields (Qiao et al., 2019). Consistent with the result of a meta-analysis focusing on global measurements (Ye et al., 2020), long-term addition of N fertilizers was found to induce decrease in fungal diversity mainly due to the decrease in soil pH under N fertilizations, as soil pH was the dominating factor in fungal communities of this study (Fig. S1, Table 2). Soil acidification caused by addition of N fertilizers had been previously reported as the primary limitation to further enhancement of crop yield (Ning et al., 2020), which might be the reason that the soybean yields under N fertilizations was slightly lower than those under P, K fertilizations. The contents of soil TC, TN and NO3–-N were directly enhanced by addition of N fertilizers also significantly influenced fungal abundances, communities, network structure and abundances of responding taxa (Table 2, Fig. S9). Although increase of crop residual from above ground could enrich soil nutrients decomposed by microbes, smaller changes in above-mentioned soil properties under P, K fertilizations suggested that addition of N fertilizers changed fungal communities by direct input of the nutrients to soils. This was possibly due to slower nutrients turnover under N fertilizations with lower microbial decomposition activity and diversity (Fig. 1B). Past studies revealed that the increase of microbial diversity could promote crop yields (Wagg et al., 2011; Tautges et al., 2016), but no direct link was observed in this study examined by the SEM model (Fig. 5).
Nevertheless, it was interesting to find that fungal community resistance positively affected soybean yields examined by the SEM model (Fig. 5). The higher community resistance in P, K fertilization treatments facilitated uptake of nutrients by crops (Fig. 2, Fan et al., 2019), and led to advantageous conditions for crop growth and less competition among fungal species relative to N fertilizations (Fig. S7). The higher resistance of microbial communities could be linked to a higher level of functional redundancy, which maintain microbial diversity and function in time of environmental disturbances (Allison and Martiny, 2008). This might be another reason why fungal taxa respond less to P, K fertilizations. Additionally, RP and RPN mostly mediated by addition of N fertilizers directly and negatively affected fungal community resistance, which suggested their indirect negative influences on yield stability. Therefore, continuous addition of chemical N fertilizations rather than P, K fertilizations probably exerted detriment to the evolution of fungal communities, calling for more attentions on to the utility of N fertilizations for crop production and sustainable agroecosystems.
5 Conclusion
In summary, our study revealed that N fertilizations significantly increased fungal abundances, decreased fungal diversity and shifted fungal communities when compared to NoF. Contrarily, P, K fertilizations did not exert significant impact on fungal communities with lower community amplitude. Higher connectedness and node degree detected in network of P, K fertilizations relative to those of N fertilizations indicated that more harmonious and stable fungal communities were moderately affected by addition of P, K fertilizers. Compared to the solitary members of responding pathogens for P, K fertilizations, N fertilizations strongly increased the number of potential plant pathogens in the responding pool and their positive interactions, which probably accelerated the propagation of fungal pathogenicity across the whole microbial network community. These responding pathogens, mostly induced by N fertilizations, indirectly affected soybean yields by exerting negative influences on fungal community resistance. Our study highlighted that N fertilizations rather than P, K fertilizations induced more changeable, disharmonious and inefficient soil fungal communities. These findings are of great importance in selecting optimal fertilization strategies for sustainable agroecosystems.