The long-term experimental field site with the wheat-rice rotation system is situated in Changshu, Jiangsu province, China (31°18′N, 120°37′E, 6 m asl), and was established in 2005. The site has a humid subtropical monsoon climate with an average annual rainfall of ≈1063 mm and mean minimum and maximum temperatures in 2013 of 2 and 36°C, respectively.

Every year, after rice is harvested, wheat is sown using 150 kg seeds ha^-1 in October or November. Rice is transplanted in June using two seedlings per hill at 13 cm × 28 cm spacing after wheat has been harvested. The rice plots are flooded with 5 cm of standing water from July to September. Crops are harvested manually at ground level by sickle, and the above-ground biomass removed. After the wheat and rice harvests, the plots are plowed to a depth of 20 cm and fertilizers are applied before sowing and planting as basal fertilizer.

The randomized block plots were 6 m × 7 m in size. Four treatments with three replicates each were established in a randomized block design: (1) NNF, in which 750 kg ha^-1 superphosphate (12% P₂O₅) and 183 kg ha^-1 muriate of potash (60% K₂O) were applied; (2) NPK CF, in which 391 kg ha^-1 urea (46% N), 750 kg ha^-1 superphosphate and 183 kg ha^-1 muriate of potash (mixed and pilled) were applied; (3) OIMF, in which 1500 kg ha^-1 OIMF (11.0% organic C, 12.0% total N, 4.1% P₂O₅ and 4.1% K₂O, mixed and pilled) was applied, with an extra application of 28.5 kg ha^-1 superphosphate and 48.5 kg ha^-1 muriate of potash to reach total quantities of N, P₂O₅, and K₂O that were equal to the CF plots; (4) OF, in which 4500 kg ha^-1 OF only (26.4% organic C, 2.5% total N, 1.6% P₂O₅, and 1.3% K₂O, from composted pig manure and rice straw by Tianniang Ltd., of Changshu, China) was applied. All the fertilizers were applied as basal fertilizers in October 25th, 2012 and June 5th, 2013 after the harvest of rice and wheat, respectively, and before sowing or transplanting of the subsequent crop. The details on the experimental design can be found in Wang et al. (2016).

Bulk soil samples at a depth of 0–20 cm were collected in 2013 in March 1st (Mar, wheat tillering), April 1st (Apr, wheat jointing), May 1st (May, wheat heading), June 4th (Jun, wheat ripening), July 7th (Jul, rice tillering), August 16th (Aug, rice jointing), September 15th (Sep, rice heading), and October 31st (Oct, rice ripening). For each plot, four cores (5 cm in diameter) were collected, which were pooled and sieved (2 mm) to remove the aboveground plant materials, roots and stones. Fresh soil was used for DNA extraction immediately after collection. Sub-samples were frozen at -80°C.

Soil pH and electrical conductivity (EC) were determined with a soil-to-water ratio of 1:2.5. Soil organic carbon (SOC) was determined by the K₂Cr₂O₇ oxidation-reduction titration method (Schollenberger, 1931). The Kjeldahl method was used for total nitrogen (TN) estimation (Page, 1982). Soil available nitrogen (AN, nitrate plus ammonium) was extracted with 2 M KCl and quantified using a BRAN+LUEBBE AutoAnalyzer 3. Soil available phosphorus (AP) was extracted with sodium bicarbonate and then measured by the molybdenum-blue method (Olsen, 1954). Soil available potassium (AK) was extracted with ammonium acetate and measured using a flame atomic absorption spectrophotometer (Brown, 1998). Soil moisture content was determined by drying the soil at 105°C for 12 h. Soil temperature (1 dm) was determined at 14:00 on each sampling day.

DNA was extracted from 0.25 fresh soil with the MoBio PowerSoil^TM DNA Isolation Kits (Mo Bio Laboratories, Carlsbad, CA, USA) according to the manufacturer’s instructions. The quantity and quality of DNA extracts were assayed using a Nanodrop ND-2000 UV–VIS Spectrophotometer (NanoDrop Technologies, Wilmington, DE, USA), and the DNA was stored at -80°C for future analyses.

