Long-Term Organic–Inorganic Fertilization Regimes Alter Bacterial and Fungal Communities and Rice Yields in Paddy Soil

Microorganisms are the most abundant and diverse organisms in soils and have important effects on soil fertility. In this study, effects of the long-term fertilization treatments no fertilizer (CK), chemical fertilizer (nitrogen–phosphorus–potassium (NPK)), and organic–inorganic fertilizer (NPK and organic fertilizer (NPKM)) on rice yield and soil bacterial and fungal community diversity, structure, composition, and interaction networks were evaluated. Of the three treatments, the highest rice yield was in NPKM. Bacterial richness was significantly higher in NPKM than in NPK. Fertilization treatment significantly altered β diversity of communities, species composition of bacterial and fungal communities, and structure of soil microbial networks. The most complex bacterial and fungal interaction co-occurrence network with the highest average degree and numbers of edges and nodes was in NPKM. Relative abundance of the plant growth-promoting fungus Trichoderma increased significantly in NPKM compared with CK and NPK. The results of the study indicate that bacterial richness and microbial community member interactions (network complexity) might be suitable indicators of soil biological fertility. This research provides new insights on the effects of different fertilization regimes on responses of soil bacterial and fungal communities and their contributions to crop yield. New indicators such as bacterial richness and complexity of microbial interaction networks are also identified that can be used to evaluate soil biological fertility.


