Organic-Inorganic Fertilization Built Higher Stability of Soil And Root Microbial Networks Than Exclusive Mineral Or Organic Fertilization

Root microbiome is critical for plant health and performance. Many studies have assessed the impact of agricultural management on soil microbiome. But a comprehensive understanding of how root microbiota is affected by soil types and fertilization is still lacking. It is clear yet whether the stability of root microbiome is affected by fertilization regimes, and whether in the same patterns as soil microbiome. Methods We conducted a long-term experiment and investigated the impact of soil type, plant type and fertilization regimes on soil and root bacterial communities using high-throughput sequencing and network analysis. Results Our results indicated that microbial network under combined organic-inorganic fertilization had higher stability than exclusive inorganic or organic fertilizer. In addition, fertilization exhibited stronger effects on root microbiome than on soil microbiome. While total nitrogen mainly contributes to the variance of root microbiome, pH and soil organic matter were responsible for the differences of soil microbiome. Bacteroidetes and Firmicutes appeared as important drivers in soil and root microbiome amended with organic fertilizer, whereas Actinobacteria was enriched in the soil microbiome under inorganic fertilizer. Our results clearly indicated the responsive shifts of soil and root microbiome to different fertilization regimes, and gave hints for developing better fertilization practices and establishing healthy root associated microbiota.


Introduction
Root-associated microbiota play a key role in plant performance and productivity [1], often referred to as the "second genome" of plants [2]. By secreting root exudates, plants selectively recruit microbes in the rhizosphere, establishing resource-rich hotspots distinct from the bulk soil [3][4][5]. Though as important players in agro-ecosystems, few studies have assessed the impacts of agricultural practices on rhizosphere. Much of previous works focused exclusively on soil microbiome [6,7], rarely considering both soil and root. How root microbiome is affected by agricultural practices, and whether to the same extent as soil microbiome, still remain elusive.
While management-induced shifts in bulk soil microbiomes affect environmental outcomes, plantregulated rhizosphere communities are more directly relevant to yield outcomes. This knowledge can contribute to managing rhizosphere interactions that promote both plant productivity and agroecosystem sustainability [8]. Clearly understanding the responsive shifts of soil and root microbiome to different fertilization regimes is vital for developing better fertilization practices and for further improving soil fertility and function.
To disentangle the complexities of soil and root microbiome, a system-level understanding of community function and structure is needed [9]. In the past decades, network-based approaches have been found useful in unravelling microbe-microbe associations in complex environments, ranging from human guts to oceans and soils [10][11][12][13]. By investigating the co-occurrence patterns, network analysis could offer new insights into potential interactions and reveal niche spaces shared by community members [14,15].
However, previous studies of agricultural management mainly focused on network complexity [16], yet few have considered network stability, which is ecologically important.
Communities are considered more stable if more limited shifts occur in response to environmental perturbation; and are more likely to return to its previous state after the perturbation [17][18][19]. The gut microbiome, for instance, is often noted for its stability and the stability is considered critical for host health and well-being [20]. In natural systems, the stability of microbial networks decreased with increasing environmental stress [21]. Although presumably linked with community functions, such information about network stability in agro-ecosystems is lacking. The effect of agricultural management, especially fertilization regimes, on microbial networks is poorly understood.
With these ideas in mind, we conducted a large-scale experiment and investigated the impact of soil type, plant type and fertilization regimes on soil and root bacterial communities using amplicon sequencing and network analysis. We speci cally asked: (1) Do soil and root microbial communities differ in their responses to fertilization practices? (2) Which microbes are the indicator taxa for speci c fertilization regimes? (3) How do fertilization practices impact network stability of soil and root microbial networks respectively?

