A total of 88 soil samples were collected from agricultural fields spanning the main cereal-producing areas of China, including 36 black soils in Northeast China, 22 fluvo-aquic soils in the North China Plain, and 30 red soils in southern China (Supplementary Figure S1). Considering the major crop in each area, the soils were planted by the major crop in each area during the sampling time, which was summer maize in Northeast China and North China Plain, while paddy rice in South China. Mean annual precipitation (MAP) and mean annual temperature (MAT) of the sampling sites were obtained from the China Meteorological Data Service Center (https://data.cma.cn/en).

Five 1 m × 1 m sampling squares were assigned to each site in a plot of 20 m × 20 m. Three topsoil cores (0 cm−20 cm) were randomly collected from each sampling square during the crop growth season (June–July 2017). All the soil cores from the same sampling site (5 × 3 = 15) were combined into a single sample and transported to the laboratory on ice within 48 h. The fresh soils were sieved through a 2 mm sieve after the rocks and plant residues were removed by hand. To minimize the effects of weather, agricultural management, and plant growth, the soils were incubated for 2 weeks at 25°C under controlled water conditions in the laboratory (upland soils were adjusted to 45% water-holding capacity and paddy soils were flooded using the deionized water with a 2 cm water layer). After incubation, soil samples were harvested and divided into two fractions: one fraction was frozen in liquid nitrogen and stored at −80°C for DNA extraction and the other was air dried for the analyses of soil physiochemical properties.

The soil texture was analyzed using the fixed pipette method (Dane and Topp, 2002). SOC was determined using the heat K₂Cr₇O₂-H₂SO₄ titration method (Bao, 2000). The soil total nitrogen content was determined using the Kjeldahl method (Bremner, 1965). Soil-available phosphorus (AP) was extracted using 0.5 M NaHCO₃ and determined using the ammonium molybdate method (Olsen, 1954). Soil pH was determined using a pH meter (Mettler Toledo, Zurich, Switzerland) at a water-to-soil ratio of 1:2.5.

Total soil DNA was extracted using a Mobio PowerSoil DNA kit (Mo Bio Laboratories, Inc., Carlsbad, CA, USA) according to the manufacturer's instructions. DNA quality and concentration were quantified using a NanoDrop spectrophotometer (NanoDrop Technologies; Wilmington, DE, USA) and via the electrophoresis on a 1% agarose gel. The isolated DNA was stored at −20°C for subsequent analyses.

Real-time quantitative PCR analysis was performed to determine the abundance of aerobic bacteria using gltA as a marker gene. The primers used were CS680F (5′AYG CCG ABC AYG ARY WSA A-3′) and CS904R (5′TAS ACS SGR TGR CCR AAG CCC AT) (Castro et al., 2012). The PCR system included 5 μL of 2 × SYBR PreMix ExTag (Takara, Japan), 1 μL each of forward and reverse primers (10 mmol/L), 1 μL of DNA template (diluted to 10 ng/μL), and 2 μL of ddH₂O. The PCR conditions were as follows: 95°C for 5 min, followed by 40 cycles of 95°C for 30 s, 56°C for 30 s, and 72°C for 30 s. The fluorescence signal was collected at the end of each cycle. A standard curve was prepared by 10-fold gradient dilution of plasmids containing the target fragment and was considered validated when R² > 0.98, with an amplification efficiency of 90%−110%.

The gltA gene from each sample was amplified using the CS680F/CS904R primer set and specific barcodes. The PCR system included 12.5 μL of 2 × PreMix Ex Taq (Takara, Japan), 2 μL of each forward and reverse primer, 1 μL of DNA template (diluted to 25 ng/μL), and 7.5 μL of ddH₂O. The PCR conditions were as follows: 95°C for 5 min; 30 cycles of 95°C for 30 s, 56°C for 30 s, and 72°C for 30 s; and a final extension at 72°C for 10 min. PCR products were purified using a Tiangen gel extraction kit (Tiangen, Beijing, China) and sequenced on an Illumina HiSeq PE2500 platform (Illumina, San Diego, CA, USA).

Pair-end raw reads were assembled, screened, and trimmed using QIIME (Caporaso et al., 2010). Low-quality tags were removed as well as those that were < 180 bp or more than 250 bp. The clean sequences were exposed to a chimera filter and frameshift correction in the RDP Fungene pipeline (Fish et al., 2013), and subsequently clustered into OTUs at 97% similarity. Only OTUs affiliated with more than two sequences observed in more than two samples were retained. A representative sequence was selected for each OTU, and BLAST against the GenBank non-redundant nucleotide database was used for taxonomic assignment. All datasets were normalized to the same sequencing depth (20,000 sequences per sample) for subsequent analyses. Raw reads were deposited in the NCBI Sequence Read Archive under the accession number PRJNA1067286.

All statistical analyses were performed using R (v4.0.2; http://www.r-project.org/), except for the kriging interpolation algorithm. Ordinary kriging interpolation in ArcGIS 10.7 software was used to draw the spatial distributions of the abundance and richness of gltA-harboring bacteria. OTUs richness calculation and PCA analysis were performed using the estimateR and rda function “vegan” package, respectively (Oksanen et al., 2022). Differences in abundance and richness of gltA-harboring bacteria communities in the different group sequences were calculated using the Kruskal-Wallis test in the “EasyStat” package (Wen, 2020). Permutational multivariate analysis of variance (PERMANOVA) was used to evaluate significant differences between the groups. Both the dissimilarity matrix of the gltA-harboring bacteria community and soil texture were calculated using the vegdist function in the “vegan” package (Oksanen et al., 2022), and linearized to calculate the relationship between the two indices.

To evaluate the relative importance of environmental factors on the gltA abundance, richness, and community composition of gltA-harboring bacteria, a random forest model was performed using the “rfPermute” packages (Archer, 2022). Linear discriminant analysis (LDA) effect size (LEfSe; Segata et al., 2011) was performed using the Galaxy online pipeline (http://huttenhower.sph.harvard.edu/galaxy/) to discover marker taxa in different soil types using default parameters.

The zonal distribution of the abundance (indicated by gltA copy number) and diversity (indicated by gltA richness) of gltA-harboring bacteria was revealed using Kriging interpolation (Supplementary Figure S2). The gltA abundance and diversity in the samples ranged from 1.3 × 10⁷-1.2 × 10¹⁰ g⁻¹ and 276–563, respectively (Supplementary Figures S2A, B), and unimodally changed with increasing latitude, with both peaking at the middle latitude of ~35° (N) (Figures 1a1, b1). A similar trend was observed for gltA abundance with increasing longitude but not for gltA diversity (Figures 1a2, b2). Species annotation detected 1144 gltA-harboring OTUs affiliated with 36 phyla (Supplementary Figure S3). On average, Proteobacteria, Acidobacteria, Actinobacteria, and Verrucomicrobia, accounted for 68% of the total reads in all samples. However, these phyla were less dominant in the central region than in the northern or southern regions (Supplementary Figure S3). As shown in Figures 1c1–f1, c2–f2, Proteobacteria and Actinobacteria presented a latitude-associated zonal distribution with the highest relative abundance at 31°-32° (N), whereas a longitude-associated zonal distribution was found for Actinobacteria and Verrucomicrobia, with the highest relative abundance at ~113° (E).

Soil type remarkably affected gltA abundance, with the highest content in fluvo-aquic soil (3.3 × 10⁹ copies g⁻¹), followed by red soils (1.5 × 10⁹ copies g⁻¹), and the lowest content in black (0.8 × 10⁹ copies g⁻¹) (Figure 2a1, p < 0.05). Similar results were observed for gltA diversity, with no significant differences between the red and black soils (Figure 2b1). However, a comparison between paddy and upland soils revealed no significant difference in either the abundance or diversity of gltA (Figures 2a2, b2). Bray-Cutis dissimilarity indicated significant differences in the overall community composition of gltA-harboring bacteria among the soil types, as well as between upland and paddy soils (Figures 2c1, c2). Nevertheless, PCoA and PERMANOVA analyses revealed a significant effect of soil type on gltA-harboring bacterial communities (Supplementary Figure S4, R² = 0.432, p = 0.001). Additionally, soil type affected the relative abundance of the dominant phyla more strongly than land use. As shown in Figures 2c1–f1, significant differences among the soil types were observed for all four dominant phyla, except Acidobacteria. Generally, Proteobacteria and Actinobacteria were more abundant in the black soil, whereas Verrucomibia was more abundant in the red soil. In contrast, difference in the relative abundance of the dominant phyla between the upland and paddy soils, except more abundant Actinobacteria in paddy soils than in upland soils (Figures 2c2–f2).

gltA copy number and richness increased linearly along the PC1 axis of soil texture (Figures 3a, b). Similarly, a significant positive correlation was observed between the Bray and Curtis distance of gltA-harboring bacterial community and soil texture (Figure 3c). These results indicated that soil texture was tightly correlated with the community of gltA-harboring bacteria. Further analysis revealed that the abundance, diversity, and community composition of gltA-harboring bacteria were most closely related to clay content, followed by sand, but not to silt content (Figures 3d–i). In particular, both the copy number and richness of gltA significantly increased with the relative clay content (R² = 0.176–0.303, p < 0.0001), but decreased with that of sand (R² = 0.069–0.087, p < 0.05). The random forest model explained 34.0%, 36.1%, and 81.5% of the variation in gltA abundance (Figure 4a), richness (Figure 4b), and community composition of gltA-harboring bacteria (Figure 4c), respectively. Among all climatic and soil properties, clay content was the most important predictor of gltA abundance, followed by annual precipitation and sand content (p < 0.01). Similarly, clay content was a significant predictor of the richness and community composition of gltA-harboring bacteria (Figures 4b, c, p < 0.01).

LEfSe analysis revealed that 39, 13, and 13 biomarkers were enriched in black, fluvo-aquic, and red soils, respectively (Figure 5). Most of the gltA-harboring bacteria enrichen in black soil and were affiliated with Proteobacteria (30.8%) and Actinobacteria (28.2%), followed by Acidobacteria (15.4%) and Gemmatimonadetes (10.3%). Planctomycetota and Bateroidetes were uniquely enriched in fluvo-aquic and red soils, respectively. The enrichment of Verrucomicrobiota was also found in red (23.1%) and black (5.1%) soils, but not in fluvo-aquic soil. The relative abundances of the biomarkers were strongly correlated with the contents of the different soil aggregate fractions (Figure 5). Generally, the biomarkers in black soil were positively correlated with the contents of sand and clay but negatively correlated with silt, whereas the biomarkers in fluvo-aquic soil were negatively correlated with clay and silt but positively correlated with sand, and the biomarkers in red soil were positively and negatively correlated with silt and sand.