Succession of Composition and Function of Soil Bacterial Communities During Key Rice Growth Stages

Elucidating the succession of soil microbial communities and microbial functions at key plant growth stages is a major goal of microbial ecology research. In this study, we investigated the succession of soil bacteria during four fertilizer treatments (control, NPK, NPK + pig manure, and NPK + straw) and at three crucial rice growth stages (tillering, heading, and ripening) in paddy soil from a rice-wheat cropping system over a 10-year period. The results showed that the bacterial community and function composition of the control treatment was significantly different from that of the other treatments with NPK fertilizers, and S1 from others stages (ANOSIM, p < 0.05). The application of pig manure could reduce the effects of applying NPK fertilizers on bacterial communities in heading and ripening stages, but the effects of straw returning is not obvious. Variance partitioning analyses (VPA) suggested that pH, OM, and AK appeared to be key factors responsible for the microbial community changes observed in all the treatments or stages. The correlation results showed the bacterial families different between S1 and other stages such as Micromonosporaceae, Nocardioidaceae, Gaiellaceae, and Anaerolineaceae etc., were correlated with bacterial KEGG metabolic pathways. In addition, the topological of the soil bacterial community network with more nodes, links and higher Maximal degree at the heading stage and maintained relatively similar topological structures at the heading and ripening stages. However, the topological of the functional networks at the ripening stage were a small yet complicated co-occurring network with 209 nodes, 789 links, higher Average connectivity (avgK), and Maximal degree. These results suggest an obvious succession of soil bacteria and bacterial function at the key rice growth stages, but the topological of functional network structure of bacteria changes a little in the early and middle stages of rice, while its changes significantly in the ripening stage of rice growth.

Introduction

The rice-wheat cropping system is the main grain production system in the Yangtze River Basin in China (Wang et al., 2016a), and soil acidification, continuous cropping obstacles and soil-borne diseases are the major issues affecting rice-wheat production in this area. Previous case studies showed that the application of organic amendments, such as pig manure and straw, could either improve the soil community structure (Blanchet et al., 2016) or alleviate soil acidification (Sun et al., 2015) and continuous cropping obstacles (Zhao et al., 2016). Fertilization can also significantly change the community structure of soil microorganisms. A proper understanding of the mutual relationship between the soil microbial structure and fertilizer management at different phenological stages of rice growth will aid the selection of appropriate and efficient fertilizers for different phenological stages of rice growth.

Due to the rapid responses of bacteria to changes in the soil environment, bacterial communities were selected as early bio-indicators of soil quality (Chen et al., 2017). Previous studies revealed that root exudates have a significant impact on the soil microbial community (Shi et al., 2013). In addition, the interactions between the rice root and soil microorganisms usually induce growth promotion and disease suppression (Edwards et al., 2015). But at different growth stages, crop root exudates commonly undergo dynamic changes (Li et al., 2016). Several studies using different methods for testing diversity reported that the soil bacterial and functional gene communities show obvious differences among different plant growth stages (Noll et al., 2005; Rui et al., 2009; Krause et al., 2010; Edwards et al., 2015; Li et al., 2016; Wang et al., 2016a). Soil microbial communities also play an important role in the maintenance of ecosystem function and sustainability (Paul and Clark, 1989). Therefore, it is essential to explore the internal dynamics among soil bacteria, the functional gene community, and the fertilizer and plant growth stage, which are important factors that influence the ecosystem function and sustainability of the cropping system.

Molecular ecology network analysis has been used as a powerful tool to understand soil microbial ecosystems and complex microbial assemblages (Zhou et al., 2011; Deng et al., 2012). The co-occurrence (positive or negative interaction) among OTUs/genera can be defined by a molecular ecological network obtained using random matrix theory-based methods (Zhou et al., 2011; Deng et al., 2012). These visual network structures provide detailed perspectives on microbial assemblages (Shi et al., 2016). The potential microbe-microbe interactions display significant changes among different growth stages of upland plants, and the complexity of the microbial network increases during the growth of the Avena fatua (Shi et al., 2016). In addition, based on the 16S rRNA gene sequences, phylogenetic investigation of communities by reconstruction of unobserved states (PICRUSt) has been widely used to predict the functional capabilities of microbial communities (Langille et al., 2013). In this study, we explored the soil microbial and microbial function networks of a rice-wheat cropping system at different rice growth stages.

Four fertilizer treatments (control, NPK, NPK + manure, and NPK + straw) were applied to a rice field and paddy soil with a distribution of three rice growth stages (tillering, heading, and ripening), and samples were collected for the analysis of the bacterial community structure. The objectives of this study were to (1) understand the effects of fertilizers and rice growth stages on the soil bacterial community structure and bacterial metabolic functions, (2) explore the key soil environmental variables that strongly affected the succession of the bacterial community during the paddy stages in a rice-wheat cropping system, and (3) understand the variations between the topological shift of the soil bacterial networks and that of the metabolic functional networks during the rice growth stages.

Materials and Methods

Field Description and Experimental Design

The experiments were performed in a paddy field in Changzhou, Jintan, China (31°39′N, 119°28′E), during rice-wheat rotation in 2005. The local climatic conditions is a northern subtropical monsoon with an annual precipitation of 1063 mm and an average annual temperature of 15.3°C. Based on the FAO soil classification system, the field soil is classified as stagnic anthrosols. A randomized block design with three replicates (each plot consisted of 8 m × 5 m) was used in the field experiment. Four common fertilization treatments were designed: A, control treatment without fertilizers; B, 100% NPK chemical fertilizer treatment (based on the conventional dosage for this region, namely, 300 kg/ha urea, 120 kg/ha P₂O₅, and 100 kg/ha K₂O); C, 50% NPK fertilizer treatment plus 6,000 kg/ha pig manure; and D, 100% NPK fertilizer treatment plus crop straw. All other fertilizers and crop straw were applied as basal fertilizers before planting, whereas N fertilizers were used as basal and supplementary fertilizers (basal fertilizers: tillering supplementary fertilizers: panicle supplementary fertilizers = 4:3:3).

Soil Sampling and Analysis of Soil Properties

