We sampled soils in the Xilingol Grassland, Inner Mongolia, China (45°50′7.74 ″N, 116°33′21.19 ″E). The average annual rainfall and temperature are ∼300 mm and 1.6°C. Past overgrazing resulted in severe grassland degradation in this region. Based on the grassland degradation index (Dong et al., 2019), non-degraded (ungrazed for the last 10 years, GDI = 0.93), lightly degraded (GDI = 0.75), moderately degraded (GDI = 0.63) and severely degraded (GDI = 0.52) sites were chosen for study. All sites occurred within a 5 km × 5 km area to minimize differences in climate and original vegetation composition among them. Vegetation characteristics including plant species, biomass, abundance, vegetation cover, and soil properties including total nitrogen, total carbon, total phosphorus, total potassium, pH and NO³⁻ -N, were also determined and listed in Supplementary Table S1.

At each site we established three plots (1 × 1 m), and from each plot we randomly sampled five soil cores (diameter 4 cm) after removing the humus layer. Two soil depths (0–10 cm and 10–20-cm) were sampled separately but for each depth the five soil cores were pooled to form a single independent sample per plot during the growing season (August) of 2021. Twenty-four samples were obtained (4 sites × 2 depths × 3 replicate plots). Soil samples (200 g) were sieved using a 2-mm mesh to remove debris and stones. Sieved samples were stored in a −80°C freezer prior to DNA extraction and at 4°C prior to chemical analyses. We also recorded the richness and abundance of each plant species from the 1 × 1 m sampling plots.

A FastDNA SPIN Kit (MP Biomedicals, Santa Ana, CA, United States) was used to obtain soil genomic DNA from 500 mg soil samples based on the manufacturer’s protocol. We diluted the extracted DNA 10 ngμL⁻¹ using a spectrophotometer (NanoDrop ND-1000, NanoDrop Technologies), which was then fragmented to ∼350 bp using an ultrasonicator (Covaris M220, Gene Company Limited, China) to construct paired-end libraries using a TruSeq DNA Sample Prep Kit (Illumina, San Diego, CA, United States). Then, blunt-end fragments were linked with adapters possessing the full complement of sequencing primer hybridization sites. The Illumina HiSeq 4000 PE150 platform was employed for pair-end sequencing (Illumina Inc., San Diego, CA, United States) based on the manufacturer instructions of the HiSeq 3000/4000 PE Cluster Kit and HiSeq 3000/4000 SBS Kits.

