Soil Autotrophic Bacterial Community Structure and Carbon Utilization Are Regulated by Soil Disturbance—The Case of a 19-Year Field Study

Abstract

The roles of bacterial communities in the health of soil microenvironments can be more adequately defined through longer-term soil management options. Carbon dioxide (CO₂) fixation by autotrophic bacteria is a principal factor in soil carbon cycles. However, the information is limited to how conservation tillage practices alter soil physiochemical properties, autotrophic bacterial communities, and microbial catabolic diversity. In this study, we determined the changes in autotrophic bacterial communities and carbon substrate utilization in response to different soil management practices. A replicated field study was established in 2001, with the following soil treatments arranged in a randomized complete block: conventional tillage with crop residue removed (T), conventional tillage with residue incorporated into the soil (TS), no tillage with crop residue removed (NT), and no tillage with residue remaining on the soil surface (NTS). Soils were sampled in 2019 and microbial DNA was analyzed using high-throughput sequencing. After the 19-year (2001–2019) treatments, the soils with conservation tillage (NTS and NT) increased the soil’s microbial biomass carbon by 13%, organic carbon by 5%, and total nitrogen by 16% compared to conventional tillage (T and TS). The NTS treatment increased the abundance of the cbbL gene by 53% in the soil compared with the other soil treatments. The cbbL-carrying bacterial community was mainly affiliated with the phylum Proteobacteria and Actinobacteria, accounting for 56–85% of the community. Retaining crop residue in the field (NTS and TS) enhanced community-level physiological profiles by 31% and carbon substrate utilization by 32% compared to those without residue retention (T and NT). The 19 years of soil management lead to the conclusion that minimal soil disturbance, coupled with crop residue retention, shaped autotrophic bacterial phylogenetics, modified soil physicochemical properties, and created a microenvironment that favored CO₂-fixing activity and increased soil productivity.

1. Introduction

Soil microbial communities serve as the reservoir of biological processes and play a critical role in converting organic carbon into plant nutrients [1,2]. It is estimated that approximately 80 to 90% of the soil biogeochemical processes are mediated by microbial activities, including organic carbon degradation [3]. An essential biosynthetic process in soil is carbon dioxide (CO₂) fixation by autotrophic bacteria, a key process of carbon cycling [4] which is responsible for the sequestration of atmospheric CO₂ into the soil [5]. In agroecosystems, autotrophic CO₂-fixing bacteria are present in various types, including distinct physiological groups, such as cyanobacteria and purple photosynthetic and chemoautotrophic bacteria [6]. The groups of bacteria can incorporate ¹⁴C into their microbial biomass through a series of biochemical reactions [7,8]. Those reactions are related to the community composition of autotrophic bacteria [9], influencing the intensity of CO₂ fixation.

The Calvin–Benson–Bassham cycle (i.e., the CBB cycle) is the most important autotrophic CO₂ fixation pathway [10,11] which influences the primary carbon source in an ecosystem [12]. Ribulose-1,5-bisphosphate carboxylase/oxygenase (known as RubisCO), a key enzyme in the CBB cycle, regulates the process of CO₂ reduction and the oxygenolysis of ribulose-1,5-bisphosphate [10,11]. RubisCO exists in four forms (referred as I, II, III, and IV) with different structures, catalytic activities, and O₂ sensitivities, of which form I’s RubisCO mainly occurs in plant and photo- and chemo-autotrophic bacteria [13,14]. In the plant–soil–microbiome environment, autotrophic microbial diversity can be reflected by the functional biomarkers (i.e., the cbbL genes). These genes encode the large subunit of RubisCO form I that has a large scale of sequences for phylogenetic analysis [15,16].

Many anthropogenic activities, such as tillage and crop rotation, affect the composition and diversity of the microbial community in soil [17,18]. Soil tillage affects soil water conservation [19], fertility enrichment [20], and carbon sequestration [21,22]. Reduced or zero tillage (i.e., minimal to no disturbance to the soil profile) coupled with crop residue retention can increase the nutrient supply potential [23] and enhance carbon cycling [24] while decreasing the carbon footprint [25,26], which helps alleviate the challenge of global climate change [27]. Diversified crop rotation increases soil carbon [28,29], improves microbial diversity [30], and enhances the system’s productivity and resiliency. However, little is known about the long-term impact of multiple years of soil tillage and cereal–legume alternate rotation on the function of autotrophic CO₂-fixing bacteria. It is unclear how the nexus of autotrophic bacterial community composition—CO₂ fixation—tillage practices impact soil microenvironments.

To elucidate those effects, we established the field experiment in 2001 on the western Loess Plateau of China. Different tillage treatments have since been applied to soil that has a cereal–legume crop rotation system. After the 19-year treatments (in 2019), soil was sampled from each plot in three replicates and analyzed for its biophysiochemical properties. The abundance and community composition of autotrophic CO₂-fixing bacterial communities were analyzed using real-time quantitative PCR and high-throughput sequencing. Soil carbon substrate utilization and its relation to autotrophic bacterial community were also determined. In a previous paper [31], we reported the diversity of bacterial communities and the characterization of the phylogenic composition in relation to soil management. In the present paper, we present the findings on (1) the changes in soil physicochemical parameters, autotrophic bacterial community composition, and soil catabolic diversity in response to the 19-year tillage treatments, and (2) the possible mechanisms responsible for shifts in the soil’s autotrophic bacterial communities in relation to tillage. These findings from the long-term field study bring us to a new level of understanding of the interaction among autotrophic bacterial community, carbon fixation, and soil management practices.

