Bacterial community composition is an important predictor of surface soil fertility across different land use types: a case study in the Three Gorges Reservoir area

Experimental design and soil sampling

The study site was located at Zhangjiachong watershed (110°57′E, 30°46′N;133 to 632 m above sea level, 1.62 km² area) which is part of the Three Gorges Reservoir area and about 12 km southwest of Three Gorges dam. According to the WorldClim database (Fick & Hijmans, 2017), annual mean temperature is 10 °C and annual precipitation is 1,200 mm in this area. The area is characterized by yellow-brown soils that originate from the weathering of granite and/or sandstone. According to China’s soil classification system, these soils are known as Alfisols, while in the United States, they are categorized as both Alfisols and Ultisols (Wang, Xu & Liu, 1988). On sloping terrains, typical landscape in this area, the granite-derived soils often have a sandy structure. The predominant land uses in this region include agricultural land (26.7%), forests (60.6%), orchards for economic crops (4.6%), grasslands (2.0%), and barren hillsides (4.9%). Soil erosion primarily happens on sloping croplands and in areas with sparse vegetation, especially hills with slopes over 25° (Meng et al., 2024).

Eleven sampling plots were selected according to land use types. The land use types were characterized by a gradient intensity of human activity (Table 1). We scored the human activity according to a previous study (Cheng et al., 2022), and whether one type of human activity existed or not in a given type of land use. We considered sloping lands to be the land use type influenced mostly by human activity (scored 6, Table 1), as the crop fields with Brassica campestris and Arachis hypogaea might be disturbed by agricultural activities, such as plowing, fertilizing, crop rotation, trampling, weeding, pest and disease controlling and harvesting. In addition, trampling, weeding, pests and disease controlling and harvesting were considered to exist in flat lands as well, because orange and tea pickings would occur there. We treated forests as the land use type with the lowest intensity of human activity, such as occasional logging.

To avoid sampling bias of plant attributes (as those plants are assumed to affect soil properties and hence soil microorganisms), a series of plant species were included within each type of land use. The three types of land use were: sloping croplands (sloping lands), flat croplands (flat lands), and secondary forests (forests). The predominate plants in sloping lands were orange (Citrus reticulata) and tea trees (Camellia sinensis), and crop species, e.g., B. campestris (in spring) and A. hypogaea (in autumn). In total, six plots of sloping land were sampled, among which three were planted with hedgerows to prevent soil erosion (plot No. 1-6 in Table 2). The major plant species in flat lands included orange and tea trees, and there were no other trees nor crops. Two sampling plots were set up from flat land (plot No. 7, 8 in Table 2). We selected three kinds of forests, in which the constructive species were camphor tree (Cinnamomum camphora), camphor and tea trees, and Fagus longipetiolata. Three sampling plots were selected, one in each of the forest types (plot No. 9-11 in Table 2).

Soil samples were collected in spring (April) and autumn (September) 2016. The five-point sampling method was applied during soil collection (Han et al., 2024). Five soil cores were collected from each quadrat: one from each of the four corners and one from the center. These cores were then pooled together to form a composite sample. Due to the plot size and condition variation, the soil was sampled in three 2 m × 2 m quadrats for crop lands, and three 5 m × 5 m quadrats for forests. In the meantime, plant richness was recorded in each quadrat. To avoid high levels of plant residue in the shallow soil layer and to minimize the effect of soil depth on the microbial community, soil was sampled at a depth of 5–10 cm instead of the conventional 0–10 cm. In total, 66 soil samples were collected, with three replicates for each plot at each season. An overview of the plot information is available in Table 2. Photographs demonstrating the natural condition of each land use type sampled are in Fig. S1.

Soil properties

