Vertical and temporal variations of soil bacterial and archaeal communities in wheat-soybean rotation agroecosystem

Field description

Our field experiment was carried out at the Chinese Academy of Agricultural Sciences, located in the North China Plain (39.97N, 116.33E), an important food-producing region in China. This region has a warm temperate monsoon climate, with an average annual temperature of 11–13 °C and average annual precipitation of 500–700 mm (80% rainfall occurs from July to August). The soil has a loamy texture and belongs to the fluvo-aquic soil class. The field was fertilized before wheat planting (September 10) under the wheat-soybean rotation system; other agricultural practices, including weeding and tillage, were identical to those done by the local agricultural management.

Soil sampling and physicochemical analysis

Soil sampling was conducted on October 18 (wheat season) and December 18 (wheat season) of 2019, and on April 17 (wheat season) and August 17 of 2020 (soybean season). In brief, six plots (1 × 1 m) were randomly selected within an area of 20 × 20 m in our study site. For each plot, five soil cores were collected between two plant individuals from a 100 cm long vertical profile that corresponds to depths of 0–20 cm (surface soil layer), 20–50 cm (middle soil layer), and 50–100 cm (deep soil layer). Therefore, each soil sample was a mixture of five soil cores (3.8 cm diameter) for a given soil layer. 72 soil samples (6 plots × 4 sampling times × 3 soil depth) were collected from the six sampling sites. One part of the as-collected soil samples was filtered through a sieve of less than two mm and kept at −80 °C for subsequent extraction of DNA. The other part was sieved with less than one mm and less than 0.25 mm then stored at 4 °C or room temperature to determine the soil physicochemical properties. Soil physicochemical characteristics including soil total carbon (TC), total nitrogen (TN), dissolved organic carbon (DOC), dissolved organic nitrogen (DOC), nitrate (NO${}_3^{-}$-N) and ammonium (NH${}_4^{+}$-N), moisture, pH, clay, silt, and sand content were determined as described by Bao (2000) and Yang et al. (2017).

Miseq sequencing and bioinformatics

The soil DNA was extracted using a PowerSoil DNA Isolation Kit (MO BIO Laboratories, USA). The V4 region of the 16s rDNA was amplified using primer pairs 515F/806R (Caporaso et al., 2012). Primer 515F contained a 12 bp barcode unique to each sample for Miseq sequencing detection. All PCR reactions followed Caporaso et al. (2011) within a 25 µL reaction system and were amplified in triplicate. The PCR products were then pooled, purified, and sequenced on an Illumina MiSeq platform at Majorbio Biotech Co., Ltd. (Shanghai, China). The raw sequence data have been deposited on the NCBI SRA (accession No. PRJNA766099).

The raw sequences were trimmed to the shortest sequence length using QIIME 2 Pipeline Version 1.8.0 (Bolyen et al., 2018). Subsequently, the sequences were dereplicated, with all singletons discarded. Sequences were then error-filtered and grouped into amplicon sequence variants (ASVs) using the Deblur software (Amir et al., 2017), with ASVs containing less than two reads removed. Further, potential chimeras were discarded using the Vsearch software. The number of sequences per sample was rarefied to 1,000,000 using the “vegan” package in R. To annotate their taxonomy, the ASVs were blasted against the silva 16s database. Then, the ASV abundance tables were rarefied at 845,591 for bacteria and 10,413 for archaea to ensure even sampling depth within each prokaryote group.

Data analysis

ASVs with relative abundances below 0.01% were defined as “rare” taxa, while those with relative abundances above 0.1% were defined as “abundant” taxa (Jiao & Lu, 2020). The prokaryotic functional profiles were predicted using the FAPROTAX (Louca, Parfrey & Doebeli, 2016). Soil fertility index is a synthetic variable calculated from the sum of z-score transformation of TN, TC, DON, DOC, NH${}_4^{+}$-N, and NO${}_3^{-}$-N. Bacterial and archaeal alpha-diversity indices, such as richness and Pielou evenness index, were calculated using the “vegan” package (Oksanen et al., 2013). Two-way analyses of variance (ANOVAs) were then performed to examine the effects of soil depth, sampling time, and their interaction on alpha-diversity indices and the proportion of rare ASVs. All data were tested for normality and homogeneity of variance before two-way ANOVAs. A post-hoc test was used to determine paired comparisons among the treatments at a 5% significance. Random forest analysis (Breiman, 2001) was used to explore the soil physicochemical drivers of bacterial and archaeal richness and evenness using the “randomForest” package (Liaw & Wiener, 2002). The “rfPermute” package (Archer, 2016) was then utilized to estimate the significance of important metrics for a random forest model by permuting the response variable.

