Enhanced response of soil respiration to experimental warming upon thermokarst formation

As global temperatures continue to rise, a key uncertainty of terrestrial carbon (C)–climate feedback is the rate of C loss upon abrupt permafrost thaw. This type of thawing—termed thermokarst—may in turn accelerate or dampen the response of microbial degradation of soil organic matter and carbon dioxide (CO 2 ) release to climate warming. However, such impacts have not yet been explored in experimental studies. Here, by experimentally warming three thermo-erosion gullies in an upland thermokarst site combined with incubating soils from five additional thermokarst-impacted sites on the Tibetan Plateau, we investigate how warming responses of soil CO 2 release would change upon upland thermokarst formation. Our results show that warming-induced increase in soil CO 2 release is ~5.5 times higher in thermokarst features than the adjacent non-thermokarst landforms. This larger warming response is associated with the lower substrate quality and higher abundance of microbial functional genes for recalcitrant C degradation in thermokarst-affected soils. Taken together, our study provides experimental evidence that warming-associated soil CO 2 loss becomes stronger upon abrupt permafrost thaw, which could exacerbate the positive soil C–climate feedback in permafrost-affected regions.

https://doi.org/10.1038/s41561-024-01440-2heterotrophic component were separated by soil collars with different lengths (that is, the 5-cm-long and 42-cm-long collars for measuring Rs and heterotrophic respiration (Rh), respectively) 27,28 .Given that Rh measurement would create an artificial situation without root exudates in deep collars (see Methods for details), the measured Rh is not the actual Rh, but a potential Rh under our specific experimental conditions (hereafter referred to as Rhp).The two consecutive years' measurements revealed that both Rs and Rhp displayed hump-shaped trends over the growing season regardless of treatments (Fig. 1a,c and Extended Data Fig. 3a,c), with the Rhp contributing to a major portion of Rs (~67%; Supplementary Fig. 1).Considering the above-mentioned methodological limitation of Rhp and its major contribution to Rs, our following analyses concentrated on Rs while the Rhp data were only used to supplement the Rs results.We showed that, in unwarmed plots, the Rs significantly increased by an average of 55.4 ± 8.0% (mean ± s.e.m.) upon permafrost collapse (Fig. 1).The higher Rs upon thermokarst formation may be explained by the following reasons: (1) after permafrost collapse, the increase in soil permeability and thinning of the diffusive boundary layer 29 would open up previously obstructed vertical and lateral flow paths and thereby enhance transportation of organic matter, which favours microbial CO 2 emission 29,30 ; (2) the above environmental changes and decrease in volumetric soil moisture (from 57.1 ± 0.4% to 51.5 ± 0.7%; Extended Data Fig. 4c,d) would also allow more oxygen to fill the voids of the organic layer, thus exacerbating the SOM decomposition 29,30 ; and (3) permafrost collapse would disrupt soil aggregates 31 , weakening SOM stability and then accelerating SOM decomposition 32 .Overall, these changes would lead to the higher Rs upon permafrost collapse observed in this study.
about the responses of the ecosystem C cycle to abrupt permafrost thaw, termed thermokarst 8 .Thermokarst occurs owing to the melting of ice-rich permafrost, which can alter the trajectory of the permafrost C cycle 9 .Given that ongoing climate warming would affect many ecosystem processes in both thermokarst and non-thermokarst areas, an interesting question is whether warming effects on soil CO 2 flux could differ between these two distinct landforms.If warming effects are amplified upon thermokarst formation, the positive permafrost C-climate feedback will be intensified.It is thus imperative to test whether the warming response of soil CO 2 release will become stronger where thermokarst landscape occurs.
Among the generally identified thermokarst landforms (lake, wetland and upland thermokarst), upland thermokarst (estimate of 0.91 × 10 6 km 2 ) often occurs on moderate slopes or along water courses 8,9 .Upland thermokarst formation can restructure land surface morphology, causing abrupt changes in soil environments, organic matter composition and microbial properties 10,11 .For example, observational studies demonstrate a substantial decline in soil organic C quality 12 and soil microbial diversity 11,13 following upland permafrost collapse.These changes could not only alter the rate of soil-respired CO 2 (refs.8,14), but also affect its response to the ongoing warming.On the one hand, compared with the high-quality soil substrates, low-quality soil substrates in the collapsed area require higher activation energy to be decomposed 15,16 , and their decomposition is thus more sensitive to elevated temperature.On the other hand, thermokarst-induced decrease in microbial diversity 11 probably exerts dual effects on the temperature response of SOM decomposition: the reduction of microbial diversity may lead to a decline in warming effect because of the reduced microbial metabolic potential 17 , but it may also elevate the warming effect because higher biodiversity is associated with greater functional resistance 18,19 , and thus the declined microbial diversity would render microbial mediated functions (for example, SOM decomposition) more sensitive to environmental change, such as climate warming.However, current in situ warming experiments have all been conducted in non-thermokarst landscapes [20][21][22] .Whether and how warming effects on soil CO 2 fluxes would change upon thermokarst formation, to date, has never been explicitly explored in any experimental studies in permafrost areas.
Here we identified three thermo-erosion gullies at a large thermokarst feature in Shaliuhe on the northeastern Tibetan Plateau (37° 28′ N, 100° 17′ E, ~3,900 m above sea level; Extended Data Fig. 1a), which has the largest distribution of alpine permafrost on Earth and stores a large amount of C (~15. 3 -46.2Pg C at a depth of 3 m; refs.23,24).In each gully, we conducted in situ warming experiments inside (collapsed area) and outside the gully (non-collapsed area) using open-top chambers (OTCs; Extended Data Fig. 2).Based on these well-replicated warming experiments, we compared warming effects on soil CO 2 flux (hereafter soil respiration; Rs) between collapsed and non-collapsed areas.To verify the broader representativeness of single-site findings, we collected soils from the Shaliuhe site and five additional upland thermokarst-impacted sites (Extended Data Fig. 1 and Supplementary Table 1) to examine how temperature sensitivity of CO 2 release (the factor by which respiration rises as a result of a 10 °C increase, Q 10 ) responds to permafrost collapse using a laboratory incubation experiment.Further, based on in situ warming experiments, combining soil physiochemical analyses, solid-state 13 C nuclear magnetic resonance 25 and metagenomic sequencing 26 , we explored the associations of Rs response with various biotic and abiotic variables to identify the key factors that were most strongly correlated with the differences in warming effects.