In 2015, the soil was sampled at three rice-growing stages, namely, at the tillering stage in July (14 days after seed sowing), the booting stage in September (56 days after seed sowing) and the ripening stage in October (98 days after seed sowing). A randomized block design with three replicates (each plot consisted of 8 m × 5 m) was used in the field experiment. For each plot, Soil samples (diameter of 2.5 cm and depth of 0–20 cm) from seven to 11 randomly selected locations were pooled to obtain one biological replicate. After the removal of rice roots and stones, the pooled soil sample was placed in a sterile plastic bag in an ice box until transport to a laboratory for further sample storage or chemical analyses (Wang et al., 2016c).

Soil DNA Extraction and High-Throughput Sequencing

The total soil DNA from 0.5 g of soil was extracted using a fast PowerSoil DNA Isolation Kit (MP Biomedicals, Santa Ana, CA, United States) according to the manufacturer’s instructions. The bacterial V4-V5 hypervariable region of the 16S rRNA gene was amplified by PCR using the primers 515F (5′-CCTACGGGAGGCAGCAG-3′) and 907R (5′-TTACCGCGGCTGCTGGC-3′). PCR reactions were performed in triplicate 20 μL mixture containing 4 μL of 5 × FastPfu Buffer, 2 μL of 2.5 mM dNTPs, 0.8 μL of each primer (5 μM), 0.4 μL of FastPfu Polymerase (Transgen, China), and 10 ng of template DNA. The PCR products from the soil samples were sent out for pyrosequencing using an Illumina MiSeq PE300 (Illumina Inc., San Diego, CA, United States) at the Majorbio Bioinformatics Technology Co., Ltd. (Shanghai, China) (Wang et al., 2016c). The Trimmomatic (Bolger et al., 2014) and flash (Magoč and Salzberg, 2011) programmes were used to process the V4-V5 tagged sequences. The 16S rRNA operational taxonomic units (OTUs) with an identical cut-off of 0.97 were processed using the Usearch software platform (version 7.1)¹ (Edgar et al., 2011), and chimeric OTUs were then removed using the Uchime (version 4.2.40) method (Edgar, 2013). A subset of 47,847 sequences per sample was selected for downstream analyses. The raw sequences are available through the NCBI Sequence Read Archive (Accession No. SUB3055082).

Statistical Analyses

The values of richness estimators [coverage estimator (Ace) and abundance-based (Chao1)] and diversity indices (Shannon and Simpson) were determined using Mothur (Schloss et al., 2011). Principal coordinate analysis (PCoA), analysis of similarities (ANOSIM), redundancy analysis (RDA), Monte Carlo permutations and variance partitioning analyses (VPA) were performed using R (Version 3.1.2, vegan package). The bacterial metabolic function profiles (KEGG and COG) were generated using PICRUSt. Comparisons of the top 100 bacterial families and bacterial metabolic functions among the different groups (fertilizer treatments and growth stages) were performed using STAMP (Parks et al., 2014). The correlations between bacterial abundance (at the family level, forward selected by STAMP) and PICRUSt-generated KEGG metabolic pathway profiles were plotted using R (version 3.1.2, corrplot2 package). A network analysis was performed using the Molecular Ecological Network Analyses Pipeline² (Zhou et al., 2010, 2011; Deng et al., 2012). The “output for Cytoscape visualization” was run in the “greedy modularity optimization mode,” and the files for network graph visualization were generated using a Gephi interactive platform (Bastian et al., 2009). The other statistical analyses were conducted using SPSS 13.0.

Results

Soil Physicochemical Characteristics

The soil organic material (OM), total K (TK), available N (AN), available P (AP), and pH were significantly different (p < 0.05) among the four fertilizer treatments at all three rice growth stages (Tables 1, 2). In addition, the soil pH, available K, and OM showed obvious increase or decrease (p < 0.05) at the three different stages (Table 1). The soil pH, TP, and AN were more likely to show changes among different fertilization treatments and rice growth stages than the other soil physicochemical characteristics (Table 2, Two-way ANOVA, p < 0.001). At all stages, particularly the ripening stage (S3), the control without fertilizers (A) had a higher soil pH but lower contents of OM, TP, TK, and AP compared with the other three fertilization treatments. Additionally, the addition of pig manure to the NPK treatment (C) significantly increased (p < 0.05) the OM and AP contents.

TABLE 1

Table 1. Chemical properties at different soil stages.

TABLE 2

Table 2. Effects of the sampling time and fertilization treatment on soil physicochemical properties (two-way ANOVA, ^∗p < 0.05. ^∗∗p < 0.01. ^∗∗∗p < 0.001).

Bacterial α-Diversity

Bacterial richness and diversity indices were calculated based on randomly selected subset of 47,847 sequences per soil sample. Significant differences (p < 0.05) in the Ace, Chao1, Shannon, and Simpson indices were found among the different sampling times (Supplementary Table S1), but no significant difference in Ace and Chao1 was found among the different fertilization treatments and rice growth stages. The Shannon index was also not significantly different between the four fertilizer treatments and the three stages. Significant differences in the Simpson index were found among the S2C, S1B, S1C, and S1D samples, but no significant difference was found among the four fertilizer treatments at the same growth stage. The bacterial α-diversity indices were not correlated with the four fertilization treatments and the three rice growth stages.

Effects of the Rice Growth Stage and Fertilizer Treatment on Soil Bacterial Communities and Bacterial Functions

The PCoA and ANOSIM results showed that the bacterial community similarity distance was influenced by the fertilization treatment and rice growth stage (Figure 1 and Tables 3, 4). The samples from the four fertilizer treatments were divided by the PCoA axis1 (Figure 1). In particular, the control treatment without fertilizers (A) was quite different from the other three fertilization groups, as also shown by ANOSIM analysis (p < 0.05, Table 3). In addition, based on the PCoA axis1, the C treatment (NPK plus pig manure) differed from the B (NPK) and D (NPK plus straw) treatments at the S2 and S3 stages (Figure 1). However, no distinction was found between the B (NPK) and D (NPK plus straw) treatments (p > 0.05, Table 3). In contrast, the influence of the rice growth stages was mainly reflected in the PCoA axis2. The samples collected at S1 were distinguishable from those collected at S2 and S3, and ANOSIM analysis showed that the separation was still significant (Table 4, p < 0.01). In addition, the COG_L2 (Supplementary Figure S1A) and KEGG_L3 (Supplementary Figure S1B) profiles that were generated with the bacterial 16S rRNA gene results were used for the PCoA analysis of soil bacterial function. The results showed that the soil bacterial function could be separated among the S1 and other two stages but not the four fertilization treatments (Supplementary Figure S1).