The effects of soil depth, sampling time, and their interaction on soil bacterial and archaeal community compositions were evaluated using a permutational multivariate analysis of variance (PERMANOVA) with 999 permutations in the “vegan” package (Oksanen et al., 2013). The bacterial and archaeal community compositions were subsequently ordinated using principal coordinates analysis (PCoA) based on the Bray–Curtis dissimilarity matrices in the “vegan” package. A ternary plot was used to demonstrate the distribution of bacterial and archaeal ASVs along the vertical soil profile.

Soil prokaryotic co-occurrence networks for the surface, middle, and deep soil layers were built. Each network was based upon 24 prokaryotic soil communities. ASVs with relative abundances > 0.5% that also occurred in > 50% of the communities were included in the networks to focus solely on abundant ASVs. Spearman’s correlation coefficients were calculated between ASVs using the “Psych” package (Revelle, 2015). P-values for multiple tests were calculated using the false discovery rate (FDR), as described by Benjamini & Hochberg (1995). The correlations with a Spearman’s coefficient < 0.8 and a P-value > 0.001 were eliminated (Widder et al., 2014). Subsequently, the number of positive and negative correlations, average degree, connectedness, and modularity was calculated in each network using the “igraph” packages. Then, threshold values of Zi and Pi were computed following the reference (Guimera & Amaral, 2005) for each node to determine their topological role. All analyses above were carried out in R (v.3.6.2).

Sequencing data analysis and prokaryotic diversity

A total of 185,093,268 reads were obtained after quality control and chimera checks, assigned as 225,608 ASVs. The read numbers were normalized to 1,000,000 for each sample, resulting in a normalized dataset containing 163,862 prokaryotic ASVs. Then, the read numbers of bacteria and archaea were normalized to 845,591 and 10,413 after species annotation, resulting in two normalized datasets containing 162,802 bacterial ASVs and 213 archaeal ASVs. Among these ASVs, 7,382 ASVs (4.5%) occurred in at least half of all samples.

The bacterial and archaeal diversities were assessed using the richness and Pielou evenness indices. One of the most interesting findings was that the vertical pattern of soil bacterial and soil archaeal diversity differed. Two-way ANOVA analysis indicated that both bacterial richness and evenness were significantly affected by soil depth (Table 1). Generally, bacterial richness and evenness were lowest in the deep soil layer and peaked in the surface or middle soil layer (Figs. 1A and 1C). In contrast, archaeal richness was unaffected by soil depth, while its evenness showed a decreasing trend along with the soil profile (Figs. 1B and 1D). We did not find an apparent temporal variation of archaeal richness or evenness (Table 1). Although bacterial richness did not show temporal variation (Table 1), its evenness index was significantly higher in August than in December (Fig. 1C). Random forest analysis indicated that TN, TC, and sand content were the key determinants for bacterial and archaeal diversities (except for archaeal richness, Figs. 1E–1H). The soil physicochemical variables are presented in Table S1.

Soil prokaryotic community composition

PERMANOVA analysis indicated that both bacterial and archaeal community compositions were significantly impacted by soil depth (Bacteria: r² = 0.38, P < 0.001; Archaea: r² = 0.70, P < 0.001), but unaffected by sampling time and their interaction. This pattern was further evidenced by PCoA ordination based on Bray–Curtis dissimilarity, which indicated that the bacterial and archaeal community compositions were separated by soil depth (Figs. 2A and 2B). We then tested for the temporal variations of the bacterial and archaeal community compositions at each soil depth. We observed that the archaeal soil community did not show any temporal variations at any soil depth (Fig. S1). In contrast, the bacterial soil community composition exhibited strong temporal variation in the middle soil layer (PERMANOVA: r² = 0.13, P = 0.006, Fig. S1).

