For the hereby work, six representative eco-regions of China were sampled for study, specifically: (i) the karst area of southwestern China (KS), (ii) the agro-pastoral ecotone of southwestern China (AS), (iii) the Qinghai–Tibet Plateau (QT), (iv) the Loess Plateau (LP), the (v) forest-grassland ecotone (FG), and the (vi) deserts of Hexi Corridors (HC) (Supplementary Figure S1). The selected areas have a mean annual temperature spanning from −22.9°C to 28.6°C, a mean annual precipitation of 44.9–1815 mm, and an elevation that ranges from 27 m to 8,305 m a.s.l. The region exhibits distinct climatic and vegetation zones from east to west, shaped by climatic conditions and soil characteristics: subtropical monsoon climate (primarily encompassing KS and parts of AS) and plateau mountain climate (spanning parts of AS and QT) dominate eastern areas. KS and AS feature expansive yellow-brown, red, and cinnamon soils, supporting subtropical evergreen-deciduous broad-leaved forests and alpine meadows. QT’s black felty soils and chernozem foster alpine grasslands. Transitioning westward, temperate continental monsoon climate (LP, FG, HC) emerges. LP’s brown calcic, aeolian sandy, and loessal soils sustain temperate grasslands, sandy semi-shrub grasslands, and warm temperate deciduous broad-leaved forests. FG’s boggy/meadow soils and dark-brown/aeolian sandy soils host similar warm temperate deciduous broad-leaved forests. The westernmost Calcic-Orthic Aridosols in HC create desert grassland ecosystems (Duan et al., 2025).

Soil sampling took place from June to July 2019 and during July 2020. Based on the mean annual precipitation gradient, a total of n = 41 locations (6 locations for KS, and likewise 6 for AS, 9 for QT, 5 for LP, 6 for FG, and 9 for HC) were chosen at intervals of approximately 100–110 km along four transects (Supplementary Figure S1). All samples collected from the minimally disturbed natural soils were collected during the peak of the growing season. At each location, surface litter was removed within a 10 × 10 m plot, then along the plot’s diagonal line five 1 × 1 m quadrats were established. From each quadrat, along that diagonal line, three replicate soil samples (0–20 cm depth) were obtained and homogenized to provide a single composite soil sample per plot. All composite soil samples were individually packed in sterilized polyethylene bags, and taken rapidly to the lab using portable refrigerators. Each composite soil sample was then split into two subsamples: one that was stored at 4°C for later biochemical analysis, the other was stored at −80°C prior DNA extraction and molecular analysis. All relevant variables and location information of our study’s soil samples are detailed in Supplementary Table S1.

The MAP (mean annual precipitation), MAT (mean annual temperature), and NDVI (normalized difference vegetation index) data for the wider experimental area were obtained from the Chinese Meteorological Database.¹ Soil pH was measured using an E20-FiveEasy pH meter (Mettler Toledo, Giessen, Germany). EC (electrical conductivity), an indicator of the soluble salt content of soil, was measured using an electric conductometer. Both soil parameters were measured using a soil-water suspension (5:1, v/v mixture of deionized water and fresh soil) after the samples were subjected to shaking for 30 min. Both the content of SOC (soil organic carbon) and TN (soil total nitrogen) were quantified using a carbon-hydrogen-nitrogen elemental analyzer (2400 II CHN Elemental Analyzer, Perkin Elmer, Boston, MA, United States). Climatic factors and soil properties datasets are summarized in Supplementary Table S1.

Total genomic DNA was extracted from each soil sample using the E.Z.N.A.^® Soil DNA Kit (Omega Bio-Tek, Norcross, GA, United States) following the manufacturer’s instructions. The concentration and purity of the extracted DNA was determined using TBS-380 and NanoDrop2000 spectrophotometers. DNA extract quality was checked on 1% agarose gel.