Soil samples were air dried, and soil physicochemical parameters were measured as described previously (Xu et al., 2020; Xu et al., 2021c). Briefly, soil pH was determined by a pH meter (UB-7, Denver, USA) after 1 min of sonication and shaking for 30 min of soil solution (1:5 weight (soil)/volume (distilled water)). The air-dried moisture content (MC) was measured gravimetrically. Total N (TN), Total P and total potassium (TK) was determined by the semimicro-Kjeldahl methods (Bremner, 1996), Mo-Sb colorimetric method (Murphy & Riley, 1962) and flame atomic emission spectrophotometry, respectively. The available N in the soil was measured by the alkaline hydrolysis diffusion method (Mulvaney & Khan, 2001); available P by Olsen’s method (Olsen, Cole & Watanabe, 1954); available potassium (AK) was extracted with a NH₄OAc leaching-flaming luminosity (Zhou, Zhang & Niklas, 2014), and soil organic carbon (SOC) was measured by the soil organic matter by K₂Cr₂O₇- H₂SO₄ oxidation method (Nelson & Sommers, 1996).

Assessment of soil fertility

To construct a soil fertility index, we selected seven key items: the SOC, total and available N, P, K contents. These screened items consisted of the key nutrient parameters which indicate the provision of nutrients for plant growth and microbial metabolism. The z-score approach (averaging the z-score values of the above-mentioned soil fertility parameters) was used to obtain a quantitative soil fertility index (Maestre et al., 2012). The data related to these soil environmental traits used in this study are offered in Data S1.

Soil microbial Biolog incubation

The metabolic activity of the microbial community was measured using the Biolog ECO-plates (Biolog Inc., Hayward, CA, USA) (Garland & Mills, 1991). Soil samples were stored at 4 °C for 24 h prior to measurement. Before analysis, the samples were incubated at 25 °C for 12 h to stimulate microbial metabolic activity. Soil microorganisms were extracted with an aliquot of 0.5 g of soil sample mixed with 24.5 mL of saline solution (0.85% w/v). Then the solution was diluted 1000 times and moved into Biolog ECO-plates. The incubation was conducted in a light-free box (to stimulate the dark environment in soil) at 25 °C for 240 h (Manjunath et al., 2018). Average well color development (AWCD) was calculated to describe metabolic activity of the microbial community (Equation (1)). (1)${\displaystyle AWCD=\sum \left(Ci-Ri\right)/31}$

where Ci is the optical density (at 590 nm) of each well with sample, while Ri is the average optical density of the control wells. Plates with an incubation time of 96 h were chosen for further analysis, because this time point represented the logarithmic phase of incubation, which could be a good proxy for microbial metabolic activity (Fig. S2).

Bacterial 16S rRNA gene sequencing

Soil samples for sequencing analyses were stored at −20 °C in a refrigerator until measurement. Soil metagenomic DNA was extracted using the D5625-01 (Omega, Norcross, GA, USA) DNA extraction kit according to the manufacturer’s instruction. To ensure that contamination was avoided, we included negative H₂O, no negative PCR mix and primer controls in the PCR reactions. For each sample, triplicate PCR reactions were performed. The V3-V4 hypervariable regions of the bacterial 16S rRNA gene were amplified using a PCR Amplifier (2720, ABI, USA), with the primer set 515F (5-GTGCCAGCMGCCGCGGTAA-3) and 806R (5-GGACTACVSGGGTATCTAAT-3) (Caporaso et al., 2011). The PCR reactions were performed in a 25 µL mixture and included 0.25 µL of Q5 high-fidelity DNA polymerase (M0491L; NEB, Ipswich, MA, USA), five µL of 5 × reaction buffer, five µL of 5 × GC buffer, five µL of dNTP (10 mM), one µL of forward primer (10 µM), one µL of reverse primer (10 µM), one µL of DNA template, and 11.25 µL of ddH₂O. The thermal program for amplification was 5 min of initial denaturation at 98 °C, 10 s at 98 °C, 30 s at 50 °C, and 30 s at 72 °C for 30 cycles before a final extension of 5 min at 72 °C.