At the phylum level, the bacterial community was mainly dominated by Actinobacteria, Proteobacteria, Firmicutes, Chloroflexi, and Acidobacteria, while the archaeal community was mainly dominated by Crenarchaeota in our study site, accounting for 69.19% of the total prokaryotic sequences (Fig. 2C). Next, we investigated the changes in relative abundance patterns of prokaryotic phyla associated with the soil profile. The relative abundance of nearly all phyla was significantly different among soil layers for all sampling times, except for Gemmatimonadetes and Chloroflexi. We observed that Proteobacteria, Firmicutes, Acidobacteria, Planctomycete, and Myxococcota, especially, decreased along with the soil profile, while Actinobacteria and Crenarchaeota exhibited the opposite trend (Fig. 2C, Table S2). Most strikingly, the proportion of unclassified ASVs showed an increasing trend along with the soil profile (Fig. 2C).

The shifts in bacterial and archaeal communities were also reflected in the abundant ASVs, with 50%–72.49% of total abundant ASVs unevenly distributed along the soil profile among sampling times (Figs. 3A–3D). ASVs enriched in the surface soil layer, mainly Solirubrobacterales and Propionibacterales, accounted for a more significant proportion than those in the middle and deep soil layers (Fig. 3E). ASVs enriched in the middle and deep soil layer, primarily classified as Rokubacterales and Gaiellales, differed from those in the surface soil layer (Fig. 3E). Meanwhile, abundant ASVs also exhibited temporal variations: there were less enriched ASVs (50%) in October than in other sampling times (71.5%–72.5%, Figs. 3A–3D). Interestingly, the proportion of unclassified ASVs was less abundant in the surface soil than in the middle and deep soil layers (Fig. 3E).

Predicted prokaryotic function

FAPROTAX was adopted to annotate and then screen the key ecological functions of the bacterial communities that contributed to soil biogeochemical cycling. A total of 48 functional groups were identified. Although a minor proportion of all ASVs (5.04%) was assigned to at least one predicted functional group, the predicted functional composition was significantly correlated with the soil prokaryotic community composition as confirmed by the Mantel test (r = 0.64, P < 0.001). Among these functional groups, we observed that 41 predicted functional groups exhibited significant differences among depths, while 21 functional groups displayed apparent temporal variations (Table S3).

Rare bacterial and archaeal ASVs

The proportion of rare bacterial richness and abundance was highest in the surface soil layer and showed a decreasing trend along with the soil profile (Figs. 4A and 4C). On the other hand, the proportion of rare archaeal richness and abundance were not impacted by soil depth (Figs. 4B and 4D). Random forest analysis further indicated that the rare bacterial community was the key determinant of the soil fertility index and functional diversity (Fig. S2). Moreover, the potential role of 43 archaeal and bacterial phyla on soil fertility was assessed. Among these phyla, Latescibacterota, Zixibacteria, Dadabacteria, Dependentiae, Firmicutes, Micrarchaeota, and Calditrichota were the main predictors of the soil fertility index (Fig. S3); these phyla, except Firmicutes, were relatively rare taxa, accounting for 0.00034%–0.31% of the total reads. We further explored the potential functions of rare phyla through FAPROTAX: Desulfobacterota accounted for 90.65% of sulfate respiration and 46.26% of iron respiration; Cyanobacteria accounted for 6.10% of nitrogen fixation and 97.50% of phototrophy; Halanaerobiaeota accounted for 62.90% of cellulolysis; and Euryarchaeota and Halobacterota accounted for 14.33% and 0.12% of dark hydrogen oxidation, respectively.

Prokaryotic co-occurrence network

We constructed a soil prokaryotic co-occurrence network for each soil layer. The network size generally became larger along with increasing soil depth, with 214, 247, and 274 nodes in the surface, middle, and deep soil layers, respectively (Figs. 5A–5C). Additionally, we observed that only a small proportion of nodes (54) were shared among these three networks. The network of the surface soil layer was analogous to the network of the middle soil layer (153 nodes shared) but was quite dissimilar to that of the deep soil layer (65 nodes shared). The surface and middle soil layers’ networks were more connected than that of the deep soil layer (Figs. 5A–5C). This pattern was affirmed by their topological characteristics (Table S4). Furthermore, the number of positive/negative links gradually decreased along with the soil profile (Table S4).

