Long-term continuously monocropped peanut significantly changed the abundance and composition of soil bacterial communities

Field experiment and soil sampling

The field experimental site was set up in Laixi experimental farm at the Shandong Peanut Research Institute, Qingdao, China (36°50′N, 120°31′E). A gradient of consecutive monocropped peanut experiment plots that we previously developed was used. Independent pools were applied, with the size of each pool being 4 m long, 1.5 m wide, and 1 m high. We collected soil from the plow layer of the cultivated land that had previously been planted with a wheat–maize rotation, thoroughly mixed, and then added to the pools. In order to construct the gradient of continuous monocropping plots, sweet potatoes were rotated with peanut in some plots during the experiment. Peanuts were planted in May each year and harvested in October. All field management, including the planting pattern and the use of water and fertilizer were consistent among different pools. After harvest, the plots lay fallow until the next planting season. The distance between rows was 16.5 cm and the ridge width was 70 cm, with two peanut seeds per hole. Before planting, about 300 kg ha⁻¹ of urea, 750 kg ha⁻¹ of calcium magnesium phosphate, and 225 kg ha⁻¹ of potassium chlorate were applied as fertilizer. The weeds were controlled using glyphosate (41% active ingredient in 3 L ha⁻¹) before planting. Plastic film mulching was applied during cultivation. In addition, to minimize the influence of other factors including variety on the measured indexes, the planting position of the same peanut variety was fixed during the continuous cropping period. The longest-running pools have now been continuously monocropped with peanuts for 16 years.

According to our previous study (Jiao et al., 2015), we selected four peanut varieties with distinct responses to monocropping to analyze their soil bacterial community structure in plots monocropped for 1, 2, 11, and 12 years. The four varieties included two large fruit varieties, Huayu 917 and Huayu 50, and two small fruit varieties, Huayu 26 and Huayu 20. Based on previous testing of their yield-related indexes, Huayu 20 was the most sensitive to long-term monocropping, Huayu 917 was intermediate, and Huayu 26 and Huayu 50 were tolerant. There were three replicates for each treatment. The bulk soils were collected at the full-bloom stage and five randomly selected replicate test plants were used for each sample. The soils around individual plants were sampled using a soil probe (1.5 cm diameter) at 5–10 cm soil depth, and at distances of 3–5, 8–10, and 13–15 cm from the main root. The root zone soils were then mixed together. These distances were chosen because the roots at this depth and these distances were relatively abundant and the soil microbial community around the peanut root system could be well characterized (Chen et al., 2012). Physico-chemical properties of pH, organic matter, and available nitrogen (N), phosphorus (P), and potassium (K) in the soil samples were determined using routine methods (Lu, 1999). Samples Y1, Y2, Y11, and Y12 were for monocropping of 1, 2, 11, and 12 years, respectively; H20, H26, H50, and H917 were samples for four peanut varieties with distinct responses to monocropping.

DNA extraction, PCR, and high-throughput sequencing

For each soil sample, total DNA was extracted with the Power Soil DNA Isolation Kit (MoBio Laboratories Inc., Carlsbad, CA, USA) following the manufacturer’s instructions. The concentration and purity of the DNA was checked by electrophoresis on 1.0% (w/v) agarose gels. Soil bacterial communities were analyzed with amplicon sequencing on an Illumina HiSeq platform. Two primers, 515F (5′-GTGCCAGCMGCCGCGGTAA-3′) and 907R (5′-CCGTCAATTCCTTTGAGTTT-3′), for the V4–V5 region of the 16S rRNA gene were applied. Both the forward and reverse primers had a 6-bp barcode unique to each sample. Each soil sample was independently amplified in triplicate, and the triplicate PCR reaction products for each sample were pooled and then all products were purified using a GeneJET™ gel extraction kit (Thermo Scientific, Waltham, USA). The purified amplicons from each sample were pooled in equimolar concentrations and the mixture was then sequenced on an Illumina HiSeq platform at Novogene Co. Ltd, Beijing, China.

Sequencing data processing and statistical analysis