Amplification of 16S rRNA gene V4 hypervariable region was performed using PCR mixers containing 15 μL of Master Mix (Qiagen Inc., Valencia, CA, USA), 0.5 μM 515F primer (5′-GTGCCAGCMGCCGCGGTAA-3′), 0.5 μM 806R primer (5′-GGACTACHVGGGTWTCTAAT-3′; Caporaso et al., 2012), 10 ng DNA template, and ddH₂O to a final volume of 30 μL. The Illumina sequencing adapter and the reverse primer contain a 6-bp barcode specific for each example. PCR was performed in triplicate for each sample using the following program: 1 min at 98°C for initial denaturing, followed by 35 cycles of 98°C for 10 s, 50°C for 30 s and 72°C for 30 s with the final extension for 5 min at 72°C. After amplification, the triplicate PCRs were pooled and purified using the PCR cleanup Kit (Axygen Biosciences, Union City, CA, USA). Sequencing was performed using the Illumina MiSeq platform (Caporaso et al., 2012) at Novogene Bioinformatics Technology Co., Ltd., Beijing, China.

The UPARSE standard pipeline was applied to analyze the sequence data (Edgar, 2013). Briefly, short reads (<250 bp) sequences were filtered out before downstream analysis. Sequences with ≥ 97% similarity were clustered into operational taxonomic units (OTUs). All sequences were assigned using the RDP classifier to identify taxonomies with a confidence threshold of 0.8 (Wang et al., 2007). An average of 1.9% of the sequences was identified as being of archaeal origin and these sequences were also included in further data analyses.

As the number of sequences per sample ranged from 4,241 to 32,286, the data from each sample was resampled to a depth of 4,241 sequences for each sample using MOTHUR (Schloss et al., 2009), and the resulting new OTU table used for the further analyses. MOTHUR was then used to calculate α-diversity, expressed as Shannon and Simpson indexes, for each sample. Multivariate ANOVAs based on the Shannon and Simpson indexes and 16S rRNA gene copies were performed to determine the effects of stage and fertilization regime on the bacterial diversities and abundance using SPSS 18.0. Hierarchical cluster analysis was carried out based on the Bray–Curtis distances of OTU composition using MOTHUR. Non-metric multidimensional scaling (NMDS), based on the Bray–Curtis distances of phyla relative abundance (RA). To compare the effects of stage, fertilization regime and their interactions on bacterial community structures, permutational multivariate analysis was performed based on OTU RA. Correlations between the Euclidean distance of soil properties and bacterial OTUs RA were calculated using the Mantel test. Redundancy analysis (RDA) was performed to analyze relationships between soil properties and bacterial dominant taxa, and followed by 999 permutations to test significance. All these analyses were performed using the ‘vegan’ package in R (R Core Team, 2014). One-way ANOVAs of soil properties were performed to compare the differences among the four fertilization regimes within the same sampling stage, followed by Tukey’s HSD as a post hoc test using SPSS 18.0. Pearson correlation analyses were performed using SPSS 18.0. Linear discriminant analysis effect size (LEfSe) was performed to identify significant differences in bacterial taxa between groups. The Kruskal–Wallis (KW) sum-rank test was used in LEfSe analysis to detect the features with significantly different abundances between assigned classes, and then liner discriminant analysis (LDA) was performed to estimate the effect size of each feature (Segata et al., 2011).