FIGURE 1

Figure 1. PCoA (Bray-Curtis distance index) plots allowing visualization of the differences in the bacterial community structure between samples (based on OTU information). The background color of the sample names indicates the sampling time: purple indicates the S1 stage, nattier blue indicates the S2 stage, and grass green indicates the S3 stage. The color of the sample points indicates the four fertilizers treatments: red indicate the control treatment (A); brown indicate the NPK treatment (B); blue indicate the NPK + pig manure treatment (C), and orange indicate the NPK + straw treatment (D). The soil samples were arranged left-to-right based on the four fertilizers treatments and top-to-bottom based on the three rice growth stages (the later stages are found at the top).

TABLE 3

Table 3. ANOSIM analysis between four fertilizer treatments.

TABLE 4

Table 4. ANOSIM analysis between rice growth stages.

Correlations of Soil Properties With Soil Bacterial Communities and Functions

We performed a redundancy analysis (RDA) to understand the OTU-level bacterial community structure (Figure 2). Six parameters, namely, pH, OM, TN, AN, AP, and AK (Supplementary Table S2), were selected based on variance inflation factors with 999 Monte Carlo permutations. The p-values obtained for OM, AK, and pH were significantly lower than those of TP and TK. These six environmental variables together explained 31.13% of the total variance. In addition, the first axis of the RDA accounted for a higher variation of 14.76%, whereas the second axis accounted for a lower variation of 8.20% (Figure 2A). The A treatment (control without fertilizers) could be distinguished from the other treatments based on the direction of change in the pH, OM and TN. However, S1 was distinguished from S2 and S3 by the direction of changes in AK and AN. Variance partitioning analyses (VPA) were further performed to assess the contributions of soil pH and nutrients to the bacterial community variances. The results showed that the pH alone could explain the higher variation of 6.42%, whereas the OM, AK, TN, and AP contents explained lower variations of 3.05, 2.97, 2.38, 2.17, and 2.39%, respectively. These results suggested that pH, OM, and AK appeared to be key factors responsible for the microbial community changes observed in all the treatments. In addition, the RDA of bacterial function (Figure 2B, COG L2) also showed that S1 was distinct from S2 and S3 based on the different directions of change in the AK and AN. The control treatment could be distinguished from the others by the direction of change in the pH, OM, and TN (Figure 2B).

FIGURE 2

Figure 2. RDAs between soil samples at different stages based on OTU (A) and COG_L2 (B) levels. The colors of the sample names indicate the sampling time: red indicates the S1 stage, brown indicates the S2 stage, and blue indicates the S3 stage. The hollow squares indicate the control treatment (A); the circles indicate the NPK treatment (B); the triangles indicate the NPK + pig manure treatment (C), and the solid squares indicate the NPK + straw treatment (D). OM, organic matter; TN, total N; TP, total P; TK, total K; AN, available N; AP, available P; AK, available K.

Correlation Between the Variation in the Abundance of Soil Bacteria and Bacterial Function During the S1 and S2/S3 Stages

The top 100 bacterial families found at the highest abundances were used to evaluate the effects of the rice phonological stage and fertilizer treatment on the bacterial community. The STAMP results (average abundance, n = 3, p < 0.05) showed that the top 100 families exhibited significant differences between S1 and the other stages (Figure 3). From S1 to S3, the abundances of Anaerolineaceae, Desulfobulbaceae, Planctomycetaceae, Rhodospirillaceae, and Geobacteraceae increased, whereas those of Comamonadaceae, Opitutaceae, Gemmatimonadaceae, Bacteriovoracaceae, Nocardioidaceae, Gaiellaceae, and Cystobacteraceae decreased. In addition, compared with the control without fertilizers (A), the application of NPK fertilizer significantly influenced the bacterial family community (Supplementary Figure S2). Nitrospinaceae, Nocardioidaceae, Nitrosomonadaceae, Cytophagaceae, and Gaiellaceae were found at high abundances in the control treatment (A), whereas Acidobacteriaceae, Xanthomonadaceae, Cystobacteraceae, Gemmatimonadaceae, and Geobacteraceae were enriched in the NPK fertilizer treatments (B, C, and D).

FIGURE 3

Figure 3. Top 100 microorganism families during the different sampling times identified using STAMP. The different colors of the oblong blocks indicate the different rice growth stages. The blue blocks indicate the rice tillering stage (S1), and the yellow blocks indicate all the other samples (S2 and S3 stages) (p < 0.05, average proportion, n = 3).

The KEGG (Kyoto Encyclopaedia of Genes and Genomes) and COG metabolic pathways of soil bacteria were analyzed using PICRUSt. The results for metabolic functions were similar to the bacterial abundance results. Specifically, both the KEGG and COG metabolic pathways exhibited significant differences among the rice growth stages in the control without fertilizers (A) and the chemical fertilizer treatments (Supplementary Figures S3, S4). The results of the COG_L2 functional difference showed that there was a significant functional difference between S1 and S2 and S3 (Supplementary Figure S3, P < 0.05). Such as, [A] RNA processing and modification, [C] Energy production and conversion, [I] Lipid transport and metabolism, [K] Transcription, etc. For the fertilizer treatment, the control without fertilizers (A) was significantly different from the other three treatments in 15 functional classifications (Supplementary Figure S3, p < 0.05). Such as, [I] Lipid transport and metabolism, [E] Amino acid transport, and metabolism, [H] Coenzyme transport and metabolism, etc. (Supplementary Figure S3). There were no significant functional metabolic differences between B and D treatments, and only [A] RNA processing and modification differed between C and D treatments (Supplementary Figure S3, P < 0.05).

