Responses of rhizosphere soil bacteria to 2-year tillage rotation treatments during fallow period in semiarid southeastern Loess Plateau

Description of experimental site

A 5-year field study was established at the Wenxi Dryland Wheat Agriculture Station (35°20′N, 111°17′E), Shanxi Province, China to observe the effects of the different tillage practices of subsoiling (SS), deep-tillage (DT), and no-tillage (NT) on wheat yield and water use efficiency between 2009 and 2013. The experiment was conducted from July 2013 to June 2015. The site represents the typical semiarid climate of the southeastern Loess Plateau and the experimental field soil was classified as silty clay loam (Chinese soil taxonomy). Winter wheat (Triticum aestivum L.) was sown in early October of every year and harvested in mid-June of the following year. Rain-fed agriculture is popular in this area due to scarce rainfall and no irrigation; the average annual precipitation of the site is 450–630 mm with no irrigation performed at the site of the study. More than 60% of the total rainfall fell during the fallow period (July to September). A brief description of the experimental site is provided in Table 1.

Experimental design

Three tillage methods were applied during the fallow season. These were: (1) subsoiling to 30–40 cm (SS); (2) deep tillage to 25–30 cm (DT) and (3) no-tillage (NT) during two years with four tillage combinations SS-SS, DT-DT, SS-DT, and NT-NT (control). The wheat crop was harvested on June 8th, 2013 and June 10th, 2014 and wheat stubble with a height of 30–40 cm was left in the field to preserve water in the soil. The three tillage methods were introduced ten to fifteen days after the harvest using two different ploughing machines. The experiment was laid out in a randomized block design with a factorial arrangement and three replications. The total area of the experimental field was 12 × 666.7 m².

Soil sampling

Winter wheat was harvested on June 14th, 2015 and rhizosphere soil samples were collected from depths of 0–20 cm, 20–40 cm, and 40–60 cm from four tillage rotation treatments. Three subsamples were collected from each sampling site. The whole plant was collected, including additional soil from a 15 cm radius around the base of the plant. Soil samples were placed into sterile petri plates and the roots were removed. Non-rhizosphere soil was dislodged by shaking the roots gently. The rhizosphere soil that remained on the roots was collected using a flame-sterilized tweezer and sieved through a 2-mm mesh to eliminate large rocks and roots (Yang et al., 2016). Each composite soil sample was homogenized and stored at −20  °C for less than 24 h before DNA extraction. Thirty-six soil samples (four tillage rotation treatments, three replicate samples per tillage rotation treatment) were analyzed.

DNA extraction of the soil samples

Total bacterial genomic DNA samples were extracted from the soil using Fast DNA SPIN extraction kits (MP Biomedicals, Santa Ana, CA, USA), following the manufacturer’s instructions. The samples were stored at −20 °C until further analysis. The quantity and quality of extracted DNA was measured using a NanoDrop ND-1000 spectrophotometer (Thermo Fisher Scientific, Waltham, MA, USA) and agarose gel electrophoresis, respectively.

16S rDNA amplicon pyrosequencing of soil samples

PCR amplification of the V3-V4 regions of bacterial 16S rDNA was performed using forward primer 338F (5′-ACTCCTACGGGAGGCAGCA-3′) and reverse primer 806R (5′-GGACTACHVGGGTWTCTAAT-3′), with sample-specific 7-bp barcodes incorporated into the primers for multiplex sequencing. The PCR components contained 5 µl of Q5 reaction buffer (5 ×), 5 µl of Q5 High-Fidelity GC buffer (5 ×), 0.25 µl of Q5 High-Fidelity DNA Polymerase (5U/ µl), 2 µl (2.5 mM) of dNTPs, 1 µl (10 uM) each of Forward and Reverse primers, 2 µl of DNA Template, and 8.75 µl of ddH₂O. Thermal cycling was conducted with the initial denaturation at 98 °C for 2 min, followed by 25 cycles of denaturation at 98 °C for 15 s, annellation at 55 °C for 30 s, an extension at 72 °C for 30 s, and a final extension at 72 °C for 5 min. The PCR products were then purified with Agencourt AMPure Beads (Beckman Coulter, Indianapolis, IN) and quantified using the PicoGreen dsDNA Assay Kit (Invitrogen, Carlsbad, CA, USA). Amplicons were pooled in equal amounts following the individual quantification step and pair-end 2 ×300 bp sequencing was performed using the Illumina MiSeq platform with the MiSeq Reagent Kit v3 at Shanghai Personal Biotechnology Co. Ltd (Shanghai, China).

