Soil Properties
Table 1 presents the chemical characteristics of the soil samples. Although the chemical properties showed no statistical difference among the three treatments under the tobacco-rice rotation system; the SOC, TN, and AP contents increased in the soils that undergone OM and NPK fertilizer treatment to different degrees than in those that did not receive fertilizer treatment (NF). However, the soil chemical properties were significantly different among the three treatments under the continuous tobacco system, with significantly lower pH and C/N values and higher SOC, TN, AP, and AK contents in the OM-treated soils than in NF- and NPK-treated soils. Additionally, the chemical properties of the soil under the same fertilizer treatments and two different tillage systems were also analyzed. Compared with those in the soils under the tobacco-rice rotation system, the pH and AK in the soils of the three treatments under the continuous tobacco system were all higher (P < 0.05). Compared with RNPK fertilization, CNPK fertilization resulted in reduced SOC and TN contents in the soil (P < 0.05), whereas COM led to an increased AP (P < 0.05). PERMANOVA showed that the variations in the soil chemical properties could be explained by the tillage system (R2 = 0.25, P < 0.001) and fertilization practice (R2 = 0.35, P < 0.001). Among these chemical properties, the soil pH and AK were more significantly affected by the tillage system than by fertilization (R2 = 0.59, P < 0.001; R2 = 0.17, P = 0.007) (R2 = 0.58, P < 0.001; R2 = 0.17, P < 0.001) (Table S3). In conclusion, the soil nutrients were influenced by the tillage system and fertilization. The application of compost pig manure could improve soil fertility to varying degrees under a continuous cropping system.
Table 1
Soil biochemical properties among different fertilization regimes.
Tillage mode | Treatment | pH | SOC (g·kg–1) | TN (g·kg–1) | AP (mg·kg–1) | AK (mg·kg–1) | C/N |
Rotation System (R) | NF | 6.49 ± 0.25A | 14.72 ± 0.63A | 1.42 ± 0.07A | 63.26 ± 0.52A | 59.91 ± 5.90A | 10.38 ± 0.25A |
NPK | 6.42 ± 0.09A | 15.69 ± 0.70A* | 1.53 ± 0.08A* | 64.32 ± 1.96A | 51.73 ± 2.65A | 10.26 ± 0.07A |
OM | 6.40 ± 0.11A | 15.40 ± 0.36A | 1.53 ± 0.06A | 67.48 ± 2.39A | 55.18 ± 9.72A | 10.09 ± 0.26A |
Continuous System (C) | NF | 7.42 ± 0.05a* | 13.14 ± 0.84b | 1.25 ± 0.06b | 57.36 ± 7.03b | 72.53 ± 5.90c** | 10.53 ± 0.17a |
NPK | 7.08 ± 0.18a* | 12.51 ± 0.21b | 1.20 ± 0.06b | 75.29 ± 8.81b | 124.56 ± 11.19b** | 10.46 ± 0.42a |
OM | 6.68 ± 0.07b* | 16.90 ± 0.81a | 1.76 ± 0.10a | 108.82 ± 1.79a** | 162.40 ± 8.92a** | 9.58 ± 0.14b |
Note: Mean ± standard deviation (n = 3). Uppercase letters represent significant (P < 0.05) differences among rotation systems; lowercase letters represent significant (P < 0.05) differences among continuous systems. Asterisks indicate significant differences between fertilizer treatments. *P < 0.05, **P < 0.01. SOC: soil organic matter; TN: total N; AP: available P; AK: available K; C: N: SOC/TN; R: rotation system; C: continuous system; NF: no fertilizer; NPK: chemical NPK fertilizer; and OM: compost pig manure only. |
Bacterial community diversity
A total of 1 443 469 high-quality 16S rRNA sequences were generated in this study (sequences ranged from 80 051–80 381) (Table S4). Supplementary Fig. S1 shows the 16S rRNA rarefaction curves. This indicated that the sequencing depth was sufficiently deep to cover most of the bacterial groups, including some rare species. Five indices (observed species, Shannon, Simpson, good coverage, and Pielou evenness indices) of the α-diversity showed no statistical difference among the three treatments in the two cropping tillage systems, whereas the observed species and Shannon index decreased under continuous cropping (Table S5). The nested PERMANOVA analysis suggested that the total variation in the bacterial community was mainly explained by the tillage system (R2 = 0.439, P = 0.002), whereas the fertilization practice did not show a significant influence on the bacterial community (R2 = 0.173, P = 0.089). Significantly higher bacterial diversity was observed in the soil of rotation system than in that of continuous tobacco system (P < 0.05, Wilcoxon rank-sum test) (Fig. 1a, b). Measures of the richness and diversity for the ASV definition of 99% similarity for 18 soil samples revealed that the tobacco-rice rotation system exhibited stronger responses than the continuous tobacco system in bacterial community diversity.
The β-diversity of the bacterial communities was also analyzed using principal coordinate analysis (PCoA) with unweighted UniFrac distance matrices (Fig. 1c, PERMANOVA: F = 3.72, P = 0.001). The results demonstrated distinct community separation of the two tillage patterns, suggesting that the tobacco soil bacterial community composition under long-term fertilization was mainly driven by the cropping system. Soil bacterial communities from the three treatments under the continuous system were separated along the Axis.2 principal coordinate, which indicated that fertilizer application influenced the bacterial community, especially with OM, compared to no-fertilizer treatment (NF) (Fig. 1c).
Microbial community composition and sensitive bacterial taxa
To compare the similarity and dissimilarity in the bacterial community compositions among all the samples in the two tillage patterns, Venn diagrams with unique and shared ASVs were generated. Here, 1,498 and 1,019 ASVs were shared among the three treatments in the rotation and continuous systems, accounting for 63.15 and 49.25% of the total observed ASVs, respectively (Fig. 2). The 16S rRNA gene sequences were affiliated with 35 phyla and 380 genera in the rotation system and 25 phyla and 316 genera in the continuous system. At the phylum level, the dominant phylum in all soil samples was Actinobacteriota, representing an average of 35.8% of all sequences, followed by Proteobacteria, Firmicutes, Chloroflexi, and Acidobacteriota. The predominant 10 phyla accounted for approximately 85.75-93.00% of the shared ASVs in all samples (Table S6). The relative abundance of the predominant phylum was similar among the RNF, RNPK, and ROM. However, phyla Firmicutes, Bacteroidetes, Actinobacteria, and Verrucomicrobia were significantly (P < 0.05) different among CNF, CNPK, and COM. The LEfSe (LDA > 3 or LDA < − 3, P < 0.05) was used to identify Micromonosporaceae and actinobacterium_WWH12 as different cultivar level taxa in RNPK and COM, respectively (Fig. 2). Additionally, 894, 941, and 882 ASVs were shared among the NF, NPK, and OM treatments, respectively (Fig. S2). Notably, compared with those under the continuous system, the relative abundances of Chloroflexi, Desulfobacterota, MBNT15, and Nitrospinota were significantly enriched under the rotation system, except for Methylomirabilot (Fig. 3).
At the genus level, the top five genera in all soil samples were Bacillus, Gaiella, Sphingomonas, Nocardioides, and MB-A2-108, representing an average of 12.0% of all the sequences (Table S6). The top 35 genera (average relative abundance > 1.0%) were analyzed. Compared to that in RNF, KD4-96 was enriched in RNPK and ROM. Microvirga, Rubrobacter, and Bryobacter were significantly different among the three treatments in the continuous system (P < 0.05). In-depth analyses were then performed at the genus level; the dominant (> 1%), common (0.1–0.1%), and rare (< 0.1%) genera were classified based on the relative abundance of the respective sequences (Lyons et al., 2015; Li et al., 2018). The total detectable microbial genera were not significantly different between fertilization treatments, with those under continuous system being approximately 20% lesser than those under rotation system, especially the common and rare species (Table S7). Moreover, a significant number of bacterial genera showed clear differences across the tillage systems (Fig. S4). Eighty-four bacterial genera showed significant differences between RNF and CNF, including 3 dominant, 36 common, and 45 rare bacterial genera. Similarly, 119 genera showed significant differences between RNPK and CNPK, including 8 dominant, 62 common, and 89 rare bacterial genera. A total of 131 genera showed significant differences between ROM and COM, including 8 dominant, 60 common, and 105 rare bacterial genera. Interestingly, 31 genera, including 1 dominant, 18 common, and 29 rare bacterial genera, belonging to 24 classes and 13 phyla, were significantly enriched under both rotation and continuous systems (Fig. 4). Twelve genera identified under rotation system were not detected in continuous system. Of the 19 genera with relative abundances greater than 0.1%, BIrii41 and Pirellula were significantly enriched under continuous system, whereas the other 17 genera were significantly enriched under rotation system (P < 0.05). In conclusion, compared with fertilization, the tillage system had a greater impact on soil bacterial abundances, especially for common and rare genera.
Functional gene prediction
Each ASV in the dataset was analyzed through the PICRUSt2 to predict the functional map of the soil bacterial community among different treatments. Heatmap calculating expression abundance multiples showed that more complex functions changed with the fertilizers and tillage treatment. The results showed that rotation and continuous systems presented distinct diversity patterns (Fig. 5a). The relative abundance of the top 20 metabolic pathways was not significantly different among the three treatments under the rotation system; however, two metabolic pathways (L-serine and glycine biosynthesis I and cis-vaccenate biosynthesis) were significantly different among CNF, CNPK, and COM (Table S8). These 20 metabolic pathways could be roughly divided into four categories, including the amino acid biosynthesis pathway, tricarboxylic acid (TCA) cycle, nucleic acid metabolism, and diacylglycerol biosynthesis, which were significantly higher under the rotation system than under the continuous system (Fig. 5). Principal component analysis (PCA) showed distinct metabolic pathways among all treatments, distinct community separation of the two tillage patterns, and separation of the three treatments along the PC2 principal coordinate under a continuous system. This suggests that compared with fertilization, the tillage system had a greater impact on soil bacterial function; rotation could improve the metabolic function of bacteria, which may contribute to the formation of a more stable and healthy soil ecosystem.
Abundance of N-cycling genes
The bacterial function prediction showed that the nitrogen cycle was the main function. To further investigate the genetic differences in the bacteria involved in the nitrogen cycling pathway, we determined the effect of the fertilizer treatments on the abundance of soil bacterial nifH, archaeal amoA, bacterial amoA, nirS, nirK, and nosZ genes via real-time PCR (Fig. 6). The results revealed that all treatments could be arranged into two groups. The abundance of the functional genes under the continuous system was generally lower than that under the rotation system. Compared with those under the continuous system, the nifH, nirS, nirK, and nosZ gene copy numbers were higher under the rotation system (P < 0.05) and showed similar patterns. There were no significant differences among the six nitrogen cycling genes in the RNF, RNPK, and ROM treatments. However, notable changes were observed among the different fertilization treatments in the C system; COM had a higher abundance of the six N-cycling genes than CNPK and CNF. The PERMANOVA analysis showed that variations in the six N-cycling genes were mainly explained by the tillage system (R2 = 0.48, P < 0.001) rather than the fertilization practice (R2 = 0.16, P = 0.02). Therefore, compared to fertilization, tillage systems had a greater impact on soil nitrogen cycling gene activities. The application of compost pig manure could significantly increase the absolute abundance of N-cycling genes under continuous cropping systems.
Correlations among environmental variables, bacterial community, and nitrogen functional genes
Redundancy analysis was conducted to assess the relationships between the environmental variables and bacterial communities at the phylum level. The results demonstrated that, in different tillage systems, the soil chemical properties (i.e., the pH, SOC, TN, AP, AK, and C/N) had different effects on the bacterial communities (Fig. 7). The contribution of the six indicators to rotation system was not significant (P = 0.525). However, critical environmental factors, i.e., AP and AK, were the vital factors that significantly (P = 0.006) explained 76.2 and 72.7% of the variance (adjusted R2) under the continuous system, respectively. The influence of each soil variable on the community structure between rotation and continuous systems under the same fertilization was further analyzed via RDA (Fig. S5). The entire community structure variance in the NF treatment was significantly explained by the pH (R2 = 0.922, P = 0.039), whereas AK (R2 = 0.967, P = 0.025) and AP (R2 = 0.876, P = 0.006) were the main factors in soils fertilized with NPK and OM, respectively (Table S9).
Based on the nonrandom aggregation patterns of bacterial communities influenced by tillage systems, network models were constructed to further explore the topological and taxonomic characteristics of the bacterial co-occurrence patterns. The network of rotation system (consisting of 372 edges with an average degree of 7.440) was smaller than that of continuous system (662 edges with an average degree of 11.414) (Table S10), indicating that the network of continuous system (Fig. 8b) displayed greater complexity and connectivity than that of rotation system (Fig. 8a). The nodes in the network were classified into 19 bacterial phyla; most of the main bacteria belonged to five phyla: Proteobacteria, Actinobacteria, Acidobacteriota, Chloroflexi, and Firmicutes. Keystone species play an important role in maintaining the structure and function of microbial communities. These tags were located in the center of the network and in a highly connected location. Therefore, the Proteobacteria and Firmicutes phyla were identified as the keystone taxa in rotation and continuous systems, respectively.
To explore possible interactions between the bacterial communities and functional genes involved in bacterial nitrogen metabolism, co-occurrence networks were established to identify the different associations among the microbial taxa (relative abundance > 0.1%) and functional genes in the tobacco soil. Genera alphaI_cluster and Dongia (phylum Proteobacteria) had strong positive correlations with nirK and nosZ (R > 0.8 and P < 0.01), which indicated that these two genera played important roles in the processes of denitrification in the rotation system. Genera Clostridium_sensu_stricto_1, Pedomicrobium, SC-I-84 (phylum Proteobacteria), SWB02, and Turicibacter (phylum Firmicutes) had strong positive correlations with nifH, amoB, nirK, and nosZ (R > 0.8, P < 0.01), which indicated that these genera played important roles in the processes of nitrogen fixation and denitrification in the continuous system. Additionally, nifH, amoB, and nirK were the hubs of the network, indicating that nitrogen fixation and nitrite reduction were likely the most important processes in the continuous system.
In summary, the RDA showed that fertilization under different tillage modes had different effects on the soil microorganisms. AP and AK were critical environmental factors in the continuous tobacco cropping systems. Otherwise, we found that different cropping system has the potential to fundamentally alter the co-occurrence network’s properties and structure. Which could alter the key taxa of soil microbiomes in tobacco agroecosystems.