The abundance of total bacterial 16S rRNA gene copies was quantified by quantitative PCR as well the abundances of Alphaproteobacteria, Betaproteobacteria, and Acidobacteria to complement the results of the high-throughput sequencing. Quantitative PCR (QPCR) assays were performed using the SYBR^® Premix Ex TaqTM (Perfect Real Time) kit (TaKaRa Biotechnology Co., Dalian, China) using the 7500 Fast Real-Time PCR System (Applied Biosystems, Branchburg, New Jersey, USA). Each reaction was performed in a 25 μL volume containing 12.5 μL of SYBR^® Premix Ex TaqTM (2×, Takara), 0.5 μL ROX Reference dye II (50×, TaKaRa Biotechnology Co., Dalian, China), 10 μL dd H₂O, 1 μL (∼10–30 ng) DNA template and 0.5 μL (5 μM) of each primer. Bacterial 16S rRNA primers were Eub338 (5′-ACT CCT ACG GGA GGC AGC AG-3′) and Eub518 (5′-ATT ACC GCG GCT GCT GG-3′), Alphaproteobacteria class primers were Eub338 (5′-ACT CCT ACG GGA GGC AGC AG-3′) and Alfa685 (5′-TCT ACG RAT TTC ACC YC TAC-3′); Betaproteobacteria class primers were Eub338 and Beta510 (5′-TCA CTG CTA CAC GYG-3′); Acidobacteria primers were Acid31 (5′-GAT CCT GGC TCA GAA TC-3′); and Eub518 (5′-ATT ACC GCG GCT GCT GG-3′) (Fierer et al., 2005). The PCR protocols for Alphaproteobacteria, Betaproteobacteria, and Acidobacteria were initiated with an activation step at 95°C for 30 s, followed by 40 cycles with 5 s at 95°C for denaturation and 34 s at 60°C for annealing and extension. Reactions were tested for possible PCR inhibitors contained in DNA extractions by comparing results across serial dilutions. No significant inhibitory effects of DNA extracts was detected. Standard curves were developed by serially diluting plasmids with known positive inserts to final concentrations of 10⁸ to 10² gene copies μL^-1. QPCR efficiencies were 111, 115, 95, and 98% for bacteria, Alphaproteobacteria, Betaproteobacteria, and Acidobacteria reactions, respectively, and R² values for all four assays were >0.99.

All sequence data has been deposited in the NCBI Sequence Read Archive (SRA) database under the accession number SRX1117202.

Soil organic carbon content was significantly higher in the OF treatment as compared to the other three fertilizer treatments (Supplementary Table S1). OF plots also tended to have higher soil TN content, pH and EC values relative to the other three treatments. However, the CF and OIMF treatments appeared to decrease the soil pH and EC values and increase soil NO₃^- content. Soil NH₄⁺ content tended to be higher in the CF and OIMF treatments. Soil AK contents were significantly higher for the NNF treatment compared to the other three treatments; however, CF treatment had the lowest AK contents. Soil AP content tended to be higher in OF treatment in August, September, and October.

Across the total dataset, 1,330,350 sequences from 96 soil samples were grouped into 13,379 OTUs at the 97% similarity cut-off level. An average of 98.1% of the sequences were classified as belonging to the bacteria, while 1.9% of them were designated as being of archaea. The most dominant phylum across all samples was the Proteobacteria, which accounted for an average of 34.0% of the total sequences, followed by Acidobacteria, Chloroflexi, and Actinobacteria, which represented 8.8, 7.3, and 5.9% of sequences, respectively (Supplementary Table S2).

As shown in the Table 1, plant growth stage significantly affected bacterial abundance, which generally increased from 4.56 × 10⁹ g^-1 dry soil in March to peak level of 1.03 × 10¹⁰ g^-1 dry soil in May (Figure 1). There after bacterial 16S rRNA gene copies dropped back to a minimum 3.70 × 10⁹ g^-1 dry soil on the first sampling stage of the rice rotation (July), after which, bacterial abundance again increased and became stable. Fertilization regime had a relative weak influence on the bacterial abundance across the eight stages (Table 1), but the effect of the interaction of stage and fertilization was significant.

Permutational multivariate analyses based on bacterial RA of phyla and OTUs revealed that stage and fertilization regime both had significant effects on the bacterial community structures (Table 1). The plant growth stage yielded higher R² values than fertilization regime at both the phylum and OTU levels, indicating that the influence of crops growth stage on the bacterial structure was stronger compared to fertilization regime. NMDS analysis based on the RA of phyla showed that samples of bacterial community structure were grouped by their sampling stages, but not by different fertilization regime, and that adjacent stages had more similar bacterial community structures than more temporally separated stages (Figure 2). Similar results were also observed by the cluster analysis based on OTU data (Supplementary Figure S2).