Sample collection
Samples were collected from 9 sites in China (Fig. 1, Additional le 1). The sampling sites spanned from the Northeast Semi-humid Plain to the Huanghuaihai Semi-humid Plain, representing 3 typical soil types in China, namely black soil, cinnamon soil, and uvo-aquic soil. The sampling sites were distributed in 3 provinces, i.e. Jilin (black soil), Beijing (cinnamon soil), and Henan ( uvo-aquic soil). In black soil (Northeast Semi-humid Plain), two plant types, spring maize (in open elds) and eggplants (in the greenhouse), were selected. All black soil sites were treated with inorganic N, P and K fertilizers (NPK). In uo-aquic soil (Huanghuaihai Semi-humid Plain), the elds were cropped with winter wheat (October to June) and summer maize (June to September) (i.e. wheat maize rotation). All uvo-aquic soil sites were fertilized with mineral NPK. In cinnamon soil (Huanghuaihai Semi-humid Plain), spring maize, wheat maize rotation, as well as fallow elds were also sampled. For spring maize, ve different kinds of fertilization were applied, including inorganic P and K fertilizers without N (PK), inorganic N, P, and K fertilizers (NPK), half substitution of the inorganic fertilizer by manure (1/2NPK + 1/2M), and equal substitution of the inorganic fertilizer with manure (M), respectively. For wheat maize rotation, ve fertilization regimes were applied, namely no fertilization (CK), NPK, inorganic fertilization plus straw Five soil cores (at 0-20 cm depth) were collected in each plot between plant rows and were pooled as one replicate of bulk soil. Each replicate contained around 80 g of bulk soil. About 10 g of each wellmixed bulk soil replicate were put into sterile falcons on ice, transferred to the lab immediately and stored at -80°C until DNA extraction. The rest bulk soil was used for physiochemical analysis.
In each sampled subplot, ve plants were selected. Plant roots were taken out from the soil and shacked to remove the loosely-attached bulk soil. The remaining soils that attached to plant roots including that need to be brushed off from plant rhizoplane were sampled and mixed as rhizosphere. The rhizosphere samples were treated in the same way as the bulk soil described above and were stored at -80°C until DNA extraction. The rest rhizosphere was used for physiochemical analysis.

Physiochemical analysis
The samples were sieved at 2 mm and kept at 4°C until analysis. Samples were analyzed for water content, pH, total nitrogen (TN), NO 3 -N, NH 4 -N and soil organic matter (SOM). The properties were determined according to the methods described in previous studies [22].
DNA extraction, PCR, library preparation, and sequencing 0.3 g from each soil/rhizosphere sample was used for DNA extraction using Gri ths' protocol [23]. DNA was ampli ed using the PCR primer pair 515F and 806R targeting the variable region 4 of 16S rRNA gene [24]. PCR reactions were carried out in 30 µL reactions containing 15 µL of Phusion® High-Fidelity PCR Master Mix (New England Biolabs, UK), 0.2 µM of forward and reverse primers, about 10 ng template DNA and the remaining volume sterile distilled water. Thermal cycling consisted of initial denaturation at 98 ℃ for 1 min, followed by 30 cycles of denaturation at 98 ℃ for 10 s, annealing at 50 ℃ for 30 s, and elongation at 72 ℃ for 30 s, followed by 72 ℃ for 5 min. PCR products were checked on 2% agarose gel. Triplicate PCR products were then puri ed with GeneJET TM Gel Extraction Kit (Thermo Scienti c, US) according to the manufacturer's instructions. Libraries were generated using Ion Plus Fragment Library Kit 48 rxns (Thermo Scienti c, US) following the manufacturer's recommendations. The library was quality-assessed and quanti ed on the Qubit@ 2.0 Fluorometer (Thermo Scienti c, US). All libraries were pooled into equal concentrations. The library was then sequenced on the Ion S5 TM XL platform (Thermo Scienti c, US) and single-end reads were generated.

Data analysis
All statistical analyses were conducted in R (v4.0.2). A work ow of the analysis steps presented below and the gures generated from each step is given in Fig. S1 (Additional le 3: Fig.S1). Alpha diversity Estimates of α-diversity (Shannon index) were calculated at each rarefaction level in usearch. We tested the effects of compartment, soil type, plant type and fertilization in overall samples and in subset samples. The normality of the dataset was checked using Shapiro-Wilk test and the homogeneity of variance across groups was computed using Levene's test. For the two-group comparison, the differences were tested using Student's t test if the dataset is normally distributed, or Wilcoxon test otherwise. For comparison of more than two groups, the differences were tested using one-way ANOVA if the samples have equal variance, or Kruskal-Wallis test otherwise. Tukey's Honest Signi cant Differences test was carried out for pair-wise comparison using the R package TukeyC [33] if applicable. Beta diversity We conducted a general analysis of β-diversity on the bacterial communities with all the samples together and then performed more speci c hypothesis testing. For the general analysis, we normalized the ltered OTU sequence counts using the "trimmed means of M" (TMM) method with the BioConductor package edgeR [34] and expressed the normalized counts as relative abundance counts per million (CPM). We then carried out unconstrained principle coordinates analysis (PCoA) on Bray-Curtis dissimilarities to quantify the major variance components of β-diversity. For in-depth analysis, we performed constrained analysis of principal coordinates (CAP). All ordination analyses were performed using the R package phyloseq [35].
The community dissimilarity was tested with permutational analysis of variance (PERMANOVA) and permutational analysis of multivariate dispersions (BETADISP) using the functions adonis and betadisp, respectively, in the vegan package [36] with 10 4 permutations. Where applicable, pairwise differences between the groups were assessed with the function pairwise.perm.manova from the package RVAideMemoire [37]. Statistical signi cance of the CAP was assessed using the permutest function in the vegan package with 10 4 permutations.

