Microbial functional changes mark irreversible course of Tibetan grassland degradation

Literature study

Literature considering the effect of pasture degradation on SOC, N, and clay content, as well as bulk density (BD), was assembled by searching (i) Web of Science V.5.22.1, (ii) ScienceDirect (Elsevier B.V.) (iii) Google Scholar, and (iv) the China Knowledge Resource Integrated Database (CNKI). Search terms were “degradation gradient”, “degradation stages”, “alpine meadow”, “Tibetan Plateau”, “soil”, “soil organic carbon”, and “soil organic matter” in different combinations. The criteria for including a study in the analysis were: (i) a clear and comprehensible classification of degradation stages was presented, (ii) data on SOC, N, and/or BD were reported, (iii) a non-degraded pasture site was included as a reference to enable an effect size analysis and the calculation of SOC and N losses, (iv) sampling depths and study location were clearly presented. (v) Studies were only considered that took samples in 10 cm depth intervals, to maintain comparability to the analyses from our own study site. The degradation stages in the literature studies were regrouped into the six successive stages (S0–S5) according to the respective degradation descriptions. In total, we compiled the results of 49 publications published between 2002 and 2020.

When SOM content was presented, this was converted to SOC content using a conversion factor of 2.0³². SOC and N stocks were calculated using the following equation:

where elemental stock is SOC or N stock [kg ha⁻¹]; content is SOC or N content [g kg⁻¹]; BD is soil bulk density [g cm⁻³] and depth is the soil sampling depth [cm].

The effect sizes of individual variables (i.e., SOC and N stocks as well as BD) were quantified as follows:

where ES is the effect size in %, D is the value of the corresponding variable in the relevant degradation stage and R is the value of each variable in the non-degraded stage (reference site). When ES is positive, zero, or negative, this indicates an increase, no change, or decrease, respectively, of the parameter compared to the non-degraded stage.

Experimental design of the field study

Large areas in the study region are impacted by grassland degradation. In total, 45% of the surface area of the Kobresia pasture ecosystem on the TP is already degraded². The experiment was designed to differentiate and quantify SOC losses by erosion vs. net decomposition and identify underlying shifts in microbial community composition and link these to changes in key microbial functions in the soil C cycle. We categorized the range of Kobresia root-mat degradation from non-degraded to bare soils into six successive degradation stages (S0–S5). Stage S0 represented non-degraded root mats, while stages S1–S4 represented increasing degrees of surface cracks, and bare soil patches without root mats defined stage S5 (Supplementary Fig. 1). All six degradation stages were selected within an area of about 4 ha to ensure equal environmental conditions and each stage was sampled in four field replicates. However, the studied degradation patterns are common for the entire Kobresia ecosystem (Supplementary Fig. 1).

Site description

The field study was conducted near Nagqu (Tibet, China) in the late summer 2013 and 2015. The study site of about 4 ha (NW: 31.274748°N, 92.108963°E; NE: 31.274995°N, 92.111482°E; SW: 31.273488°N, 92.108906°E; SE: 31.273421°N, 92.112025°E) was located on gentle slopes (2–5%) at 4,484 m a.s.l. in the core area of the Kobresia pygmaea ecosystem according to Miehe et al.⁸. The vegetation consists mainly of K. pygmaea, which covers up to 61% of the surface. Other grasses, sedges, or dwarf rosette plants (Carex ivanoviae, Carex spp., Festuca spp., Kobresia pusilla, Poa spp., Stipa purpurea, Trisetum spp.) rarely cover more than 40%. The growing season is strongly restricted by temperature and water availability. At most, it lasts from mid-May to mid-September, but varies strongly depending on the onset and duration of the summer monsoon. Mean annual precipitation is 431 mm, with roughly 80% falling as summer rains. The mean annual temperature is −1.2 °C, while the mean maximum temperature of the warmest month (July) is +9.0 °C².

A characteristic feature of Kobresia pastures is their very compact root mats, with an average thickness of 15 cm at the study site. These consist mainly of living and dead K. pygmaea roots and rhizomes, leaf bases, large amounts of plant residue, and mineral particles. Intact soil is a Stagnic Eutric Cambisol (Humic), developed on a loess layer overlying glacial sediments and containing 50% sand, 33% silt, and 17% clay in the topsoil (0–25 cm). The topsoil is free of carbonates and is of neutral pH (pH in H₂O: 6.8)⁵. Total soil depth was on average 35 cm.

The site is used as a winter pasture for yaks, sheep, and goats from January to April. Besides livestock, large numbers of plateau pikas (Ochotona) are found on the sites. These animals have a considerable impact on the plant cover through their burrowing activity, in particular the soil thrown out of their burrows, which can cover and destroy the Kobresia turf.