The results of the KEGG_L2 functional difference showed that S1 and S2 and S3 has a significant (P < 0.05) difference of nearly 21 functional classifications, including Energy Metabolism, Amino Acid Metabolism, Carbohydrate Metabolism, Membrane Transport like, etc. (Supplementary Figure S4, p < 0.05). But there were only 5 significant functional classifications between the S2 and S3 groups (Supplementary Figure S4, P < 0.05). For the fertilizer treatment, the control without fertilizers (A) was significantly different from the other three treatments in 14 functional classifications (Supplementary Figure S4, p < 0.05). Such as, Lipid Metabolism, Amino Acid Metabolism and Glycan Biosynthesis and Metabolism like etc. There were 4 functional classifications differences between B and C treatments, while there was no significant difference between the others treatments (Supplementary Figure S4).

The correlations among the differentiated bacterial families, soil properties and KEGG metabolic pathways are illustrated in Supplementary Table S3 and plotted in Supplementary Figure S5 (Ma et al., 2014). The abundances of bacterial families, such as Micromonosporaceae, Nocardioidaceae, Gaiellaceae, Anaerolineaceae, Syntrophobacteraceae, and Syntrophaceae, were correlated with bacterial KEGG metabolic pathways. Specifically, 70–80% of the identified bacterial families were significantly correlated with amino acid metabolism (Such as Alanine aspartate and glutamate metabolism, Tryptophan metabolism, Phenylalanine tyrosine, and tryptophan biosynthesis, Tyrosine metabolism, etc.), the biosynthesis of other secondary metabolites (Isoquinoline alkaloid biosynthesis, beta-Lactam resistance, Flavone and flavonol biosynthesis, etc.), carbohydrate metabolism (Citrate cycle (TCA cycle), Butanoate metabolism, Propanoate metabolism, etc.), glycan biosynthesis and lipid metabolic pathways (Protein kinases, Peptidoglycan biosynthesis, Glycosaminoglycan degradation, etc.) (Supplementary Figure S5). We also used the Spearman’s rank-order correlation method to evaluate the relationships between KEGG metabolic pathways and soil properties (Supplementary Figure S5) and found that the soil TK, AN, pH, and AK were significantly correlated with most metabolic pathways or differentiated bacterial families. In addition, the differentiated families were also significantly correlated with the primary soil properties, such as pH, AK, and AN (Supplementary Figure S6). These results showed that the soil properties, differentiated bacterial families and KEGG metabolic pathways were closely related to each other.

Molecular Ecology Network Analysis of Soil Bacterial Communities and Functions at the Different Rice Growth Stages

The potential bacteria-bacteria interactions were analyzed using bacterial co-occurrence networks (Deng et al., 2012). Both the network size and the degree of connectivity showed differences among the rice growth stages. At the S2 and S3 stages, the bacterial assemblages formed larger and more complex network topological structures, with increasing numbers of links and nodes, compared with the S1 stage networks (Figure 4 and Table 5). The S1 networks were smaller than the others. The increased numbers of nodes in the S2 network were connected and formed larger modules, whereas the size of each module in the network of the S3 stage was larger than that of the S1 network and smaller than that of the S2 networks. The increased complexity of the S2 networks displayed an average degree of increase and shorter harmonic geodesic distances (HD) (Table 5; Deng et al., 2012).

FIGURE 4

Figure 4. Succession of soil bacteria (A, OTU level) and bacterial metabolic pathway (B, KEGG L3 produced by PICRUSt) networks during the different stages of rice growth (tillering, heading and ripening). The networks represent the random matrix theory co-occurrence models derived from 12 biological replicates at each time point. In the models, the nodes represent the OTUs or KEGG L3 metabolic pathways. The red links indicate significant positive correlations between two individual nodes, whereas the blue links indicate negative covariations. The modules are randomly colored at each growth stage, and the nodes in the modules with less than six nodes are shown in black.

TABLE 5

Table 5. Global properties of the networks based on OTUs and KEGG_L3.

We also used the KEGG L3 metabolic pathway profile to construct bacterial functional networks for three rice growth stages. The bacterial functions exhibited denser connections at the KEGG L3 level. The S1 and S2 stages exhibited more similar network topological structures, but analyses of the average degree (avgK), average clustering coefficient (avgCC), and average path distance (GD) in all three networks revealed that the S2 network contained 20% more links than the S1 (Table 5). The structure of the S2 network is more complex than that of the S1 network (with higher avgK, a higher avgK indicates a more complex network). The number of nodes and links and the values of other network indices obtained for the S3 network were clearly higher than those found for the S1 and S2 networks. At the S3 stage, a greater number of bacterial functional pathways were connected and formed an intimate yet complicated co-occurrence network.

Discussion

Effects of Treatments and Rice Stage on the Soil Chemical Properties

The soil bacterial communities are commonly influenced by the availability of soil nutrients and their ability to absorb these limited nutrients (Noll et al., 2005). The soil nutrient status could be significantly affected by different chemical and organic fertilizer treatments. In this study, the soil nutrient status showed distinct changes among different stages of rice phonology. Similarly, the soil physical and chemical gradients are dependent on plant development over time (Krause et al., 2010). Some key soil properties, such as pH, AK, OM, and AN, were significantly affected by the fertilizer treatments and rice growth stages, and previous studies showed that pH, AK, OM, and AN are critical factors that influence the soil bacterial communities (Anderson et al., 2014; Wu et al., 2014; Sun et al., 2015; Zhou et al., 2015; Ding et al., 2016; Li et al., 2017). The changes in soil properties among different rice growth stages is likely due to the effects of the slow release of fertilizer, the uptake of nutrients by the plant and various root secretions (an effect on bulk soil mediated by secretions looks like a longshot, e.g., organic acids, amino acids, sugars, and growth substances) (Huang et al., 2014; Reinhold-Hurek et al., 2015).

Effects of Treatments and Rice Stages on the Soil Bacterial Community Composition