Intensified Rs response to warming in thermokarst area
In the field experiments, we simultaneously measured Rs and its heterotrophic component, soil temperature and moisture three times per month during the growing season in 2019 and 2020.Rs and its The OTC deployment increased soil temperature in the upper 10 cm by an average of 1.3 °C (from 6.0 to 7.3 °C) across the thermokarst landforms (non-collapsed and collapsed areas) and growing seasons (2019 and 2020; Extended Data Fig. 4a,b).However, volumetric soil moisture at various depths was insensitive to experimental warming (control versus warming: 55.4 versus 57.8%, 61.3 versus 63.4% and 56.1 versus 51.4% across non-collapsed and collapsed areas at 10, 20 and 30 cm depth, respectively; Extended Data Fig. 5d-f).Warming treatment substantially increased Rs by an average of 45.6 ± 1.8%, which should mostly result from the top 10 cm soil due to the limited warming effects in deeper soils (Extended Data Fig. 5a-c).By taking advantage of the simultaneously measured Rs in warmed and control plots between the collapsed and non-collapsed areas, we explored whether and how warming effects would change upon thermokarst formation.To this end, we first computed the normalized Rs response because the rising magnitude of topsoil (0-10 cm) temperature by OTC differed between the two land features, as indicated by the significant interaction between permafrost collapse and experimental warming (Extended Data Fig. 4a,b).The normalized warming effects were calculated using the following formula 33 : where Rs Warming and Rs Control are Rs (μmol m −2 s −1 ) in warmed and control plots, respectively, and T Warming and T Control are the topsoil temperatures (°C) in warmed and control plots, respectively.After normalization, we found that the Rs response was much stronger upon thermokarst formation (increased 4.2-and 6.7-fold in 2019 and 2020, respectively, both P < 0.001; Fig. 1b,d).To further validate these results, we estimated the apparent Q 10 of Rs by fitting Rs with the corresponding topsoil temperature in unwarmed plots using an exponential function 27 (see Supplementary Note 1 for detailed calculations).Results showed that the apparent Q 10 was significantly higher upon permafrost collapse (Extended Data Fig. 6a,b), supporting the results from the in situ warming experiments.Notably, similar patterns of the normalized warming effects (Extended Data Fig. 3) and the apparent Q 10 (Extended Data Fig. 6c,d) were observed for Rhp, providing additional evidence for a stronger soil CO 2 response to experimental warming in thermokarst landscapes.
To test the universality of the single-site findings, we collected soil samples from five additional thermo-erosion gullies across a 550 km permafrost transect on the Tibetan Plateau (Extended Data Fig. 1).We then conducted a laboratory incubation using soils from the three gullies at the Shaliuhe site plus five additional gullies from Ebo, Mole, Reshui, Huashixia and Huanghe sites.The Q 10 was calculated from the soil CO 2 release rate at two incubation temperatures 34 (5 °C and 15 °C; Methods).During most of the measurement occasions and throughout the incubation period, significant increases in Q 10 values were observed in collapsed relative to non-collapsed areas (ranging from 12% to 53%) at all sites except for Huashixia (Fig. 2 and Extended Data Fig. 7).Although we observed a consistent response of warming-induced soil CO 2 release to thermokarst formation both in field and laboratory experiments, the field-based Q 10 values (non-collapsed versus collapsed: 5.1 ± 0.1 versus 7.2 ± 0.1; Extended Data Fig. 6) were much higher than those from the laboratory incubation (non-collapsed versus collapsed: 2.2 ± 0.1 versus 2.7 ± 0.2; Fig. 2).Such a difference is probably induced by methodological discrepancy.For the incubation experiment, temperature is the major factor linked to changes in soil CO 2 emissions.While in the field, other factors may covary (for example, elevated substrate availability and microbial biomass) with temperature increase over the growing season 29 , which could then interact with temperature to promote Rs and yield larger apparent Q 10 values (see Supplementary Note 2 for detailed discussions).