The amplicons with about 300 bp size were selected by running 2% agarose gel electrophoresis in 1.0× TAE buffer and were purified using the AxyPrep DNA Gel Extraction Kit (AP-GX-250; Axygen). After quantified on a microplate reader (FLX800; BioTek) with the Quant-iT PicoGreen dsDNA Assay Kit (P7589; Invitrogen, Waltham, MA, USA), the purified amplicons from one sample were pooled together and normalized in equimolar amounts.

DNA libraries were constructed with equimolar mixed amplicons (100 ng) according to the manufacturer’s instruction (Illumina). Finally, target DNA was sequenced on an Illumina MiSeq platform. The sequencing data is deposited in NCBI database. The accession number is PRJNA511782. The data also can be accessed via Figshare (https://doi.org/10.6084/m9.figshare.27367386).

Bioinformatics

Paired-end reads were merged using FLASH software (Magoč & Salzberg, 2011). The QIIME2 pipeline was then employed for subsequent analyses (Bolyen et al., 2019). Quality controls were processed using “quality-filter q-score” with default parameters. The amplicon sequence variants (ASVs) were obtained using the DADA2 denoise algorithm (Callahan et al., 2016). Using the “qiime feature-table filter-features” command, ASVs were further filtered with the following criteria: ASVs present less than two times, less than two soil samples and those from mitochondria and chloroplasts. Using “qiime feature-table rarefy” command, all samples were rarefied to nine, 860 sequences per sample, based on the sample with the fewest reads (see details in the Qiime2 scripts in the Code S1). For taxonomic classification of the 16S rRNA gene, the SILVA ribosomal RNA gene database (Release 138) was employed (Quast et al., 2012). The original and rarefied datasets can be found at Data S2. The scripts and these datasets are also deposited at Figshare (https://doi.org/10.6084/m9.figshare.27367386).

Statistical analysis

The normality of data distribution was assessed using the ‘shapiro.test’ function in R. To ensure comparability across different parameters, tests for differences in soil fertility, microbial metabolic activity, and bacterial diversity were conducted between seasons (Spring and Autumn) and across three land use types using consistent methods. As the Shannon diversity index did not meet the normality assumption, the non-parametric Kruskal-Wallis test was applied, utilizing the ‘kruskal’ function from the agricolae package v1.3-7 (De Mendiburu, 2023). Using the ‘vegan’ package v2.6.4 (Oksanen et al., 2011), non-metric multidimensional scaling (NMDS) and permutational multivariate analysis of variance (PERMANOVA) were used complementary to check the differences in bacterial β diversity between seasons and land use types. Random forest analysis was used to identify the biomarkers governing soil fertility, using the ‘randomForest’ package v4.7.1.1 (Liaw & Wiener, 2002). The significance of each biomarker’s importance with respect to the soil fertility index was assessed using the ‘rfPermute’ package v2.5.2 (Archer, 2013).

Evaluation of the direct and indirect effects of human activity, plant diversity, pH, and microorganisms on soil fertility was conducted by partial least squares path modeling (PLS-PM), using the ‘plspm’package v0.5.1 (Sanchez, 2013). In this model, the biomarkers (bacterial taxa) identified by random forest modeling were treated as the proxies of microbial composition, which was set as a latent variable, while the other factors were manifest variables. These latent variables included two or more manifest variables. The relative contribution of the manifest variables was estimated by their loadings in the model. The models were conducted for 999 permutations using bootstraps to obtain reliable estimates of path coefficients (representing the direction and strength of the linear relationships between variables and loadings of the indicators). All statistical analyses were conducted in R 4.3.3 (R Core Development Team, 2010). The R scripts used in this study can be accessed in the Code S2 or via Figshare (https://doi.org/10.6084/m9.figshare.27367386).

Soil properties of different land use types

All soil samples used in our study showed relatively low pH values (between five and six) and low organic C contents (from 10 to 20 mg kg⁻¹, Table S1). Soil fertility index, calculated from total and available N, P, and K, and SOC, was significantly and positively correlated with soil organic C, total N, available N and available K (Fig. S3). In addition, SOC, AN and AK showed relatively higher correlation strength (Fig. S3). The soil fertility index significantly differed among different land use types (KW test: χ²=11.135, df = 5, P = 0.049), which was lower on sloping lands compared with other land use types (Fig. 1). However, the influence of season on soil fertility was not significant regarding each type of land use (Fig. 1). In addition, when land use was controlled, soil fertility was higher in orange plots (considering sloping lands only) and Fagus plots (considering forests only) than in plots with other plant species (Fig. S4).

Variation of soil microbial communities

It is important to note that the rarefaction curve of each sample gradually reached a steady state with a sequencing depth higher than 5,000 sequences per sample (Fig. S5). Thus, resampling each sample to 9,860 sequences was enough to detect soil bacterial diversity. In total 29,294 ASVs were obtained; after filtration and resampling, the final dataset includes 12,683 ASVs (Data S2), which can be affiliated to 44 phyla 506 genera (Data S3–Data S4).

The top 10 soil bacterial phyla were Proteobacteria, Acidobacteriota, Actinobacteriota, Gemmatimonadota, Chloroflexi, Myxococcota, Verrucomicrobiota, Planctomycetota, Methylomirabilota, and Firmicutes (Fig. 2A). Proteobacteria was with higher amplicon relative abundances in samples collected in spring than those collected in autumn (paired Wilcox test: W = 256, P < 0.001); Proteobacteria (KW test: χ²=11.867, df = 2, P = 0.003) and Acidobacteriota (KW test: χ²= 45.134, df = 2, P < 0.001) had higher amplicon relative abundances in forests than croplands (sloping lands and flat lands). The highest amplicon relative abundance of Actinobacteriota was observed in sloping croplands (KW test: χ²= 39.94317, df = 2, P < 0.001, Fig. 2A).

The representative bacterial genera were Subgroup2 of Acidobacteriota, MND1, Acidibacter, Acidothermus, Haliangium, Pseudolabrys, Paraburkholderia, SC.I.84, Rokubacteriales, and Bradyrhizobium (Fig. 2B). Amplicon relative abundance of Haliangium (paired Wilcox test: W = 201.5, P < 0.001) and Rokubacteriales (paired Wilcox test: W = 705, P = 0.040) were higher in autumn than spring; while Subgroup2 (KW test: χ²= 43.854, df = 2, P < 0.001), Acidibacter (KW test: χ²= 25.655, df = 2, P < 0.001) and Acidothermus (KW test: χ²= 21.769, df = 2, P < 0.001) were with greater dominance in forests than croplands (sloping lands and flat lands, Fig. 2B).

According to the 16S rRNA gene sequencing data, bacterial diversity was significantly differed among two seasons and land use types (KW test: χ² = 51.832, df = 5, P < 0.001). The bacterial Shannon diversity was significantly higher in spring than autumn, and in cropland than forests (Fig. 3A). Regarding the same season and land use types, bacterial Shannon diversity showed minor difference among different plants, which was high in tea than orange plots in spring samples from sloping land and was lower in the Fagus forest than other plants in spring (Fig. S6).

Notably, differences in AWCD (which indicates microbial metabolic activity) among seasons were dependent on specific land use types (KW test: χ²=13.622, df = 5, P < 0.05), which was statistically similar between two seasons regarding sloping land and was higher in spring than autumn regarding flat land and forests (Fig. 3B). Additionally, in spring, sloping land demonstrated lower AWCD value than flatland and forests. A comparison of the samples from the same land use type indicated that AWCD was higher in orange plots (among all sloping land plots) and Camphor trees (among all forest plots) but was lowest in Fagus plots (among all the forest plots Fig. S7). In addition, AWCD differed in plots with or without hedgerows in the sloping lands (Fig. S7A).

According to the NMDS plot and PERMANOVA results, bacterial community composition was significantly different between autumn and spring and across different land use types (Fig. 4 and Table 3). This suggested that human-induced changes in land use significantly altered soil bacterial compositions, regardless of sampling season considered.

Environmental factors influencing soil microbial communities and soil fertility

The relative abundances of the 16S amplicon sequences related to key bacterial phyla were regressed against the soil fertility index in our random forest model. Only phyla with amplicon relative abundances higher than 0.01% were considered. The model accounted for 23.27% of the variance in soil fertility. Amplicon relative abundance of various phyla, such as Proteobacteria, Chloroflexi, Verrucomicrobiota and Gemmatinonadota and Acidobacteriaota were marginal (P < 0.1) or significant (P < 0.05) indicators of soil fertility in the study area (Fig. 5).

Path analysis demonstrated that bacterial community composition was the dominant driver of soil fertility. Moreover, microbial metabolic activity was directly and negatively affected by human activity. The negative effect of human activity on soil fertility was mitigated by the bacterial community composition (Fig. 6A). While amplicon relative abundances of Gemmatimonadota and Choroflexi were negatively correlated with soil fertility index (Fig. S8), the amplicon relative abundances of Acidobacteriota and Proteobacteria were the key drivers of soil fertility (Figs. S6B and S8).

Land use types and above—ground plant species affect soil microbial communities

We recognize potential biases introduced during the amplification process that could affect the accuracy of taxonomic predictions. Nonetheless, we highlight that despite these limitations, our study provides valuable insights into the impact of micro-habitat factors on soil prokaryotic taxa. This is supported by numerous prior studies that have successfully used amplicon sequencing techniques to gather detailed community information on soil bacteria (Zhang et al., 2016b; Xu et al., 2020; Xu et al., 2021a; Xu et al., 2021b; Xu et al., 2024). Considering the disadvantages of Biolog ECO-plates and 16S rRNA gene sequencing, we combined both techniques to explore the soil microbial community. Previous studies suggested that Biolog ECO plate-based microbial metabolic activity in soils varied with land use (Rutgers et al., 2016), dominant vegetation (C. camphor, Pinus massoniana, or Lespedeza bicolor), and root proximity (Tam et al., 2001). This study further demonstrated that both microbial metabolic activity and bacterial diversity varied among seasons and land use types. However, the patterns of microbial metabolic activity and bacterial diversity with respect to land use type differed (Fig. 2). These differences may be a consequence of the method used: Biolog ECO-plates and 16S rRNA gene sequencing. Firstly, Biolog ECO-plates are limited to measuring the microorganisms (e.g., bacteria and fungi) which could utilize the 31 types of carbon sources, while 16S rRNA gene sequencing identifies bacteria, regardless of the carbon source used. Secondly, the Biolog ECO-plate incubation conditions favor fast growing microorganisms. Consequently, ECO-plates can identify microorganisms present in soil in small quantities or even in spore form, which become active after obtaining an easy carbon source. These organisms may not be identifiable using sequencing alone, due to low DNA concentrations in the soil sample. For example, the sloping croplands harbored a more diverse bacterial community than the forests, but most of these bacteria were in low amplicon relative abundance and did not survive Biolog incubation. The higher nutrient content and amplicon relative abundances of core bacterial taxa (e.g., Acidobacteria and Proteobacteria) in forests, compared to sloping croplands, may lead to higher microbial metabolic activity and reduced bacterial diversity.

Within this watershed, human activity has led to distinct land use types. With the same soil type, different land use types showed different soil nutrient availability, bacterial diversity, and microbial metabolic activity (Figs. 1–4). Such differences were even observed among plots with different plant species. This is not surprising, since soil characteristics (Dybzinski et al., 2008; Huang et al., 2019) and microbial community are influenced by plant species (Burns et al., 2015), vegetation types, plant diversity and plant growth stages (Chen et al., 2013; Lozano et al., 2014; Phazna Devi et al., 2020; Zheng et al., 2021a; Zheng et al., 2021b; Zhu et al., 2021). In addition, we observed significant differences in soil fertility, microbial metabolic activity between soil with or without hedgerows in the sloping lands (Figs. S3–S5), highlighting the essential role of aboveground plants in controlling soil fertility and soil microbial communities. Our results suggested that activity and diversity of soil microorganisms varied among seasons and land use types. This further reflected the complex connections between land use, soil fertility, plant attributes, and soil microbial community.

The higher vegetation coverage in forests may leave more plant-derived organic matter in soils. Plant litter and root exudates may change the soil environment (Cotrufo et al., 2015). In addition, the decomposition of fallen oranges may also provide a driving force for higher soil organic matter. The accumulated organic matter may in turn retain moisture through the complex metabolic pathways of soil bacteria and therefore sustain higher soil fertility compared with other sampling plots. Moreover, the organic matter provided by root-associated processes could offer more substrates for microbial metabolism; consequently, the microorganisms selected by the Biolog ECO-plates were able to survive, because more than nine types of the carbon sources in the ECO-plates were consistent with root exudates (Choi & Dobbs, 1999). Therefore, it is reasonable that both high microbial metabolic activity and high soil fertility were observed in soils from forests and orange sloping lands.

Human disturbances decrease surface soil fertility and alter soil microbial communities

Soil microbial communities are sensitive to disturbances (Shade et al., 2012). Sloping croplands underwent the greatest level of human activity, such as ploughing, rotation, fertilizing, and trampling. Disturbances in sloping croplands and crop rotation could introduce new microbial species into the soil and decrease the abundance of previously existing ones. Consequently, the number of microbial species in sloping croplands increased, while having an insignificant impact on soil function due to low abundance. Although previous studies have suggested that rare taxa can exhibit unique metabolic activities (Bickel & Or, 2021; Campbell et al., 2011) and are involved in particular metabolic functions, dominant species play a greater role in ecosystem functioning than subdominant and rare species (Treplin, Pennings & Zimmer, 2013; Xu et al., 2021b; Xu et al., 2021c). The flat croplands were agricultural lands that experienced trampling when oranges and tea were harvested, and the intensity of this type of disturbance was lower than that in sloping croplands. The forests should have the lowest level of disturbance caused by human activity. Thus, the introduction of non-native bacterial species would be lower in forest soils compared to flat croplands. Hence, soils in sloping croplands harbored more diverse bacterial species than flat croplands and forests, where the core amplicons were in higher relative abundances. These findings indicated that with a growing intensity of disturbance, soil microbial diversity and community composition might vary significantly.

The random forest and path modeling results showed that human activity had negative effects on soil fertility, and such negative effects were mitigated by the bacterial community composition. Similarly, the bacterial community composition was found to mitigate global change impacts (nitrogen addition, altered rainfall frequency) on ecosystem multifunctionality in a biocrust system (Liu et al., 2017). In addition, the results of path modeling signified that among all the factors we studied here, microbial community composition was the key factor mediating soil fertility. Analogously, the microbial community composition is also suggested to be the dominant factor controlling soil respiration rates in a cultivated land (Liu et al., 2018). Our study further indicated that the amplicon relative abundances of the core bacterial taxa (i.e., Proteobacteria and Acidobacteriota) were important factors driving soil fertility, and these two phyla were more abundant in soils with higher fertility. Acidobacteriota are positively correlated with soil organic carbon (Vasileiadis et al., 2013). These taxa are well-known copiotrophs that thrive in nutrient-rich soil conditions (Yao et al., 2014; Yao et al., 2017; Xu et al., 2020). In addition, Proteobacteria are predicted to be copiotrophs (Goldfarb et al., 2011; Tong et al., 2021), which grow quickly and increase soil fertility through plant root-associated metabolic processes (Oh et al., 2012). The high amplicon relative abundance of Acidobacteria and Proteobacteria, along with their importance in mediating soil fertility, further supported the notions that abundant taxa are active and important in mediating biogeochemical cycling (Cottrell & Kirchman, 2003; Treplin, Pennings & Zimmer, 2013; Xu et al., 2021b; Xu et al., 2021c).