Although the bacterial α-diversity indices, such as Ace, Chao1, Shannon, and Simpson, were not significantly different (p < 0.05) among different fertilizer treatments and rice growth stages, the bacterial beta-diversity (PCoA and RDA) results showed that the bacterial community showed significant differences among the fertilization treatments and rice stages. At all the stages, the bacterial community composition obtained with the A treatment (control) was quite different from that obtained with the other three fertilization groups, and the B (NPK) and D (NPK + straw) treatments were grouped together. The inclusion of pig manure in the C treatment improved the bacterial community composition at the S2 and S3 stages. These results strengthened the notion that the amendment of livestock manure could prevent the negative effects of the long-term application of chemical fertilizers on the bacterial community in a field (Sun et al., 2015), such as the decrease in the number of microorganisms, the decrease in the diversity index, and the decrease in the abundance of beneficial microorganisms such as rhizobium and myxobacteria (Wang et al., 2016c). In this study, pig manure exerted stable effects on the bacterial community that appeared to persist throughout all rice stages.

Yu suggested that studies of the bacterial community should pay more attention to different stages other than fertilization (Yu et al., 2015). In this study, we selected three critical rice stages and found significant differences in the bacterial community. The succession of the bacterial community from the early growth stage (S1) to S2 and S3 could be clearly distinguished, whereas bacterial succession became relatively stable at the S2 and S3 stages. The finding that the soil bacterial community was relatively stable at the middle and later stages was also confirmed by PLFA analyses of a rice field (Noll et al., 2005). The researchers noted that aerobic and anaerobic conversion caused significant changes in the bacterial community at the early rice growth stages; specifically, the bacterial community gradually transitioned from the early to the late community structure at the middle stage but not at the late stage (Noll et al., 2005). Simultaneously, the differences between the bacterial community and functional gene succession processes at different stages (or sampling times) were also found in rice (Rui et al., 2009), wheat (Jing et al., 2015), maize (Li et al., 2014), nitrogen-fixing organisms (Wang et al., 2016b), and methanotrophs (Krause et al., 2010).

Growth-Stage Related Dynamics of the Soil Bacterial Function Community Structures

Soil bacterial functional genes show changes among the treatments and rice stages (Langille et al., 2013). The PCoA, RDA, and STAMP results indicated that the bacterial functional gene communities were affected by both the fertilization treatments and the rice developmental stages. We found an obvious succession of the bacterial functional community during rice development. In addition, the NPK fertilizers had significant effects on bacterial metabolic functions. Previous studies have revealed that higher-quality organic fertilizer management could improve the soil bacterial functional community (Guan et al., 2011; Yu-Hong et al., 2011). The bacterial functional structures at the S1 stage were significantly different from those at the S2 and S3 stages, which presented relatively stable bacterial functional structures. Other studies revealed a similar effect in the nifH gene (Yuan et al., 2015; Wang et al., 2016b), carboxylic acids (Yu et al., 2015), and anammox rate (Li et al., 2016). It is likely that the combination of physical (root structure characteristics) and chemical (root exudation) changes contributes to the succession of the bacterial functional structure during rice development (Yuan et al., 2015). In addition, the RDA results also proved that the soil available nutrients AN and AK appeared to be key factors in the succession of the bacterial function community.

In addition, based on the correlations between the most abundant differentiated bacterial families and KEGG metabolic pathways and between soil properties and KEGG metabolic pathways, we concluded that there is a significant correlation among bacterial composition, KEGG metabolic pathways and soil properties. The succession of the bacterial compositions, KEGG metabolic pathways and soil properties occurred nearly concurrently with the progression through the rice developmental stages.

Molecular Ecology Network Analysis of Soil Bacterial Communities and Functions at Different Rice Growth Stages

Molecular ecology network analysis has been widely used to explore the interactions between soil microorganisms and function genes (Zhou et al., 2010, 2011; Ren et al., 2015), and the results of these analyses have helped us understand the dynamic changes in microbial niches or function genes (Shi et al., 2016). In addition, recent studies have indicated that a network analysis of bacterial communities might allow the assessment of soil ecology or their contribution to habitat niches (Weiss et al., 2016). Compared with that at the early stage, the network structures at the S2 and S3 stages generated from the bacterial and KEGG function assemblages exhibited increased complexity. The interactions between bacteria or their KEGG functions varied during different plant growth stages. The same phenomenon has also been detected in vegetative rhizosphere soil (Shi et al., 2016), and changes in network complexity have even been found between fallow and mature periods in tobacco fields (Niu et al., 2016), between different seasons in eutrophic Lake Mendota (Kara et al., 2013), and between different years in salt marsh chronosequence (Dini-Andreote et al., 2014). The succession of networks during rice growth was reflected in the soil chemical properties and bacterial community changes.

The key question addressed in this study is the primary factor driving the succession of the network, and the quality and quantity of root exudates might be this factor because organic acids show changes among different rice phonological stages (Li et al., 2016). In addition, significantly lower root exudates are found in bulk soil (Li et al., 2016). Previous research suggests that roots promote the development of soil microbial niches, which results in interactions and covariations among the bacteria (Shi et al., 2016). The differences in temperatures, precipitation and sunlight hours among different stages might be responsible for the differences in the growth characteristics of cultivated bacteria. Climate change has altered plant phenology and microbial communities (Classen et al., 2015).

However, the network of functional genes, unlike the bacterial network, at the S3 stage exhibits a more complex (higher avgK) network organization than those at the S1 and S2 stages. The succession of the soil bacterial community network structure was only partially consistent with the changes in function. A more complex soil bacterial community network was formed at the middle stage of rice development, and the topological of the community network maintained a relatively similar structure (the parameters of the community network are nearly consistent) during the middle and late stages. It is worth noting that the topological of the functional network structure in the early and middle stages is more similar, but formed a complicated co-occurrence network at the S3 stage. The topological of functional network structure of bacteria changes a little in the early and middle stages of rice, while its changes significantly in the ripening stage of rice growth. It is possible that the topological of functional networks are closely related to the rice maturity stage. In addition, during rice growth, the formation of complex bacterial functional networks occurred later than the formation of complex bacterial networks, although changes in the microbial structure and functions are commonly observed under different environmental conditions (Mäder et al., 2002; Graham et al., 2014; Riah-Anglet et al., 2015), and the topological of changes in soil ecology function were not consistent with those in the community composition in a rice paddy (Chen et al., 2015). As described by Fry, the function of a grassland ecosystem is influenced by the experimental treatments but is not related to associated changes in the microbial community (Fry et al., 2016).

Author Contributions