2. Materials and Methods

2.1. Site Description

A field experiment was started in 2001 at the Rainfed Agricultural Education Center (35°28′ N, 104°44′ E; 1971 m a.s.l.) of Gansu Agricultural University in China. The site is on the western Loess Plateau and has a temperate, semiarid, continental monsoon climate [32]. Its solar radiation is 5.68 KJ m⁻² and sunshine duration is 2476 h annually [33]. The mean temperature is 6.5 °C and the 10 °C-based accumulative temperature is 2339 °C. The average precipitation was 392 mm per year, with 54% occurring from June to September in most years. The experimental plots were on level terrain with dark loessal Calcaric Cambisol soil (FAO/UNESCO soil classification system, 1990). Prior to the start of the experiment in 2001, the site had a long history of the continuous cropping of flax (Linum usitatissimum L.) under conventional tillage practices [34].

2.2. Experimental Design

A randomized complete block design was used to accommodate the four soil management treatments in each of the three replicates. The description of the treatments and the implementation of the tillage practices and crop residue retention are summarized in

Table 1. Treatment structure and the description of the tillage practices and crop residue implementation in the long-term field experiment started in 2001.

Treatment Abbreviation a	Tillage	Crop Residue
T	Conventional tillage, plowed to the 20 cm depth in the fall and harrowed the following spring	Crop residue was removed out of the field at harvest
NT	No tillage	Crop residue was removed out of the field at harvest
TS	Conventional tillage, plowed to the 20 cm depth in the fall and harrowed the following spring	Crop straw was chopped and incorporated in the soil via the fall plowing
NTS	No tillage	Crop straw was chopped, spread, and remained on the soil’s surface

^a The treatment descriptions and abbreviations will be used in the other Tables and Figures throughout the article.

. More details on the tillage and field operations of the treatments can be viewed in a previous publication [34]. In brief, all the treatments were implemented in both the wheat (Triticum aestivum L.) and field pea (Pisum arvense L.) phases of the rotation each year. For the TS and T treatments, the soil was plowed to a depth of 0.20 m in the fall and harrowed in the following spring’s seeding time. The experiment included four tillage and crop residue management treatments per spring. In each of the three replicates, the following treatments were arranged in plot sizes of 80 m² (4 m wide × 20 m long). The crops were managed according to the recommended agronomic practices. The crops were planted in late March to early April, varying slightly from one year to the other due to weather conditions. The seeding rate was 187.5 kg ha⁻¹ for wheat and 180 kg ha⁻¹ for pea, on average, varying slightly each year due to seed size and percent germination. The seeding rates were aimed to achieve a plant population of 250 plants per m⁻² for wheat and 80 plants per m⁻² for field pea. The wheat was fertilized with urea nitrogen (46% N) at 105 kg N ha⁻¹ and the pea at 20 kg N ha⁻¹. Both crops received a phosphorus fertilizer (Ca(H₂PO₄)2.2H₂O, with 16% P₂O₅) at 46 kg P₂O₅ ha⁻¹. The fertilizers were broadcast at sowing. No potassium (K) was applied as the soil contained exchangeable K greater than 300 mg kg⁻¹. No herbicides were applied to any crop because of the dry weather and little weed pressure. No irrigation was applied. At maturity (usually in mid-July to early August), the entire plot was harvested using a plot harvester and the biomass and yields were recorded.

2.3. Soil Sampling and Physicochemical Property Measurement

Bulk soils were sampled in June 2019 when pea plants were at the mid-flowering stage and wheat plants were at the late stage of flowering. During flowering, crop plants are most active, with roots interacting with soil microbes under semiarid environments [35,36]. Sampling at this stage can provide researchers with invaluable information about rhizospheres and their interactions with environmental factors [37]. From each plot, five soil samples were randomly taken from a 0 to 0.2 m depth using a 50 mm diameter iron soil corer, and the five samples were bulked for each plot. After passing through a 2 mm mesh sieve, two subsamples were taken from each plot: one was for DNA extraction and carbon substrate utilization assessment and the other was for soil property analysis.

The soil pH was determined with a pH meter (Mettler Toledo FE20, Shanghai, China) using a deionized suspension with 1:2.5 soil: water ratio (mass: volume). Soil organic carbon (SOC) was analyzed using the Walkley–Black wet oxidation method and total N (TN) was measured with the Kjeldahl method. Olsen phosphorus was examined using the colorimetric method with 0.5 M NaHCO₃. A UV-1800 spectrophotometer (Mapada Instruments, Shanghai, China) was used to determine soil NH₄⁺-N and NO₃⁻-N. Soil microbial biomass carbon and nitrogen (SMBC and SMBN, respectively) were measured using the Chloroform method. In each plot, soil moisture was measured gravimetrically and converted into volumetric water content using soil bulk density, which was measured by sampling the intact core of the soil using known-volume metal rings, drying the soil, and weighing the dried soil.