Sampling design

The vertical and horizontal extent of the surface cracks was measured for each plot (Supplementary Table 2). Vegetation cover was measured and the aboveground biomass was collected in the cracks (Supplementary Table 2). In general, intact Kobresia turf (S0) provided high resistance to penetration as measured by a penetrologger (Eijkelkamp Soil and Water, Giesbeek, NL) in 1 cm increments and four replicates per plot.

Soil sampling was conducted using soil pits (30 cm length × 30 cm width × 40 cm depth). Horizons were classified and then soil and roots were sampled for each horizon directly below the cracks. Bulk density and root biomass were determined in undisturbed soil samples, using soil cores (10 cm height and 10 cm diameter). Living roots were separated from dead roots and root debris by their bright color and soft texture using tweezers under magnification, and the roots were subsequently washed with distilled water to remove the remaining soil. Because over 95% of the roots occurred in the upper 25 cm⁵, we did not sample for root biomass below this depth.

Additional soil samples were taken from each horizon for further analysis. Microbial community and functional characterization were performed on samples from the same pits but with a fixed depth classification (0–5 cm, 5–15 cm, 15–35 cm) to reduce the number of samples.

Plant and soil analyses

Soil and roots were separated by sieving (2 mm) and the roots subsequently washed with distilled water. Bulk density and root density were determined by dividing the dry soil mass (dried at 105 °C for 24 h) and the dry root biomass (60 °C) by the volume of the sampling core. To reflect the root biomass, root density was expressed per soil volume (mg cm⁻³). Soil and root samples were milled for subsequent analysis.

Elemental concentrations and SOC characteristics

Total SOC and total N contents and stable isotope signatures (δ¹³C and δ¹⁵N) were analyzed using an isotope ratio mass spectrometer (Delta plus, Conflo III, Thermo Electron Cooperation, Bremen, Germany) coupled to an elemental analyzer (NA 1500, Fisons Instruments, Milano, Italy). Measurements were conducted at the Centre for Stable Isotope Research and Analysis (KOSI) of the University of Göttingen. The δ¹³C and δ¹⁵N values were calculated by relating the isotope ratio of each sample (R_sample = ¹³C/¹²C or ¹⁵N/¹⁴N) to the international standards (Pee Dee Belemnite ¹³C/¹²C ratio for δ¹³C; the atmospheric ¹⁵N/¹⁴N composition for δ¹⁵N).

Soil pH of air-dried soil was measured potentiometrically at a ratio (v/v) of 1.0:2.5 in distilled water.

Lignin phenols were depolymerized using the CuO oxidation method²⁵ and analyzed with a gas chromatography-mass spectrometry (GC–MS) system (GC 7820 A, MS 5977B, Agilent Technologies, Waldbronn, Germany). Vanillyl and syringyl units were calculated from the corresponding aldehydes, ketones, and carboxylic acids. Cinnamyl units were derived from the sum of p-coumaric acid and ferulic acid. The sum of the three structural units (VSC = V + S + C) was considered to reflect the lignin phenol content in a sample.

DNA extraction and PCR