Microbiome network construction and analysis
OTUs with the relative abundance no less than 0.05% in at least one third samples were selected. Cooccurrence networks were constructed using Spearman rank correlations from R package psych [38]. Signi cant correlations (r > 0.7 and FDR adjusted p < 0.001) were visualized using Gephi [39]with the Fruchterman-Reingold layout.
We then quanti ed two network properties that have been associated with stability of ecological communities in perturbation studies: (1) how compartmentalized the network is, and (2) the number and strength of positive/negative correlations, via cohesion [40] and modularity [21] analyses, respectively. Modularity Modules (groups of taxa whose abundances are more correlated/anti-correlated with each other than the rest of the community) were identi ed using the Clauset-Newman-Moore algorithm (greedy_modularity_communities) from the Python package networkx [41]. We then calculated modularity, a measure of whether connections tend to occur within or between modules, using the quality function in the Python package networkx [41]. We calculated one value of modularity for each microbial network under different fertilization regimes.
Large positive modularity values (i.e., close to 1) indicate that more connections occur within, rather than between, modules compared to random chance. Communities with high modularity tend to be more stable, as the impact of losing a taxon is restricted to its own module, thus preventing the effects of that taxon's extinction from propagating to affect the rest of the network. Cohesion Cohesion is an abundance-weighted, null model-corrected metric based on pairwise correlations across taxa [40]. We used the author-recommended 'taxa shu e' null model with provided R code to calculate both positive and negative cohesions for each microbial community. The proportion of negative to positive cohesion was calculated as the absolute value of negative : positive cohesion.
By characterizing positive and negative co-occurrences separately, cohesion provides insights into associations among taxa caused by both positive and negative species interactions and/or by both similarity and differences in the niches of microbial taxa [42]. The ratio of absolute value of negative : positive cohesion indicates whether negative interaction or positive interaction dominant in the co cooccurrence networks. Identi cation of fertilization sensitive OTUs (fsOTUs) Complementary approaches were adopted to identify the OTUs differ under varied fertilization regimes.
We rst carried out correlation based indicator species analysis with the R package indicspecies [43] to calculate the point-biserial correlation coe cient (r) of an OTU's positive association to one or a combination of fertilization practice. The analysis was conducted with 10 4 permutations and considered signi cant at p < 0.05. Additionally, we tested for differential OTUs among fertilization regimes using likelihood ratio tests (LRT) with the R package edgeR [34]. OTUs differed in abundance among fertilization regimes (false discovery rate (FDR) corrected p value < 0.05) were identi ed. We then de ned OTUs that were con rmed by both indicator species analysis and LRT as fertilization sensitive OTUs (fsOTUs).

Bipartite networks
The fsOTUs were visualized using bipartite networks. The networks were constructed using the Fruchterman-Reingold layout with 10 4 permutations as implemented in the R package igraph [44].