Bacterial community structures were apparently changed varied somewhat during the course of the year as environmental conditions and plant growth stages changed, but the dominant taxa remained Proteobacteria, Acidobacteria, Chloroflexi, and Actinobacteria (Supplementary Table S2). We further divided the eight sampling stages into two groups, representing the wheat rotation (March–June) and the rice rotation (July–October), and performed LEfSe analyses to examine which taxa differed most in RA between the two rotations. The cladogram (Figure 3) shows that many taxa were common between the wheat and rice stages of the rotation (shown in yellow), but there were also some specific differences. For example, the RA of Actinobacteria, Gammaproteobacteria, and Alphaproteobacteria in soils under wheat rotation were higher than in soils under rice rotation, while RA of Verrucomicrobia was significantly higher in the rice rotation soils (Figure 3; Supplementary Table S2).

With respect to bacterial diversity, multivariate ANOVA revealed that plant growth stage and fertilization regime had no significant effects on either Shannon or Simpson indexes (Table 1), indicating that bacterial diversity levels were stable both across the plant growth stages and under different fertilization regimes.

Although showing small effect than crop growth stage, fertilization did significantly affect bacterial community structures, as revealed by permutational multivariate analyses. Persistent differences in bacterial taxa patterns between fertilization regimes, regardless of the different growth stages, were revealed by LEfSe analyses (Figure 4).

The CF and NNF treatments were compared pairwise to identify the influence of the nitrogen application on the bacterial community compositions (Figure 4A). The results revealed that across all stages of the wheat-rice rotation, the NNF treatment resulted in significantly higher (P< 0.05) RAs of the Cyanobacteria phylum, Cyanobacteria class and Cyanobacteria order, while the CF treatment had a higher RA of Acidobacteria_GP1.

The CF and OIMF treatments were compared pairwise to determine how bacterial composition was affected by the partial substitution of OF for CF (Figure 4B). The results showed that the RA of Archaea was significantly higher in the CF plots across the whole wheat-rice rotation seasons. Within the Archaea, the phyla Chloroflexi and Thaumarchaeota, and the classes Anaerolineae and Nitrososphaerales were also significantly higher in the CF treatment.

Chemical fertilizer and OF treatments were compared pairwise to evaluate the difference in bacterial communities when the soil nutrition was switched from CF to OF alone (Figure 4C). The OF treatment exhibited a higher RA of Proteobacteria (3.5% higher than CF, Supplementary Table S2) and especially the sub-class Betaproteobacteria. Input of CFs increased the RA of the Gemmatimonadetes phylum, the Gemmatimonadetes class and Gemmatimonadales order, as well as the Acidobacteria phylum (7.6% higher than OF, Supplementary Table S2).

The OF and OIMF treatments were compared pairwise to detect the difference in bacterial communities when the nutrition from pure organic matter was partially substituted with CF (Figure 4D). Similar to the CF treatment, OIMF also enhanced the RA of the Gemmatimonadetes phylum, the class Gemmatimonadetes and the order Gemmatimonadales, and the RA of the Acidobacteria phylum. The application of pure OF increased the RA of Betaproteobacteria.

Quantitative PCR was performed to determine the abundance of the Acidobacteria phylum and the α and β sub-classes of the Proteobacteria affording comparison with results obtained by high-throughput sequencing (Figure 5). Compared with the CF and OIMF treatments, the OF treatment had more abundance of Betaproteobacteria (Figure 5B), which was in agreemtent with the LEfSe results. Acidobacteria tended to be higher in the NNF treatment relative to the other three treatments, but consistantly higher Acidobacteria levels were not found in CF and OIMF treatments compared with the OF treatments by QPCR analysis (Figure 5C), which was inconsistent with the results from the LEfSe analysis. Proteobacteria/Acidobacteria (P/A) ratio is suggested be indicative of higher soil nutrient availability, since the specific primer for the phylum Proteobacteria could not be confirmed, only α- and β- Proteobacteria were determined by QPCR to caculate the P/A ratio. Results showed that the OF treatment tended to have higher P/A ratios than CF and OIMF treatments across all the eight growth stages (Figure 5D). This was consistent with the high throughput sequencing that indicated that the OF treatment had the highest P/A ratio (4.05), while the CF treatment had the lowest ratio (3.64; Figure 5E). Although the RAs of Acdiobacteria by QPCR and by high throughput sequencing were significant (P< 0.001, Pearson correlation analysis), the RAs of α- and β- Proteobacteria by these two methods were not well correlated (P> 0.05).