WW performed all the experiments, coordinated the data analyses, and prepared the manuscript. XL, YC, and XY contributed to the preparation of the manuscript and data analyses. ZheC and WR conceived the research study. ZLC and HW supervised the entire study.

Funding

This work was financially supported by the “973” Program (2015CB150502), the National Key R&D Program (2016YFD0200309), the National Natural Science Foundation of China (No. 41371319), Mechanisms of Marginal Land Productivity Expansion and Technologies for Storing Grain in Land (KFZD-SW-112), and Genetically Modified Organisms Breeding Major Projects of China (Grant No. 2014ZX08011-003).

Conflict of Interest Statement

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Acknowledgments

We thank Biomics (Beijing) Biotech Co., Ltd. for assistance in bioinformatics analysis.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fmicb.2019.00421/full#supplementary-material

Footnotes

1. ^http://qiime.org/
2. ^http://ieg2.ou.edu/MENA/

References

Anderson, C., Hamonts, K., Clough, T., and Condron, L. (2014). Biochar does not affect soil N-transformations or microbial community structure under ruminant urine patches but does alter relative proportions of nitrogen cycling bacteria. Agric. Ecosyst. Environ. 191, 63–72. doi: 10.1016/j.agee.2014.02.021

CrossRef Full Text | Google Scholar

Bastian, M., Heymann, S., and Jacomy, M. (2009). Gephi: an open source software for exploring and manipulating networks. ICWSM 8, 361–362.

Google Scholar

Blanchet, G., Gavazov, K., Bragazza, L., and Sinaj, S. (2016). Responses of soil properties and crop yields to different inorganic and organic amendments in a Swiss conventional farming system. Agric. Ecosyst. Environ. 230, 116–126. doi: 10.1016/j.agee.2016.05.032

CrossRef Full Text | Google Scholar

Bolger, A. M., Lohse, M., and Usadel, B. (2014). Trimmomatic: a flexible trimmer for Illumina sequence data. Bioinformatics 30, 2114–2120. doi: 10.1093/bioinformatics/btu170

PubMed Abstract | CrossRef Full Text | Google Scholar

Chen, J., Liu, X., Li, L., Zheng, J., Qu, J., Zheng, J., et al. (2015). Consistent increase in abundance and diversity but variable change in community composition of bacteria in topsoil of rice paddy under short term biochar treatment across three sites from South China. Appl. Soil Ecol. 91, 68–79. doi: 10.1016/j.apsoil.2015.02.012

CrossRef Full Text | Google Scholar

Chen, Z., Wang, H., Liu, X., Zhao, X., Lu, D., Zhou, J., et al. (2017). Changes in soil microbial community and organic carbon fractions under short-term straw return in a rice–wheat cropping system. Soil Tillage Res. 165, 121–127. doi: 10.1016/j.still.2016.07.018

CrossRef Full Text | Google Scholar

Classen, A. T., Sundqvist, M. K., Henning, J. A., Newman, G. S., Moore, J. A., Cregger, M. A., et al. (2015). Direct and indirect effects of climate change on soil microbial and soil microbial-plant interactions: what lies ahead? Ecosphere 6, 1–21. doi: 10.1890/ES15-00217.1

CrossRef Full Text | Google Scholar

Deng, Y., Jiang, Y.-H., Yang, Y., He, Z., Luo, F., and Zhou, J. (2012). Molecular ecological network analyses. BMC Bioinformatics 13:113. doi: 10.1186/1471-2105-13-113

PubMed Abstract | CrossRef Full Text | Google Scholar

Ding, J., Jiang, X., Ma, M., Zhou, B., Guan, D., Zhao, B., et al. (2016). Effect of 35 years inorganic fertilizer and manure amendment on structure of bacterial and archaeal communities in black soil of northeast China. Appl. Soil Ecol. 105, 187–195. doi: 10.1016/j.apsoil.2016.04.010

CrossRef Full Text | Google Scholar

Dini-Andreote, F., Silva, M. D. C. P., Triado-Margarit, X., Casamayor, E. O., Van Elsas, J. D., and Salles, J. F. (2014). Dynamics of bacterial community succession in a salt marsh chronosequence: evidences for temporal niche partitioning. ISME J. 8, 1989–2001. doi: 10.1038/ismej.2014.54

PubMed Abstract | CrossRef Full Text | Google Scholar

Edgar, R. C. (2013). UPARSE: highly accurate OTU sequences from microbial amplicon reads. Nat. Methods 10, 996–998. doi: 10.1038/nmeth.2604

PubMed Abstract | CrossRef Full Text | Google Scholar

Edgar, R. C., Haas, B. J., Clemente, J. C., Quince, C., and Knight, R. (2011). UCHIME improves sensitivity and speed of chimera detection. Bioinformatics 27, 2194–2200. doi: 10.1093/bioinformatics/btr381

PubMed Abstract | CrossRef Full Text | Google Scholar

Edwards, J., Johnson, C., Santos-Medellín, C., Lurie, E., Podishetty, N. K., Bhatnagar, S., et al. (2015). Structure, variation, and assembly of the root-associated microbiomes of rice. Proc. Natl. Acad. Sci. 112, E911–E920. doi: 10.1073/pnas.1414592112

PubMed Abstract | CrossRef Full Text | Google Scholar

Fry, E. L., Manning, P., Macdonald, C., Hasegawa, S., De Palma, A., Power, S. A., et al. (2016). Shifts in microbial communities do not explain the response of grassland ecosystem function to plant functional composition and rainfall change. Soil Biol. Biochem. 92, 199–210. doi: 10.1016/j.soilbio.2015.10.006

CrossRef Full Text | Google Scholar

Graham, E. B., Wieder, W. R., Leff, J. W., Weintraub, S. R., Townsend, A. R., Cleveland, C. C., et al. (2014). Do we need to understand microbial communities to predict ecosystem function? A comparison of statistical models of nitrogen cycling processes. Soil Biol. Biochem. 68, 279–282. doi: 10.1016/j.soilbio.2013.08.023

CrossRef Full Text | Google Scholar

Guan, G., Tu, S.-X., Yang, J.-C., Zhang, J.-F., and Li, Y. (2011). A field study on effects of nitrogen fertilization modes on nutrient uptake, crop yield and soil biological properties in rice-wheat rotation system. Agric. Sci. Chin. 10, 1254–1261. doi: 10.1016/S1671-2927(11)60117-X