The DNA extracts were fragmented to an average size of about 400 bp, using the Covaris M220 (Gene Company Limited, China). For the paired-end library construction, NEXTFLEX Rapid DNA-Seq (Bio Scientific, Austin, TX, United States) was used. Adapters containing the full complement of sequencing primer hybridization sites were ligated to the blunt-end of each fragment. Next, the paired-end sequencing was carried out on an Illumina HiSeq Xten system (Illumina Inc., San Diego, CA, United States) at the Majorbio Bio-Pharm Technology Co., Ltd. (Shanghai, China) by using HiSeq X Reagent Kits and following the manufacturer’s instructions.²

The generated data from sequencing were analyzed on the free Majorbio Cloud Platform online³ (Ren et al., 2022). Briefly, the paired-end Illumina reads were trimmed of their adaptors, and any low-quality reads (i.e., having a length <50 bp, or a quality value <20, or N bases) were removed by the fastp (v 0.20.0) tool⁴ (Chen et al., 2018). The resulting high-quality read pairs from the same sample were then assembled into contigs by using MEGAHIT (v1.1.2) software⁵ (Li et al., 2015) with kmer values ranging from 47 to 97 (step = 10). Contigs having a length ≥300 bp were selected for gene prediction and functional annotation.

An open reading frame (ORF) for each contig was predicted by MetaGene⁶ (Noguchi et al., 2006). Those predicted ORFs with a length ≥100 bp were retrieved and translated into amino acid sequences using the NCBI translation table. Next, a non-redundant gene catalog was constructed using CD-HIT (v4.6.1)⁷ (Fu et al., 2012) based on a minimal 90% sequence identity and 90% coverage. To calculate the abundance of genes at a 95% identity threshold, the obtained high-quality reads were aligned to non-redundant gene catalogs via a SOAP aligner (v 2.21)⁸ (Li et al., 2008). For their taxonomic identification, Diamond (v0.8.35)⁹ (Buchfink et al., 2015) was used with an e-value <1 × 10⁻⁵ and the alignments searched against the NCBI microbial NR database. The predicted gene fragments were searched against KEGG,¹⁰ NCycDB (Tu et al., 2019) and SCycDB (Yu et al., 2021) reference databases using Diamond (e-value <1 × 10⁻⁵) for functional annotation. Hmmscan¹¹ was used to search against the Carbohydrate-Active Enzymes (CAZy) database¹² (e-value cutoff = <1 × 10⁻⁵).

To minimize the effects of sequencing depth on statistical analysis, the calculation method of species and gene abundance is Reads Per Kilobase Million (RPKM) (Ye et al., 2024):

where R_i represents the abundance value of Gene _i in a given sample, i.e., the number of Reads compared to Gene _i in that sample; L_i means the nucleotide length of Gene _i ; and ∑ 1 n R j represents the sum of reads corresponding to all genes in that sample.

Microsoft Excel 2019 and R (v 3.2.1) were used for statistical data analyses. The differences in Chao1 index for the α-diversity of taxonomic (species) and functional traits [KOs (KEGG Orthology) and CAZyme-encoding genes] across the six eco-regions were performed by using the nonparametric Kruskal–Wallis test and the post-hoc pairwise Wilcoxon rank-sum test, with p-values adjusted through the Benjamini–Hochberg method. The α-diversity¹³ estimates were calculated using the diversity function of the “vegan” package¹⁴ in the R computing platform (v 3.2.1)¹⁵ (Oksanen, 2017). Kruskal–Wallis test and Wilcoxon rank-sum test were performed through functions “kruskal.test” and “wilcox.test” in package “stats” in R (v 3.2.1). To identify taxa and metabolic pathways as biomarkers, the linear discriminant analysis (LDA) effect size (LEfSe) algorithm (LDA >2.5, p < 0.05) (Segata et al., 2011) was applied. The abundance (i.e., RPKM, reads per kilobase per million mapped reads) of all microbial groups (bacteria, archaea, viruses, and eukaryota) at phylum level was visualized in histogram plot (s). RPKM abundances of differentially abundant (DA) KEGG pathways, normalized by z-score across all data sets, visualized using heatmap, and the sample-based dendrogram was performed using Euclidean distance metric in package “pheatmap”¹⁶ in R (v 3.2.1). Principal co-ordinates analysis (PCoA) of taxonomic traits (species) and functional traits (KOs, KEGG Orthology database), based on their Bray–Curtis (BC) distances, was carried out using the “vegan” package (see text footnote 14) in R (v 3.2.1) (Oksanen, 2017). The differences of taxonomic traits (species) and functional traits among six eco-regions were tested using permutational multivariate analysis of variance (PERMANOVA) with 9,999 permutations. Distance-based redundancy analysis (db-RDA) was implemented to explore the effects of environmental factors on the SOC’s decomposition and biosynthesis, N cycling, and S cycling across six eco-regions using the Bray–Curtis dissimilarities distance matrix. The individual effects of each explanatory variable on response variables from db-RDA were estimated using the hierarchical partitioning method (“rdacca.hp” function from rdacca.hp. package) (Lai et al., 2022).