2.4. Soil Microbial Catabolic Diversity

For each soil sample from each plot (i.e., three replicates), community-level physiological profiles (CLPP) were determined using Biolog Eco-PlatesTM (Biolog Inc., Hayward, CA, USA) which contained 31 widely used C substrates along with a control well. The 31 C sources included nine carboxylic acids, seven carbohydrates, three other substrates, six amino acids, and two amines/amides [38]. Each well of the Eco-plates was inoculated with 150 μL of soil suspensions before being incubated at 25 °C. The carbon substrate utilization based on well color development was measured with optical density (OD) values at a wavelength of 590 nm (color development plus turbidity) and 750 nm (turbidity only) every 24 h for 7 d using the Biolog method (Biolog Inc., Hayward, CA, USA). Average well-color development (AWCD) was used to assess the microbial metabolic activity, as follows [39]:

$$\mathrm{AWCD}={\displaystyle \sum }\frac{C-R}{31},$$

where C is the OD value of each well and R is the absorbance reading of the control well. Substrate richness and evenness, defined as the Shannon–Weaver index (H’), were calculated as follows:

$$H^{\prime }=-{{\displaystyle \sum }}_{i=n}{p}_i(ln{p}_i),$$

where pi represents the ratio of color development of a certain well to the sum of the color development of all wells in a micro-plate.

2.5. DNA Extraction and Quantification of cbbL Gene Abundance

Soil microbial DNA was extracted from each approximately 0.5 g of fresh soil sample, in triplicate, using a DNeasy PowerSoil Kit (QIAGEN, Inc., Hilden, Germany ) following the manufacturer’s protocol. The concentrations of the extracted DNA were determined with a NanoDrop ND-1000 spectrophotometer (Thermo Fisher Scientific, Waltham, MA, USA) and the quality of the extracted DNA was determined using agarose gel electrophoresis.

The cbbL gene copy number was assessed using real-time quantitative PCR (qPCR) on a LightCycler 480 II real-time PCR system (Roche Diagnostics, Mannheim, Germany). The cbbL genes were amplified with the forward primer K2f (5′-ACCAYCAAGCCSAAGCTSGG-3′) and the reverse primer V2r (5′-GCCTTCSAGCTTGCCSACCRC-3′) (reference). The reaction mixture contained 1 μL of DNA template, 5 μL of 2 × LightCycler^® 480 SYBR Green I Master, 0.2 μL of forward primer, 0.2 μL of reverse primer, and 3.6 μL of nuclease-free water, and it made up a final volume of 10 μL. The thermal protocol was 10 min of incubation at 95 °C, followed by 40 cycles of 10 s at 95 °C and 30 s at 60 °C. Melting curve analysis was performed to validate the specificity of the amplified PCR products. The standard curves were produced using triplicate 10-fold dilutions of plasmid DNA with inserted cbbL gene. PCR amplification efficiencies were 94% with an R² value of 0.99.

2.6. High-Throughput Sequencing of the cbbL Gene

The primer set (K2f-V2r) for cbbL gene real-time qPCR was used for the cbbL gene amplification for each sample, respectively. 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 (5 U/μL), 2 μL (2.5 mM) of dNTPs, 1 μL (10 uM) of each forward and reverse primer, 2 μL of DNA template, and 8.75 μL of ddH₂O. The PCRs were performed with the following procedures: 2 min of initial denaturation at 98 °C, followed by 25 cycles of 15 s at 98 °C, 30 s for annealing at 55 °C, 30 s for elongation at 72 °C, and a final extension at 72 °C for 5 min. PCR amplicons were extracted from 1% agarose gels and purified with Agencourt AMpure beads (Beckham Coulter, Indianapolis, IN, USA) and quantified using the PicoGreen dsDNA assay kit (Invitrogen, Carlsbad, CA, USA). High-throughput sequencing data were obtained through commercial laboratory services (Personal Biotechnology Co., Shanghai, China) and deposited in the NCBI database with accession PRJNA689959.

2.7. Data Processing and Bioinformatics Analysis