The paired-end reads from the original DNA fragments were merged with FLASH software (Edgar, 2013). The QIIME 1.7.0 software package was used to quality-filter and process the raw sequencing reads (Caporaso et al., 2010). The UCLUST method was used to delineate the operational taxonomic units (OTUs) at a threshold of 97% identity (Caporaso et al., 2010). A representative sequence for each OTU was taxonomically classified using the Silva Database (https://www.arb-silva.de/) based on Mothur algorithm (Quast et al., 2013). The OTU abundance data were normalized using a standard sequence number corresponding to the sample with the least sequences. The output normalized data were used for the subsequent analyses.

To assess the bacterial alpha-diversity of each sample, analyses including rarefaction curves, Shannon index, and Chao1 index were performed using QIIME (Version 1.7.0) and displayed with R software (Version 2.15.3). Beta diversity analyses were used to estimate sample differences in bacterial community compositions. In order to analyze the influence of peanut monocropping on soil bacterial community composition, heatmap analysis was performed with R software. Cluster analysis was preceded by principal component analysis (PCA), which was applied to reduce the dimensions of the original variables using the FactoMineR package and ggplot2 package in R software (Version 2.15.3). Principal coordinate analysis (PCoA) was performed to obtain principal coordinates and visualize the complex multidimensional data. The PCoA analysis was performed using the WGCNA, stat, and ggplot2 packages in R software (Version 2.15.3). Unweighted pair-group method with arithmetic means (UPGMA) clustering was performed as a type of hierarchical clustering method to interpret the distance matrix using average linkage and was conducted using QIIME software (Version 1.7.0). Taxa with significantly diverse abundance with long-term monocropping at different levels were further investigated using the Metastats method and t-test using R software (Version 2.15.3). Tax4Fun functional prediction was achieved by the nearest neighbor method based on the minimum 16S rRNA sequence similarity. The data set was deposited in the NCBI-Sequence Read Archive with the submission accession number PRJNA559575.

Physico-chemical properties of soil

Physico-chemical properties of pH, organic matter, and available N, P, and K in the peanut monocropping soil were tested and analyzed (Table S1). The soil pHs in samples from 2-year monocropped plots (8.30–8.59) were slightly higher than those from 1-year monocropped plots (7.20–8.05), but were obviously lower in samples from 11- and 12-year monocropped plots (6.46–6.92). The contents of available P were lower in long-term (15.12–21.12 mg kg⁻¹) than in short-term monocropping samples (39.37–48.58 mg kg⁻¹). The contents of available K were also lower in long-term (59.98–86.13 mg kg⁻¹) than in short-term monocropping samples (99.42–120.67 mg kg⁻¹). However, contents of organic matter and available N showed no obvious changes with long-term peanut monocropping.

Characteristics of sequencing data

In order to compare the soil bacterial community structure and composition of the four peanut varieties to short- and long-term monocropping periods, a total of 48 16S rRNA gene libraries were analyzed using high-throughput sequencing. The sampling variables, including four peanut varieties and four monocropping periods, are shown in Table 1. A total of 2,541,421 quality-filtered sequences obtained from the 48 samples, ranging within 37,605–68,117 with an average length of 373 bp, resulted in identification of a total of 9028 OTUs applying a 3% sequence dissimilarity cutoff. After data normalization, 1,805,040 quality-filtered sequences affiliated with 8856 OTUs were obtained. The bacterial complexity of the 48 samples was estimated on the basis of alpha-diversity (OTU number, Chao1 index, and Shannon index) and showed relatively higher bacterial diversity in samples from the 1- and 2-year compared to the 11- and 12-year monocropped plots (Table 1). Rarefaction curves (Fig. S1) and indices of richness and diversity including all samples’ Chao1 and Shannon indexes (Table 1) tended to approach a saturation plateau, indicating that the majority of bacterial diversity was recovered by the surveying effort.

Bacterial community structure and composition

In total, 99.73% of the identified 8856 OTUs were related to bacteria, and the other 24 OTUs accounting for 0.15% of the total quality-filtered sequences were affiliated with Archaea. There were 41 bacteria phyla detected across all samples and 0.73% of bacterial sequences were unclassified at the phylum level. The top 10 phyla accounted for 95.49% of the total sequences and 79.67% of the total OTUs. The dominant phyla across all samples were Proteobacteria, Actinobacteria, Acidobacteria, and Chloroflexi, accounting for 30.71, 21.59, 18.04, and 7.69% of the total sequences, respectively, and 26.61, 7.43, 6.49, and 10.77% of the total 8856 OTUs, respectively. The phyla Gemmatimonadetes, Firmicutes, Planctomycetes, Nitrospirae, Bacteroidetes, and Thermomicrobia were relatively abundant and diverse, accounting for 4.66, 3.75, 3.13, 2.96, 1.62, and 1.36% of the total sequences, respectively, and 3.42, 4.41, 14.76, 0.61, 3.97, and 1.19% of the total OTUs, respectively (Fig. S2).

At the class level, 92 bacteria taxa were identified. The dominant 19 taxa (relative abundance > 1%) accounted for > 88.43% of the sequences in each sample. The most abundant taxa (relative abundance > 5%) were Acidobacteria (16.19%), Alphaproteobacteria (10.11%), Betaproteobacteria (9.48%), Thermoleophilia (9.22%), Actinobacteria (7.31%), Deltaproteobacteria (5.77%), and Gammaproteobacteria (5.29%). The other relatively abundant taxa included Gemmatimonadetes, Bacilli, Nitrospira, MB-A2-108, Planctomycetacia, Acidimicrobiia, Anaerolineae, KD4-96, Holophagae, Thermomicrobia, Sphingobacteriia, and TK10, each accounting for 1.05–4.66% of the total sequences (Fig. 1A).

At the genus level, 494 taxa were detected, accounting for 31.41% of the total sequences; and the top 50 taxa accounted for 75.04% of the total identified sequences at the genus level. The dominant genera (relative abundance > 1%) were Bacillus, Gaiella, Mizugakiibacter, and Streptomyces, accounting for 2.22, 1.79, 1.35, and 1.01% of the total sequences, respectively. Relatively abundant genera were Sphingomonas (0.90%), Haliangium (0.89%), Gemmatimonas (0.80%), Bryobacter (0.76%), Nocardioides (0.74%), Arthrobacter (0.71%), Steroidobacter (0.68%), Solirubrobacter (0.66%), and Microvirga (0.57%) (Fig. 1B).

Bacterial community variation at different levels with monocropping time

In order to illuminate the bacterial community variation across long-term monocropping time of peanut, UPGMA, PCA, and PCoA were used to cluster the bacterial communities within soil samples. Overall, the UPGMA tree showed that replicate samples were grouped together, and soil samples with similar monocropping time but different peanut varieties were clustered into a group (Fig. 2). Both PCA and PCoA also showed obvious clustering of bacterial communities based on monocropping time (Fig. 3, Fig. S3). Samples from 1- and 2-year monocropping plots were grouped into one class, and samples from 11- and 12-year plots were grouped into another. However, the cluster distance was shorter among samples from 1- and 2-year than among samples from 11- and 12-year monocropping plots. These analyses indicated that bacterial communities were significantly influenced by the monocropping period of peanut, and showed more obvious changes under long-term monocropping.

The heatmap analyses showed that soil bacterial community composition was significantly diverse at different levels between short- and long-term monocropping. At the phylum level, among the 41 bacteria taxa identified, nine taxa showed an obvious decrease in abundance with monocropping time and five taxa showed an obvious increase in abundance (Fig. S4, Fig. 4). The average abundance of Acidobacteria accounted for 20.30 and 19.83% in the 1- and 2-year libraries, respectively, and decreased to 15.74 and 16.28% in the 11- and 12-year libraries. The phyla Planctomycetes, Nitrospirae, and Bacteroidetes were relatively abundant in the 1- and 2-year libraries, accounting for 3.12, 2.96, and 1.62% of the total sequences, respectively, and correspondingly decreased in the 11- and 12-year libraries by 20.23, 17.47, and 27.81%. Additionally, the taxa Armatimonadetes, Latescibacteria, JL-ETNP-Z39, Thermotogae, and Caldiserica, for which relative abundances were <1%, also had clear decreasing tendencies. In general, decreases in their abundance were accompanied by decreases in their diversity (Table S2). In contrast, the abundance of some taxa increased with monocropping time. The average abundance of Gemmatimonadetes, Firmicutes, and Elusimicrobia accounted for 3.91, 2.65, and 0.16% in the short-term monocropping libraries, respectively, and correspondingly increased to 5.40, 4.84, and 0.21% for the long-term. The other taxa with increasing abundance all had relatively low abundance, accounting for 0.01–0.04% of the total sequences, and included Parcubacteria and Chlamydiae. The diversity of the three relatively abundant taxa (Gemmatimonadetes, Firmicutes, and Elusimicrobia) did not show obvious changes with monocropping time, but that of the other low-abundance taxa tended to increase with monocropping time (Table S2).

Analysis of the microbial community at lower taxonomic levels (genus or species) can provide better phylogenetic resolution than at higher taxonomic levels (Ho et al., 2016). At the genera level, the top 50 most-abundant taxa accounted for 75.04% of the total identified sequences that were affiliated with 494 taxa. The abundance distributions of the top 50 genera in each sample were presented in a heatmap (Fig. 5). Some taxa showed obvious changes with monocropping time (Figs. S5 and S6). The dominant genera Bacillus (2.22%) and Mizugakiibacter (1.35%) increased their average abundance from 1.49 and 0.67%, respectively, in the 1- and 2-year libraries to 2.96 and 2.02% in the 11- and 12-year libraries. However, the number of OTUs affiliated with the genera Bacillus (13–14 OTUs) and Mizugakiibacter (four OTUs) in each sample was relatively low and changed little. Seven genera (Gemmatimonas, Marmoricola, Candidatus Solibacter, Tumebacillus, Jatrophihabitans, Sporosarcina, and Pseudolabrys) accounted for 0.17–0.80% of the total sequences, and increased their average abundance by >90%. The genera Nocardioides, Solirubrobacter, and Reyranella accounting for 0.21–0.74% of the total sequences also showed increased abundance with monocropping time, with increase rates of 26.44–57.20%. Some taxa showed an opposite trend, with abundance decreasing with monocropping time. They were Steroidobacter, Candidatus Entotheonella, Nitrosospira, Pirellula, Piscinibacter, Ramlibacter, Lysobacter, Skermanella, and Planctomyces, which accounted for 0.20–0.68% of the total sequences, and the abundance decrease rate ranged within 12.86–51.13%. Overall, abundance of these genera began to decline in the 2-year monocropping libraries, and decreased more for the long-term libraries.

Most sequences, accounting for 91.54% of the total sequences, were not classified at the species level, but some taxa among the identified species also obviously changed with monocropping time (Fig. S7, Fig. 6). Among the top 10 most-abundant species, Bacillus aryabhattai and Bacillus funiculus accounting for 0.58 and 0.42% of the total sequences, respectively, showed increased abundance in long- compared to short-term monocropping libraries. The increase rate of these two taxa reached 116.23 and 177.70%, respectively. The relatively abundant taxa of Bacillus luciferensis, Tumebacillus ginsengisoli, Bacillus decolorationis, and Streptacidiphilus luteoalbus, accounting for 0.01–0.10% of the total sequences, increased their relative abundances by > 118%. The taxa Paenibacillus alginolyticus and Streptosporangium violaceochromogenes also presented increased abundance with increase rates ranging within 52.07–84.49%, but the change trends differed among peanut varieties. The species Paenibacillus alginolyticus and Streptosporangium violaceochromogenes showed no obvious changes in abundance in Huayu 917 and Huayu 50 libraries, respectively. In contrast to the above taxa, Chitinophaga ginsengihumi, Lysobacter yangpyeongensis, and Phyllobacterium myrsinacearum, accounting for 0.03–0.05% of the total sequences, showed obvious decreases in abundance and the abundance decrease rate ranged within 32.57–61.59%.

Prediction of bacteria community functions

The predominant category of the predicted functional genes was affiliated with metabolism (48.71%), followed by genetic information (21.03%), environmental information processing (12.89%), and cellular processes (7.62%) (Fig. S8). The PCA was used to cluster the predicted functional pathway within soil samples at three KEGG levels. The analyses at KEGG levels 2 and 3 and KEGG ortholog (KO) level all showed similar clusters to that of bacterial community structure (Fig. 7, Fig. S9). Samples from 1- and 2-year monocropping plots were grouped into one class, and samples from 11- and 12-year plots were grouped into another.

Bacterial community functions also presented significant changes at different pathway levels across long-term monocropping. For KEGG pathways at KEGG level 2 that are involved in carbohydrate metabolism, endocrine system, excretory system, nucleotide metabolism, transport and catabolism, transcription, biosynthesis of other secondary metabolites, and aging showed decreasing trends with long-term monocropping, but pathways affiliated with membrane transport and cellular community prokaryotes showed increased abundance (Table S4). There were 390 KEGG pathways at KEGG level 3 predicted in our libraries. The abundant pathways (>1.0%) that presented obvious changes in abundance with monocropping period were further examined. The pathways involved in DNA repair and recombination proteins, purine metabolism, transfer RNA biogenesis, exosome, amino acid related enzymes, pyrimidine metabolism, mitochondrial biogenesis, ribosome, oxidative phosphorylation, carbon fixation pathways in prokaryotes, aminoacyl tRNA biosynthesis as well as alanine, aspartate, and glutamate metabolism showed decreased abundance in long- compared to short-term monocropping samples. In contrast, pathways involved in ribosome biogenesis, messenger RNA biogenesis, RNA degradation, quorum sensing, glyoxylate and dicarboxylate metabolism, and lipid biosynthesis proteins were enriched in long-term monocropping samples (Table S4).