Trimmomatic 0.36, a read-trimming tool, was used to process the raw data (Bolger et al., 2014) to remove adapters and reads of moderate quality. Sequences with Q values < 20 were excluded. Short reads were assembled with De bruijn-graph-based assembler Megahit (v1.0.6) (Li et al., 2015), using a 19.7%–40.1% mapped ratio across all the samples. The open reading frames (ORFs) of each metagenomic sample were estimated with Prodigal (v2.60) (Hyatt et al., 2012), and those of lengths >100 bp were translated into amino acid sequences according to the NCBI translation table (https:// www.ncbi.nlm.nih.gov/Taxonomy/taxonomyhome.html/index.cgi?chapter = tgencodes#SG1). CD-HIT (http://www.bioinformatics.org/cd-hit/) was used to cluster sequences from gene sets with a 95% sequence identity (90% coverage) into a non-redundant gene catalog, which were then mapped with Bowtie 1.1.2. Next, sam2counts 0.91 was employed to convert the mapped results to reference sequence counts to generate a gene table for functional annotation.

CAZymes in the metagenome contigs were annotated following gene calling with pipeline dbCAN2 and non-redundant protein (Yin et al., 2012), according to the manual curation of sequence alignments from the CAZy database (July 2020). Non-redundant protein taxonomical annotations were determined using the easy-taxonomy mode and nr database of MMseqs2 (Steinegger and S¨oding, 2017). Next, we mapped the high-quality sequences of the samples on the predicted gene sequences to assess the abundance of these genes using Salmon (Patro et al., 2017). We then used the transcripts per kilobase million (TPM) mapped reads to normalize the abundance values in the metagenomes, and normalized the total sequencing number to one million reads per sample to avoid the sequencing depth bias in the comparative analysis.

We determine soil total nitrogen (TN), total carbon (TC), total potassium, total phosphorus, extractable phosphorus, pH, water content, NH₄ ⁺-N, NO₃ ⁻-N, organic carbon, microbial biomass carbon (MBC), microbial biomass nitrogen (MBN), and microbial biomass phosphorus (MBP). Soil TC and TN were determined by performing elemental analyses (Elementar Vario MACRO, Germany). An automated discrete analysis was used to quantify ammonium nitrogen (NH₄ ⁺-N) and nitrate nitrogen (NO₃ ⁻-N) (CleverChem 380, Germany). Total K was determined by flame photometry (Lu, 1999). Total and extractable P concentrations were determined using a spectrophotometer (UV-1600 spectrophotometer, Beijing) (Lu, 1999). The standard Walkley-Black potassium dichromate oxidation method (Nelson and Sommers, 1982) was used for organic carbon. A 1:1 ratio of soil to water was used to measure soil pH (HANNA, Padova, Italy). Additionally, we gravimetrically measured soil water content (SWC) by oven drying to constant mass at 105°C, and MBC, MBN, and MBP were determined using the chloroform fumigation extraction method.

The effect of grassland degradation severity and soil depth on CAZyme gene family read abundances were determined using two-way analyses of variance (ANOVAs). If ANOVA indicated significant effects, then the means of treatments were compared using Tukey’s honestly significant difference test (p < 0.05). Two-way ANOVA and multiple comparisons were also performed on each of the glycoside hydrolase (GH), glycosyltransferase (GT), polysaccharide lyases (PL), carbohydrate esterase (CE), carbohydrate-binding module (CBM), and auxiliary activity (AA) families.

Non-metric multidimensional scaling (NMDS) analysis was performed using Bray-Curtis distances to assess variation in the CAZyme gene community structure. Significant differences in CAZyme gene communities among different degradation severities and soil depths were determined using the PerMANOVA test (Vegan package in R, Oksanen et al., 2013).

After confirming that the structure of a CAZyme gene family was significantly affected by treatment, indicator species analyses were conducted to investigate genes with significantly different abundances using the ‘signassoc’ function in the ‘indicspecies’ package (DeCáceres and Legendre 2009) according to each gene’s abundance in a sample. The mode zero (site-based) was used and p-values were obtained following Sidak’s correction for multiple testing.

We assessed the association of CAZyme gene composition with plant community composition or soil characteristics using the Mantel tests. Partial mantel tests were performed between CAZyme gene composition and plant community composition or soil characteristics after controlling for other factors.

Variance partitioning of Bray-Curtis distances of CAZyme gene compositions was performed using a multiple regression for distance matrices (Swenson, 2014) to assess the significance of plant community compositions and soil characteristics in determining CAZyme gene composition.

We built structural equation models (SEM) with Mantel R values as inputs in AMOS 20.0 (Arbuckle, 2011) to determine the direct and indirect effects of plant communities on soil characteristics of the CAZyme gene families. Basing on the Mantel test and variance partitioning, which showed a significant interactive effect between plant communities and soil characteristics, a conceptual model was built whereby plant communities affected soil characteristics and both types of variables jointly affected CAZyme gene composition. SEM models with observations were compared using a maximum likelihood estimation method. Based on the methodology of Hu and Bentler (1999), we assessed the goodness of fit using a root mean squared error of approximation (RMSEA) <0.05 and TLI (Tucker-Lewis Index) >0.95.

There were a total of 1,676,180 CAZyme gene sequences (minimum 67,091; maximum 72,168; mean 69,840 per sample). GT and GH were identified as the two dominant families of genes, representing an average of 32.29% and 25.33%, respectively.

Degradation severity and soil depth, and their interaction, were significant with respect to CAZyme gene sequence abundance (p < 0.05; Figure 1; Supplementary Table S2). CAZyme gene abundance was significantly different between the two soil depths (p < 0.05). For the lightly and moderately degraded grassland, the abundance of the CAZyme genes were significantly higher in the 10–20 cm depth soil than in 0–10 cm soil depth. For the non- and severely degraded grassland, CAZyme gene abundances were not different between the two soil layers. Among the different degradation severities, CAZyme gene abundance significantly differed between non-degraded and severely degraded, between lightly and severely degraded, and between moderately and severely degraded, but not among non-degraded, lightly degraded and moderately degraded grassland (Figure 1; Supplementary Table S3).

The NMDS together with the two-way PerMANOVA (p < 0.05 for degradation severity, soil depth, and their interaction) indicated a significant difference in the structure of CAZyme gene families across degradation severity and soil depth (Figure 2; Table 1).

Indicator species analysis identified 201 gene families that differed significantly in abundance across treatments (Figure 3). Of these, 100 gene families were significantly higher in severely degradation grassland, of which 52 had a soil depth exceeding 10–20 cm and 48 exceeding 0–10 cm (Figure 3). Further, 49 out of the 100 gene families enriched in severely degraded sites were GHs, of which 20 had a soil depth exceeding 10–20 cm, and 29 exceeded 0–10 cm (Figure 3).

Mantel tests showed a significant relationship between the structure of CAZyme gene families and plant communities (p < 0.05; Table 2) and soil characteristics (p < 0.05; Table 2), at both of the soil depths. In the partial Mantel test, after controlling for soil characteristics, plant communities showed a significant relationship with structure of CAZyme gene families (p < 0.05; Table 2). However, soil characteristics were not significantly related to CAZymes after controlling for plant community structure (p > 0.05; Table 2).

Using the distance matrix multiple regression, plant communities and soil characteristics variance partitioning was used to examine their relative influence on CAZyme gene structure. These two factors together explained 30% and 22% for CAZyme composition at soil depths of 0–10 cm and 10–20 cm, respectively (Figure 4). The interpretation strength is mainly due to the interaction of plant communities and soil characteristics, with 23% and 18% for CAZyme structure at soil depths of 0–10 cm and 10–20 cm.

We found that plant community structure and soil characteristics explained much of CAZyme structure in the SEM model for both of the soil layers (for 0–10 cm, Chi-square = 0.375, df = 3, probability level = 0.945, RSMEA = 0.001, TLI = 0.999; for 10–20 cm, Chi-square = 0.004, df = 2, probability level = 0.998, RSMEA = 0.001, TLI = 0.999) (Figure 5). The final model explained 33% and 37% of the variation in structure of CAZyme gene families for soil depths of 0–10 cm and 10–20 cm, respectively. Total carbon, organic carbon and MBC explained significant variation in the structure of CAZyme composition. Furthermore, the SEM model also identified the interactive effect of plant community structure and soil characteristics. Plant community structure had an indirect effect on CAZyme composition by influencing soil total carbon or soil organic carbon. Soil total carbon and organic carbon also exerted an indirect effect on CAZyme composition by influencing MBC.

This study provided empirical evidence that abundance and composition structure of the soil C cycling functional genes, CAZyme gene families, varied across a gradient of grassland degradation. Some CAZyme families were overrepresented in the severely degraded grassland. The result of high CAZyme abundance and low litter input may lead to reduced C sequestration. The decrease in SOC in the severely degraded grassland suggests that unstabilized organic matter of plant origin (plant litter and root exudates) is rapidly mineralized. Moreover, a recent study reported that microbial functional changes result in an irreversible degradation of Tibetan Plateau pastures (Breidenbach et al., 2022). Thus, conventional strategies excluding grazing to mitigate severe degradation is unlikely to lead to ecosystem recovery without intervention. We urgently need to seek for ecological restoration approaches for severely degraded grassland. Reseeding may be one approach to restoring severely degraded grassland, but one must be aware that conventional N fertilization often carries the risk of N leaching and pollution (Mejias et al., 2021). For the non-degraded, lightly degraded and moderately degraded grasslands, more reasonable levels of grazing (e.g., rotational grazing; Dong et al., 2022) compared to the current grazing pattern could be explored to prevent further ecosystem degradation.

An important limitation associated with this present study was that we only performed sampling on a single occasion with only three replicate of each degradation severity. Thus, it is difficult for us to generalize our findings to other grassland systems or to our grassland at different times of year. Thus, more studies comprising larger temporal and spatial scales across various degraded and restored grassland ecosystems are needed. Furthermore, future studies might fruitfully consider microbial activity with respect to C sequestration using meta-transcriptome and metaproteome approaches in the prediction of C accumulation.