Illumina MiSeq of the bacterial communities from soil samples

Species of the bacterial communities were assayed with a genotypic fingerprinting approach using the Illumina MiSeq technique. Total DNA concentrations that met the requirements of amplification were used for subsequent analysis (Illumina MiSeq high-throughput sequencing, Enterprise Group Personal Biotechnology Co., Ltd., Shanghai, China).

Sequence analysis of the soil samples

Raw sequencing data were filtered and saved in the paired-end FASTQ form. Sets of four lines were regarded as sequencing reads in the FASTQ file. Low-quality sequences were filtered out according to the following criteria: sequence lengths of < 150 bp, average Phred scores of < 20, and sequences that contained ambiguous bases, mismatched bases of 5′-terminal > 1 or mononucleotide repeats of > 8 bp. The remaining high-quality sequence reads were then assembled using FLASH software (version 1.2.7) (Magoč & Salzberg, 2011). The raw sequencing reads matching the 7-bp barcodes were assigned to individual samples and identified as valid sequences which were further filtered using the Quantitative Insights Into Microbial Ecology method (Qiime, version 1.9.0). Chimera sequences were removed using the UCHIME method of Mothur (version 1.31.2) (Schloss et al., 2009; Caporaso et al., 2010; Edgar, 2010).

The remaining high-quality sequences were clustered into operational taxonomic units (OTUs) at 97% sequence identity using the UCLUST method of Qiime (Blaxter et al., 2005; Edgar, 2010), with a representative sequence selected from each OTU using default parameters. Taxonomic information of each representative OTU sequence was then obtained using the UCLUST method of Qiime by comparing sequence databases (Greengene, Release 13.8) (Altschul et al., 1990; De Santis et al., 2006). A streamlined OTU list was obtained for subsequent analysis by removing OTUs with an abundance value less than 0.001% of the total sequence number (Bokulich et al., 2013). A phylogeny tree for representative sequences of OTU was constructed using the FastTree method (Price, Dehal & Arkin, 2009).

Bioinformatic analysis of soil samples

Sequence data bioinformatics analyses of soil samples were performed using Qiime and R packages (Version 3.2.0). According to the species abundance of each soil sample in the OTU list, four commonly used biodiversity indices were calculated utilizing the summary single command of the Mothur software, including two community richness indices (Chao1 and ACE) and two community diversity indices (Shannon index and Simpson index). The ecological indices, Chao1 and ACE, were used to estimate the total number of species (Pitta et al., 2010); higher values indicated greater community richness. The Shannon and Simpson indices were used to reflect alpha diversity (Shannon & Weaver, 1998; Mahaffee & Kloepper, 1997). A higher Shannon index but lower Simpson index value indicated greater community diversity. According to the results from the OTU list, the species proportion of each sample were obtained at different taxonomic levels (phylum, class, order, family, and genus, respectively) using Qiime (De Santis et al., 2006). The species proportion was visualized on an intuitive bar chart reflecting the community structure of the sample at different taxonomic levels, for example, abundance bar charts at the phylum and genus levels. Taxon abundances at both levels were compared statistically among samples or groups by Qiime (Oberauner et al., 2013). Principal component analysis (PCA) was conducted based on the genus-level compositional profiles using the vegan package of R (Ramette, 2007). The classification information on the genus level was clustered and drawn in Heatmap with different colors according to the abundance similarity and variation of samples and groups. The prediction of functional species and relevant statistical analyses were carried out using PICRUSt according to the abundance at the OTU level (Langille et al., 2013).