Topsoil samples were collected from n = 41 locations, covering six distinct geographical regions (Supplementary Figure S1 and Supplementary Table S1). 579.3 gigabases (Gb) of paired-end sequence data were generated averaging 93.6 million paired reads per sample. De novo assembly of sequencing data yielded a non-redundant gene catalog for all 41 locations. The total length of this non-redundant assembly was 24.6 Gb (means contig N50 length of 543 bp), from which 54.7 million partial genes >100 bp were predicted. After removing redundancy by clustering the genes by identity (>90%), and by shortening gene coverage (>90%), a total of 25.4 million non-redundant genes were deposited in the metagenomic libraries. Bacterial genes were the most predominant among all species comprising 97.69% of all sequences, followed by 2.22% for archaea, 0.07% for eukaryota, and only 0.02% for viruses (Figure 1a and Supplementary Table S2). Our results showed that the number of genes sequenced had their lowest abundance in the Hexi Corridor deserts (region HC) (1632676). Regarding the five other eco-regions, gene abundance ranged as FG, LP, QT, AS, and KS from lowest to highest, respectively. The low abundance at the desert ecosystem could be attributed to the harsh conditions present, characterized by extreme drought, extreme temperature variation (e.g., great thermal difference between day and night), and low soil fertility (Coleine et al., 2024; D’Odorico et al., 2013; Duan et al., 2022).

The Chao1 index values for taxonomic (species) and functional (KOs, KEGG Orthology) traits for all microbial groups from the trialed eco-regions are shown in Figures 1b,c respectively. For taxonomic traits, the Chao1 richness of the entire soil microbiome ranged from 12383.95 (at AS group) to 13569.53 (at HC group) with an overall mean (±SD) and median of 12961.03 ± 551.81 and 12938.09, respectively (Supplementary Tables S3, S5). For functional traits, whole microbiome’s Chao1 richness based on the KOs ranged from 6540.04 (at LP group) to 7505.37 (at HC group) having an overall mean of 6875.65 ± 476.89 and a median of 6771.20 across the six regions (Supplementary Tables S4, S6). Significant differences on soil microbiome’s Chao1 indices (species or KOs) between the selected region samples (Supplementary Tables S5, S6). Overall, the Chao1 index for the α-diversity of taxonomic and functional traits at HC and FG was higher than that from other regions. Interestingly, although soil microbial abundance and diversity in desert ecosystems are generally lower than that in less stressful environments, Chao1 richness of taxonomic traits (species) and their functional traits (KOs) in the HC region exceeded that from the other five regions (KS, AS, QT, LP, and FG) (Figures 1b,c).

Principal coordinates analysis (PCoA) revealed significant differences at the β-diversity of both taxonomic (species) (PERMANOVA, R² = 0.415, p = 0.001; Figure 1d) and functional traits (KOs) (PERMANOVA, R² = 0.432, p = 0.001; Figure 1e) between different regions. Evidently, whether taxonomic or functional, the microbiome traits from the five regions (KS, AS, QT, LP, and FG) were geographically separated along the first principal coordinate from those traits of HC. Overall, the β-diversity of taxonomic and functional traits at HC was obviously different from other regions (Figures 1d,e). Again, these pronounced differences are mainly attributed to the unique conditions of the desert habitat (D’Odorico et al., 2013; Hu et al., 2019).