Redundancy analysis analysis revealed that variation in bacterial community composition was significantly correlated with specific soil properties (Figure 6; for additonal detail see Supplementary Table S3). For instance, Betaproteobacteria and Bacteroidetes were positively correlated with soil SOC content and negatively correlated with soil moisture, while Acidobacteria, Thaumarchaeota, and Chloroflexi tended to prefer high moisture and low SOC content soil conditions. Additionally, as revealed by the RDA biplot and the further permutations test, soil temperature and TN content was the key predictors of bacterial community structure, and soil moisture, SOC, AK, and AP contents also played significant roles. Mantel tests also demonstrated that soil TN content was the factor having greatest influence on bacterial community structure in the wheat rotation, while soil temperature was dominant in the rice rotation (Table 2).

Since soil temperature and TN contents were the strongest drivers of bacterial community variation. These two soil properties were further correlated with the RA of each phylum as examined by correlation models (Figure 7). Most of the bacterial RAs were positively correlated with soil temperature (P< 0.01), including Verrucomicrobia, Acidobacteria, and seven other phyla, while Actinobacteria and Gemmatimonadetes were negatively correlated with soil temperature (Figure 7A). Soil TN was positively correlated (P< 0.01) with the RAs of Proteobacteria and Fusobacteria, while TN contents and the RAs of Acidobacteria and Planctomycetes were negatively correlated (Figure 7B).

In present study, plant growth stage was determined to be the most important factor for explaining the observed variation in soil bacterial community structure. Interestingly, bacterial diversity remained steady over the different stages, which did not support our hypothesis mentioned before. Thus, although bacterial community composition fluctuated over time, changes in RA and species replacements were driving these differences as opposed to species loss or gain.

It is difficult to determine the exact cause of the observed shifts in bacterial communities in response to rotation stage in the wheat-rice cultivation system. Several factors were strongly correlated with plant growth stages, including soil temperature, soil moisture, crop species and growth. Soil temperature varied greatly over the time of this study (from 8.9 to 31.1°C), and this factor was one of the key predictors in the RDA analysis. The effect of temperature could be both direct and indirect. It will have a direct effect on microbial metabolic activity (Schindlbacher et al., 2011) and growth (Bärcenas-Moerno et al., 2009). While numerous indirect effects through crop growth and the decomposition of organic matter (Glanville et al., 2012). The higher 16S rRNA gene copy numbers in May and September corresponding well with the wheat and rice heading stages when the plants had their most vigorous growth (Yoshida, 1981). Soil microbes would be expected to respond differently to fertilizer treatments as plant activity varies through growth stages, and in line with this notion, we observed a significant interactive effect between stage and fertilization treatment on bacterial abundance.

The difference between wheat and rice in the quality and quantity of root exudates could also induce changes in the bacterial community, especially for heterotrophic species (Buyer et al., 2002). Although we did not measure root exudates, previous researchers have reported that root exudates are strongly affected by crop growth stage, and that this relationship in turn influences the succession of bacterial communities (Yang and Crowley, 2000; Duineveld et al., 2001; Xu et al., 2009). In addition, plant species has been acknowledged as an important factor influencing the diversity and structure of microbial communities (Smalla et al., 2001; Kowalchuk et al., 2002), and we also observed difference in soil bacterial communities under wheat versus rice crops, although this impact might be less than pure temporal dynamics.

Flooding of soil between July and September during the rice season had a great influence on the soil bacterial population due to the onset of anoxic conditions (Wang et al., 2015). During this period, oxygen was rapidly depleted in the bulk soil (Liesack et al., 2000), which changed the soil conditions drastically from the wheat ripening (June) to rice tilling (July). These changes would suppress aerobic bacteria while stimulating many previously dormant anaerobic bacteria. Such stimulation was exemplified by the increase in obligate anaerobic bacteria such as Verrucomicrobia (Bergmann et al., 2011).

The fertilizer treatments significantly impacted the soil chemical properties in our study, thereby affecting the resident bacterial communities. The impacts of different fertilizations on bacterial communities have been reported by previous studies on paddy soil (Yuan et al., 2013; Dong et al., 2014). However, the degree to which such effects persist throughout the whole wheat and rice rotation had not previously been explored, and we observed consistent impact of fertilizer regime despite the backdrop of large temporal changes in microbial community structure.