Rs response associated with SOM quality and microbial genes
The stronger response of Rs upon thermokarst formation may be associated with many biotic and abiotic factors, such as plant, edaphic and microbial properties, and so on.To identify the key factors that were most strongly correlated with the different warming effects on Rs, we measured a number of potential factors and compared their differences between collapsed and non-collapsed areas.Specifically, we first compared plant attributes and observed that these variables remained stable upon permafrost collapse (Supplementary Fig. 2a-d), indicating that the warming responses of Rs may be minimally related to plant attributes.We then explored changes in soil characteristics and microbial properties, and found that there were significant changes in these variables (Fig. 3a-i, Extended Data Figs.8 and 9, and Supplementary Figs.2-4; see Supplementary Note 3 for detailed descriptions).Finally, based on the above comparisons together with the single variable relationships between Rs response and these factors, we constructed the structural equation modelling (SEM) to tease apart the pathways and relative importance of edaphic and microbial factors (Fig. 3j and Supplementary Figs.2q and 5).The SEM analysis revealed that the normalized Rs change was directly linked to substrate quality and the abundances of microbial functional genes (Fig. 4a), with the standardized path coefficients of 0.31 and 0.29, respectively (Fig. 4b).In addition, soil moisture and clay+silt were indirectly correlated with the normalized Rs change via substrate quality and the abundances of microbial functional genes (Fig. 4a).Similar results were also found https://doi.org/10.1038/s41561-024-01440-2 for Rhp response (Extended Data Fig. 10).These results highlighted that substrate quality and microbial functional genes were probably the key factors related to the differences in warming effects between collapsed and non-collapsed areas.
While the statistical correlations may not infer the cause-effect relationships, the results allow for a causal speculation on these associations by prior ecological principles.We observed a positive correlation between normalized Rs change and substrate quality (Fig. 3j).In our study site, the recalcitrant C components were found to increase upon thermokarst formation (Fig. 3b).In such a case, we may expect that, based on the 'C quality-temperature' hypothesis 15,35 , the SOM decomposition is more responsive to experimental warming in collapsed than that in non-collapsed areas because the decomposition of recalcitrant substrates requires higher activation energy, which is therefore more sensitive to warming 35 .Then an interesting question is why there was an increase in the recalcitrant C components.Such a change may be explained by the preferential usage of the labile C by microbial decomposition and hydrologic transport/lateral flux 12,36 .On the one hand, the easily degradable compounds can be rapidly accessed by soil microorganisms after permafrost collapse 37,38 , while the more recalcitrant alkyl C would accumulate and eventually become more abundant within soils of collapsed vegetation patches 12 .On the other hand, the weakened mineral absorption (as indicated by reduced clay+silt content) and increased hydrological transport in collapsed soils (as reflected by reduced soil moisture; Fig. 4a) may lead to the leaching of labile SOM 39 , while the recalcitrant compounds are more readily adsorbed to the leftover mineral soils 40 , thus reducing the SOM quality.
Although there were no obvious relationships between soil CO 2 flux response and bacterial diversity changes (Extended Data Fig. 9a,b), the normalized Rs change was positively correlated with abundances of microbial functional genes involved in SOM decomposition (Fig. 4), particularly those involved in recalcitrant substrate decomposition (that is, those C degradation genes for pectin and aromatic compounds; Fig. 3j).Increases in microbial functional genes were associated with reduced soil moisture and clay+silt content upon permafrost collapse (Fig. 4a), which allowed more oxygen as for a greater microbial functioning.The increased abundance of these functional genes after thermokarst formation implies an elevated potential for the decomposition of recalcitrant SOM 11,41 , which may enhance the Rs response at higher temperatures based on the 'C quality-temperature' hypothesis 15,35 .We also found a pronounced increase in the normalized response of oxidase activity in the collapsed soils (Extended Data Fig. 8), implying an enhanced degradation of recalcitrant SOM and further supporting the higher warming response of Rs upon thermokarst formation.
Apart from the above pathways supported by our statistical results, there may be other processes underlying the different warming effects that were not measured or excluded from the SEM analysis, such as physiochemical protection and other environmental changes.First, the larger Rs response may be associated with destabilized aggregate protection.Permafrost collapse as a disturbance together with the decreased clay+silt (Supplementary Fig. 2i) is supposed to weaken physical occlusion within aggregates 31 and thus increase substrate accessibility, which would then alleviate constraints to https://doi.org/10.1038/s41561-024-01440-2microbial decomposition and further favour larger warming response of Rs.Despite the vital role of reactive minerals in attenuating the temperature response of microbial decomposition 35 , it might be less important in our study because organically complexed Fe increased upon permafrost collapse (Supplementary Fig. 2l), which would attenuate Rs response from the theoretical perspective but contrasted from the stronger response detected here.Second, the larger Rs response may also be related to interaction effects of temperature with other environmental changes.Following permafrost collapse, decrease in soil moisture (Extended Data Fig. 4c,d) and associated increase in soil permeability 42 and thinning of the diffusive boundary layer 43 would allow as much as substrates and oxygen for microbial assess.The enhanced supplies of substrates and oxygen would then amplify the temperature effect on soil CO 2 flux, favouring stronger normalized Rs response.In support of this argument, the observed difference of Q 10 between collapsed and non-collapsed areas was smaller in laboratory incubation compared with that in the field (0.5 ± 0.2 versus 2.1 ± 0.2; Fig. 2 and Extended Data Fig. 6a,b), indicating that the interaction effects would amplify the Rs response detected in situ.

Implications for soil CO 2 release under climate warming
By experimentally warming three thermo-erosion gullies in an upland thermokarst site combined with incubating soils from another five thermokarst-impacted sites on the Tibetan Plateau, this study provided multiple lines of evidence that the warming-associated soil CO 2 loss became stronger under upland thermokarst formation.In such a case, thermokarst-affected soils were estimated to respire an additional 422 ± 51 g C m −2 yr −1 compared with the non-thermokarst landforms under a 2 °C warming (the warming magnitude in permafrost regions by the end of the twenty-first century under a moderate scenario 44 ).As a preliminary exploration of global importance, we extrapolated the warming response of Rs to all upland thermokarst regions in the Northern Hemisphere (0.91 × 10 6 km 2 ) 9 and found 0.4 ± 0.05 Pg C yr −1 of the additional soil CO 2 loss (see Supplementary Note 4 and Supplementary Fig. 6 for detailed calculations).This warming-induced additional CO 2 release upon thermokarst formation corresponds to roughly a quarter of the projected permafrost soil CO 2 losses (2 Pg C yr −1 ) by the end of the twenty-first century 45 .Our estimate indeed has uncertainties owing to methodological limitations.For instance, the OTC warming method would cause artefact chamber effects that are not associated with warming (for example, local wind shielding, rainfall interception and so on) 46 , which may lead to an overstatement of the warming effects and thus bias of the estimated additional CO 2 emission (see Supplementary Note 4 for details); the spatial heterogeneity of climate, vegetation and soil properties may also influence the Rs responses and lead to uncertainties of the global extrapolation.Despite this, our findings highlight that the intensified warming response of soil CO 2 efflux upon thermokarst formation cannot be neglected when projecting the future trajectory of permafrost C-climate feedback.More broadly, our findings may be extrapolated to other disturbed versus non-disturbed sites, where disturbances such as tillage 47 and fire 48 led to the loss of labile compounds, while soil erosion 49 , deforestation and land use change 50 altered community composition.Consequently, the strength of soil CO 2 release may be underestimated if the permafrost collapse and other analogous disturbances are ignored.Extending studies that focus on warming effects simultaneously in disturbed versus non-disturbed sites are thus required to be performed in diverse regions and ecosystems for a more accurate estimation of the soil-respired CO 2 under future warming.