From the plot of Zi (a value measuring within-module connectivity) against Pi (a value measuring among-module connectivity), the different roles of each node in the network were identified (Figs. 5D–5F). Notably, there were more keystone species observed in the deep soil layer (six) than in the surface (one) and middle (three) soil layers (Figs. 5D–5F). These keystone species were not shared among the three networks. The classification of these keystone species is listed in Table S5.

Different vertical patterns between bacterial and archaeal diversities

We observed that the deep soil layer harbors the lowest bacterial diversity, which is consistent with previous studies from the alpine meadow (Xu et al., 2021), soybean cropland (Hao et al., 2020), and deciduous birch forest (Mundra et al., 2021). This finding confirmed the ecological principle that more extreme environments are expected to be inhabited by a less diverse community of microbes (Gaston, 2000). In contrast to previous studies, we observed that bacterial richness and evenness in the middle soil layer were roughly equal to that in the surface. This observation may have the following two explanations. First, this is likely attributable to the soil carbon content which did not differ between the surface and middle soil layer. Soil organic carbon content is the most important factor that influences soil microorganisms in addition to pH (Fierer, 2017; Kukreti et al., 2020) and was observed to be the best predictor of bacterial diversity, as indicated by random forest analysis in this study. Second, on a depth-weighted basis, the microbial biomass in the middle soil layer may not be less than that in the surface soil layer, thus, containing an equal diversity to the surface soil layer.

For archaea, its richness was unaffected by soil depth and cannot be explained by any measured soil variable. Archaea are reported to be highly adaptable to environmental stress (Cao et al., 2012), so may survive in the nutrient-poor and oxygen-deficient deep soil layer. In addition, we detected that archaeal communities are highly uneven in deep soil, meaning that there is extreme dominance by one or a few archaeal taxa. Therefore, the archaeal communities in deep soil would be less resistant to environmental stress (Wittebolle et al., 2009).

Vertical and temporal patterns of bacterial and archaeal communities

Both bacterial and archaeal community compositions exhibited distinct distribution patterns with the vertical soil profile, mainly due to the drastic differences in soil carbon and nitrogen content and soil texture (Jobbágy & Jackson, 2000; Eilers et al., 2012; Gu et al., 2017). We also detected different ecological preferences for certain phyla. For instance, Actinobacteria and Crenarchaeota showed increasing abundance along with the vertical soil profile. This result corresponds with the conclusions of previous studies that archaea occupy a larger proportion in the deep soil layer due to their anaerobic characteristic. However, some bacterial phyla, including Proteobacteria, Firmicutes, and Acidobacteria, displayed a sharply decreasing trend in the deep soil layer. This trend may be owing to the copiotrophic-oligotrophic trade-off theory (Ramirez, Craine & Fierer, 2012) and/or the deficiency in oxygen availability. For instance, Proteobacteria are typically copiotrophic and capable of degrading various organic materials (Eilers et al., 2012). Thus, the higher abundances of Proteobacteria in the surface layer may, in part, reflect increased organic carbon availability.

We only detected a slight temporal variation in bacterial and archaeal community compositions. In a meta-analysis, Shade et al. (2013) found that soil microbial communities were consistently less temporally variable than other ecosystems. For time-series field studies, it is impossible to collect samples in the exact same location due to destructive soil sampling. Therefore, soil microbial community dynamics may be masked by the high small-scale spatial heterogeneity (Docherty et al., 2015). Although the use of the soybean in a crop rotation would increase soil total nitrogen, there is no apparent difference in soil TN and TC among seasons (Table S6) in the present study, perhaps due to the removal of crop residues after harvest. As proposed by previous studies, temporal variations of soil microbial communities are driven by soil physicochemical characteristics (Yang et al., 2017). Consequently, the minimal temporal variations might be explained, in part, by the unapparent difference in soil TC and TN. Alternatively, soil contains many dormant microbes (Lennon & Jones, 2011), which were indiscriminately captured using the 16s rRNA sequencing approach (Shade et al., 2013). For the reasons above, the soil prokaryotic communities only fluctuated within a narrow range over the given time period.

Vertical pattern of rare ASVs and their potential function

Deep sequencing enabled us to explore reliable rare taxa in soil. For the first time, the vertical patterns of rare bacterial and archaeal taxa were evaluated. We observed that the proportion of rare bacterial ASVs decreased along with the soil profile. The root system in the surface soil layer was more abundant than that in the middle or deep soil layers (Chaparro, Badri & Vivanco, 2014) and may provide more niches for rare species with specialized functions. Alternatively, it was proposed that rare species may occur from stochastic processes (Jousset et al., 2017), whereas deterministic processes would reduce rare species. Since soil nutrient availability generally decreased along with the vertical soil profile, stochastic processes would play a less predominant role in the deep soil layer.