CrossRef Full Text | Google Scholar

Huang, X.-F., Chaparro, J. M., Reardon, K. F., Zhang, R., Shen, Q., and Vivanco, J. M. (2014). Rhizosphere interactions: root exudates, microbes, and microbial communities. Botany 92, 267–275. doi: 10.1139/cjb-2013-0225

CrossRef Full Text | Google Scholar

Jing, H., Xia, X., Liu, H., Zhou, Z., Wu, C., and Nagarajan, S. (2015). Anthropogenic impact on diazotrophic diversity in the mangrove rhizosphere revealed by nifH pyrosequencing. Front. Microbiol. 6:1172. doi: 10.3389/fmicb.2015.01172

PubMed Abstract | CrossRef Full Text | Google Scholar

Kara, E. L., Hanson, P. C., Hu, Y. H., Winslow, L., and McMahon, K. D. (2013). A decade of seasonal dynamics and co-occurrences within freshwater bacterioplankton communities from eutrophic Lake Mendota, WI, USA. ISME J. 7, 680–684. doi: 10.1038/ismej.2012.118

PubMed Abstract | CrossRef Full Text | Google Scholar

Krause, S., Lüke, C., and Frenzel, P. (2010). Succession of methanotrophs in oxygen–methane counter-gradients of flooded rice paddies. ISME J. 4, 1603–1607. doi: 10.1038/ismej.2010.82

PubMed Abstract | CrossRef Full Text | Google Scholar

Krych, L., Hansen, C. H., Hansen, A. K., Van Den Berg, F. W., and Nielsen, D. S. (2013). Quantitatively different, yet qualitatively alike: a meta-analysis of the mouse core gut microbiome with a view towards the human gut microbiome. PLoS One 8:e62578. doi: 10.1371/journal.pone.0062578

PubMed Abstract | CrossRef Full Text | Google Scholar

Langille, M. G., Zaneveld, J., Caporaso, J. G., McDonald, D., Knights, D., Reyes, J. A., et al. (2013). Predictive functional profiling of microbial communities using 16S rRNA marker gene sequences. Nat. Biotechnol. 31, 814–821. doi: 10.1038/nbt.2676

PubMed Abstract | CrossRef Full Text | Google Scholar

Li, F., Chen, L., Zhang, J., Yin, J., and Huang, S. (2017). Bacterial community structure after long-term organic and inorganic fertilization reveals important associations between soil nutrients and specific taxa involved in nutrient transformations. Front. Microbiol. 8:187. doi: 10.3389/fmicb.2017.00187

PubMed Abstract | CrossRef Full Text | Google Scholar

Li, H., Yang, X., Weng, B., Su, J., Nie, S. A., Gilbert, J. A., et al. (2016). The phenological stage of rice growth determines anaerobic ammonium oxidation activity in rhizosphere soil. Soil Biol. Biochem. 100, 59–65. doi: 10.1016/j.soilbio.2016.05.015

CrossRef Full Text | Google Scholar

Li, X., Rui, J., Mao, Y., Yannarell, A., and Mackie, R. (2014). Dynamics of the bacterial community structure in the rhizosphere of a maize cultivar. Soil Biol. Biochem. 68, 392–401. doi: 10.1371/journal.pone.0178425

PubMed Abstract | CrossRef Full Text | Google Scholar

Ma, J., Prince, A. L., Bader, D., Hu, M., Ganu, R., Baquero, K., et al. (2014). High-fat maternal diet during pregnancy persistently alters the offspring microbiome in a primate model. Nat. Commun. 5:3889. doi: 10.1038/ncomms4889

PubMed Abstract | CrossRef Full Text | Google Scholar

Mäder, P., Fliessbach, A., Dubois, D., Gunst, L., Fried, P., and Niggli, U. (2002). Soil fertility and biodiversity in organic farming. Science 296, 1694–1697. doi: 10.1126/science.1071148

PubMed Abstract | CrossRef Full Text | Google Scholar

Magoè, T., and Salzberg, S. L. (2011). FLASH: fast length adjustment of short reads to improve genome assemblies. Bioinformatics 27, 2957–2963. doi: 10.1093/bioinformatics/btr507

PubMed Abstract | CrossRef Full Text | Google Scholar

Niu, J., Rang, Z., Zhang, C., Chen, W., Tian, F., Yin, H., et al. (2016). The succession pattern of soil microbial communities and its relationship with tobacco bacterial wilt. BMC Microbiol. 16:233. doi: 10.1186/s12866-016-0845-x

PubMed Abstract | CrossRef Full Text | Google Scholar

Noll, M., Matthies, D., Frenzel, P., Derakshani, M., and Liesack, W. (2005). Succession of bacterial community structure and diversity in a paddy soil oxygen gradient. Environ. Microbiol. 7, 382–395. doi: 10.1111/j.1462-2920.2005.00700.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Parks, D. H., Tyson, G. W., Hugenholtz, P., and Beiko, R. G. (2014). STAMP: statistical analysis of taxonomic and functional profiles. Bioinformatics 30, 3123–3124. doi: 10.1093/bioinformatics/btu494

PubMed Abstract | CrossRef Full Text | Google Scholar

Paul, E., and Clark, F. (1989). Occurrences and Distribution of Soil Organics. Soil Microbiology and Biochemistry. San Diego, CA: Academic Press, 81–84.

Reinhold-Hurek, B., Bünger, W., Burbano, C. S., Sabale, M., and Hurek, T. (2015). Roots shaping their microbiome: global hotspots for microbial activity. Ann. Rev. Phytopathol. 53, 403–424. doi: 10.1146/annurev-phyto-082712-102342

PubMed Abstract | CrossRef Full Text | Google Scholar

Ren, D., Madsen, J. S., Sørensen, S. J., and Burmølle, M. (2015). High prevalence of biofilm synergy among bacterial soil isolates in cocultures indicates bacterial interspecific cooperation. ISME J. 9, 81–89. doi: 10.1038/ismej.2014.96

PubMed Abstract | CrossRef Full Text | Google Scholar