Site description
Climate warming has intensified thermokarst formation on the Tibetan Plateau 51 .The area affected by thermokarst formation has increased 40-fold during the past five decades 52 .To explore how warming responses of soil CO 2 release will change upon upland thermokarst formation, we established an in situ warming experiment at an upland thermokarst site and incubated soils from five additional thermokarst-impacted sites across a permafrost transect of 550 km (Extended Data Fig. 1).The elevation of the six study sites ranged from 3,515 to 4,707 m above sea level.The mean annual temperature and precipitation varied from −3.1 to 2.6 °C and from 353 to 436 mm, respectively.All sites were characterized as swamp meadows dominated by sedges including Kobresia tibetica, Kobresia royleana and Carex atrofuscoides.The soils were classified as cryosols based on the World Reference Base for Soil Resources 53 .All sites were located in discontinuous permafrost zones, with the active layer thickness being ~0.7-1.1 m.There was no obvious sign of cryoturbation in the six selected sites.See Supplementary Table 1 for detailed information.
Among the six study sites, we selected three thermo-erosion gullies at the Shaliuhe site (Extended Data Fig. 2) and set up the warming experiments in collapsed and non-collapsed areas for each thermo-erosion gully.In the field warming experiments, we measured a series of biotic and abiotic factors such as vegetation and microbial properties, soil physicochemical properties, and substrate quality.To evaluate the generalizability of the response pattern obtained from the field warming experiments, we conducted laboratory incubations using soil samples collected from the Shaliuhe and five additional thermokarst-impacted sites.The specific methods are described below.

In situ warming experiment
Experimental design.The field warming experiment was established during 2017-2018 at the Shaliuhe site (37° 28′ N, 100° 17′ E).The location is ~30 km away from Gangca County, Qinghai Province, China (Extended Data Fig. 1a).The site has a large SOC density (C amount per area; Supplementary Fig. 7), which is much higher than the average on the Tibetan Plateau 54 .We identified three thermo-erosion gullies at a large thermokarst feature on a hillslope with a gentle south-facing incline of 9° (ref.55).The three gullies were approximately 100 m from each other.A 20 × 20 m 2 fenced area was set in a relatively flat area at the end of each gully.Meanwhile, an area of the same size was set in the adjacent non-collapsed grassland (Extended Data Fig. 2).Ten 2 × 2 m 2 blocks were randomly arranged in each enclosure, and two plots within each block were randomly assigned to control and warming treatments.An OTC was employed to warm the near-surface air and the topsoil.The OTC (60° inclination of the panels, 0.4 m tall and 1.0 m basal diameter) was made of transparent polymethyl methacrylate and placed in the field to achieve year-round warming.The OTC warming method indeed has some limitations (see details in Supplementary Note 5), but has been widely used in remote areas without electricity 46,56,57 .In total, six warming experiments (three thermo-erosion gullies × two areas (thermokarst and non-thermokarst areas)) each with ten replicates were established at the Shaliuhe site.It should be noted that all blocks were laid on vegetated patches rather than on bare soil within the thermokarst landscape (Extended Data Fig. 2f).Despite the inevitable soil translocation after permafrost collapse, topsoils in vegetated patches have retained their original shape due to the protection by mattic epipedon 58,59 , which allows more direct comparison between the collapsed areas and non-collapsed areas.Considering that ignoring bare soil within the thermokarst landscape might lead to bias in extrapolation of the thermokarst effect globally, we conducted a supplementary experiment to explore whether and how warming effects on Rs would differ between bare soil and vegetated patches.Our additional measurements showed that there was a similar response of Rs to experimental warming between bare soil and vegetated patches, indicating no significant influence of ignoring bare soil on the extrapolation (see Supplementary Note 6 for details).
Rs measurements.The Rs and its heterotrophic component were monitored using a Li-8100A Automated Soil CO 2 Flux System (LI-COR Inc.).The trenching method was used to separate Rs components 27,28 .Specifically, we installed two PVC collars (20 cm in diameter) in each plot, one 5-cm-long collar was pushed 3 cm into the soil for measuring Rs, and the other 42-cm-long collar was inserted to a depth of 40 cm for measuring Rh.Given that most roots (~90%) were distributed within 0-40 cm at our study site 13 , the 42-cm-long collar was able to block the majority of roots and thus the respiration from roots.To minimize the effects of soil disturbance, PVC collars were installed in June 2018, one year before the first measurement.Rs and its heterotrophic component were measured three times a month during the growing seasons (May to October) in 2019 and 2020.Living plant materials within the collars were removed at least one day before flux measurement to eliminate aboveground plant respiration 60 .All measurements were conducted at 9:00 a.m.-12:00 p.m. local time.During the period of soil CO 2 flux measurement, we did not observe any obvious loosening contact of soils with the inner edges of the collars (Supplementary Fig. 8).
Although the trenching method has been widely used for separation of Rs components 27,28,61 , the exclusion of active roots in the deep collars would create an artificial situation without root exudates, which fuel a large share of the microbial respiration.To avoid misleading, we used the concept 'potential heterotrophic respiration (Rhp)' to signal that it was not the actual Rh, but a potential Rh under our specific experimental conditions.Considering this point, the Rhp data obtained in this study were only used to support major findings derived from Rs measurements.At the same time as flux measurement, soil temperature and moisture in the 0-10 cm layer were measured adjacent to collars manually using a TP101 portable thermometer (Shanghai Automation Instrument Factory) and a TDR100 portable soil moisture probe (Spectrum Technologies Inc.), respectively.Additionally, we determined topsoil temperature and moisture inside the collars four times during the experimental period, and observed no significant differences between inside and outside collars (Supplementary Fig. 9).