Raw fastq files were quality-filtered and merged using VSEARCH pipeline (version 2.13.4) on the platform of QIIME 2. The sequences with a quality score <20, having ambiguous bases, or those containing mononucleotides repeats of >8 bp were removed. The operational taxonomic units (OTUs) were clustered with a 97% sequence similarity cut-off using VSEARCH pipeline (v. 2.13.4). All representative reads in the OTUs were taxonomically classified and annotated based on the GenBank^® nucleotide sequence database (http://ncbi.nlm.nih.gov; accessed on 15 March 2020).

2.8. Statistical Analysis

Data were analyzed using SPSS software (SPSS Inc., Chicago, IL, USA). Significant differences between treatments in soil properties, cbbL gene abundance, and catabolic diversity were determined at the probability of <0.05.

The QIIME2 package was employed to determine alpha diversity, and the community evenness and richness were expressed using the Shannon and Simpson indexes. The variation of cbbL-harboring bacterial community composition across the samples was determined using principal co-ordinates analysis (PCoA) based on the genus-level compositional profiles. Hierarchical clustering analysis was performed with the unweighted pair-group method of arithmetic means (UPGMA) based on Bray–Curtis distance matrices. The differences in the taxonomic and phylogenetic communities were compared by clustering them using the R “vegan” package (R Foundation for Statistical Computing, Vienna, Austria). The significant differences in microbial structure among treatments was determined by permutational multivariate analysis of variance (PERMANOVA) and analysis of similarities (ANOSIM).

Network analysis was performed to investigate the co-occurrence patterns among the cbbL-harboring bacterial communities. The OTUs with a relative abundance of higher than 0.01% were retained and the “psych” package of R version 4.0.2 was used to analyze the preprocessed data and calculate the spearman correlation coefficient matrix. Statistical correlations with a cut-off at an absolute r value of higher than 0.6 and a p value of lower than 0.05 were retained for further analysis. The software Gephi (version 0.9.2) was employed for network visualization and network properties measurement. To identify the keystone bacterial taxa driving community turnover under different treatments, the abundance of bacterial taxonomic OTUs was regressed against the four field treatments based on random forest classification analysis using the “randomForest” package of R version 4.0.2. Lists of biomarkers taxa were ranked to compare the importance of OTUs using 10-fold cross-validation, implemented using R ‘rfcv’ function. The top 15 bacterial biomarkers were chosen according to the stabilized cross-validation error curve. The relative abundance of biomarker taxa was illustrated using the “heatmap” package of R version 4.0.2. The Spearman’s correlation coefficients and the Mantel test with 999 permutations in the “vegan” package of R were used to assess the correlations of soil properties and the distribution of the dominant genera of cbbL-harboring bacteria based on the distance matrix.

3. Results

3.1. Soil Physiochemical and Biological Properties

The soil’s biophysiochemical properties significantly differed between treatments (

Table 2. Biophysiochemical properties of the soils under different management practices.

Soil Parameter a	Treatment b
	NT	NTS	T	TS
pH (H2O)	8.13 ± 0.03 a c	8.10 ± 0.03 a	8.13 ± 0.02 a	8.14 ± 0.07 a
SOC (g kg−1)	12.77 ± 0.61 ab	14.03 ± 0.40 a	12.37 ± 0.81 b	13.13 ± 0.60 ab
TN (g kg−1)	0.86 ± 0.07 ab	1.07 ± 0.13 a	0.78 ± 0.08 b	0.88 ± 0.09 ab
TP (g kg−1)	0.56 ± 0.02 a	0.66 ± 0.16 a	0.54 ± 0.06 a	0.55 ± 0.05 a
TK (g kg−1)	18.61 ± 0.90 ab	19.29 ± 0.44 a	17.39 ± 0.77 b	18.86 ± 0.67 ab
Olsen p (mg kg−1)	13.43 ± 2.71 a	15.89 ± 1.95 a	14.04 ± 4.30 a	15.38 ± 1.03 a
NO3−-N (mg kg−1)	37.54 ± 3.38 a	35.79 ± 1.02 a	36.26 ± 3.95 a	33.84 ± 3.20 a
NH4+–N (mg kg−1)	1.28 ± 0.29 a	1.48 ± 0.09 a	1.40 ± 0.13 a	1.20 ± 0.13 a
Moisture (%)	13.57 ± 0.43 a	14.76 ± 0.97 a	13.56 ± 0.35 a	13.87 ± 0.44 a
Bulk density (g cm−3)	1.25 ± 0.03 b	1.19 ± 0.02 c	1.31 ± 0.02 a	1.23 ± 0.03 bc
Total porosity (%)	53.05 ± 0.28 a	53.52 ± 0.79 a	50.43 ± 0.72 b	51.25 ± 0.32 b
SC (mm h−1)	77.16 ± 2.24 b	86.28 ± 2.09 a	73.05 ± 0.53 c	86.21 ± 1.45 a
Temperature (°C)	23.20 ± 0.96 a	23.63 ± 0.70 a	23.10 ± 0.66 a	24.03 ± 0.55 a
SMBC (mg C kg−1)	185.51 ± 8.82 ab	217.12 ± 11.35 a	164.69 ± 12.31 b	192.75 ± 16.92 ab
SMBN (mg N kg−1)	19.27 ± 0.66 b	22.01 ± 1.61 a	18.22 ± 0.61 b	18.20 ± 0.53 b

^a SOC, soil organic carbon; TN, total nitrogen; TP, total phosphorus; TK, total potassium; SC, saturation conductivity; SMBC, soil microbial biomass carbon; SMBN, soil microbial biomass nitrogen. ^b The treatment descriptions are summarized in Table 1. ^c Data (means ± SD, n = 3) in each row followed by different letters are significantly different at p < 0.05.

). The contents of SOC, TN, TK, SMBC, and SMBN were significantly greater for the NT and NTS treatments compared to the T and TS treatments (p < 0.05). Soil organic carbon was 3.2%, 13.4%, and 6.1% greater for the NT, NTS, and TS treatments, respectively, compared with the T treatment. The highest soil TN concentration (1.07 g kg⁻¹) was obtained in the NTS soil. The TK for the NTS treatment was 10.9% (p < 0.05) greater than that for the NT, TS, and T treatments. The soil’s bulk density was 4.6%, 9.2%, and 6.1% lower for the NT, NTS, and TS treatments, respectively, compared with the T treatment. The NT, NTS, and TS treatments increased SMBC by 12.6%, 31.8%, and 17.0%, respectively, compared to the T treatment. In addition, SMBN for the NTS treatment was greater (p < 0.05) than that for the T treatment (p < 0.05), with no difference among NT, TS, and T. There was no difference in soil pH, TP, NO₃⁻-N, NH₄⁺–N, soil temperature, and moisture among the treatments (p > 0.05).

3.2. cbbL Gene Abundance and the Relation to Soil Properties

The soils under the NTS treatment had more copy numbers of cbbL genes than the soils under the other treatments, and the increment was 52% when NTS was compared with the T treatment (p < 0.01;

Table 3. The cbbL gene copy numbers of the soil samples under the different soil management practices.

Treatment Abbreviation a	cbbL Abundance (108 Copies g−1 Soil)
NT	2.17 ± 0.42 b b
NTS	3.51 ± 0.36 a
T	2.31 ± 0.30 b
TS	2.39 ± 0.27 b

^a The treatment descriptions are summarized in Table 1. ^b Values in each column followed by different letters are significantly different at p < 0.05.

). Correlation analysis across the treatments showed that the soil cbbL gene was positively correlated with soil moisture, saturation conductivity, SMBC, and SMBN, but it was negatively correlated with bulk density (p < 0.05; Table S1).

3.3. Autotrophic Bacterial Community Diversity and OTUs Richness

After quality control, a total of 468,535 filtered sequences with 32,250 to 43,646 sequences per sample remained (Table S2). The rarefaction curves were close to the saturation phase with the increase in sample size. There was sufficient sequencing depth and the OTUs were representative of the overall bacterial community libraries. To compare soil carbon fixation and microbial community diversity, the same survey effort level of 30,637 sequences was randomly selected from each soil sample (95% of the smallest number of reads). In total, 4988 OTUs were observed across all samples, with the number of OTUs ranging from 503 to 1692, varying with the soil samples (Table S2). The four tillage and crop residue management treatments shared 451 OTUs, accounting for 9.04% of the core OTUs (Figure 1a). There were differences in the OTU-sharing between treatments: the NTS and NT treatments shared the highest number of OTUs (1107 OTUs, 22.19%), whereas the TS and T treatments shared the least number of OTUs (683 OTUs, 13.69%) and the NTS and T treatments shared 18.04% of the OTUs (900 OTUs). Moreover, the NTS treatment harbored the highest number of unique OTUs (1458 OTUs, 29.23%) and the T treatment harbored the lowest number of unique OTUs (308 OTUs, 6.17%).

The alpha diversity of bacterial populations differed significantly among soil treatments, as shown by the Chao1, Simpson, Shannon, and observed OTUs indexes (p < 0.05; Figure 1b). The Shannon index for the NTS treatment was significantly higher compared to the T treatment and the diversity was ranked NTS > NT > TS > T. No differences in the Chao1 and observed OTUs indexes were observed among the four treatments (p > 0.05; Figure 1b). The PCoA results showed similarities and differences in CO₂-fixing bacterial community structure across the different treatments (Figure 1c), where total variance explained by PCoA1 and PCoA2 was 35.5 and 13.6%, respectively. The bacterial communities of the soil under the NT and TS treatments were clustered closely, whereas the community compositions in the soil under the NTS and T treatments differed from those of the NT and TS treatments. Permutational multivariate analysis of variance (PERMANOVA) and analysis of similarities (ANOSIM) confirmed the significant differences in the compositions of the soil cbbL-carrying bacterial communities under different tillage practices and crop residue managements (p < 0.05; Table S3).

3.4. Composition and Networks of Soil Autotrophic Bacteria Communities

The composition of the bacterial communities was determined for soils from each of the three field replicates (e.g., T1, T2, and T3 for the T treatments in Figure 2 and Table S4). At the phylum level, the bacterial communities were dominated primarily by the phyla Proteobacteria (38.52–80.94%) and Actinobacteria (5.72–22.65%), both accounting for an average of 56–85% of the relative abundance of the bacterial communities (Figure 2a; Table S4). The phyla Euryarchaeota and Cyanobacteria had low relative abundances. The tillage treatment had a significant (p < 0.05) impact on the relative abundance of a bacterial community. The relative abundance of the phylum Proteobacteria was 34.2, 41.6, and 33.5% lower in the soil under NT, NTS, and TS, respectively, while the relative abundance of the phylum Actinobacteria was 158.4, 245.7, and 142.5% higher, respectively, than that measured in the soil under the T treatment.

At the class level, alpha-Proteobacteria and beta-Proteobacteria dominated the bacterial communities, containing 14.89–30.71% and 5.71–22.65% of the total cbbL gene sequences, respectively (Figure 2b; Table S4). The class Rubrobacteria was the most dominant within the phylum Actinobacteria, which accounted for 0.57–4.00% of the total relative abundance. At the genus level, the cbbL-carrying bacteria were mainly affiliated with Bradyrhizobium, Variovorax, Gaiella, and Steroidobacter, accounting for 20–35% of the identified communities (Figure 2c; Table S4). The tillage treatments had an impact on the community composition, as shown in the top 10 genera; compared to the T treatment, the NTS treatment increased the relative abundance of genus Bradyrhizobium by 59.3% and reduced the relative abundance of genus Variovorax by 60.7%. However, no significant differences in the taxonomic composition were observed among the NT, TS, and T treatments. The hierarchical cluster analysis (Figure 2c) had results similar to those obtained from PERMANOVA and ANOSIM (Table S4), where the bacterial communities formed two main clusters (T and TS vs. NT and NTS), and the largest differences were between the NTS and T treatments.

Co-occurrence network analyses were conducted to evaluate the different co-occurrence patterns of the cbbL-harboring bacterial communities under the different treatments (Figure 3). In total, there were 35 nodes and 71 links, with an average degree of 4.057 and a graph density of 0.119 in the autotrophic bacterial network (Table S5). The ratio of positive edges was higher than the ratio of negative edges in the cbbL-harboring bacterial communities, with the ratio of positive edges accounting for 85.92%. The genera Methylocella and Advenella were the most connected keystone genera in the cbbL-harboring bacterial communities, and both were from the phylum Proteobacteria, with significant correlations with other genera (Figure 3; Table S5).

3.5. Biomarkers in the Autotrophic Bacterial Communities

To determine the key biomarkers discriminating the cbbL-harboring bacterial communities under different treatments, we performed a regression of the abundances of bacterial OTUs against the tillage treatments based on the random forests machine learning algorithm. The top 15 keystone OTUs were chosen as the respective bio-marker taxa because the cross-validation error curve was stabilized when using these OTUs for both bacterial communities (Figure 4). The biomarker taxa were mainly affiliated with the classes of alpha-Proteobacteria and beta-Proteobacteria in the cbbL-harboring bacterial communities (Table S6). The heatmaps in Figure 4 further illustrate the scaled abundances of these biomarker taxa in response to soil tillage management. It was found that the phylum Proteobacteria played more crucial roles than the other phyla across soil samples for the soil’s cbbL-harboring bacterial communities.

3.6. Soil Microbial Catabolic Diversity and Carbon Utilization Pattern

The soil’s microbial metabolic activity was expressed as the average well-color development (AWCD), an indicative measurement of the physiological profiles at the community level. The AWCD differed significantly between the treatments (Table S7). The NTS and TS treatments increased the AWCD index by 21.6 and 19.1%, respectively, while the NT treatment reduced the AWCD index by 15.4% in comparison with the T treatment. However, no significant differences for the H’ (Shannon–Wiener) index were observed among the four treatments.

The soil’s microbial carbon utilization patterns based on six fundamental groups differed among the treatments (Figure 5a). Compared to the T treatment, NT and TS reduced carbon resource utilization based on amines/amides by 63.6 and 35.2%, respectively, while NTS and TS increased carboxylic acid utilization by 37.1 and 44.9%, respectively. The overall carbon utilization capabilities differed significantly among the six carbon substrate groups (Figure 5b). The carbon utilization for carbohydrates, carboxylic acids, and amino acids were significantly higher than that for the polymers, miscellaneous, and amines/amides groups. The highest carbon utilization was obtained in the carbohydrates group, with a mean value of 8.808, while the lowest carbon utilization was found in the amines/amides group, with a mean value of 0.815.

The carbon utilization capability of soil microbiomes, based on seven specifical substrates (L-Arginine, L-Phenylalanine, α-D-Lactose, N-Acetyl-D-Glucosamine, D-Mannitol, D-Malic Acid, and α-Cyclodextrin), differed significantly among the tillage treatments (Figure 5c). The utilization of L-Arginine substrate in the TS treatment was 73.9% higher (p < 0.05) than that in the T treatment. The NTS treatment increased the utilization of L-Phenylalanine substrate by 85.5%, compared with the T treatment (p < 0.05). The NT treatment reduced the utilization of α-D-Lactose and N-Acetyl-D-Glucosamine substrates by 28.9 and 34.4%, respectively, compared with the T treatment (p < 0.05). The NT treatment increased the utilization of D-Mannitol substrate by 18.67% (p < 0.05) compared with the NTS, TS, and T treatments, with the latter three having a similar utilization. For the utilization of D-Malic Acid substrate, the NTS and TS treatments respectively had 38.8 and 45.6% greater (p < 0.05) utilization than the T treatment. The utilization of α-Cyclodextrin substrate was 68.9% greater in the NTS and 102.7% greater in the TS treatments compared to the T treatment, while the NT treatment reduced the α-Cyclodextrin substrate utilization by 64.8% compared with the T treatment (p < 0.05). Across the four tillage treatments, the mean carbon utilization capabilities of the soil microbiomes differed significantly among the 31 carbon substrates. D-Mannitol substrate had the highest carbon utilization, averaging 1.624, and 2-Hydroxy benzoic Acid substrate had the lowest utilization, averaging 0.023.

3.7. Correlation among Physiochemical Traits, Catabolic Diversity, and the Autotrophic Bacterial Community

There were significant (p < 0.05) correlations among the soil physicochemical properties and cbbL-harboring bacterial community compositions (Figure 6a). Spearman correlation analysis across the tillage treatments revealed that SOC was positively correlated with TN, saturation conductivity, SMBC, and SMBN, but it was negatively correlated with bulk density (Table S8). Total soil N was correlated positively with TK, saturation conductivity, and SMBN, while soil bulk density was negatively correlated with saturation conductivity, SMBC, and SMBN. Further, the Mantel test revealed that the structure of the autotrophic bacterial communities was highly correlated with SOC (r = 0.315, p < 0.05), TK (r = 0.328, p < 0.05), bulk density (r = 0.724, p < 0.05), total porosity (r = 0.415, p < 0.05), saturation conductivity (r = 0.349, p < 0.05), and SMBC (r =0.377, p < 0.05) (Figure 6a).

Redundancy analysis, showing the relationship between carbon source utilization and the structure of cbbL-harboring communities at the genus level, revealed that the community composition varied from 14.20 to 20.68% (Figure 6b). The carbon utilization of polymers, amino acids, carbohydrates, and carboxylic acid were clustered to the edge of the genus Bradyrhizobium and Methyloversatilis, and they were negatively correlated to the genus Rhodopila. Among the carbon substrates, polymers explained 16.6% of the variance (F = 2.0, p = 0.178) and amino acids explained 8.1% of the variance (F = 0.9, p = 0.414) across the soil samples, and they were the key indicators of carbon utilization capacity (with 999 permutations) (Table S9). The abundance and diversity of the genus Gaiella were closely related to miscellaneous. The 19 years of tillage treatments had an impact on bacterial community structure, with the genera Bradyrhizobium and Gaiella presenting highly in the soil under the NTS treatment, while the genus Variovorax was high in the soil with continuous T practices.

4. Discussion

4.1. Diversity of Soil cbbL-Carrying Bacterial Communities Shaped by Crop Residue Retention

The biomarker cbbL genes, which encode a large subunit of ribulose-1,5-bisphosphate (RUbP) carboxylase/oxygenase in the Calvin–Benson–Bassham cycle, are considered important indicators of CO₂-fixing autotrophic microbes in soil. We found a significant shift in the phylogenetic diversity of cbbL-carrying bacteria in response to soil management, and the largest differences in diversity existed between the NTS and the T treatments. The NTS treatment reduced the relative abundance of Proteobacteria while increasing the abundance of Actinobacteria at the phylum level, and at the genus level, the NTS treatment increased the relative abundance of Bradyrhizobium while reducing the relative abundance of Variovorax. The mechanisms for the shift in the bacterial community composition due to different soil management treatments are unclear, but we suggest the following possible mechanisms: (a) the minimal physical disturbance to the soil, coupled with continuous crop residue input (in each growing season, year-after-year), provides cumulative feedback to the soil’s microenvironment, forming a balanced C to N ratio over the years, which favors microbial activities. A high N availability in soil is found to promote Proteobacteria while prohibiting Actinobacteria. A high C availability in soil is found to stimulate Proteobacteria while prohibiting other microbes compared to low N or low C soils. (b) The 19 years of continuous input of crop residue as organic material in the NTS treatment increased C sources for energy, essential for the cell synthesis of organisms [17], leading to the increased abundance of cbbL genes encoding ribulose-1,5-bisphosphate carboxylase oxygenase for CO₂ fixation. The decomposition of the input material increased the substrate sources, leading to the great shift in the taxonomic composition of cbbL-carrying bacterial communities. (c) The shift in the microbial community composition due to soil management led to the formation of a microenvironment favorable to the humification process of SOC, increasing the size of the organic C pools. Additionally, (d) the facultative CO₂-fixing bacterial members in soil may grow heterotrophically with the addition of crop residue (besides autotrophic activities) [40,41]. Nevertheless, our results, in combination with the findings reported by others [42,43,44], bring us to a new level of understanding that long-term crop residue input with minimal disturbance can shape microbial community composition, benefiting CO₂ fixation processes and soil C cycling.

4.2. cbbL-Carrying Bacterial Catabolic Diversity in Relation to Soil Management

Bacterial OTU richness and catabolic diversity were highest in the soil under no till with crop residue retention among the soil treatments evaluated, suggesting that soil management does alter the diversity of soil microbiota, provided that the soil is treated for a longer period of time (19 years in the present study). Other researchers have found that the catabolic diversity is reflected by the shift in microbial community composition accompanied by the different capacities of carbon substrate utilization [45]. However, the catabolic diversity found in bulk soils may not be detectable in rhizosphere soils [46], since the expression of genes for soil nutrient reduction varies with soil microenvironments. Further, we found that the 19 year soil treatments did not affect soil pH (Table 1). The cbbL-carrying taxa were mainly affiliated with the types of bacteria, regardless of the level of pH tolerance. The diversity of autotrophic bacterial communities was highest in soil under the NTS treatment and lowest in soil under the T treatment, mainly due to the continuous lack of disturbance and crop residue input that created an ideal microenvironment that favored microbial activities. Our results support the findings of some researchers that soil pH does not have a direct role in impacting bacterial community diversity [47,48], and they do not agree with some others [49,50]. We suggest that the relationship between soil pH and bacterial community diversity is likely driven by pH-related soil nutrients. For example, the heavy use of N fertilizers can reduce the diversity and abundance of cbbL gene-containing microorganisms [51], mainly due to soil N decreasing pH, which prohibits some beneficial microbial activities [17,52].

4.3. Relationships among Physiochemical Properties, CO₂ Fixation, and Bacterial Diversity

Soil is a comprehensive ecosystem in which microbial communities have a strong correlation with soil management options, host plants, and edaphic factors [53,54], including nutrient availability, ion exchange capacity, and SOC turnover [55,56]. In our study, the 19 year soil management practices led to substantial changes in physiochemical and biological properties. No till coupled with crop residue retention significantly increased cbbL gene abundance, SOC, SMBC, and SMBN while decreasing bulk density compared to conventional tillage. Our results add new knowledge that in a no-till soil microenvironment, bulk density, SMBC, and SMBN impact the cbbL gene abundance, and that the increased cbbL gene abundance, in turn, increases CO₂ fixation, affirming the positive role of cbbL-carrying bacteria in soil’s C cycling. The cbbL-carrying facultative autotrophs played a more prominent role than the obligate autotrophs in soil CO₂ fixation, as we identified both facultative and obligate bacteria through sequencing, and we found that the relative abundance of Proteobacteria at the phylum level was highest in all the soil treatments evaluated. Proteobacteria is a well-known facultative aerobic autotroph that assimilates CO₂ via the Calvin–Benson cycle. However, the CO₂-fixing capability of autotrophs may vary with other factors, such as plant species [57], crop fertilization [51], and soil amendments [58].

4.4. Implication of Long-Term Field Studies

Bacteria play an important role in soil C cycling, some of which serve as the decomposers of root exudates, plant litter, and other carbon sources which break down the organisms and convert organic energy into forms useful to the rest of the organisms in the food web. The conversion process releases nutrients and thus improves soil fertility. However, limited information is available regarding how the bacterial community composition and diversity are influenced by long-term soil management options, such as tillage practices and crop residue input. In the scientific literature, some aspects of the relationship between soil bacterial communities and anthropogenic activities are documented, for example, soil tillage practices shape microbial community composition [59,60], change microbial diversity [46,61], and provide feedback to soil physiochemical properties [62]. However, most published findings are primarily drawn from short-term (< 10 years) experiments, such that the conclusions are often inconclusive, inconsistent, or even controversial from one study to another. Running a large field experiment continuously, year after year, is undoubtedly costly, but the experience of conducting the 19-year field experiment has provided us with the high confidence for making concrete conclusions on the subject area.

The improved soil microenvironmental conditions due to no till and crop straw retention for 19 years has led to increased soil productivity. Pea yield ranged from 1030 to 1980 kg ha⁻¹ and aboveground plant biomass ranged from 4300 to 6400 kg ha⁻¹ (Figure S1). Spring wheat followed a similar trend of treatment effect as that obtained in field pea. The effect of soil management on crop yields in 2019 was similar to that averaged cross the 19 study years; the NTS, TS, and NT treatments increased (p < 0.01) yields by 32.5, 25.1, and 10.0%, on average, of wheat and pea, respectively, compared to the T treatment. The soils receiving residue retention treatments (i.e., NTS and TS) had significantly higher productivity (22.7% more) compared to the soil without crop residue retention (i.e., NT and T).

Knowledge of the complex nexuses among soil’s physiochemical properties, microbiomes, and anthropogenic activities (such as tillage practices and crop residue input) is highly demanded in the development of strategies to improve the health of soil microenvironments. Our findings are well positioned to meet the needs. One of the important aspects in the assessment of the impact of long-term soil management practices is to quantify the relative weights or percent impact between the two subfactors—tillage practice and crop straw retention. There is a need to understand whether the combination of the two subfactors actually generates additional synergy between the two. For example, the NTS treatment in our study had the greatest benefits in terms of soil microbial diversity, C-cycling, and nutrient enrichment compared to the rest of the treatments; however, we were unable to differentiate the relational weight of the impact between tillage practice and straw retention. Did the benefits stem from the NTS treatment due to straw input (x% of the benefits), no tillage (y% of the benefits), or the combination of the two? Nevertheless, our long-term experiment has been running since 2001 at a site representative of the 1.56 million hectares of the China Loess Plateau. The results of the study can provide some guidelines for improving soil health and productivity in local regions, and possibly extended to areas with similar climatic conditions to the experimental site.

5. Conclusions

The results of the 19-year field experiment show that continuous disturbance to the soil profile elevated the selective pressure to the soil CO₂-fixing bacteria community, disturbed microbial catabolic activity, and created hardship for the microenvironment. Conversely, no disturbance to the soil, coupled with residue-retaining on the soil surface or crop residue being incorporated into the soil, increased the soil’s organic carbon and nutrients and lowered the soil’s bulk density. The underlying micro-driven functional changes led to enhanced soil carbon utilization. We conclude that bacterial community structure, composition, and CO₂-fixing capability are highly regulated by soil management practices, and that a minimal disturbance to the soil’s microenvironment, coupled with the retention of crop residues in the soil, will improve bacteria-involved biological activities and increase nutrient cycling and soil productivity.