Statistical analysis

The SAS package (SAS Institute Inc., North Carolina, USA) was used to conduct analysis of variance tests (ANOVA). Differences between the means of treatments were compared using the least significant difference (LSD) approach. Differences were considered statistically significant at the probability level of 0.05.

Changes in diversity of the rhizosphere soil bacterial community under different tillage rotations

The richness (Chao1 and ACE) and diversity (Shannon and Simpson) indices were examined to evaluate the abundance and diversity of rhizosphere bacterial communities in different soil layers under various tillage rotations (Table 2). The results showed that four indices decreased with an increase in the soil depth, with the most significant decrease found in the 40–60 cm soil layer. However, the overall effect of the soil layer on bacterial richness and diversity was found to be non-significant (p > 0.05), whereas the effect of tillage rotations on richness and diversity was significant (p < 0.05), indicating that the bacterial communities were altered under different tillage rotations.

In the soil layer of 0–20 cm, the Shannon and Simpson indices were significantly higher under all tillage rotations than the control, whereas, no significant difference was observed for the richness (Chao1 and ACE) indices. In the 20–40 cm soil layer, the richness indices (Chao1 and ACE) of NT-NT and SS-SS were significantly higher than that of SS-DT and were not significantly higher than DT-DT. This indicated that tillage rotations inhibited bacterial diversity and richness at depths of 20–40 cm. There was no statistical difference in the richness and diversity indices by different tillage rotations in the 40–60 cm soil layer. Thus, tillage rotation has the most significant effect on the diversity of the bacterial community in the 0–40 cm soil layer.

Impact of tillage rotations on rhizosphere soil bacterial community structure and composition

The dominant bacterial phyla found in the rhizosphere soil of three soil layers under all tillage rotations were Proteobacteria, Actinobacteria, Acidobacteria, Planctomycetes, Bacteroidetes, Gemmatimonadetes, Frimicutes, Chloroflexi, Nitrospirae, and Verrucomicrobia (Fig. 1). The relative abundances of Proteobacteria and Chloroflexi were highest in the SS-DT treatment compared with NT-NT; Acidobacteria, and Gemmatimonadetes were highest in the SS-SS treatment; Planctomycetes and Verrucomicrobia were relatively more abundant in the DT-DT treatment; Actinobacteria, Bacteroidetes, Frimicutes, and Nitrospiraewere highest in the NT-NT treatment. The abundance of other bacteria were extremely low (close to 0%) under all four tillage rotations.

Tillage rotation improved the abundances of Proteobacteria, Actinobacteria and Frimicutes and inhibited Planctomycetes, Acidobacteria, Gemmatimonadetes, Chloroflexi and Verrucomicrobia as the tillage depth increased. The smallest abundance value of Bacteroidetes was in the 20–40 cm soil layer. The largest abundance value of Nitrospirae was in the 20–40 cm soil layer. This may be due to aerobic nature of these bacteria phyla and the better soil aeration in the 0–40 cm soil layer than in 40–60 cm.

The top fifty bacteria genera were selected and listed for their compositions and abundances for the four tillage patterns under 0–20, 20–40 and 40–60 cm soil layers (Fig. 2). The diversity, richness, and evenness of the bacteria community structure gradually decreased at the genus level with an increased depth in the soil layer. The relative abundance values of Phenylobacterium, Lactococcus, and Cohnella were significantly affected by the tillage pattern and soil layer. For example, the abundance value of Phenylobacterium was less than 3% in the 0–20 cm plow soil layer and these values increased to 26.66% (SS-DT), 15.64% (DT-DT) and 4.88% (NT-NT) in the 20–40 cm plow soil layer, respectively. Its abundance values increased to 22.31% (SS-SS), 36.45% (SS-DT), 36.75% (DT-DT) and 17.66% (NT-NT) in the 40–60 cm plow pan soil layer, respectively. These results were also reflected in Lactococcus and Cohnella. Similarities among the bacteria communities were found at the genus level for SS-SS, NT-NT, DT-DT and SS-DT.