Taxonomic analysis identified a total of 105 phyla of soil microorganisms, specifically: 82 for bacteria, 12 for archaea, 10 for eukaryota, plus 1 virus (Figure 2a and Supplementary Table S7). Most of the metagenomes were dominated by Actinobacteria (18.01–62.09%) and Proteobacteria (14.66–50.31%), followed by Acidobacteria (0.61–13.10%). In addition, Chloroflexi (1.37–8.01%), Firmicutes (1.19–3.78%), Bacteria_unclassified (0.96–5.42%), Gemmatimonadetes (0.65–5.49%), Cyanobacteria (0.89–21.54%), Planctomycetes (0.72–2.54%), Euryarchaeota (0.19–27.94%), Verrucomicrobia (0.23–3.85%), Bacteroidetes (0.47–6.38%), Thaumarchaeota (0.05–3.90%), Candidatus_Rokubacteria (0.05–3.45%), Candidatus_Tectomicrobia (0.14–2.36%), Nitrospirae (0.15–1.70%), and Deinococcus-Thermus (0.22–1.51%), together accounted for ca. 98% of all metagenomic sequences derived from the topsoil samples (Supplementary Table S7). Overall, these bacterial phyla exhibit no site-specificity across China and are considered common in terrestrial ecosystems on a global scale (Delgado-Baquerizo et al., 2018).

The RPKM abundance data were screened for the KEGG pathways enriched in the six eco-regions, observation based on the linear discriminant analysis (LDA) effect size (LEfSe) method (LDA >2.5, p < 0.05). Collectively, 76 KEGG pathways (18.5% of total) were enriched in the six regions, with 4, 27, 34, 4, 7, and 0 pathways enriched in LP, KS, HC, FG, AS, and QT regions, respectively (Figure 2b and Supplementary Table S8). These 76 KEGG pathways could be classified into 23 classes that are mainly associated with seven biological processes (Figure 2c and Supplementary Table S8): carbohydrate metabolism amino acid metabolism xenobiotics biodegradation and metabolism, energy metabolism, metabolism of cofactors and vitamins, replication and repair lipid metabolism (Figure 2c and Supplementary Table S8). The RPKM-normalized read counts notably differed (more than five-fold) across the six regions, being 7.29 × 10⁶, 6.27 × 10⁶, 1.22 × 10⁷, 8.63 × 10⁶, 6.44 × 10⁶, and 6.37 × 10⁶ in LP, KS, HC, FG, AS, and QT, respectively (Supplementary Table S8). Additionally, the hereby results reveal that different eco-regions may form different functional metabolic niches. The heatmap analysis reveals that the three KS, AS, and QT regions in southern China are clustered together; likewise, the HC, LP, and FG in northern China are clustered together (Figure 2c and Supplementary Table S8).

To identify the functional genes involved in C, N, and S cycling, the obtained soil metagenomic reads were annotated using the databases of CAZy (Cantarel et al., 2008), KEGG (Kanehisa et al., 2008), NCycDB (Tu et al., 2019) and SCycDB (Yu et al., 2021).

For carbohydrate metabolism, the 41 soil samples contained six critical CAZyme-encoding genes involved in SOC decomposition and biosynthesis (Figure 3; Supplementary Figure S2 and Supplementary Table S9). The genes had uneven distribution among the samples. Specifically, as gene and location in brackets: glycoside hydrolases (GHs), glycosyl transferases (GTs), polysaccharide lyases (PLs), carbohydrate esterases (CEs), auxiliary activities (AAs), and carbohydrate-binding modules (CBMs). This differentiation generally reflects the microbial substrate affinity, and as a proxy, the decomposition and biosynthesis potential of each community against SOC (Kanehisa and Goto, 2000). Evidently, the genes encoding the GHs (organic carbon decomposition) and GTs (organic carbon biosynthesis) enzymes were the most abundant in all samples across the six regions (Figure 3; Supplementary Figure S2), followed by those coding for CEs and AAs, and last but not least for CBMs and PLs. In addition, there has been considerable variation in the abundance of genes participating in SOC decomposition and biosynthesis across all 41 samples (ranging from 61267.39 to 77190.28) (Figure 3; Supplementary Figure S2 and Supplementary Table S9).