Concerning the influence of nitrogen fertilizer input, NNF and CF treatments clearly differed, with higher RA of Cyanobacteria in the plots that only received phosphorus and potassium fertilizers (NNF). This trend is consistent with the known behavior of Cyanobacteria which have the ability to fix nitrogen biologically, and this activity is inhibited by inorganic nitrogen fertilizers inputs (Irisarri et al., 2001; Ramirez et al., 2010; Grandy et al., 2013).

Compared with OIMF treatment, CF improved the soil ammonium and nitrate content and resulted in higher RA of members of the ammonia-oxidizing archaeal Thaumarchaeota (Spang et al., 2010; Pester et al., 2011). This group has previously been correlated with soil ammonium contents (Hu et al., 2013), corresponding with our result.

Additions of OFs tend to increase soil organic C and N contents whereas CFs have a much smaller impacts (Lal, 2004). In addition, applications of N as urea, as in the current experiment, lead to acidification of soil (Guo et al., 2010) through the action of ammonium oxidizing bacteria and archaea that produce protons. Hence differences in the microbial populations between organic and inorganic fertilizer treatments are to be expected (Chu et al., 2007; Ge et al., 2008; Zhang et al., 2012). In the present study, higher RAs of Acidobacteria but lower RAs of Proteobacteria were observed in chemical fertilization treatments as compared to organic fertilization treatments. Acidobacteria generally are associated with oligotrophic, slow-growing properties as typical of a K-selected life strategy. Thus, this group is often abundant in soils with low resource availability and low organic C conditions, and it is negatively correlated with soil pH (Fierer et al., 2007). In contrast, Proteobacteria, including the α- and β- subclasses, generally exhibit copiotrophic attributes, with higher abundances nutrient-rich high C soils (Fierer et al., 2007; Newton and McMahon, 2011; Zhou et al., 2015). Correlation analyses (Supplementary Table S3) supported this ecological distinction, with Acidobacteria negatively and α-/β- Proteobacteria positively correlated with SOC and TN content. The ratio between Proteobacteria and Acidobacteria is thought to be an indicator of the nutritional status of the soil ecosystem (Smit et al., 2001; Torsvik and Øvreås, 2002), so it is interesting that this ratio was higher in the treatments receiving pure OF. QPCR also showed that the abundance of Alphaproteobacteria and Betaproteobacteria classes tended to be higher in the OF treatment. But higher Acidobacteria abundance in CF and OIMF treatments, shown in LEfSe, was not detected by the QPCR. This discrepancy with the LEfSe results may be due to the specificity of the QPCR Acid31F primer, which may not amplify all members of the Acidobacteria (George et al., 2009; Lee and Cho, 2011). In addition, the low correlation coefficient when comparing QPCR and high throughput sequencing results for the α- and β- Proteobacteria, highlights the fact that caution should be exercised when interpreting such PCR-derived results.

It is difficult to know exactly why measures of the bacterial α-diversity were unaffected by fertilization regimes in our study. Clearly, community composition is altered, which alters the RA of different taxa within our experimental study site. Apparently, a range of different selection pressures is at play, i.e., there is not a simple strong selection for a limited number of taxa that thereby lowers diversity.

The grain yields in OF treatment were as low as in the no nitrogen input treatment (Supplementary Figure S1) even though OF treatment tended to host more copiotrophic bacterial groups. OF and NNF treatments had a similar extremely low nitrogen availability, thought as the direct factor influencing the crop yield (De Ponti et al., 2012), explaining the low yields in these two treatments. However, one would not expect to identify direct impacts of soil bacterial community size or composition with crop yields. But it is valuable to build up a picture of variations in soil microbial community structure so that, in due course, this can be assessed as a part of sustainability considerations.

It should also be noted that we sampled bulk soil which is inhabited by free-living bacteria such as Cyanobacteria and Proteobacteria were inhabited (Pii et al., 2015), so the results will be biased toward such organisms rather than those concentrated in the rhizosphere which may show different trends in response to fertilization or plant growth (Ai et al., 2012; Geisseler and Scow, 2014). Thus, future studies comparing rhizosphere and bulk soils will be necessary to gain full understanding of the impacts of fertilizer regime and rotation stage on soil microbial communities in such wheat-rice systems.