Identi cation of key drivers in networks
The soil and root communities under each fertilization regime were combined to construct metanetworks. Four meta-networks were constructed consequently, in accordance with the fertilization regime NPK, M, NPK + M, and NPK + 1.5M.
We used a pair of parameters (i.e., within-module connectivity (Zi) and connectivity among modules (Pi)) [45] to describe the topological roles of individual nodes (OTUs We also used the method NetShift [48] to identify important microbial taxa which serve as "drivers" between two networks (https://web.rniapps.net/netshift). This method allows one to quantify the changes in association of a single taxon and compute the relative increase in importance of a node and thereby predict key microbial taxa.

Results
The main drivers of soil and root microbiota We conducted bacterial community pro ling of 51 soil and 44 root samples from 9 sites with 3 different soil types, 4 plant types and 8 fertilization regimes (Fig. 1). A total of 2 758 622 high-quality sequences was yielded (range 8 165 − 57 939; median 28 701; Additional le 2). In sum, we identi ed 18 097 bacterial zOTUs (zero radius OTU) across all samples.
The major phyla of the bacterial community were Proteobacteria, Bacteroidetes, Actinobacteria, Acidobacteria, and Firmicutes (Fig. S2). Soil and root, as different microbial habitats, were found inhabited by speci c sets of microbes (Fig. 2a). Principal coordinate analysis (PCoA) indicated that microbial communities were clearly separated by soil types (Fig. 2b; Table S1). The discrete outlier in the bacterial communities was consistent with TN ( Fig. S3b) and SOM (Fig. S3c). Soils supported higher species richness than roots ( Fig. S4 and Table S2).
For in-depth analysis, we employed canonical analysis of principal coordinates (CAP). Soil types could explain 13% of the variance in the soil microbiome and 15% of the root microbiome, both of which were con rmed by pairwise PERMANOVA test ( Fig. S5; Table S3). A higher diversity of soil and root microbiome was found in cinnamon soil than in uvo-aquic soil by the comparison of shannon index. Unexpectedly, no differences were observed in the α-diversity between black and cinnamon soil, though black soil is considered more fertile in general( Fig. S5; Table S4).
Soil and root microbiome were clearly separated by plant types, con rmed by both PCoA plots and PERMOVA tests ( Fig. S6; Table S5). However, no statistical differences of α-diversity were found between different plant types ( Fig. S7; Table S6).
We further investigated the fertilization impacts on soil and root bacterial communities. Clear differences of beta diversity among fertilization regimes were indicated by CAP and PERMANOVA tests ( Fig. S8; Fig.  3). Notably, fertilization explained a higher proportion of variation in root microbiome (74%) than in soil microbiome (53%), indicating a stronger impact of fertilization exerted on root bacterial community.
The impact of fertilization on α-diversity was only observed in root microbiome but not in soil microbiome, regardless with crop types (Fig. S8; Table S7). Not surprisingly, the lowest α-diversity of root microbiome was observed under fertilization regime PK, where N is missing. Unexpectedly, root microbiome without fertilization and fertilized with NPK did not differ in the diversity. Interestingly, the addition of organic materials (straw and manure) signi cantly lowered the α-diversity. In particular, the addition of straw (NPK + S) showed the lowest α-diversity, much lower than the addition of manure (NPK + M and NPK + 1.5M).

Fertilizatiton sensitive OTUs
To identify OTUs varied in abundance among different fertilization regimes, we employed indicator species analysis based on correlation. We further validated them using likelihood ratio tests implemented in edgeR. Finally, we de ned the OTUs supported by both methods as fertilization sensitive OTUs (hereafter: fsOTUs) and summarized them in bipartite networks ( Fig. 3; Fig. S9). The patterns were reminiscent of the effects seen in the beta diversity analyses. Each fertilization regime supports a specialized subset of soil and root bacteria. Particularly, we noted that Actinobacteria was enriched in the soil microbiome under inorganic NPK fertilizer. Instead, Acidobacteria and Bacteroidetes largely dominated the root microbiome, and were enriched with the addition of organic fertilizer. Firmicutes was enriched in under organic fertilization as well. Around 28% fsOTUs identi ed in root belonged to Acidobacteria, while only half of fsOTUs in soil (14%) were assigned to Acidobacteria. As approximation for an "effect size" of fertilization on microbial communities, the fsOTUs accounted for 9.8% and 14.1% of the total soil and root bacterial OTUs, respectively.

Network properties under different fertilization regimes
We constructed co-occurrence networks of bacterial community under each fertilization regime. The overall community taxonomy changed by fertilization practices (Fig. 4). We further quanti ed the network stability via modularity, a re ection of how compartmentalized the network is, and cohesion, a metric quantifying the degree of community complexity, respectively.
In both soil and root microbiome, combined fertilization leads to the highest modularity, indicating higher community stability (Fig. 5). Organic fertilizer resulted in higher soil network modularity than inorganic fertilizer. However, the opposite trend was observed in the root microbiome, where bacterial networks showed higher modularity with inorganic fertilizer than organic fertilizer.
Similarly, both highest negative and positive cohesion metrics were observed with combined fertilization (p < 0.01, ANOVA), indicating higher network stability (Fig. 5). However, no statistical differences of cohesion were found between inorganic and organic fertilization. As to the soil microbiome, the differences among fertilization regimes were only found with negative cohesion but not positive cohesion (p < 0.05, ANOVA). As like in the root microbiome, combined fertilization showed higher negative cohesion than inorganic NPK or organic manure alone, but no differences were found between organic and inorganic fertilization.
We further carried out canonical correspondence analysis (CCA) to investigate the environmental variables corresponding with fertilization practices. SOM and pH were found signi cantly correlated with soil microbiome, while TN was found signi cantly correlated with root microbiome (Fig. S11).

Key drivers in network shifting
The soil and root microbiota under each fertilization regime were combined to construct meta-networks.
Consequently, we obtained four meta-networks in accordance with the fertilization regime NPK, M, NPK + M, and NPK + 1.5M. On the basis of their within-module connectivity (Zi) and among-module connectivity (Pi), we identi ed a series of module hubs (nodes highly connected to other members in a module) and connectors (nodes linking different modules), which could be regarded as keystone nodes that play key roles in shaping network structure (Fig. 6). The number of module hubs were highest under combined fertilization and lowest under inorganic fertilization, which is in line with the results of modularity analysis. Under inorganic fertilization, nearly one fourth of the connector OTUs with NPK were assigned to Acidobacteria, which was less abundant in connector OTUs of other networks. Instead, Firmicutes and Candidatus_Saccharibacteria became prominent as connectors under organic and combined fertilization, whereas they were absent in the connector OTUs with inorganic fertilizer.
We further explored the potential "driver microbes" in shaping microbial networks under different fertilization regimes using the newly-developed method "NetShift". A taxon with an altered set of associations (identi ed by a high Neighbor shift (NESH) score), while being increasingly important for the whole network (identi ed by a positive delta betweenness (ΔB) score) is predicted as a "driver". Accordingly, we selected top 30 taxa of highest NESH score with positive ΔB values (Fig. 7). In the shift from inorganic (NPK) to organic (M) fertilization, Bacteroidates and Verrumicrobiota stand out as the most prominent drivers, as both of their NESH and ΔB score are high. In comparison with pure organic fertilization (M), Bacteroidates, Acidobacteria, Firmicutes, BRC1, and Gammaproteobacteria (Pseudomonadales and Xanthomonadales) contributed as important members in driving network changes under combined fertilization (NPK + M). In shift of fertilization regime from NPK + M to NPK + 1.5M, Turicibacter and Bacillus from Firmicutes were identi ed as the key drivers with the highest NESH and ΔB score. Besides, increased number of Proteobacteria, particularly Rhizobiales, were found among the driver taxa.

Discussion
In this study, we characterized soil and root microbiome from 3 different soil types, 4 plant types and 8 fertilization regimes in a long-term eld experiment. While soil types could largely determine microbial communities, fertilization practices were found as a primary factor in shaping soil and root microbiota under the same soil type. Our results suggested that combined organic-inorganic fertilization built higher stability of both soil and root microbial networks than exclusive inorganic or organic fertilization based on the analysis of network properties.
Network properties have been used to successfully predict the stability of microbial networks [11,21,42].
In particular, communities with greater modularity, reduced positive associations among taxa, and greater negative associations among taxa are more stable. Modularity could re ect biological processes such as shared ecological functions among taxa in a module [49][50][51], spatial compartmentalization [13], or similar niche requirement [52,53]. High modularity could stabilize communities by restricting the impact of losing a taxon to its own module [54,55]. Positive cohesion (positive relationships) represent high niche overlap and/or positive interactions between taxa, while negative cohesion (negative relationships) indicate divergent niches and/or negative interactions [40,42]. It is argued that positive associations can create dependency and mutual downfall [17]. In contrast, negative co-occurrences/interactions could dampen positive feedbacks and thus improve stability [21]. In this study, we found higher modularity and connectivity (i.e. cohesion) as well as a dominance of negative correlations in microbial networks under combined fertilization, indicating that combined fertilization leads to microbial community with higher network stability.
A recent study provided evidence that naturally-occurring microbiome display properties characteristic of unstable communities when under persistent stress [21]. Networks with higher stability are more robust to environmental perturbations [56]. In this sense, our ndings may indicate that the microbiota under combined fertilization is more resilient to environmental stresses. In addition, the success of pathogen invasion in the rhizosphere was reported to depend on the network structure of resident bacterial communities [57]. Therefore, the structure and stability of root community are highly important for plant health and tness. Indeed, our previous results showed combined fertilization resulted higher crop yields than exclusive manure application than solely mineral fertilization [22]. The yield increase by combined fertilization was also con rmed in other long-term experiments [58], with enhanced soil nutrient availability, microbial biomass, enzymatic activities and soil nitrogen processes [59]. In a four-decade nutrient fertilization experiment, the application of combined inorganic fertilizers and cow manure led to the most resistant microbial community, which was associated with the lowest relative abundance of potential fungal plant pathogens after 35 years of nutrient fertilization [60]. In brief, our results are in line with the notion that host can bene t from increased microbiome stability.
Interestingly, the in uence of fertilization is stronger on root microbiome, but less signi cant on soil microbiome, indicating compartment-speci c responses of bacterial community. Our CAA analysis revealed that TN is the environmental factor responsible for the community variation of root microbiome, whereas pH and SOM explained the soil community differences. While pH is well-known for its decisive role in selecting bacterial community, it seems that C and N factors drive the soil and root microbiome  [70]. Cultured members of these phyla have genomic pathways for the breakdown of complex, plant-derived polysaccharides. Considered together, they may play an important role in decomposing complex organic material, and thereby contribute to the community shift from inorganic to organic fertilization. With extra manure applied (i.e. NPK + 1.5M), a dominance of Firmicutes was observed as the potential drivers. Firmicutes are likely to increase in nutrient-rich conditions [71]. For instance, the increase in Firmicutes in gut microbiota was often correlated with obesity [72]. In our case, the large dominance of Firmicutes under the fertilization regime NPK + M might be an indication of overfertilization. Fertilization regimes that relies exclusively based on inorganic inputs may result in root selection of microbial communities more dependent on easily accessed C and disrupt the plant's ability to select for a prokaryotic community that mineralizes nutrients from existing organic matter.

Conclusions
Overall, we found that fertilization regimes had strong impacts on soil and root microbiome. Network analysis with modularity and cohesion indicated that microbial network under combined fertilization had higher stability than inorganic or organic fertilizer alone. In addition, the response of root microbiome to fertilization is stronger than soil microbiome and exhibited different patterns. While TN contributes mostly to the variance of root microbiome, pH and SOM could largely explain the differences in soil microbiome. Bacteroidetes and Firmicutes appeared as important drivers in soil and root microbiome amended with organic fertilizer, while Actinobacteria was enriched in the soil microbiome under inorganic NPK fertilizer alone. Our study imply that combined organic-inorganic fertilization might be a sound practice better than exclusive mineral or organic fertilization. However, the risk of over-fertilization still need to be taken care of.

Declarations
Availability of data and materials The clean reads were deposited under the accession number PRJCA004095 in the GSA database (https://bigd.big.ac.cn/gsa/). The scripts were deposited at https://github.com/yysmile2014/Fertilization.

Figure 1
Geographical distributions of the 9 sampling sites in China. The green shaded area represents croplands.
The samples were collected in 2 typical agro-climatic area with 3 typical soils in China, namely black soil, cinnamon soil, and uvo-aquic soil, represented by different colors. The sampling sites were distributed in 3 provinces, Jilin (a), Beijing (b), and Henan (c), which were illustrated in detail in the right panels. The plant types were indicated by different shapes. Fertilization were indicated with I and II in Beijing, and other sites were under NPK treatment.      The top 30 taxa of "driver microbes" in microbial networks during the shift of fertilization regime. Each dot represents a taxa in the microbial networks. The X axis denotes the delta between score (ΔB), implying the changes of importance of each taxa in the network in comparison with the former network.
The size of the dot corresponds to the NESH score, indicating the changes of node associations. Taxa with high NESH score and positive delta between value were predicted as "driver microbes".