Study sites, experimental design, and field sampling
The Tibetan alpine permafrost region, the largest area of permafrost in the middle and low latitudes of the Northern Hemisphere64, stores substantial soil C (15.3–46.2 Pg C within 3 m depth)65,66,67. With continuous climate warming, permafrost thaw has triggered the formation of widespread thermokarst landscapes across this permafrost area13,68. To explore the impacts of thermokarst formation and development on soil C dynamics, we collected topsoil samples (0–15 cm) from a thaw sequence in 2014 and from five additional sites spread across the region in 2020. The thermokarst landscape was characterized as thermo-erosion gullies (Supplementary Table 2). The elevation of these six sites is between 3515 and 4707 m. The mean annual temperature across this area ranges from −3.1 to 2.6 °C, and the average annual precipitation varies from 353 to 436 mm. The vegetation type across these sites is swamp meadow, with the dominant species being Kobresia tibetica, Kobresia royleana and Carex atrofuscoides. Although the dominant species did not change after permafrost collapse, the forb coverage increased along the thaw sequence and across the five additional thermokarst-impacted sites. The main soil type is Cryosols on the basis of the World Reference Base for Soil Resources69, with soil pH ranging from 5.6 to 7.3 (Supplementary Fig. 1e). The active layer thickness varies between 0.7 and 1.1 m across the six study sites and the underlying soil parent material is either siliciclastic sedimentary or unconsolidated sediments (Supplementary Table 2).
To evaluate the dynamics of soil C fractions after permafrost collapse, we collected soil samples across the Tibetan alpine permafrost region based on the following two steps (Supplementary Fig. 5). In the first step, we established six collapsed plots (~15 × 10 m) along a thaw sequence (located in Shaliuhe close to Qinghai Lake, Qinghai Province, China), which had been collapsed for 1, 3, 7, 10, 13, and 16 years before the sampling year of 201413. The collapse time of each plot was estimated by dividing the distance between the collapsed plot and the gully head by the retreat rate (~8.0 m year−1; the rate of the head-wall retreat was determined by Google Earth satellite images and in situ monitoring)13. Then, we set up six paired control (non-collapsed) plots adjacent to these collapsed plots. To limit experimental costs, we selected three paired control and collapsed plots (collapsed for 1, 10, and 16 years, representing the early, middle, and late stages of collapse) to examine the responses of POC, MAOC and OC-Fe to permafrost collapse (Supplementary Fig. 5). Within each collapsed plot, we collected topsoil (0–15 cm) samples from all vegetated patches (Supplementary Fig. 6), and then evenly selected 10 vegetated patches for this study considering the heavy workload and high cost. In each selected vegetated patch, 5–8 soil cores were sampled and completely mixed as one replicate. Within each control plot, topsoil samples were randomly collected from five quadrats at the center and four corners of the plot. In each quadrat, 15–20 soil cores were sampled and mixed as one replicate. Thereby, ten replicates were acquired in each collapsed plot (n = 10), and five replicates were obtained in each control plot (n = 5). In total, we acquired 45 soil samples, including 30 samples from the three collapsed plots and 15 samples from the non-collapsed control for subsequent analysis.
In the second step, to further verify the universality of collapse effects on SOC fractions, we collected topsoil (0–15 cm) samples from an additional five similar sites located near the towns of Ebo, Mole, Huashixia, and Huanghe across a 550 km permafrost transect in August 2020 (Fig. 1). Specifically, paired collapsed and control plots (15 × 10 m) were established at the end of a gully and in adjacent non-collapsed areas in each site (Supplementary Fig. 5). In the collapsed plot, we set five 5 × 3 m quadrats at the center and four corners of the plot, and then collected topsoil samples within all the vegetated patches in these quadrats. In each quadrat, all the collected soil cores (15–20 cores) were completely mixed as one replicate, and finally, five replicates were acquired in each collapsed plot (n = 5). Similarly, five replicates were obtained from the five quadrats in each control plot (n = 5). In total, we collected 50 topsoil samples across these five thermokarst-impacted sites. After transportation to the laboratory, all the soil samples were handpicked to remove surface vegetation, roots and gravels, and sieved (2 mm) for subsequent analysis.
It should be noted that the space for time approach was only used for the permafrost thaw sequence, not for the other five sites over the regional scale. Across these five sites, we focused on the impact of permafrost collapse on POC, MAOC as well as OC-Fe by comparing soil C fractions inside and outside the gully in each site rather than among the study sites. Given the low coefficient of variation of parameters (i.e., edaphic variables and soil minerals) in the control plot of each site (Supplementary Table 3), the pristine soils in each site could also be regarded as homogeneous70, and the differences in parameters inside and outside the gully could be attributed to the effects of permafrost collapse. Along the permafrost thaw sequence, to verify whether the plots with different collapse times (1, 10, and 16 years) were comparable, we analyzed a series of parameters (i.e., vegetation biomass, edaphic variables, and soil minerals) for the three control plots which were located outside the gully but adjacent to three collapsed plots within the gully (Supplementary Fig. 5). By comparing aboveground biomass, belowground biomass, SOC, soil moisture, pH, bulk density, soil texture, and soil minerals (see below for details of the analytical method), we observed that the above parameters were not significantly different among the three control plots along the thaw sequence (all P > 0.05; Supplementary Fig. 7). These comparisons demonstrated that the study area was homogeneous before permafrost thaw and thus it was reasonable to adopt the space for time approach along the permafrost thaw sequence.
It should also be noted that the collected topsoil samples used in this study were less affected by physical mixing and translocation due to thaw phenomena at the thermokarst-impacted sites. Specifically, to examine changes in soil properties upon permafrost thaw, we chose to collect topsoil within the vegetated patches rather than from the exposed soil areas in the collapsed plots (Supplementary Fig. 6). These vegetated patches (40–60 cm thickness) are formed during the landscape fragmentation after permafrost collapse13. Although permafrost collapse inevitably led to soil translocation, these vegetated patches maintained their original shapes, especially for the topsoil because it is protected by mattic epipedon in this swamp meadow ecosystem on the Tibetan Plateau (which has an intensive root network protecting soils against interference)71,72. Moreover, we collected 0–15 cm of topsoil within the vegetated patches, in which soil cores were at least 10 cm away from the edge of the patch. Due to these two points, topsoil should not be mixed with the subsoil in our case. To test this deduction, we compared the non-collapsed (control) plot with the collapsed plot occurring for 1 year (the early stage of the permafrost thaw sequence), and observed no significant differences in soil properties such as bulk density, SOC, pH, soil texture and soil minerals (all P > 0.05; Supplementary Fig. 8). These comparisons illustrated that permafrost collapse did not cause soil physical mixing for the topsoil samples involved in this study, and soil layers were comparable between the collapsed and control plots.
SOC fractionation
We separated POC and MAOC from bulk soils using a fractionation method based on a combination of density and particle size18 using the following three steps. First, 10 g of soil was put into a 100 mL centrifuge tube, and added with 50 mL of 1.6 g cm−3 NaI. After being completely mixed, the mixture was sonicated and then centrifuged at 1800 × g. The floating particulate organic matter, together with the supernatant, was poured into a GF/C filter membrane for filtration, completely washed with deionized water, and then dried at 60 °C to constant weight. Then, the C content of the particulate organic matter was determined as POC. Second, deionized water were added to the remaining soils in the tube to wash out any residual NaI. The washed soils were then separated with a 53-μm sieve. The residues on the sieve (>53 μm) were dried and determined as heavy POC. Third, the organic matter that passed through the sieve (<53 μm) was oven-dried and grinded, and the C content was determined as MAOC. The POC, heavy POC, and MAOC contents were determined using an elemental analyzer (Vario EL III, Elementar, Hanau, Germany) after these fractions were fully grinded.
We also quantified the content of OC-Fe using the citrate-bicarbonate-dithionite approach25,28. Specifically, each soil sample was evenly divided into two parts: one subsample was extracted using 0.27 M trisodium citrate, 0.11 M sodium bicarbonate, and 0.1 M sodium dithionite as the treatment group, and the other subsample was treated with 1.85 M sodium chloride and 0.11 M sodium bicarbonate as the control group28. We then measured the C content in these two subsamples with an elemental analyzer (Vario EL III, Elementar, Hanau, Germany). The difference in the C content of the soil residues between the treatment and control groups is the OC-Fe content.
Measurements of plant and soil physicochemical properties
To examine the influence of vegetation and soil properties on soil C fractions, we measured aboveground and belowground biomass, topsoil temperature and moisture, pH, bulk density, and soil texture. Briefly, both aboveground and belowground biomass was determined by the harvesting method. All aboveground vegetation within the frame quadrat (25 × 25 cm, the number of replicates being the same as for the soil samples) was cut off and dried to constant weight at 65 °C. The aboveground biomass was determined based on the dried biomass. Roots were washed free of attached soils and separated into live and dead roots according to their color and tensile strength73. Live roots were then dried at 65 °C and weighed to calculate belowground biomass (the number of replicates being the same as for the aboveground biomass).
Soil temperature in the top 15 cm was determined with a digital thermometer (DS 1922L, Wdsen Electronic Technology Co., Shanghai, China), and soil moisture was measured by drying 20 g fresh soil sample at 105 °C to constant weight. Bulk density was determined by the oven-dried soil mass divided by the container volume. Soil texture was examined using a particle size analyzer (Malvern Masterizer 2000, Malvern, Worcestershire, UK) after eliminating organic matter and carbonates by utilizing hydrogen peroxide (30%) and hydrochloric acid (3 M), respectively74,75. Soil pH was analyzed by using a pH probe (PB-10, Sartorius, Göttingen, Germany) for a soil-water mixture (soil: water = 1:5). Soil C concentration was measured with an elemental analyzer (Vario EL III, Elementar, Hanau, Germany). Given that inorganic C was not detected in soil samples with a carbonate content analyzer (Eijkelkamp 08.53, Eijkelkamp, Giesbeek, Netherlands), soil C was equal to SOC. It is worthy to note that biotic and abiotic parameters along the thaw sequence, such as aboveground and belowground biomass, soil temperature and moisture, pH, bulk density, soil texture, and SOC, were reanalyzed from published data in refs. 13, 68, while the parameters across the five additional sites at the regional scale were measured in this study.
Microbial necromass C determination
We determined soil microbial necromass C, an important source of MAOC, based on the amino sugars76. Amino sugar analyses were conducted using the following procedure77. Specifically, ~0.4 g of soil was mixed with 10 mL HCl (6 M) at 105 °C for 8 h. After cooling, the hydrolysate was added with 100 μL of myo-inositol, filtered with glass fiber filters, and dried via a rotary evaporator at 52 °C. We redissolved the residues and adjusted the pH to 6.6–6.8. The precipitates were removed after centrifuging, and the supernatant was dried. Amino sugars were dissolved with methanol and centrifuged (1000 × g for 10 min) to dislodge salts. After the addition of 100 μL N-methyl-glucamine (internal standard), the residues were then transformed to aldononitrile derivatives through heating with 300 μL of derivatization reagent (hydroxylamine hydrochloride and 4-dimethylamino pyridine mixing with pyridine and methanol [4:1; v-v]) at 80 °C for 30 min. The derivatives were then acetylated through acetic anhydride (1 mL) under 80 °C for 20 min, and after cooling to room temperature dichloromethane (1.5 mL) was added. We removed excessive derivatization reagents as completely as possible by extracting them with 1 M HCl and MilliQ water. The organic phase owning to containing amino sugar derivatives was quantified by employing a gas chromatography (Agilent 6890 A, Agilent Technologies, Palo Alto, USA) coupled with a flame ionization detector and an HP-5 capillary column (25 m × 0.25 mm × 0.25 μm).
After the examination, we obtained three types of amino sugars, including muramic acid, glucosamine, and galactosamine. Bacterial necromass C was determined based on the content of muramic acid, and fungal necromass C was calculated from the content of glucosamine and muramic acid. Specifically, bacterial necromass C = muramic acid × 45, where 45 represents the conversion ratio of muramic acid to bacterial necromass C78,79; fungal necromass C = (mmol glucosamine – 2 × mmol muramic acid) × 179.17 × 9, among which 179.17 is the glucosamine molecular weight, and 9 indicates the conversion ratio from fungal glucosamine to fungal necromass C78,79. Microbial necromass C is the sum of bacterial and fungal necromass C. For a detailed description involving conversion formulas between amino sugars and microbial necromass C, see ref. 76.
Soil mineral analysis
Given the critical role of soil minerals in protecting organic C, we measured two types of secondary minerals, including Fe oxides and phyllosilicates (known as clay minerals). For Fe oxides, we analyzed the amount of pedogenic Fe oxides, poorly crystalline Fe oxides, and organically complexed Fe oxides. Briefly, pedogenic Fe oxides were extracted using the citrate-bicarbonate-dithionite method44, and poorly crystalline Fe oxides were extracted with acid-ammonium oxalate19. Sodium pyrophosphate was used to extract organically complexed Fe oxides48. The concentration of elements in the solution was examined using an inductively coupled plasma optical emission spectrometer (iCAP 6300, Thermo Fisher Scientific, Waltham, USA).
With respect to phyllosilicates, we employed X-ray diffraction (XRD) analysis to identify clay-sized minerals (<2 µm) following the modified method described by ref. 80. (i) Organic matter and carbonates were removed by hydrogen peroxide (30%) and acetic acid (1 M), respectively. These samples were then washed with deionized water. Based on Stoke’s law, the clay fraction was isolated through sedimentation81. (ii) We prepared three oriented clay specimens for each sample by air-drying glass slide mounts coated with clay suspension. The first one was an air-dried clay specimen, the second clay slide was saturated with ethylene glycol, and the third was treated with a combination of ethylene glycol saturation and heating to 550 °C80. (iii) These three clay slides were scanned between 3 and 30° (2θ) with a 0.02° (2θ) increment using a D8 Advance X-ray diffractometer (Bruker, Karlsruhe, Germany)80,81.
Additionally, we examined the specific surface area of mineral-associated organic matter. Before the analysis, we removed organic matter by hypochlorite oxidation (five times with 1 M NaClO)82. After that, these treated samples were first degassed with helium at 325 °C for 4 h83. Then, nitrogen was dosed on the surfaces at −196 °C under a partial pressure (P/P0) ranging from 0.04 to 0.30 in a surface area and porosity analyzer (BSD-PM2, Beishide Instrument, Beijing, China). The specific surface area was calculated based on the multipoint Brunauer-Emmett-Teller approach84.
Statistical analyses
All data were tested and transformed to ensure the normality of variance, and then analyzed with the following three steps. First, we used one-way ANOVAs with Least Significant Difference (LSD) multiple comparisons to analyze the differences in SOC and its associated fractions (POC, heavy POC, MAOC, and OC-Fe) as well as biotic (aboveground and belowground biomass, microbial residual C) and abiotic (soil temperature and moisture, pH, texture, Fe oxides, clay minerals, and specific surface area) variables along the thaw sequence. In these analyses, permafrost thaw stages (non-collapse, collapsed for 1 year, 10 years, and 16 years) were treated as the between-subject effect. Based on the same statistical methods, we also analyzed the differences in soil C fractions (or biotic and abiotic variables) between collapsed and control plots at each of the regional thermokarst-impacted sites.
Second, we performed ordinary least squares regression analyses to explore the relationships of POC and OC-Fe with potential factors, including vegetation parameters (aboveground and belowground biomass), soil properties (soil temperature and moisture, pH, and texture), and Fe minerals (pedogenic Fe oxides, poorly crystalline Fe oxides, and organically complexed Fe oxides) along the thaw sequence. We then conducted a random forest model to explore the relative importance of explanatory variables in affecting POC and OC-Fe using the randomForest package.
Third, we performed linear regression analyses to assess the relationships of response ratio (RR, the ratio of two variables between collapsed and control plots [Variablecollapse/Variablecontrol]) of POC and OC-Fe with RR of explanatory variables across the thermokarst-impacted sites at the regional scale. Statistical differences were considered to be significant at the level of P < 0.05. All the statistical analyses were performed in R 4.0.4 (R Core Team, 2020).
Source: Ecology - nature.com