Microbial rare taxa constitute an important “genomic reservoir” and could be used for a variety of functions (Jousset et al., 2017). It was evidenced that rare taxa contributed to sulfate reduction (Pester et al., 2012) and have the potential for bioremediation (Pascoal, Magalhães & Costa, 2020). Recently, the ecological importance of rare taxa in terrestrial and aquatic ecosystems has been increasingly recognized (Jousset et al., 2017; Chen et al., 2020). In the present study, the community of rare taxa was observed to be a key determinant of soil fertility and ecosystem multifunctionality, indicating that rare taxa contribute greatly to agroecosystem functioning. Likewise, Chen et al. (2020) also observed that rare bacterial and fungal taxa (less than 3% of total reads) were the major drivers of ecosystem multifunctionality in long-term fertilized soils. As revealed by Random Forest analysis, rare prokaryotic phyla (e.g., Latescibacterota, Zixibacteria, and Dadabacteria) had an over-proportional role in the soil biogeochemical cycling process. In contrast, the dominant bacterial and archaeal phyla (e.g., Proteobacteria, Actinobacteria, and Acidobacteria) showed little control over soil fertility. Our analysis further revealed that rare phyla influence narrow-range functions, such as dark hydrogen oxidation, sulfate respiration, and cellulolysis. Other studies highlighted that rare taxa exert great influence over sulfate reduction (Pester et al., 2010; Dell’Anno et al., 2012) and organic pollutant degradations, where their removal greatly impacted ecosystem function (Jousset et al., 2017). Taken together, our knowledge of the ecological role of the rare microbial taxa is still obscure with more attention needed.

Vertical patterns of prokaryotic co-occurrence networks

Our results emphasized that the prokaryotic co-occurrence network was more complex in the surface and middle soil layers than in the deep soil layer, reflected by the higher connectedness and average degree. Likewise, Xu et al. (2021) observed a more complex and more extensive microbial network in the surface than in the deep soil layer of meadows and shrubland. Higher network complexity in the surface and middle soil layer might be explained, in part, by greater amounts of substrate and nutrient in the surface soil layer (Creamer et al., 2016), as microbial interactions feeding on these substrates would be largely strengthened. In fact, a large number of studies have found that soil microbial networks would be more connected after nutrient amendments (Yang et al., 2019; Yang et al., 2020). Moreover, the microbial richness and biomass were extremely low in the deep soil layer (Eilers et al., 2012; Xu et al., 2021), reducing the opportunity for different species to interact. It has been proposed that complex networks contribute to better environmental adaption and higher efficiency of resources transfer than simple networks (Morriën et al., 2017). In this sense, the complex prokaryotic networks in the surface and middle soil layers would lead to better biogeochemical functions and provide a stable soil environment for plant growth.

We observed a gradual increase in proportion of negative links along the soil profile by assessing the correlations among microbial species, which indicates that microbial competition for the same resources increases with soil depth. The niche differentiation would be a possible explanation for this result. The surface and middle soil layers are easily impacted by external disturbance; the roots within these layers also create heterogeneous environments (Morriën et al., 2017), leading to significant niche differentiation and reduction of competition. In contrast, the deep soil layer is deficient in roots and is less impacted by external disturbance, possessing a more homogeneous environment and weak niche differentiation (Faust & Raes, 2012). The nutrient-limited environment in the deep soil layer would also induce intense competition among species.

In conclusion, we fully evaluated the vertical and temporal patterns of soil bacterial and archaeal communities and co-occurrence networks in a wheat-soybean rotation agroecosystem. Our results showed that soil depth as a more predominant factor than sampling time in shaping soil prokaryotic communities and co-occurrence networks and revealed different vertical distribution patterns of bacteria and archaea. With deep sequencing at a depth of millions of sequences per sample, we can investigate the rare microbial taxa and reveal their potential in soil biogeochemical function. Overall, our better understanding of vertical and temporal patterns of soil prokaryotic communities should enable more strategic agricultural management to mitigate adverse environmental impacts and improve crop productivity.