For nitrogen (N) metabolism, biogeochemical N cycling between inventories is often attributed to the six common N-transformation processes: nitrogen fixation, nitrification, denitrification, nitrate reduction, nitrogen transport/nitrate assimilation, and organic nitrogen metabolism (Kelly et al., 2021; Kuypers et al., 2018) (see the N-cycling model in Figure 4 and Table 1). Hereby, 62 gene families responsible for N cycling were detected from all samples across the six regions; these gene families could be categorized into 7 individual pathways (Figure 4 and Table 1; Supplementary Tables S10, S11). (1) Nitrogen fixation: 3 key gene families (nifD, nifH, and nifK) encoding the nitrogenase complex were identified in all regional groups (KS, AS, QT, LP, FG, and HC samples), nifW gene family was not present in either the KS or AS group. The highest and lowest abundance of gene families involved in nitrogen fixation were found in HC (95.52) and AS groups (11.68) respectively; the former nearly 9 times greater than the latter. Overall, the abundance of these four gene families was measurable yet consistently low across all 41 soil sites; (2) Nitrification: 7 key gene families (amoB_A, amoB_B, amoC_A, amoC_B, nxrB, and hao) encoding ammonia monooxygenase were identified in all six regions. amoA_A gene and nxrA gene families were identified in KS, AS, and QT. The gene families involved in nitrification were most abundant in the QT group (40.86), the least for HC (10.98) with the former 3.72 times higher than the latter. Among them, amoB_A and hao were the two most abundant gene families participating in the nitrification process; (3) Denitrification: 15 marker gene families (nirK, nosZ, napA, nirS, norB, narG, narH, narZ, norC, narJ, narI, napC, napB, and narV) encoding key enzymes for denitrification were identified in all six regions. One gene family narW was only identified in KS and AS. These gene families attained their highest abundance in the QT group (2217.52), being lowest in the HC (1529.79) with the former 1.45 times higher than the latter. The three most abundant gene families involved in the denitrification process were nirK, nosZ, and napA; (4) Nitrate reduction: 5 key functional gene families (nirB, nirD, nrfC, nrfA, and nrfD) encoding key enzymes for DNRA were identified in all six regions. One gene family nrfB was not present in HC. The highest abundance of these gene families involved in the DNRA pathway was found in the LP group (1674.78) and the lowest in the FG group (1283.15), with the former 1.07 times higher than the latter. Among them, nirB was the gene family present in greatest abundance. Meanwhile, 6 marker gene families (nasA, NR, nirA, narB, narC, and nasB) encoding key enzymes for ANRA pathway were detected in all six regions, with nasA and NR being the most abundant Furthermore, the highest abundance of gene families participating in the ANRA pathway occurred in the KS group (3789.56) and the lowest in FG (2883.24), the former was 1.31 times higher than the latter; (5) Nitrogen transport: 4 marker gene families (NRT, nrtA, nrtB, and nrtC) encoding key enzymes for nitrogen transport/nitrate assimilation were identified in all six regions. These gene families were found to be the most and least abundant in the KS (699.27) and HC (324.96) groups, respectively, the former being 2.15 times higher than the latter; (6) Organic nitrogen metabolism: 17 marker gene families (glnA, gs_K00264, gs_K00265, gs_K00266, gs_K00284, nmo, asnB, gdh_K00260, gdh_K00261, gdh_K00262, gdh_K15371, glsA, ureA, ureB, ureC, ansB and nao) encoding key enzymes for organic nitrogen metabolism were identified in all six regions. The highest abundance of these gene families was found in the LP group (18472.78) and the lowest in the HC group (14146.65); the former was 1.31 times higher than the latter. Of these, glnA and gs_K00266 were the two most abundant genes involved in the organic nitrogen metabolism process. Thus, the significant differences in the total gene abundances of N metabolism among the different eco-regions were summarized, the list findings indicate that the soil microbiome of LP is the environment with the highest gene abundance for N metabolism (26120.01), followed by that of QT (25974.50), KS (25623.967), and AS (24699.21), with FG (22462.91) and HC (20375.07) having the lowest abundance (Figure 4 and Supplementary Table S10).