Determination of vegetation and edaphic properties
Vegetation survey.To characterize the associations between Rs response and plant properties, we determined the normalized difference vegetation index (NDVI), vegetation community composition and plant height during the peak growing season (late July) in 2019.Specifically, we set a 25 × 25 cm 2 quadrat near each collar, and measured the NDVI using an ADC multi-spectral digital camera (Tetracam).We also counted the number of species, estimated the relative coverage and measured the maximum height of each species in the quadrat.Based on these measurements, we obtained the NDVI, species richness, Shannon index and community-level mean plant height for each plot.
Soil sampling.We collected 60 samples (three thermo-erosion gullies × two areas (thermokarst and non-thermokarst areas) × ten replicates) of topsoil (0-10 cm; organic horizon with SOC larger than 20%; Supplementary Fig. 7) using an auger (3 cm in diameter) adjacent to the PVC collars during the peak growing season in 2019.In each replicate, three soil cores were collected and then mixed into one composite sample.Subsequently, soil samples were sieved through a 2 mm mesh, with the roots and organic debris being manually picked out.Each sample was subdivided into four subsamples.One subsample was immediately frozen in liquid nitrogen (N) and subsequently stored at −80 °C for DNA extraction, one subsample was stored at −20 °C for microbial biomass analysis and incubation experiment, one subsample was kept fresh at 4 °C to determine the inorganic N concentration and the activities of extracellular enzymes, and the remaining part was air-dried and ground for all other analyses.To determine the differences in abiotic and biotic https://doi.org/10.1038/s41561-024-01440-2factors between collapsed and non-collapsed areas, we measured a series of variables such as soil physicochemical properties, mineral protection, substrate quality and microbial properties.

Soil physicochemical and mineral analyses.
We determined a number of edaphic variables including soil total C and N concentrations, soil inorganic N (sum of NH 4 + -N and NO 3 − -N) and phosphorus concentrations, soil pH, texture, bulk density and the water-filled pore space (see Supplementary Note 7 for details).We also measured the concentration of organically complexed Fe (Fe p ) in the bulk soil.Specifically, Fe p was extracted using 0.1 M sodium pyrophosphate 62 .The soil solutions were then shaken in the dark and centrifuged to measure the amount of Fe in the supernatants.Fe concentration was determined using an iCAP6300 inductively coupled plasma optical emission spectroscopy (Thermo Scientific).

Substrate quality.
Soil substrate quality was analysed on 60 samples using the solid-state 13 C nuclear magnetic resonance technique 63 .Specifically, soil samples were treated with 10% hydrofluoric acid repeatedly, rinsed with deionized water and freeze-dried 36 .Approximately 100 mg of the treated soil samples were analysed using an AVANCE NEO 600 MHz spectrometer (Bruker BioSpin), operating at 100.6 MHz, 8 kHz MAS spinning rate, 2 ms contact time and 6 s recycle delay 25 .The spectra were integrated into the following chemical shift regions, and the relative intensities were acquired as the relative abundances of the corresponding C components: alkyl C (0-50 ppm), O-alkyl C (50-110 ppm), aromatic C (110-165 ppm) and carboxylic C (165-220 ppm) 25 .O-alkyl C is generally considered to be more easily decomposed by soil microorganisms, while alkyl C is recalcitrant 25,36 .The alkyl/O-alkyl ratio was then used to represent substrate quality.

Microbial biomass and enzyme activity.
Considering the important role of soil microorganisms in mediating the response of Rs to experimental warming, we measured soil microbial biomass and the activities of extracellular enzymes.The contents of microbial biomass C (MBC) and N (MBN) were determined by chloroform fumigation 64,65 , with the conversion factors being 0.45 and 0.54, respectively 66,67 .We also assayed the activities of four extracellular hydrolases, β-1,4-glucosidase (BG), β-1,4-N-acetyl-glucosaminidase (NAG), leucine aminopeptidase (LAP) and phosphatase (AP) using fluorescence measurements 68 , as well as one extracellular oxidase (phenol oxidase) measured by spectrophotometry according to ref. 69

(see details in Supplementary Note 8).
Amplicon sequencing and bioinformatic analyses.The bacterial and fungal communities in the topsoil were determined using amplicon sequencing.Briefly, DNA was extracted from ~0.5 g of fresh soil by the DNeasy PowerSoil Pro Kit (Qiagen), and the concentration and quality of the DNA were analysed using a NanoDrop 2000 spectrophotometer (Thermo Fisher Scientific).Primers 515 F (5′-GTGCCAGCMGCCGC GGTAA-3′)/806 R (5′-GGACTACVSGGGTATCTAAT-3′) 70 and ITS3 (5′-GCATCGATGAAGAACGCAGC-3′)/ITS4 (5′-TCCTCCGCTTATTGAT ATGC-3′) 71 were chosen to amplify the hypervariable V4 region of the bacterial 16S rRNA gene and the fungal ITS2 rRNA gene, respectively.The PCR products were purified with the EZNA Gel Extraction Kit (Omega Bio-Tek), and the sequencing library was generated using NEBNext Ultra DNA Library Prep Kit for Illumina (New England Biolabs).Finally, the sequencing library was sequenced by an Illumina Miseq sequencer (Illumina Inc.).Raw reads were quality-filtered and merged using Trimmomatic 72 and FLASH v. 1.2.11(ref.73).Merged reads were clustered into operational taxonomic units (OTUs) at 97% identity based on USE-ARCH (v.11) 74 after removing chimeras.Bacterial and fungal taxonomy of OTUs were annotated using SILVA v. 123 (ref.75) and UNITE v. 7.1 (ref.76), respectively.OTU tables were rarefied at 54,963 for bacteria and 41,219 for fungi, respectively (see Supplementary Note 9 for more details).
Metagenomic sequencing and bioinformatic analysis.Shotgun metagenomic sequencing was used for analysing the topsoil microbial functional genes.Specifically, we performed 2 × 150 bp paired-end sequencing using the Illumina HiSeq 2500 platform.A total of 5.45 billion reads with an average of 13.6 Gbp (s.e.m. = 0.2 Gbp) of sequence data were obtained for each of the 60 samples.Raw sequences were quality-controlled using Trimmomatic (v.0.39) 72 .The clean data were then de novo assembled (k-mer options: k-min 35, k-max 95, k-step 20) using MEGAHIT (v.1.2.9) 77 .After assembly, scaftigs longer than 500 bp were retained and predicted using Prodigal (v.2.6.3) 78o obtain open reading frames (ORFs).The predicted ORFs were dereplicated using MMseqs2 (ref.79), with the longest ORF used as a representative sequence for each cluster by setting the similarity threshold as 95% and coverage as 90%.Finally, functional genes were annotated by the Kyoto Encyclopedia of Genes and Genomes (KEGG) database 80 .To evaluate the gene abundance, clean reads derived from each sample were mapped to the gene clusters using Bowtie2 (v.1.0) 81 , and counts and lengths of genes and alignments were extracted using SAMtools (v.1.3.1) 82.Abundance of genes was normalized to transcripts per million.