Contrast of soil bacterial community structures under different tillage rotations using PCA and heatmap analysis

Principal component analysis (PCA) was performed to assess similarities between the bacterial communities at the genus level in the different soil samples. In the 0–20 cm soil layer, 88% of the bacterial community composition was represented by the two axes of the PCA diagram (Fig. 3). According to PC1, the composition of the soil bacterial communities in the 20–40 cm depth showed the same effect under the three tillage rotations, and was contrary to the results of NT-NT (Fig. 4). SS-DT and DT-DT were clustered into one group, while SS-SS and NT-NT were clustered into another. SS-SS and SS-DT had the same effect as NT-NT on the composition of soil bacterial communities, according to PC2. The greatest differences in the bacterial community composition for all tillage rotations were found in the 40–60 cm soil layer (Fig. 5).

Heatmap analysis was conducted at the genus level in order to aggregate the high and low abundance of the species. The tree topologies reveal a clear division of two clusters in the three soil layers, consistent with the clusters found in PCA. This implies that approximately the same bacterial communities existed in two tillage rotations in each cluster. The two SS treatments were clustered together in the 0–20 cm plow soil layer (Fig. 6), indicating that the effect of SS-SS on rhizosphere bacteria abundance is more significant than the other treatments. However, DT-DT was clustered with NT-NT, indicating that it has little effect on the abundance of bacterial species communities. The identified bacteria genus classifications were divided into two greater clusters, which could be further divided into four or five clusters. Twenty-one aerobic bacteria (including Balneimonas, Phenylbacteria, Rubrobacter, Flavisolibacter, Bradyrhizobium, Pirellula, and Kaistobacter) were clustered together in the 0–20 cm plow soil layer, and the relative abundance of aerobic bacteria in the SS-DT was higher than that in the SS-SS. In the 20–60 cm soil layer, the DT treatments (SS-DT and DT-DT) clustered together (Figs. 7 and 8), indicating a similar community structure between the two treatments. The community structure under SS-SS was similar to that of NT-NT. The aerobic bacteria (including Balneimonas, Peanibacillus, Phenylbacteria, Cohnella, Lactococcus, and Solibacillu s) were clustered together, and the relative abundance of aerobic bacteria in the SS-DT was higher than that in the DT-DT.

Function predictive bacterial functions under the tillage rotations

To determine the functions of the soil bacteria communities under the four tillage patterns, based on the OTUs abundance in the 16S rDNA sequences, the functional genes composition, including species and abundances, were predicted and statistical analysis was performed using PICRUSt (Figs. 9, 10 and 11). Soil bacteria, as decomposers in ecosystems, are involved in cellular processes, environmental information processing, genetic information processing, the promotion of human diseases, and metabolism and organismal systems. Bacteria played a prevalent role in metabolism, environmental information processing, and genetic information processing. Of the bacteria involved in metabolism, the largest group was related to amino acid, carbohydrate, and energy metabolism. The majority of the bacteria responsible for environmental information processing were involved in membrane transport and signal transduction. Among the bacteria involved in genetic information processing, there was an advantage among groups involved in genetic information replication and repair. The abundance of rhizosphere soil bacterial metabolism in DT-DT was greatest in the 0–20 cm plow soil layer but was greatest in SS-DT in the 20–60 cm plow soil layer. However, the bacterial diversity and composition structures were similar for treatments causing higher soil perturbation, (SS-DT and DT-DT) which broke down the population balance between aerobic and anaerobic bacteria; SS-SS and NT-NT were two less disruptive tillage rotation methods that preserved the soil’s biological integrity. The NT-NT treatment may lead to soil compaction, particularly in the 20–40 cm layer. Our results suggested that SS-SS was the most effective tillage rotation practice.