For sulfur (S) metabolism, biogeochemical cycling of S between inventories is often attributed to three distinct sulfate-transforming processes (Figure 5 and Table 2): sulfur assimilation, anaerobic sulfate respiration, and sulfide oxidation (Llorens-Marès et al., 2015) (see the S-cycling model in Figure 5). In this study, from all soil samples across the six regions, a total of 15 gene families responsible for S cycling were detected. These gene families could be categorized into six individual pathways according to SCycDB (Yu et al., 2021) (Figure 5 and Table 2; Supplementary Tables S12, S13). (1) Sulfur reduction: 5 marker genes (sudA, ttrB, sudB, hydG, and sreB) encoding key enzymes participating in the Sulfur reduction pathway were identified in all six eco-regions. One gene family psrC with low abundance was not present in QT and FG. The highest abundance of these genes was found in the KS group (837.24) and the lowest in the HC group (614.17). Of these, sudA was the most abundant gene family associated with the S reduction pathway, 6.07 times of the total abundance of other genes in this pathway. (2) Sulfur oxidation: 1 marker gene family (soeB) encoding for key enzymes related to sulfite oxidation was identified in all six regions. Overall, the abundance of the gene family was measurable yet consistently low across all 41 sites. (3) Sulfur disproportionation: 1 marker gene family (phsB) with low abundance across all 41 sites. (4) SOX systems: 1 marker gene family (soxC) encoding for key enzymes related to thiosulfate oxidation was identified in all six regions. The highest and lowest abundance of these was found in the LP (108.81) and AS group (85.55) group respectively, the former being 1.27 times higher than the latter. (5) Dissimilatory sulfur reduction and oxidation: 2 marker gene families (dsrO and dsrL) with low abundance across all 41 soil sites. (6) Assimilatory sulfate reduction: 4 marker gene families (cysC, cysH, cysJ, and sir) encoding for key enzymes participating in the sulfate reduction were identified in all six eco-regions. The highest abundance of these genes was found in the LP group (736.94) and the lowest in the AS group (569.23). Of these, cysC was the most abundant gene family, 1.7 times more abundant than sir and cycH, 3.06 times more abundant than cycJ. Altogether, significant differences in the total abundances of genes related to S metabolism among different regions were detected. Specifically, the LP (1587.46) and KS (1506.16) regional groups have the highest gene abundance for S metabolism, followed by QT (1496.76) and AS (1429.73), with FG (1420.98) and HC (1341.04) having the lowest abundance (Figure 5 and Supplementary Tables S12, S13).

To understand the relationship between the microbial abundance of the functional genes involved in C, N, and S cycling, across the six eco-regions in China, when subjected to various environmental factors, distance-based redundancy analysis (db-RDA) was used (Figure 6). The db-RDA showed that the first two axes accounted for 70.19% of the variability of the microbial genes composition involved in the C cycling, whereas db-RDA1 (the x-axis) and db-RDA2 (the y-axis) accounted for 57.67 and 13.52% of the variation, respectively (Figure 6a). For N cycling, the db-RDA showed that the first two axes accounted for 81.16% of the variability of the microbial genes composition, whereas db-RDA1 (the x-axis) and db-RDA2 (the y-axis) accounted for 65.97 and 15.19% of the variation, respectively (Figure 6b). As for the S cycling, the db-RDA showed that the first two axes accounted for 79.05% of the variability of the microbial genes composition, whereas db-RDA1 (the x-axis) and db-RDA2 (the y-axis) accounted for 53.33 and 25.72% of the variation, respectively (Figure 6c). It was interesting that latitude and seven soil properties (including pH, NDVI, MAP, SOC, TN, C:N ratio, and EC) were significantly correlated with the first two axes (p < 0.01) (Figure 6 and Supplementary Table S14) of the C, N, and S cycling, and these environmental factors explained 82.21, 79.13 and 81.45% (Figure 6d and Supplementary Table S14) of the variations in microbial genes composition involved in the C, N, and S cycling, respectively.