Samples were directly frozen on site at −20 °C and transported to Germany for analysis of microbial community structure. Total DNA was extracted from the soil samples with the PowerSoil DNA isolation kit (MoBio Laboratories Inc., Carlsbad, CA, USA) according to the manufacturer’s instructions, and DNA concentration was determined using a NanoDrop 1000 spectrophotometer (Thermo Fisher Scientific, Wilmington, DE, USA). The extracted DNA was amplified with forward and reverse primer sets suitable for either t-RFLP (fluorescence marked, FAM) or Illumina MiSeq sequencing (Illumina Inc., San Diego, USA): V3 (5’-CCT ACG GGN GGC WGC AG-3’) and V4 (5’-GAC TAC HVG GGT ATC TAA TCC-3’) primers were used for bacterial 16 S rRNA genes whereas ITS1 (5’-CTT GGT CAT TTA GAG GAA GTA A-3’), ITS1-F_KYO1 (5’-CTH GGT CAT TTA GAG GAA STA A-3’), ITS2 (5’-GCT GCG TTC TTC ATC GAT GC-3’) and ITS4 (5’-TCC TCC GCT TAT TGA TAT GC-3’) were used for fungi^33,34. Primers for Illumina MiSeq sequencing included adaptor sequences (forward: 5’-TCG TCG GCA GCG TCA GAT GTG TAT AAG AGA CAG-3’; reverse: 5’-GTC TCG TGG GCT CGG AGA TGT GTA TAA GAG ACA G-3’)³³. PCR was performed with the Phusion High-Fidelity PCR kit (New England Biolabs Inc., Ipswich, MA, USA) creating a 50 µl master mix with 28.8 µl H₂O_molec, 2.5 µl DMSO, 10 µl Phusion GC buffer, 1 µl of forward and reverse primer, 0.2 µl MgCl₂, 1 µl dNTPs, 0.5 µl Phusion HF DNA Polymerase, and 5 µl template DNA. PCR temperatures started with initial denaturation at 98 °C for 1 min, followed by denaturation (98 °C, 45 s), annealing (48/60 °C, 45 s), and extension (72 °C, 30 s). These steps were repeated 25 times, finalized again with a final extension (72 °C, 5 min), and cooling to 10 °C. Agarose gel electrophoresis was used to assess the success of the PCR and the amount of amplified DNA (0.8% gel:1.0 g Rotigarose, 5 µl Roti-Safe Gelstain, Carl Roth GmbH & Co. KG, Karlsruhe, Germany; and 100 ml 1× TAE-buffer). PCR product was purified after initial PCR and restriction digestion (t-RFLP) with either NucleoMag 96 PCR (16 S rRNA gene amplicons, Macherey-Nagel GmbH & Co. KG, Düren, Germany) or a modified clean-up protocol after Moreau (t-RFLP)³⁵: 3× the volume of the reaction solution as 100% ethanol and ¼x vol. 125 mM EDTA was added and mixed by inversion or vortex. After incubation at room temperature for 15 min, the product was centrifuged at 25,000 × g for 30 min at 4 °C. Afterwards the supernatant was removed, and the inverted 96-well plate was centrifuged shortly for 2 min. Seventy microliters ethanol (70%) were added and centrifuged at 25,000 × g for 30 min at 4 °C. Again, the supernatant was removed, and the pallet was dried at room temperature for 30 min. Finally, the ethanol-free pallet was resuspended in H₂O_molec.

T-RFLP fingerprinting

The purified fluorescence-labeled PCR products were digested with three different restriction enzymes (MspI and BstUI, HaeIII) according to the manufacturer’s guidelines (New England Biolabs Inc., Ipswich, MA, USA) with a 20 µl master mix: 16.75 µl H₂O_molec, 2 µl CutSmart buffer, 0.25 or 0.5 µl restriction enzyme, and 1 µl PCR product for 15 min at 37 °C (MspI) and 60 °C (BstUI, HaeIII), respectively. The digested PCR product was purified a second time³⁵, dissolved in Super-DI Formamide (MCLAB, San Francisco, CA, USA) and, along with Red DNA size standard (MCLAB, San Francisco, USA), analyzed in an ABI Prism 3130 Genetic Analyzer (Applied Biosystems, Carlsbad, CA, USA). Terminal restriction fragments shorter than 50 bp and longer than 800 bp were removed from the t-RFLP fingerprints.

16 S rRNA gene and internal transcribed spacer (ITS) sequencing and sequence processing

The 16 S rRNA gene and ITS paired-end raw reads for the bacterial and fungal community analyses were deposited in the National Center for Biotechnology Information (NCBI) Sequence Read Archive (SRA) and can be found under the BioProject accession number PRJNA626504. This BioProject contains 70 samples and 139 SRA experiments (SRR11570615–SRR11570753) which were processed using CASAVA software (Illumina, San Diego, CA, USA) for demultiplexing of MiSeq raw sequences (2 × 300 bp, MiSeq Reagent Kit v3).

Paired-end sequences were quality-filtered with fastp (version 0.19.4)³⁶ using default settings with the addition of an increased per base phred score of 20, base-pair corrections by overlap (-c), as well as 5′- and 3′-end read trimming with a sliding window of 4, a mean quality of 20 and minimum sequence size of 50 bp. Paired-end sequences were merged using PEAR v0.9.11³⁷ with default parameters. Subsequently, unclipped reverse and forward primer sequences were removed with cutadapt v1.18³⁸ with default settings. Sequences were then processed using VSEARCH (v2.9.1)³⁹. This included sorting and size-filtering (—sortbylength,—minseqlength) of the paired reads to ≥300 bp for bacteria and ≥140 bp for ITS1, dereplication (—derep_fulllength). Dereplicated sequences were denoised with UNOISE3⁴⁰ using default settings (—cluster_unoise—minsize 8) and chimeras were removed (—uchime3_denovo). An additional reference-based chimera removal was performed (—uchime_ref) against the SILVA⁴¹ SSU NR database (v132) and UNITE⁴² database (v7.2) resulting in the final set of amplicon sequence variants (ASVs)⁴³. Quality-filtered and merged reads were mapped to ASVs (—usearch_global–id 0.97). Classification of ASVs was performed with BLAST 2.7.1+ against the SILVA SSU NR (v132) and UNITE (v7.2) database with an identity of at least 90%. The ITS sequences contained unidentified fungal ASVs after UNITE classification, these sequences were checked (blastn)⁴⁴ against the “nt” database (Nov 2018) to remove non-fungal ASVs and only as fungi classified reads were kept. Sample comparisons were performed at the same surveying effort, utilizing the lowest number of sequences by random selection (total 15,800 bacteria, 20,500 fungi). Species richness, alpha and beta diversity estimates, and rarefaction curves were determined using the QIIME 1.9.1⁴⁵ script alpha_rarefaction.py.