Incubation experiment
To validate the generality of the results obtained from the field warming experiments, we collected topsoil (0-10 cm) samples from five additional thermokarst-impacted sites in Ebo, Mole, Reshui, Huashixia and Huanghe (Extended Data Fig. 1) in August 2020.The thermokarst landforms for all sites belong to the thermo-erosion gully (Supplementary Table 1).We established a collapsed plot (10 × 10 m 2 ) inside the gully and a paired control plot in the adjacent undisturbed grassland for each thermokarst-impacted site.In each plot, we sampled five quadrats (1 × 1 m 2 ) at the centre and the four corners.Three soil cores were collected in each quadrat and homogenized as one replicate, with five replicates being ultimately obtained from each plot (n = 5).In total, 50 topsoil samples were collected from the five thermokarst-impacted sites (five thermo-erosion gullies × two plots (collapsed versus non-collapsed plot) × five replicates).All soil samples were sieved through a 2 mm mesh and we manually picked out visible roots.
An incubation experiment was conducted to estimate the temperature sensitivity of soil CO 2 emission (that is, Q 10 ) for the Shaliuhe site and the five additional thermokarst-impacted sites.Samples from all sites were incubated at two temperatures (5 and 15 °C, within the range of mean monthly air temperature across study sites) during a 60-day period.Specifically, 10 g of fresh soil was weighed into a 150 ml airtight amber jar.We adjusted soil moisture to 60% water holding capacity and kept it constant during the incubation by adding demineralized water to the samples.All samples were pre-incubated for seven days at 5 °C to recover them from physical disturbances.The rate of CO 2 release during the next 60 days of incubation was measured every 2-3 days for the first week, every 4-7 days for the next three weeks and every 10 days for the remainder of the period.All incubation jars and three empty jars were flushed using CO 2 -free air before each measurement, with the empty jars used as blanks.We collected 10 ml of headspace air after 10 hours of incubation at 5 °C and 15 °C with gas-tight syringes for CO 2 analysis using an EGM-5 infrared gas analyser (PP Systems).Notably, CO 2 dissolution was not considered because this process should be negligible compared with biotic CO 2 release in our acid and organic-matter-rich soils 83,84 .

Statistical analyses
Data were processed following four steps.First, we calculated the normalized Rs change and apparent Q 10 for the field experiments as well as Q 10 for the laboratory incubation experiment.Specifically, the normalized Rs change was adopted to facilitate comparison of the warming effects between collapsed and non-collapsed areas, that is, the warming-induced changes in Rs were normalized to per unit https://doi.org/10.1038/s41561-024-01440-2temperature change using equation ( 1) described in the main text.Apparent Q 10 was calculated based on the seasonal Rs (Supplementary Note 1).The Q 10 for laboratory incubation was calculated based on the rate of CO 2 release at the two incubation temperatures as follows 34 : where R w and R c are the average rates of CO 2 release (mg CO 2 -C g −1 d −1 ) at 15 °C (T w ) and 5 °C (T c ), respectively.
Second, linear mixed-effects models were applied to compare the differences of normalized Rs change, apparent Q 10 , abiotic and biotic factors (vegetation characteristics, edaphic properties, substrate quality and microbial attributes) between collapsed and non-collapsed areas.In the mixed-effects models, collapse status was included as a fixed effect, and plot nested in thermo-erosion gully was treated as a random effect for the field warming experiments.Accordingly, for the incubation experiment, we treated collapse status as a fixed factor and plot as a random factor in the mixed-effects models.
Third, we explored the relationships of the normalized Rs change with various factors (vegetation properties, edaphic properties, substrate quality and the abundance of microbial functional genes) using mixed-effects models.Here, the explanatory variables were included as fixed factors, and plot nested in thermo-erosion gully was treated as a random factor.The normality of the residues of the models was tested, and the data were log-transformed when necessary.The linear mixed-effects models were implemented using the nlme package 85 .
Finally, SEM was constructed to elucidate the direct and indirect linkages of the biotic and abiotic factors to the normalized Rs change.Before SEM analysis, we constructed an a priori conceptual model of hypothetical relationships (Supplementary Fig. 5).Specifically, we first chose those variables that exhibited significant changes upon permafrost collapse and also were significantly correlated with the normalized Rs change.Then we hypothesized that (1) soil substrate quality was directly linked to normalized Rs change based on the 'C quality-temperature' hypothesis; (2) microbial functional genes involved in SOM decomposition were directly associated with normalized Rs change because soil microorganisms are carriers of the SOM decomposition process; moreover, microbial functional genes and substrate quality could interact with each other; and (3) soil texture and moisture were linked to normalized Rs change directly via mediating mineral absorption and oxygen diffusion, and indirectly through their effects on substrate quality and microbial functional genes.Considering that variables characterizing substrate quality or microbial functional genes were closely correlated, we performed principal component analysis (PCA) to generate comprehensive metrics to represent each group.The first axis of substrate quality and the weighted mean of the first and second components for microbial functional genes, which explained 77.6% and 86.1% of the variations respectively, were selected for constructing the final model (Supplementary Table 2).The SEM was completed using R packages piecewiseSEM 86 and nlme 85 .All above analyses were also conducted for Rhp.All analyses were conducted in R 4. which were formed by the fragmentation of the soil surface upon permafrost collapse.Given that the upper soil was protected by mattic epipedon which had an intensive root network protecting the soil against interference 58,59 , these vegetated patches maintained their original soil structure, despite the translocation of the exposed soil upon permafrost collapse 13 .Photograph credit: G.Q. Wang.