Riah-Anglet, W., Trinsoutrot-Gattin, I., Martin-Laurent, F., Laroche-Ajzenberg, E., Norini, M.-P., Latour, X., et al. (2015). Soil microbial community structure and function relationships: a heat stress experiment. Appli. Soil Ecol. 86, 121–130. doi: 10.1016/j.apsoil.2014.10.001

CrossRef Full Text | Google Scholar

Rui, J., Peng, J., and Lu, Y. (2009). Succession of bacterial populations during plant residue decomposition in rice field soil. Appl. Environ. Microbiol. 75, 4879–4886. doi: 10.1128/AEM.00702-09

PubMed Abstract | CrossRef Full Text | Google Scholar

Schloss, P. D., Gevers, D., and Westcott, S. L. (2011). Reducing the effects of PCR amplification and sequencing artifacts on 16S rRNA-based studies. PLoS One 6:e27310. doi: 10.1371/journal.pone.0027310

PubMed Abstract | CrossRef Full Text | Google Scholar

Shi, S., Nuccio, E. E., Shi, Z. J., He, Z., Zhou, J., and Firestone, M. K. (2016). The interconnected rhizosphere: High network complexity dominates rhizosphere assemblages. Ecol. Lett. 19, 926–936. doi: 10.1111/ele.12630

PubMed Abstract | CrossRef Full Text | Google Scholar

Shi, S., Richardson, A. E., O’Callaghan, M., Firestone, M., and Condron, L. (2013). Challenges in assessing links between root exudates and the structure and function of soil microbial communities. Mol. Microb. Ecol. Rhizosphere 2, 125–135. doi: 10.1002/9781118297674.ch11

CrossRef Full Text | Google Scholar

Sun, R., Zhang, X.-X., Guo, X., Wang, D., and Chu, H. (2015). Bacterial diversity in soils subjected to long-term chemical fertilization can be more stably maintained with the addition of livestock manure than wheat straw. Soil Biol. Biochem. 88, 9–18. doi: 10.1016/j.soilbio.2015.05.007

CrossRef Full Text | Google Scholar

Wang, J., Xue, C., Song, Y., Wang, L., Huang, Q., and Shen, Q. (2016a). Wheat and rice growth stages and fertilization regimes alter soil bacterial community structure, but not diversity. Front. Microbiol. 7:1207. doi: 10.3389/fmicb.2016.01207

PubMed Abstract | CrossRef Full Text | Google Scholar

Wang, J., Zhang, D., Zhang, L., Li, J., Raza, W., Huang, Q., et al. (2016b). Temporal variation of diazotrophic community abundance and structure in surface and subsoil under four fertilization regimes during a wheat growing season. Agric. Ecosyst. Environ. 216, 116–124. doi: 10.1016/j.agee.2015.09.039

CrossRef Full Text | Google Scholar

Wang, W., Wang, H., Feng, Y., Wang, L., Xiao, X., Xi, Y., et al. (2016c). Consistent responses of the microbial community structure to organic farming along the middle and lower reaches of the Yangtze River. Sci. Rep. 6:35046. doi: 10.1038/srep35046

PubMed Abstract | CrossRef Full Text | Google Scholar

Weiss, S., Van Treuren, W., Lozupone, C., Faust, K., Friedman, J., Deng, Y., et al. (2016). Correlation detection strategies in microbial data sets vary widely in sensitivity and precision. ISME J. 10, 1669–1681. doi: 10.1038/ismej.2015.235

PubMed Abstract | CrossRef Full Text | Google Scholar

Wu, K., Yuan, S., Wang, L., Shi, J., Zhao, J., Shen, B., et al. (2014). Effects of bio-organic fertilizer plus soil amendment on the control of tobacco bacterial wilt and composition of soil bacterial communities. Biol. Fertil. Soils 50, 961–971. doi: 10.1007/s00374-014-0916-9

CrossRef Full Text | Google Scholar

Yu, C., Hu, X., Deng, W., Li, Y., Xiong, C., Ye, C., et al. (2015). Changes in soil microbial community structure and functional diversity in the rhizosphere surrounding mulberry subjected to long-term fertilization. Appl. Soil Ecol. 86, 30–40. doi: 10.1016/j.apsoil.2014.09.013

CrossRef Full Text | Google Scholar

Yuan, J., Chaparro, J. M., Manter, D. K., Zhang, R., Vivanco, J. M., and Shen, Q. (2015). Roots from distinct plant developmental stages are capable of rapidly selecting their own microbiome without the influence of environmental and soil edaphic factors. Soil Biol. Biochem. 89, 206–209. doi: 10.1016/j.soilbio.2015.07.009

CrossRef Full Text | Google Scholar

Yu-Hong, Y., Dong-Mei, C., Yan, J., Hai-Bin, W., Yu-Qi, D., Xu-Kui, G., et al. (2011). Effect of different fertilizers on functional diversity of microbial flora in rhizospheric soil under tobacco monoculture. Acta Agronomica Sin. 37, 105–111.

Google Scholar

Zhao, J., Ni, T., Li, J., Lu, Q., Fang, Z., Huang, Q., et al. (2016). Effects of organic–inorganic compound fertilizer with reduced chemical fertilizer application on crop yields, soil biological activity and bacterial community structure in a rice–wheat cropping system. Appl. Soil Ecol. 99, 1–12. doi: 10.1016/j.apsoil.2015.11.006

CrossRef Full Text | Google Scholar

Zhou, J., Deng, Y., Luo, F., He, Z., Tu, Q., and Zhi, X. (2010). Functional molecular ecological networks. mBio 1:e00169–10. doi: 10.1128/mBio.00169-10

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhou, J., Deng, Y., Luo, F., He, Z., and Yang, Y. (2011). Phylogenetic molecular ecological network of soil microbial communities in response to elevated CO2. mBio 2:e122–11. doi: 10.1128/mBio.00122-11

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhou, J., Guan, D., Zhou, B., Zhao, B., Ma, M., Qin, J., et al. (2015). Influence of 34-years of fertilization on bacterial communities in an intensively cultivated black soil in northeast China. Soil Biol. Biochem. 90, 42–51. doi: 10.1016/j.soilbio.2015.07.005

CrossRef Full Text | Google Scholar