The final ASV tables were used to compute heatmaps showing the effect of degradation on the community using R (Version 3.6.1, R Foundation for Statistical Computing, Vienna, Austria) and R packages “gplots”, “vegan”, “permute” and “RColorBrewer”. Fungal community functions were obtained from the FunGuild database⁴⁶. Plant mycorrhizal association types were compiled from the literature^{38,39,40,41,47,48,49,50}. If no direct species match was available, the mycorrhizal association was assumed to remain constant within the same genus.

Enzyme activity

Enzyme activity was measured to characterize the functional activity of the soil microorganisms. The following extracellular enzymes, involved in C, N, and P transformations, were considered: two hydrolases (β-glucosidase and xylanase), phenoloxidase, urease, and alkaline phosphatase. Enzyme activities were measured directly at the sampling site according to protocols after Schinner et al.⁵¹. Beta-glucosidase was incubated with saligenin for 3 h at 37 °C, xylanase with glucose for 24 h at 50 °C, phenoloxidase with L-3,4-dihydroxy phenylalanine (DOPA) for 1 h at 25 °C, urease with urea for 2 h at 37 °C and alkaline phosphatase on P-nitrophenyl phosphate for 1 h at 37 °C. Reaction products were measured photometrically at recommended wavelengths (578, 690, 475, 660, and 400 nm, respectively).

SOC stocks and SOC loss

The SOC stocks (in kg C m⁻²) for the upper 30 cm were determined by multiplying the SOC content (g C kg⁻¹) by the BD (g cm⁻³) and the thickness of the soil horizons (m). SOC losses (%) were calculated for each degradation stage and horizon and were related to the mean C stock of the reference stage (S0). The erosion-induced SOC loss of the upper horizon was estimated by considering the topsoil removal (extent of vertical soil cracks) of all degraded soil profiles (S1–S5) and the SOC content and BD of the reference (S0). To calculate the mineralization-derived SOC loss, we accounted for the effects of SOC and root mineralization on both SOC content and BD. Thus, we used the SOC content and BD from each degradation stage (S1–S5) and multiplied it by the mean thickness of each horizon (down to 30 cm) from the reference site (S0). The disentanglement of erosion-derived SOC loss from mineralization-derived SOC loss was based on explicit assumptions that (i) erosion-derived SOC losses are mainly associated with losses from the topsoil, and (ii) the decreasing SOC contents in the erosion-unaffected horizons were mainly driven by mineralization and decreasing root C input.

Statistical analyses

Statistical analyses were performed using PASW Statistics (IBM SPSS Statistics) and R software (Version 3.6.1). Soil and plant characteristics are presented as means and standard errors (means ± SE). The significance of treatment effects (S0–S5) and depth was tested by one-way ANOVA at p < 0.05. Prior to this, we checked for normality and homogeneity of variance using the Shapiro–Wilk test and Levene’s test, respectively.

Post-hoc multiple comparisons were carried out using the LSD or Tukey HSD ANOVA, if normality was indicated. In cases of non-normal distribution, the nonparametric Kruskal–Wallis test was implemented coupled with a Bonferroni correction. To detect relationships between various plant and soil characteristics, we used linear and nonlinear regressions. Correlations were deemed significant for single regressions at p < 0.05.

Before testing for significant differences, three outliers were detected by Grubbs outlier test (p < 0.05) in the bacterial t-RFLP and MiSeq datasets: MspI: S0_B_0-5, S1_E_0-5; BstuI: S5_A_15-35; MiSeq: S3_E_5-15 and were excluded from the analysis. T-RFLP and MiSeq data of all degradation stages were compared for significant differences with MANOVA, based on the Bray–Curtis index for dissimilarity. For pairwise multilevel comparisons, “pairwiseAdonis” was used⁵². Differences in microbial community data from t-RFLP and MiSeq were displayed in non-metric multidimensional scaling (NMDS) plots and environmental factors were correlated by canonical correspondence analysis (CCA). Statistics on community data were carried out with R statistical software (Version 3.6.1).

Reporting summary

Further information on research design is available in the Nature Research Reporting Summary linked to this article.

Source: Ecology - nature.com