Fig. 1 |
Fig. 1 | Response of Rs to permafrost collapse and experimental warming at the Shaliuhe site.a,c, Seasonal dynamics and averages of Rs during the growing seasons in 2019 (a) and 2020 (c).Data are presented with mean values ± s.e.m. (n = 30).Inserts show effects of permafrost collapse, experimental warming and their interactions on seasonal averages of Rs. b,d, Normalized Rs change in 2019 (b) and 2020 (d).The ends of the boxes represent the 25th and 75th percentiles, and the whiskers refer to the 10th and 90th quartiles.The white dot in each box shows the average (n = 30).These comparisons were analysed based on linear mixed-effects models with one-sided F tests.***P < 0.001.

Fig. 2 |
Fig.2| Difference in Q 10 between collapsed and non-collapsed areas for the six thermokarst-affected sites.a, Q 10 between collapsed and non-collapsed areas for the Shaliuhe site (SLH).Samples (n = 30) were collected from three thermo-erosion gullies in 2020.The ends of the boxes represent the 25th and 75th percentiles, and the whiskers represent the 10th and 90th quartiles.The white dot in each box represents the average.b, Q 10 between collapsed and non-collapsed areas for the five additional sites.Samples (n = 5) were collected in 2020.The dashed lines denote the average.EB, Ebo site; ML, Mole site; RS, Reshui site; HSX, Huashixia site; HH, Huanghe site.These comparisons were analysed based on linear mixed-effects models with one-sided F tests.ns, non-significant; *P < 0.05; **P < 0.01; ***P < 0.001.

Fig. 3 |
Fig. 3 | Changes in substrate quality and microbial functional gene abundances upon permafrost collapse and their relationships with the normalized Rs change.a-i, Changes in O-alkyl C (a), alkyl C (b), alkyl/O-alkyl ratio (c) and the relative abundances of C degradation genes as transcripts per million reads mapped (TPM) upon thermokarst formation for cellulose (d), hemicellulose (e), chitin (f), pectin (g), aromatics (h) and lignin (i).The ends of the boxes represent the 25th and 75th percentiles, and the whiskers represent

Fig. 4 |
Fig. 4 | Effects of abiotic and biotic factors on the normalized Rs change.a, SEM of the direct and indirect pathways that link to the normalized Rs change.Solid and dotted arrows represent positive and negative relationships, respectively.Arrow width is proportional to the strength of the relationship.The numbers adjacent to the arrows denote standardized path coefficients.Substrate quality is the first component from the PCA of the variables listed inside the rectangles, with ↑ and ↓ indicating positive and negative correlations between the factors and the first component, respectively.Microbial functional genes represent the weighted average of the first and second components from the PCA of the variables inside the rectangle, with ↑ indicating a positive correlation between the factors and the weighted average of the first and second components.*P < 0.05; **P < 0.01; ***P < 0.001; AIC, Akaike information criterion.b, Standardized direct effects of substrate quality, abundance of microbial functional genes and indirect effects of edaphic variables (clay+silt and soil moisture) on the normalized Rs change.
2.2 (ref.87).Extended Data Fig. 1 | Locations and landscapes of six thermokarst-impacted sites across the Tibetan Plateau.a, Positions of six study sites shown on the background of permafrost map on the Tibetan Plateau.The pentagon refers to the site where the in situ warming experiments were conducted and the triangles denote the locations of regional thermokarst soil sampling.The insert shows regional location on the Tibetan Plateau.b-f, Individual features of the thermokarst-impacted sites at EB (b), RS (c), HSX (d), ML (e) and HH (f).EB, Ebo site; RS, Reshui site; HSX, Huashixia site; ML, Mole site; HH, Huanghe site.The background map of permafrost distribution is from ref. 89.Photograph credit: G.Q. Wang.Extended Data Fig. 2 | Landscapes of the three thermo-erosion gullies and views of warming experiments at the Shaliuhe site.a-c, Photographs of the three thermo-erosion gullies at the Shaliuhe site.d,e, Overall views of warming experiments in non-collapsed (d) and collapsed (e) areas from panel (a).f, A photo showing the warming treatment on the vegetated patch.The opentop chambers were built within the vegetated patches in the collapsed plots,