INTRODUCTION
Chemical fertilizer is often used to compensate the loss of soil nutrients during growing stages of crops in which chemical fertility is the primary focus. By contrast, organic fertilizer is typically applied as a base fertilizer before crop growth. Many studies indicate that organic fertilizer plus chemical fertilizer treatment can increase yields compared with only chemical fertilization, even when total nitrogen (N), phosphorus (P), and potassium (K) contents are the same (Pan et al., 2009;Xie et al., 2016;Zhao et al., 2016). Addition of organic fertilizer with chemical fertilizer can slow the release of nutrients in the chemical fertilizer and thus reduce N loss and increase N use efficiency (Pan et al., 2009). The combination of fertilizer types can also optimize soil aggregate structure by increasing organic matter and thus improve soil physical structure (Bronick and Lal, 2005). However, how addition of organic fertilizers influences soil biological fertility has not been determined.
Microorganisms are the most abundant and diverse soil organisms. On average, 1 g soil contains 10 10 cultivable cells and approximately 10 4 microbial species (Liu et al., 2012;Saleem et al., 2019). As the most active component of an ecosystem, microorganisms have essential roles in energy flow, particularly in decomposition of organic matter, and biogeochemical cycling of nutrient elements (Morris and Blackwood, 2015). Therefore, growth and metabolism of microbial communities can directly and indirectly affect crop growth (Ahn et al., 2012). The most significant effects of soil microorganisms are on the decomposition of organic matter and the release of minerals. Soil microorganisms regulate decomposition of organic matter and release of minerals and are essential in forming humus and improving structure and cultivability of soil, which ultimately promote crop growth. Simultaneously, crops secrete biosynthesized organic matter and other substances into soils through the root system to provide soil microorganisms with sources of nutrients such as carbon (C) and N. Those root exudates can enrich and increase the growth of bacterial flora that are beneficial to plant growth and help resist infection by pathogenic bacteria (Hou et al., 2021). Thus, soil microbial activities and crop growth and yield are interdependent, and it is therefore important to consider activities of microorganisms when conducting research on crop yields in agroecosystems.
In agroecosystems, soil microbial biomass and diversity are potential indicators of soil quality (Bending et al., 2004), which is associated with soil productivity and crop growth and yield and can reflect soil vitality, health status, and biological fertility. Soil biological fertility is the core component of soil fertility, the key to sustainable use of agricultural land, and the foundation of a future agricultural revolution (Hatfield and Walthall, 2015). As decomposers and material cyclers of soil ecosystems, soil microorganisms are the core of soil biological fertility and are essential in regulating biological fertility. However, microbial community indicators that can be used to evaluate soil biological fertility have not been determined. To date, most studies have focused on effects of different fertilization treatments on soil bacterial communities (Gu et al., 2009;Zhao et al., 2014Zhao et al., , 2016, whereas effects on soil fungal communities and their functions have largely been ignored. Many studies have examined how different fertilization regimes shape soil microbial community diversity, structure, and composition (Ge et al., 2008;Gu et al., 2009;van der Bom et al., 2018). However, few studies focus on interactions among community members. Community member interactions are largely dependent on nutrient and energy supplies in a soil (Qiu et al., 2021) and therefore might be an important indicator of soil biological fertility.
In this study, a long-term field fertilization experiment was set up in paddy field that included three treatments: no fertilization, chemical fertilization, and chemical fertilization plus organic fertilizer addition. In addition to determining rice yields, high-throughput sequencing was used to investigate diversity, structure, composition, and interaction networks of soil bacterial and fungal communities. The aims of the study were to understand how different fertilization regimes influence soil microbial communities and rice yields and to determine which microbial community parameters might be used as indicators of soil biological fertility.

Site Description
The study was conducted in a long-term experimental field site in Changshu, Jiangsu Province, China (31°18′N, 120°37′E). The site is at the center of the Tai Lake plain region, where the cropping regime is a rotation of summer paddy rice and winter wheat. The climate is humid subtropical monsoon with average annual rainfall of 1,063 mm and annual mean minimum and maximum temperatures of 3.1°C and 33°C, respectively (Chen et al., 2016;Wang et al., 2016a). The field was tilled to an average depth of 20 cm before either sowing wheat or transplanting rice seedlings. Rice plots were flooded with 5 cm of standing water from July to September (Wang et al., 2016b). Rice was transplanted in June using two seedlings per hill at 13 cm × 28 cm spacing (Wang et al., 2016a). Following rice harvest, wheat was sown with seeds at 150 kg ha −1 in October or November every year. A sickle was used to manually harvest crops at ground level, and aboveground biomass was removed from plots. Fertilizers were applied as basal fertilizer after harvest of both rice and wheat (Wang et al., 2016b).

Field Experiment and Soil Sampling
Three treatments with four replicates were established in a randomized block design in 2005, and each replicate plot was 16 m 2 (4 × 4 m). To avoid margin effects, yield was measured in a 4-m 2 (2 m × 2 m) area only in the middle of each plot, with yield converted to yield per hectare. Experimental treatments included CK (no fertilizer control), NPK (240 kg N ha −1 , 90 kg P 2 O 5 ha −1 , and 120 kg K 2 O ha −1 ), and NPKM (The NPK and NPKM treatments contained the same total amounts of nutrients. NPKM contained 4,500 kg ha −1 organic fertilizer). The organic fertilizer was derived from composted pig manure with rice straw and contained 26.4% organic C, 2.5% total N, 1.6% P (P 2 O 5 ), and 1.3% K (K 2 O).
Soil samples were collected at 0 to 20 cm after rice harvest in October 2015. Four soil cores (5-cm diameter) were randomly collected in each 4 × 4-m plot. Soil samples from different treatments were mixed separately and sieved (2 mm) to remove plant materials, roots, and stones. To minimize changes in soil communities after sampling, DNA was immediately extracted from fresh soil. Subsamples were frozen at −80°C.

DNA Extraction
Soil DNA was extracted from 0.25 g of fresh soil with a MoBio Powersoil™ DNA Isolation Kit (MoBio Laboratories, Carlsbad, CA, United States of America) using the bead-based homogenizer protocol according to the manufacturer's instructions. Quantity and quality of DNA extracts were assayed by a Nanodrop ND-2000 UV-VIS Spectrophotometer (NanoDrop Technologies, Wilmington, DE, United States of America), and DNA was stored at −80°C until analysis.
Amplification of the fungal ITS1 (Internal Transcribed Spacer I) region was performed using a PCR reaction solution containing 12.5 μl of Master Mix (Qiagen Inc.), 0.5 μl (10 mM) of ITS1F primer (5′-CTTGGTCATTTAGAGGAAGTAA-3′; Gardes and Bruns, 1993), 0.5 μl (10 mM) of ITS2 primer (5′-GCTGCG TTCTTCATC-GATGC-3′; Baldwin, 1992), 1 μl of DNA template (10 ng μl −1 ), and 10.5 μl of ddH 2 O to a final volume of 25 μl. The PCR protocol was performed in triplicate using the following conditions: 10 min at 95°C for initial denaturing, followed by 35 cycles of 95°C for 15 s, 56°C for 15 s, and 72°C for 30 s, with a final extension at 72°C for 5 min. The Illumina sequencing adapter-ligated reverse primer contained a 6-bp bar code specific to each sample for identification (Caporaso et al., 2012). After amplification, triplicate PCR products were pooled and purified using a PCR Cleanup Kit (Axygen Biosciences, Union City, CA, United States of America). Bacterial and fungal PCR products were pooled separately for sequencing. Sequencing was performed on a single lane of an Illumina MiSeq platform at Personal Biotechnology Co., Ltd. (Shanghai, China). All sequence data have been deposited in the NCBI (National Center for Biotechnology Information) Sequence Read Archive database. The accession numbers are SRP359136 for bacterial data and SRP359138 for fungal data.

Bioinformatics Analysis and Statistics of High-Throughput Sequencing Data
High-throughput sequencing data were processed using USEARCH v.10.0 (Edgar, 2010), VSEARCH v.2.13 (Rognes et al., 2016), and in-house scripts . First, the --fastq_mergepairs command was used to merge paired-end sequences of the sequencing data and rename them. Then, the --fastx_filter command was used to remove the double-ended primers and bar codes and to perform quality control to make the error rate less than 1%. The --derep_fulllength command was used to reduce sequence redundancy. Redundant sequences were clustered into operational taxonomic units (OTU) with 97% similarity by using the -cluster_otus command, and chimeras were removed simultaneously. UPASE (Edgar, 2013) was used to select representative sequences, and then an OTU table was generated by the --usearch_global command. Species annotations were conducted using the -sintax command of VSEARCH with the SILVA database (Quast et al., 2013) for bacteria and the UNITE database (Abarenkov et al., 2010) for fungi.
Bacterial and fungal sequences were flattened to 3,935 and 5,837 sequences per sample, respectively, before α diversity indices were calculated. Chao1, Pielou's evenness, and Shannon indices were calculated to evaluate community richness, evenness, and diversity (Lu et al., 2020), respectively. Principal coordinate analysis (PCoA; Zhang et al., 2019) was used to show differences in microbial community structure among samples. To test the significance of differences in microbial community structure among different treatments, adonis analysis based on permutational multivariate analysis of variance (PERMANOVA; Zhang et al., 2019) was used. One-way ANOVA followed by Tukey's post hoc test was performed to compare means of rice yield, α-diversity indices, microbial abundance, and the Proteobacteria to Acidobacteria ratio. Ecological networks of fungi and bacteria were constructed, and network topology parameters were calculated in R (https://www.r-project.org/) using the "igraph" package (Wang et al., 2018). Pearson correlation analyses were performed to examine relations between rice yield and α diversity indices. Mantel tests (Wang et al., 2016b) were conducted to test correlations between rice yield and community composition. All data visualization was performed in R using the "ggplot2" package. When considering bacteria and fungi together, OTU tables of fungi and bacteria were combined for calculations.

Effects of Different Long-Term Fertilization Treatments on Rice Yield
Rice yield varied among different fertilization treatments (Figure 1). Compared with CK, fertilization treatments significantly increased rice yield (Tukey's post hoc test, p < 0.05). Average rice yields were 5,513 kg/ha in CK, 9,329 kg/ha in NPK, and 10,316 kg/ha in NPKM. Although NKP and NPKM had equal inputs of N, P, and K, rice yield was significantly higher in NPKM than in NPK (Tukey's post hoc test, p < 0.05).

Effects of Different Fertilization Treatments on Microbial Alpha Diversity
Chao1, evenness, and Shannon indices were calculated to estimate microbial richness, evenness, and diversity, respectively (Figure 2). Compared with NPK, NPKM significantly increased bacterial richness (Tukey's post-hoc test, p < 0.05) but decreased bacterial evenness and diversity, although differences were not significant. However, different fertilization treatments did not significantly affect fungal richness, evenness, and diversity. In Pearson correlation analyses, rice yield was positively correlated with bacterial richness and fungal evenness but negatively correlated with bacterial evenness and diversity and fungal richness and diversity, although none of the correlations were significant (Supplementary Figure S1).  Frontiers in Microbiology | www.frontiersin.org

Effects of Different Fertilization Treatments on Microbial Beta Diversity
In the unconstrained PCoA of weighted UniFrac distance, soil bacterial and fungal communities formed three distinct clusters according to different fertilization treatments (PERMANOVA: p < 0.01; Figure 3). The NPKM community was significantly different from those in CK and NPK and separated along the first coordinate axis, which indicated that the greatest variation among treatments was most likely due to the addition of organic fertilizer. According to Mantel tests, rice yield was significantly positively correlated with bacterial community structure (r = 0.462, p = 0.005), whereas fungal community structure and rice yield were not significantly correlated (r = 0.033, p = 0.408).

Effects of Different Fertilization Treatments on Co-occurrence Networks of Soil Microorganisms
Nodes in the bacterial and fungal co-occurrence network included seven bacterial phyla (Proteobacteria, Actinobacteria, Chloroflexi, Firmicutes, Nitrospirae, Gemmatimonadetes, and Acidobacteria) and three fungal phyla (Ascomycota, Basidiomycota, and Mortierellomycota; Figure 6). A unique microbial co-occurrence network formed in each treatment. The most complex network was in NPKM, followed by those in NPK and CK. Compared with CK and NPK, the NPKM network had the highest connectivity (0.097), number of edges (114), number of vertices (49), and average degree (4.6; Figure 6).
To quantify differences in co-occurrence networks among treatments, numbers of common and unique edges of CK, NPK, and NPKM networks were calculated. There were few common edges shared by any two treatments, compared with many unique edges in each treatment. Networks in CK and NPK shared only four edges, whereas the CK network had 56 unique edges and the NPK network had 74. Networks in CK and NPKM shared only three edges, whereas the CK network had 57 unique edges and the NPKM network had 111. Networks in NPK and NPKM shared only six edges, but the NPK network had 72 unique edges and the NPKM network had 108 (Figure 6).
The fungal genus Trichoderma was absent in NPK but occurred in NPKM. The sequence number of Trichoderma was significantly higher in NPKM than in CK and NPK (Tukey's post hoc test, p < 0.05), with the lowest number in NPK (Figure 6).

DISCUSSION
Soil microbial biomass and diversity are potential indicators of soil quality (Bending et al., 2004), as well as important factors in maintaining integrity and stability of soil functions in agroecosystems (Kennedy and Smith, 1995;Nannipieri et al., 2003;Brussaard et al., 2007). They are associated with the level of soil biological fertility and thereby affect crop yields. In this study, in a long-term field fertilization experiment, effects of no fertilization, chemical fertilization, and chemical fertilization plus organic fertilizer treatments on rice yields and soil bacterial and fungal communities were evaluated. High-throughput sequencing was used to evaluate composition of microbial communities. The highest rice yield among the three treatments was in NPKM. Compared with NPK, combined application of organic and inorganic fertilizers significantly increased bacterial richness but decreased bacterial evenness. Those results suggested that application of organic fertilizer increased the number of species to increase bacterial richness and enriched certain groups that preferred organic fertilizer to decrease bacterial evenness. By contrast, fertilization treatments did not significantly affect indices of fungal alpha diversity. Thus, fertilizer treatments had different effects on alpha diversity of bacterial and fungal communities. The results are consistent with those of previous studies in which organic fertilizer additions increased soil microbial abundance but not evenness (Parham et al., 2003;Sun et al., 2004;Jangid et al., 2008;Yuan et al., 2008).
Fertilization can stimulate growth of soil microorganisms and thereby affect microbial community structure (Chu et al., 2007;Esperschütz et al., 2007;Gu et al., 2009). In previous studies, long-term application of different fertilizers significantly altered the structure of soil bacterial and fungal communities (Allison et al., 2007;Wu et al., 2011;Yu et al., 2013). According to the PCoA in this study, the different fertilization treatments altered bacterial and fungal species composition and therefore microbial community structure. Mantel tests indicated that rice yield was positively correlated with bacterial community structure but not with fungal community structure. Therefore, because the response of bacterial communities to fertilization treatments was greater than that of fungal communities, bacterial community parameters might be suitable indicators of soil biological fertility.
Bacteria can be divided into two life types according to life history strategy, namely copiotroph r-strategists and oligotroph K-strategists (Pianka, 1970;Fierer et al., 2007). When soil organic matter content is high, r-strategists are usually the primary decomposers of organic matter, and microorganisms in the eutrophic group are most abundant. By contrast, when soil organic matter content is low and nutrients are lacking, K-strategists in the oligotrophic group have a competitive advantage (Fontaine et al., 2003). Actinobacteria, α-Proteobacteria, β-Proteobacteria, and Bacteroidetes are generally regarded as r-strategists in the Frontiers in Microbiology | www.frontiersin.org eutrophic group and are generally more abundant in highfertility soils (Fierer et al., 2007;Newton and Mcmahon, 2011;Leff et al., 2015;Zhou et al., 2015). Acidobacteria and Verrucomicrobia are classified as K-strategists in the oligotrophic group (Ramirez et al., 2010). In a previous study, application of organic fertilizer significantly increased relative abundances of r-strategists in the eutrophic group, such as Proteobacteria, Bacteroidetes, and Actinobacteria, whereas application of inorganic fertilizers significantly increased relative abundances of K-strategists in the oligotrophic group, such as Acidobacteria (Xun et al., 2016). In this study, the three most abundant bacterial groups were Proteobacteria (39.5%), Acidobacteria (19.1%), and Actinobacteria (6.5%), which together accounted for 65.1% of the bacterial community. In NPKM, relative abundances of r-strategists in the eutrophic group (Actinobacteria and Proteobacteria) increased, but those of K-strategists in the oligotrophic group (Acidobacteria and Verrucomicrobia) decreased significantly (Figure 4). The ratio of abundances of Proteobacteria to Acidobacteria may be an indicator of nutrient status of a soil ecosystem (Smit et al., 2001;Torsvik and Øvreås, 2002). In this study, the ratio of Proteobacteria to Acidobacteria was significantly higher in NPKM than in other treatments (Figure 7), indicating that soil nutrient status and productivity were also high. However, it should be noted that considering microbial community response at finer resolutions (e.g., family, genus and species level) may be more adequate when assigning life strategies to microorganisms (Ho et al., 2017). Notably, abundance of the fungal genus Trichoderma increased significantly in NPKM (Tukey's post hoc test, p < 0.05; Figure 6E). Trichoderma is a well-known plant growthpromoting fungus (Masunaka et al., 2011) that can significantly increase crop growth and productivity (Hyakumachi, 1994). Thus, Trichoderma might be one reason for the high rice yield in the organic-inorganic combined treatment.
Co-occurrence network analysis was conducted to examine interactions among soil microorganisms in response to different fertilization treatments. The most complicated network was in NPKM, suggesting organic fertilizer addition increased interactions among microbial community members (Figure 6). The increase in complexity was likely due to the additional C input in the organic fertilizer. Soil bacteria and fungi can use those carbon sources to generate additional energy and thus boost expression of soil microbial functions, which increases the number of microbial interactions. The results suggested that organic fertilizer addition increased stability of community structure and therefore ability to resist external interference. In addition, by increasing levels of cooperation within microbial communities, addition of organic fertilizer was likely conducive to developing soil biological fertility. Thus, complexity of microbial interaction networks might be an indicator of soil biological fertility. Notably, the fungal genus Trichoderma was not in the co-occurrence network in NPK but was in that in NPKM. Enrichment of the beneficial fungus Trichoderma in the combination organic and inorganic treatment indicated there were changes in the interaction network of soil microorganisms that could increase rice productivity. Inoculation of Trichoderma could be used to improve crop productivity in agroecosystems.

CONCLUSION
Soil biological fertility is an important component of overall soil fertility. However, how to best evaluate soil biological fertility has not been determined. Soil microbial communities are important in regulating soil biological fertility. Thus, the responses of soil bacterial and fungal community diversity, structure, composition, and interactions to different fertilization treatments were analyzed by high-throughput sequencing. The inorganic and organic fertilizer treatment had the highest bacterial richness, the most unique bacterial and fungal communities due to species selection by fertilizers, and the microbial network with the highest complexity, and as a result, the highest productivity. Trichoderma was enriched in the NPKM treatment and might be a key contributor to the increase in soil fertility. The results indicate that bacterial richness and complexity of microbial interaction networks could be used as indicators of soil biological fertility. This research provides new insights on responses of soil bacterial and fungal communities to different fertilization treatments and their contributions to crop yield. The study also identifies new indicators to evaluate soil biological fertility and indicates that inoculation with Trichoderma might improve crop productivity in agroecosystems.

DATA AVAILABILITY STATEMENT
The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found in the article/Supplementary Material.

AUTHOR CONTRIBUTIONS
CX, QH, and QS designed this experiment. TM, XH, SC, and YL carried out the experiment. TM wrote this manuscript. CX and QS revised the manuscript. All authors contributed to the article and approved the submitted version.