Fig. 3 |
Responses of potential heterotrophic respiration (Rhp) to permafrost collapse and experimental warming at the Shaliuhe site.a,c, Seasonal dynamics and averages of Rhp during the growing seasons in 2019 (a) and 2020 (c).Data are presented with mean values ± s.e.m. (n = 30).Inserts show effects of permafrost collapse, experimental warming and their interactions on seasonal averages of Rhp.b,d, Normalized Rhp change between collapsed and non-collapsed areas in 2019 (b) and 2020 (d).The ends of the boxes represent the 25 th and 75 th percentiles, and the whiskers refer to the 10 th and 90 th quartiles.The white dot in each box shows the average (n = 30).These comparisons were analyzed based on linear mixed-effects models with one-sided F tests.***, P < 0.001.Extended Data Fig. 4 | Responses of topsoil (0-10 cm) temperature and moisture to permafrost collapse and experimental warming at the Shaliuhe site.a,b, Effects of permafrost collapse, experimental warming and their interactions on soil temperature in 2019 (a) and 2020 (b).c,d, Effects of permafrost collapse, experimental warming and their interactions on soil moisture in 2019 (c) and 2020 (d).The ends of the boxes represent the 25 th and 75 th percentiles, and the whiskers refer to the 10 th and 90 th quartiles.The white dot in each box represents the average (n = 30).Both soil temperature and moisture were measured manually during Rs monitoring.These comparisons were analyzed based on linear mixed-effects models with one-sided F tests.Extended Data Fig. 5 | Effects of experimental warming on temperature and moisture in different soil layers of non-collapsed areas at the Shaliuhe site.a-c, Daily soil temperature at depth of 10 cm (a), 20 cm (b) and 30 cm (c) from January to December.d-f, Daily soil moisture at depth of 10 cm (d), 20 cm (e) and 30 cm (f) from January to December.The inserts show the averages of soil temperature and moisture across the growing season (May -October).ST, soil temperature; SM, soil moisture.The ends of the boxes represent the 25 th and 75 th percentiles, and whiskers refer to the 10 th and 90 th quartiles.The white dot in each box represents the average.Significant differences between warming and control plots were determined based on two-sided t tests.ns, non-significant; ***, P < 0.001.Both soil temperature and moisture were measured by EM50 data loggers (Decagon Devices, Washington, USA) in 2019.Extended Data Fig. 6 | Comparisons of apparent temperature sensitivity (Q 10 ) of soil respiration (Rs) and potential heterotrophic respiration (Rhp) between collapsed and non-collapsed areas at the Shaliuhe site.a,b, Apparent Q 10 of Rs during the growing season in 2019 (a) and 2020 (b).c,d, Apparent Q 10 of Rhp during the growing season in 2019 (c) and 2020 (d).Apparent Q 10 was estimated for each gully by fitting the seasonal observations of soil CO 2 flux and soil temperature in the top 10 cm (see Supplementary Note 1 for detailed calculation).Data are presented with mean values ± s.e.m. (n = 3).Significant differences between collapsed and non-collapsed areas were determined based on linear mixed-effects model with one-sided F tests.**, P < 0.01.Extended Data Fig. 7 | Temperature sensitivity (Q 10 ) estimated on each measurement occasion during the 60-day laboratory incubation for the six sampling sites.a, Q 10 between collapsed and non-collapsed areas at the Shaliuhe site.Samples (n = 30) were collected from three thermo-erosion gullies in 2020.b-f, Q 10 between collapsed and non-collapsed areas for the five additional sites.Samples (n = 5) were collected in 2020.SLH, Shaliuhe site; EB, Ebo site; ML, Mole site; RS, Reshui site; HSX, Huashixia site; HH, Huanghe site.Data are presented with mean values ± s.e.m.Significant differences between collapsed and noncollapsed areas were determined based on linear mixed-effects models.†, P < 0.1; *, P < 0.05; **, P < 0.01; ***, P < 0.001.Samples were collected and analyzed in 2020.Extended Data Fig. 8 | Comparison of normalized changes in topsoil (0-10 cm) enzyme activities between collapsed and non-collapsed areas at the Shaliuhe site.a,b, Normalized changes in the hydrolase (a) and oxidase (b) activities.Normalized changes of enzyme activities (ΔEnzyme per temperature change = (En Warming − En Control )/(T Warming − T Control ), where En Warming and En Control are soil enzyme activity (nmol g -1 h -1 ) in the warmed and control plots, and T Warming and T Control are temperature ( o C) in the warmed and control plots, respectively) were calculated in collapsed and non-collapsed areas, and weighted by the temperature difference between the warming and control treatments (T Warming -T Control ) to eliminate the difference of warming effects between collapsed and non-collapsed areas.The ends of the boxes represent the 25 th and 75 th percentiles, and the whiskers refer to the 10 th and 90 th quartiles.The white dot in each box shows the average (n = 30).These comparisons were analyzed based on linear mixed-effects models with one-sided F tests.ns, non-significant; *, P < 0.05.Samples were collected from the three thermo-erosion gullies in 2019.Extended Data Fig. 9 | Thermokarst-induced differences in topsoil (0-10 cm) microbial diversity and their relationships with the normalized soil respiration (Rs) change at the Shaliuhe site.a-d, Relationships of normalized Rs change with Shannon and Simpson indices for bacterial (a, b) and fungal (c, d) communities.The inserts denote differences in Shannon and Simpson indices upon permafrost collapse.Data are presented with mean values ± s.e.m. (n = 30).Significant differences between collapsed and noncollapsed areas were determined based on linear mixed-effects models with one-sided F tests.ns, non-significant; *, P < 0.05.Samples were collected from the three thermo-erosion gullies in 2019.Extended Data Fig. 10 | Effects of abiotic and biotic factors on the normalized potential heterotrophic respiration (Rhp) change.a, Structural equation model (SEM) of the direct and indirect pathways that affect the normalized Rhp change.b, Standardized direct effects of substrate quality and the abundances of microbial functional genes and the indirect effects of edaphic variables (clay+silt and soil moisture) on the normalized Rhp change.Solid and dotted arrows represent positive and negative relationships, respectively.Arrow width is proportional to the strength of the relationship.The numbers adjacent to the arrows denote standardized path coefficients.Substrate quality is the first component from the principal component analysis (PCA) conducted with the factors listed in the rectangles, with "↑" and "↓" indicating positive and negative correlations between the factors and the first component, respectively.Microbial functional genes represent the weighted average of the first and second components from the PCA of the factors in the rectangle, with "↑" indicating a positive correlation between the factors and the weighted average of the first and second components.*, P < 0.05; **, P < 0.01; ***, P < 0.001.AIC, Akaike information criterion.