Wang et al. (2023). Responses of soil organic carbon to climate extremes under warming across global biomes

Satellite-based observations and global observation networks (e.g., FLUXNET) have shown

evidence that even a single climate extreme event can cause significant reduction in regional

ecosystem production and its effects are usually time-lagged or last for years up to decades

5-7

However, the relatively short-term nature of satellite-based and site-level records (e.g., satellite

observation was only available from 1980s) are incapable of predicting the long-term

consequences of climate extreme shifts on the global carbon cycle

3,8

Amongst terrestrial carbon pools of the global carbon cycle, soil organic carbon (SOC) is

the largest and more than twice of the carbon stored in vegetation biomass and the atmosphere

9

Although the response of SOC to mean climate changes like warming has been extensively

studied using modelling and experimental approaches

10,

, the role of accompanied shifts of

climate extremes in regulating such response has rarely been quantified based on observational

data, particularly across large spatiotemporal extents

3

. A major challenge is that SOC changes

would be only detectable or fully manifested at the time scale of SOC turnover which is usually

decadal or even centurial in most soils

,11

. For this reason, observations after one or several

extreme events occurred in a typical year or short period cannot fully capture long-term impacts

of more extreme climate on SOC dynamics. Furthermore, replicating climate extremes in

experiments is challenging, and different extremes (e.g., droughts and heat waves) may elicit

distinct responses in soil carbon fluxes, potentially altering SOC balance

12,13

Climate extremes exert influence on SOC dynamics at various temporal scales. In the short

term, events like heavy precipitation can lead to soil waterlogging, inhibiting microbial

decomposition and/or causing spatial redistribution of SOC through erosion

14,15

. Drought alter

carbon inputs to the soil in terms of their amount, quality and timing through its effect on plant

growth and soil hydrothermal regimes

8,16

. Heatwaves can accelerate SOC decomposition

17

, but

often accompanied by drought, resulting in complex effect on SOC.

Over the long term,

different climate extremes can bring about a series of ecological and environmental shifts

3

such

as: 1) changes in soil physiochemical environment (e.g., soil aggregate transformation and

stabilization, hydrophobicity), 2) alterations in plant, soil microbial and animal community

composition and functionality

3,18

, and 3) impacts on soil hydrology and nutrient cycling

Overall, shifts in climate extremes may directly or indirectly influence long-term SOC

dynamics depending on local environmental conditions and the specific type of climate

extremes

20

3

Here we quantify the response of SOC to shifts of various climate extremes which are

projected under future climate

21,22

. This quantification builds upon SOC stock (kg m

–2

)

measurements in the 0–30 cm topsoil at 113,013 sites across the globe (ref.

, Fig. S1), along

with eight indices representing the magnitude and frequency of different extreme events (Fig.

S2, Table S1). These events include heat waves, cold waves, extreme dry and wet days from a

global map product derived by the Expert Team on Climate Change Detection and Indices

24,25

Another 25 indices were also obtained from the same map product, reflecting the threshold and

probability of temperature- and precipitation-related extremes based on long-term climate

records

26

(Table S1). Using the data, a hybrid approach combining space-for-time substitution

and meta-analysis techniques

(Fig. 1) was used to estimate the responses of SOC to two

general climate change scenarios: a warmer climate without extreme shifts (W, i.e., only mean

annual temperature is increased, which is assigned to 1.5 °C in this study), and a warmer and

more extreme climate (W+E, 1.5 °C warming plus more extreme climate, Fig. 1). Uncertainties

in future climate extremes

27

were considered by using a five-level change gradient for each of

the 33 climate extreme indices (CEIs) (Table S1). Compound extreme events (i.e., concurrent

occurrence of multiple climate extremes) may also be more common in the future

28

, but are

more difficult to predict. To demonstrate their importance, here we illustrate the influence of

shifts in two compound extremes involving two widely recognized individual extremes: heat

wave + extreme dry, and heat wave + extreme wet (Methods).

Climate extremes diversify the response of SOC to warming

Averaging across the globe, SOC is reduced by 10.5% (8.1–12.3%, 95% confidence interval)

under 1.5 °C warming (i.e., W, Figs. 2 and S3), consistent with our previous estimate at a

similar warming level

. However, this SOC reduction is significantly altered by shifts in

climate extremes (Figs. 2 and S3). For 8 out of 33 CEIs (heat wave frequency, extreme wet

magnitude, extreme wet frequency, extreme dry magnitude, max daily minimum temperature,

maximum daily precipitation, maximum precipitation of five consecutive days, precipitation

intensity), more extreme condition consistently exacerbates SOC reduction (i.e., additional

SOC reduction, Figs. 2 and S3), with the additional reduction increasing with higher extreme

levels (Figs. 2 and S3). For example, SOC reduction increases from 2% with a two-day increase

in heat wave frequency to 21% with a ten-day increase (Fig. 2). On the contrary, SOC reduction

is consistently attenuated or even reversed for two extremes (cold wave magnitude and coldest

daily minimum temperature), offsetting the negative effect of warming on SOC under more

4

extreme conditions (Figs. 2 and S3).

Several CEIs exhibit threshold or turning point regarding the rate or direction of warming-

induced SOC changes (Figs. 2 and S3). Specifically, a decrease in cold wave frequency (CWF)

initially leads to less SOC reduction until reaching a turning point of a four-day decrease, after

which further decrease in CWF results in more SOC reduction (Fig. 2). A similar pattern was

also detected for the length of summer days, the number of cool days and nights, and the

number of hot (tropical) nights (Fig. S3). For some other CEIs (i.e., extreme dry frequency,

min daily maximum temperature, max daily minimum temperature, growing season length,

number of icing days, heavy precipitation, consecutive wet days), SOC changes have an

apparent peak at certain CEI changes (Figs. 2 and S3). These varied SOC responses may stem

from trade-offs between changes in carbon inputs to soil and decomposition of SOC under

more extreme climate

29-31

, or distinct effects of climate extreme shifts on different ecosystem

processes (e.g., microbial carbon use efficiency

32

) related to SOC dynamics. For instance, the

decrease of cold wave frequency may not only benefit plant growth thereby enhancing carbon

inputs to soil but also stimulate microbial decomposition. The net effects on SOC balance

depend on the proportionate changes in carbon inputs and outputs.

The responses of SOC to single extremes such as heat wave, extreme dry, and extreme wet

are not significantly different from the responses to their compound events at their low-level

changes (i.e., 20% increase in their frequency, Fig. 2b). However, at high-level changes (i.e.,

40% increase in frequency) of both individual events, compound heat wave and extreme dry

(HWF+EDF) leads to an additional 13% SOC reduction, whereas heat wave and extreme dry

events individually only lead to additional 2% and 5% SOC reduction, respectively (Fig. 2b).

Similarly, additional SOC reduction induced by compound heat wave and extreme wet event

was nearly double compared to that induced by single shifts of heat wave or extreme wet under

high-level changes (Fig. 2b). These results demonstrate that compound climate extremes may

cause more SOC reduction under warming (especially at high-level changes) and have a more

profound and complex impacts compared to single extremes.

Global average SOC changes induced by climate extremes

We estimated percentage SOC changes attributed to climate extremes (ΔSOC

E

) under the W+E

scenario by calculating the difference between the percentage responses to W+E and to W.

Then, a linear mixed-effects model with biome-dependent slopes but without intercept (i.e.,

5

additional change should be zero under zero climate extreme changes) was fitted to infer the

relationship of ΔSOC

E

with climate extreme change levels. The regression slopes can be also

explained as the percentage of additional SOC changes induced by unit changes in climate

extremes.

On average ΔSOC

E

shows distinct relationships with the 33 CEIs in terms of both

magnitude and direction (i.e., the regression slopes for the fixed effect of climate extremes, Fig.

3a). Negative relationships are dominant (22 out of 33). However, only the length of hot

(tropical) nights (TR), cold spell duration index (CSDI), extreme dry magnitude (EDM), and

extreme wet magnitude (EWM) exert significant linear effects (i.e., the regression slope is

significantly different from zero) on ΔSOC

E

(Fig. 3a). Specifically, the slope for TR is –8%,

that is, additional percentage of SOC reduction is 8% for every five-day increase of TR. SOC

increases by 22% by per five-day decrease of CSDI. For EWM, SOC decreases by 7% per 0.1

unit increase of EWM, and some other CEIs related to extreme precipitation magnitude also

show negative effects on SOC. Positive relationships between SOC and climate extremes

mainly occur for extremes related to shortening of low temperature days such as frost days and

cold spells (Fig. 3a). In addition, the increase of low temperature-related extremes represented

by day-time and night-time minimum temperature have positive effects, while CEIs related to

high temperature extremes show negative effects, especially the maximum and range of daily

temperature. For EDM, every 0.1 unit increase leads to additional 16% SOC reduction, and the

extreme dry frequency also shows a negative effect (Fig. 3a).

Biome-specific SOC changes induced by climate extremes

SOC in different ecosystems show distinct responses to climate extremes (Figs. 3b and S4). In

temperate grasslands and croplands, SOC is negatively influenced by most climate extremes

assessed (23 and 20 out of 33, respectively), followed by tropical/subtropical forests (16 out of

33, Fig. 3b). The impacts of climate extremes on SOC in croplands may be modulated by

management intervention. Irrigation, for example, may stimulate SOC decomposition via its

positive effects on soil moisture. Besides, many crop varieties are genetically selected for high

production in optimal conditions, so they would be more sensitive to extreme events

33

. In other

ecosystems except tundra, positive effects are more common (Fig. 3b). Especially, only four

CEIs show significant negative effects in boreal forests. In tundra, 21 CEIs exert neutral (i.e.,

insignificant) effects, while only does decreasing cold wave magnitude show positive influence.

The remaining 11 extremes show negative effects, and, more importantly, the magnitude of

6

most effects is much stronger in tundra than their corresponding effects in other ecosystems

(Figs. 3b and S4). Tundra has been suggested to be the most vulnerable ecosystem to climate

change

34,35

. More extreme climate may worsen the vulnerability of SOC stocks in tundra.

The same climate extreme exerts distinct effects on SOC across ecosystems (Figs, 3b and

S4). For example, increase of heat wave magnitude has significant negative effect on ΔSOC

E

in tropical/subtropical ecosystems (forests and grasslands/savannas), but positive in temperate

forest, Mediterranean/montane shrublands and croplands (Fig. S4). Increase of extreme dry

frequency negatively influences ΔSOC

E

in five ecosystems but positively in temperate forest

and deserts. In addition, for the same type of extremes, the changes of its frequency and

magnitude (i.e., the intensity or severity of one event) have different consequences on ΔSOC

E

(Fig. S4). For example, unlike the positive effects of extreme dry frequency in temperate forests

and deserts, the positive effect of extreme dry magnitude occurs in boreal forests and

tropical/subtropical grasslands/savannas. It is noteworthy that cold spell duration and extreme

wet magnitude are the only two CEIs showing consistent positive and negative effects (if the

effect is significant) on SOC in all ecosystems, respectively (Fig. S4).

Comparison with field experiments

Drought is the most experimentally explored climate extreme in terms of its effect on

ecosystem carbon balance, albeit experimental duration is relatively short

. We compared our

estimates to a recent meta-analysis synthesizing global drought experiments

(Fig. S5). We

integrated ΔSOC

E

induced by EDM into relevant three general ecosystems in order to conduct

a comparison with field drought experiments (limited observational data do not allow us to do

a biome-to-biome comparison, Fig. S5). In grasslands, our estimation of significant negative

ΔSOC

E

is comparable to the estimation of field experiments. Forest ecosystems on average are

less vulnerable to EDM than grasslands, and our estimation is not significantly different from

experimental estimations which, however, show much larger uncertainties (Fig. S5). In

shrublands, we estimate comparable negative ΔSOC

E

between grasslands and shrublands, but

field observations are generally positive albeit insignificant. Averaging across the globe, field

manipulative experiments observed smaller response of SOC to the same level of EDM

changes (Fig. S5). In the short time (e.g., the duration of field experiments ranges from 1 to 13

years

), drought usually induces immediate reduction of soil respiration, yet carbon inputs

may do not change much and thereby exerting limited effect on net SOC balance. In the long

term, however, EDM increase may substantially decrease carbon inputs to soil through

7

increased plant mortality and decreased vegetation productivity

36

, which cannot be fully

captured by relative short manipulative experiments.

Drivers of the responses of SOC to climate extremes

We trained a machine learning model - random forest (RF) - to explore controls over the

variance of ΔSOC

E

for the eight CEIs

indicating the magnitude and frequency of heatwave,

cold waves, extreme dry and wet events (Methods). For all eight CEIs, the RF model driven by

climatic attributes related to mean and extreme climate conditions, soil properties, vegetation

characteristics, and topographic attributes explains >70% of the variance of ΔSOC

E

(Figs. 4a

and S6). Baseline SOC level is the most important individual predictor, and alone explains

15~18% of the variance of ΔSOC

E

. For all eight CEIs, ΔSOC

E

is negatively associated with

baseline SOC levels and remains constant until SOC reaches a threshold of ~300 Mg C ha

–1

(Fig. 4b). However, it should be acknowledged that baseline SOC itself should not be explained

as the driving force of SOC dynamics

37

. We find that those carbon-rich soils mainly distribute

nearby or in wetlands including peatland and coastal areas, and another small fraction in tundra

which has relatively low temperature (Fig. S7). In these systems, particulate organic carbon

(POC) usually dominates SOC

38

and is more sensitive to climate change than other SOC

components like mineral-associated organic carbon (MAOC)

39

, which may result in general

strong response of SOC to climate change. However, POC has much higher C:N ratio than

MAOC

38

, its decomposition is strongly dependent of nutrient availability which would limit its

response to climate change when POC is quite high

40

. In addition, under climate change,

particularly warming and drought-relevant events, wetland soils may become less anaerobic

and suffer from subsidence

41

, thereby altering SOC density and the estimation of SOC changes

based on soil depth. Overall, wetlands involve complex interactive responses of aboveground

(e.g., vegetation shift) and belowground process (e.g., oxidation, subsidence, shift of anerobic

and aerobic conditions) to climate change, which may play a vital role in determining wetland

SOC balance

42

Baseline climatic attributes together contribute the largest fraction (42-51%) to the

explained variance of ΔSOC

E

for all CEIs (Fig. 4a), and their influence is

obviously

nonlinear

(Figs. 4 and S8). Mean annual temperature negatively influences ΔSOC

E

(Fig. 4c), while the

influence of mean annual precipitation is marginal and much weaker except for extreme dry

magnitude (Fig. 4d vs 4a), demonstrating that SOC is generally more vulnerable to extreme

dry event in warmer and wetter regions. However, extreme climate conditions together are

8

more important than mean climate, suggesting that baseline climate extreme conditions play a

vital role in determining the effect of future climate extreme shifts on SOC. Especially, we

found strong dependency of ΔSOC

E

induced by a typical CEI on the corresponding baseline

CEI (Fig. S8). For example, if the baseline magnitude of EDM is already high, more extreme

conditions have smaller effect on ΔSOC

E

(Fig. S8). The positive effect of decreasing cold wave

magnitude (CWM) is much stronger in regions with higher baseline CWM.

Although our study does not explicitly assess the underlying mechanisms, biome-specific

responses to climate extremes are likely associated with interactive effects of water and energy

on both carbon inputs and outputs over space and time. Energy-limited systems (e.g., boreal

forests and tundra) would be more sensitive to temperature-related extremes. In boreal forests,

for example, the reduction of extreme cold events substantially stimulates plant growth

superseding enhanced SOC decomposition

43,44

, resulting in net SOC gains

45,46

. This is

corroborated with the positive effects of shortening of cold spell, lengthening of warm spell,

and other temperature-related extremes (Figs. 3 and S4). In tundra, reduction of energy

limitation may accelerate permafrost thaw and thereby stimulating microbial growth and SOC

decomposition

34,35

. Water-limited systems (e.g., grasslands and savannas) would be more

sensitive to extremes that alter precipitation seasonal pattern or worsen/loosen water limitation.

Additionally, the characteristics of vegetations (e.g., canopy height, litter layer, rooting depth)

would also mediate how ecosystems respond to climate extremes. Woody systems (e.g., forests

and shrublands) may be more tolerant to extreme wet conditions than herbaceous systems (e.g.,

grasslands and croplands)

47

, and the existence of organic layer would increase soil water

holding capacity and prevent or weaken the loss of particulate and dissolved organic carbon

via erosion or runoff

48

. Indeed, our results show that extreme wet frequency shows positive

effects in forest and shrubland ecosystems, but negative in grasslands and croplands (Fig. S4).

Overall, climate extreme shifts may alter water-energy interactions depending on local climatic

conditions, which then change plant performance and SOC decomposition and stabilization

processes such as shifts in microbial community functioning and substrate utilization

strategies

18,49

and thereby inducing diverse effects of climate extremes across ecosystems.

Global pattern of SOC responses to climate extremes

We apply the cross-validated RF models (Fig. S6) to map ΔSOC

E

induced by changes in the

eight CEIs (i.e., the magnitude and frequency of heatwave, cold wave, extreme dry, and

extreme wet) across the globe (Figs. 5 and S9) and its uncertainty (Figs. S10 and S11). Global

9

average additional ΔSOC

E

is estimated to be –139 ± 7.1, –98 ± 4.5, –27 ± 6.1, –4.5 ± 6.3, –125

± 5.8, –106 ± 5.2, –130 ± 5.8, –160 ± 5.8 Pg (mean ± standard error) for the eight CEIs,

respectively (Table S3). There are hotspots of SOC response, albeit the global average SOC

change is small under some extremes (Figs. 5, 11, and S9). For example, a 0.15 unit increase

in drought magnitude (which is equivalent to a 40% increase of its global average) only leads

to additional <1% reduction of global SOC (Fig. S12e), but induces 50 Pg SOC reduction in

tropical and subtropical regions, north America and Europe (Table S3), and increases SOC in

many high latitudinal regions, especially in Siberia (Figs. 5e and S12e). Six-day increase of

annual heat wave frequency results in 6% additional SOC reduction averaging across the globe

(Fig. S12b). This reduction is mainly in tropical and subtropical regions, north America, and

north Europe (Figs. 5b and S12b). A 6% increase of drought frequency does not significantly

alter global SOC (Fig. S12f) but leads to substantial additional reduction in north European and

north American wetlands (Fig. 5f).

More intense (0.15 unit increase in magnitude) and frequent (6% increase in frequency)

precipitation extremes result in SOC reduction by 7% and 2%, respectively, across the globe

(Fig. S12g, h). The most vulnerable regions to extreme precipitation are Europe, north America,

and some tropical and subtropical regions (Figs. 5g, h). However, more frequent heavy

precipitation can offset warming-induced SOC reduction in dry and semi-arid regions, such as

central Asia, Australia, south Africa and Sahel region (Fig. 5h). Less cold waves significantly

decrease SOC reduction (i.e., positive ΔSOC

E

) in most regions of the globe (Figs. 4b, c and

S12b, c). Our results indicate that there are distinct hotspots in ΔSOC

E

depending on the type

of climate extremes. However, there are regions where SOC is vulnerable to all climate

extremes. These regions mainly locate in Europe, southeast Asia, north America, and Amazon.

Unfortunately, these regions are also hotspots experiencing more frequent and intense climate

extremes

50,51

Limitations and uncertainties

The space-for-time substitution approach has some limitations and uncertainties, e.g., without

considering other global change factors (e.g., nitrogen deposition and wildfire) and the steady

state assumption

. However, the impacts of these factors would be relatively small compared

to the impacts of climate change represented by warming and shifts in climate extremes. In

addition, other global change factor usually exerts much larger spatial variability (e.g., wildfire),

and the relevant data is relatively limited. Croplands are the land use type greatly affected by

10

human activities in terms of both cultivation length (from years to millennia) and land

management. Thus, SOC stocks are unlikely to be under steady state in croplands. However,

excluding soil profiles from croplands, the estimated response of SOC to climate extremes did

not apparently change (Fig. S13), suggesting that steady state assumption would have little

effect on our SOC change projection under projected climate extremes. A number of other

sensitivity analyses have also demonstrated the robustness of the approach (Methods, Fig. S14-

16). Another limitation is that only were two compound events involving two types of climate

extremes assessed. The natural environment may experience multiple compound events over

different time scales. We acknowledge that more analyses focusing on compound events and

their influence are needed

52

. However, potential interconnections among different types of

climate extremes involved in compound events need to be clearly defined and quantified, which

itself remains a challenge

28

Conclusions

We used a novel approach to infer SOC changes induced by a warmer climate with more

climate extremes. Our approach implicitly takes long-term transition of ecosystems into

account, and quantifies the SOC difference between steady states under contrasting climate

conditions, which is otherwise impossible by conducting field experiments or modelling

without explicit understanding of the underlying mechanisms. The results reveal distinct

(positive, neutral or negative) responses of SOC to the shifts of different types of climate

extremes. On a global scale, on average, most climate extremes may exacerbate SOC loss under

future climate. Among ecosystems, however, the magnitude and direction of the effects are

ecosystem-dependent. Temperate grasslands, croplands, and tropical/subtropical forests are

sensitive to most climate extremes under warming. Tundra showed the strongest response to a

series of climate extremes compared to other ecosystems. Studies focusing on the impacts of

climate extremes on tundra SOC and carbon-rich soils (e.g., wetland soil) are urgently needed.

Our results and mapping provide a benchmark for identifying biome-specific critical climate

extremes and sensitive areas of SOC changes as a result of climate extremes. These results

demonstrate that climate extremes play a key role in determining SOC balance under climate

change probably via modulating local water-energy interactions thereby carbon inputs and

SOC decomposition and stabilization processes, which must be considered in Earth system

models and related climate polices to gain a more reliable prediction for SOC dynamics in

response to climate change.

11

References

1 Seneviratne., S. I.

et al.

Changes in climate extremes and their impacts on the natural

physical environment.

In: Managing the Risks of Extreme Events and Disasters to Advance

Climate Change Adaptation [Field, C.B., V. Barros, T.F. Stocker, D. Qin, D.J. Dokken, K.L.

Ebi, M.D. Mastrandrea, K.J. Mach, G.-K. Plattner, S.K. Allen, M. Tignor, and P.M. Midgley

(eds.)]. A Special Report of Working Groups I and II of the Intergovernmental Panel on

Climate Change (IPCC). Cambridge University Press, Cambridge, UK, and New York, NY,

USA, pp. 109-230

, (2012).

2 Piao, S.

et al.

The impacts of climate extremes on the terrestrial carbon cycle: A review.

Sci.

China Earth Sci.

, 1-13, (2019).

3 Reichstein, M.

et al.

Climate extremes and the carbon cycle.

Nature

500

, 287-295, (2013).

4 Frank, D.

et al.

Effects of climate extremes on the terrestrial carbon cycle: concepts,

processes and potential future impacts.

Glob Chang Biol

21

, 2861-2880, (2015).

5 Gampe, D.

et al.

Increasing impact of warm droughts on northern ecosystem productivity

over recent decades.

Nat. Clim. Change

11

, 772-779, (2021).

6 Ciais, P.

et al.

Europe-wide reduction in primary productivity caused by the heat and

drought in 2003.

Nature

437

, 529-533, (2005).

7 Zscheischler, J., Mahecha, M. D., Harmeling, S. & Reichstein, M. Detection and attribution

of large spatiotemporal extreme events in Earth observation data.

Ecol Inform

15

, 66-73,

(2013).

8 Sippel, S.

et al.

Drought, heat, and the carbon cycle: a review.

Current Climate Change

Reports

4

, 266-286, (2018).

9 Friedlingstein, P.

et al.

Global Carbon Budget 2021.

Earth Syst. Sci. Data

14

, 1917-2005,

(2022).

10 Crowther, T. W.

et al.

Quantifying global soil carbon losses in response to warming.

Nature

540

, 104-108, (2016).

11 Xiao, L.

et al.

Younger carbon dominates global soil carbon efflux.

Glob Chang Biol

28

,

5587-5599, (2022).

12 Smith, M. D., Knapp, A. K. & Collins, S. L. A framework for assessing ecosystem

dynamics in response to chronic resource alterations induced by global change.

Ecology

90

,

3279-3289, (2009).

13 Sippel, S., Zscheischler, J. & Reichstein, M. Ecosystem impacts of climate extremes

crucially depend on the timing.

Proc. Natl. Acad. Sci. U.S.A.

113

, 5768-5770, (2016).

14 Lugato, E.

et al.

Soil erosion is unlikely to drive a future carbon sink in Europe.

Science

12

Advances

4

, (2018).

15 Van Oost, K.

et al.

The impact of agricultural soil erosion on the global carbon cycle.

Science

318

, 626-629, (2007).

16 Anderegg, W. R.

et al.

Pervasive drought legacies in forest ecosystems and their

implications for carbon cycle models.

Science

349

, 528-532, (2015).

17 Anjileli, H.

et al.

Extreme heat events heighten soil respiration.

Sci Rep

11

, 6632, (2021).

18 Bardgett, R. D. & Caruso, T. Soil microbial community responses to climate extremes:

resistance, resilience and transitions to alternative states.

Philos Trans R Soc Lond B Biol

Sci

375

, 20190112, (2020).

19 Deng, L.

et al.

Drought effects on soil carbon and nitrogen dynamics in global natural

ecosystems.

Earth Sci Rev

214

, (2021).

20 Sippel, S., Zscheischler, J. & Reichstein, M. Ecosystem impacts of climate extremes

crucially depend on the timing.

Proc Natl Acad Sci U S A

113

, 5768-5770, (2016).

21 Diffenbaugh, N. S.

et al.

Quantifying the influence of global warming on unprecedented

extreme climate events.

Proc. Natl. Acad. Sci. U.S.A.

114

, 4881-4886, (2017).

22 Fischer, E. M., Sippel, S. & Knutti, R. Increasing probability of record-shattering climate

extremes.

Nat. Clim. Change

11

, 689-695, (2021).

23 Wang, M.

et al.

Global soil profiles indicate depth-dependent soil carbon losses under a

warmer climate.

Nat. Commun.

13

, 5514, (2022).

24 Karl, T. R., Nicholls, N. & Ghazi, A. in

Weather and Climate Extremes: Changes,

Variations and a Perspective from the Insurance Industry

(eds Thomas R. Karl, Neville

Nicholls, & Anver Ghazi) 3-7 (Springer Netherlands, 1999).

25 Mistry, M. A High-Resolution Global Gridded Historical Dataset of Climate Extreme

Indices.

Data

4

, (2019).

26 Karl, T. R., Nicholls, N. & Ghazi, A. Clivar/GCOS/WMO Workshop on Indices and

Indicators for Climate Extremes Workshop Summary.

Weather. Clim. Extremes

42

, 3-7,

(1999).

27 Yeon-Hee Kim, S.-K. M., Xuebin Zhang, Jana Sillmann, Marit Sandstad. Evaluation of the

CMIP6 multi-model ensemble for climate extreme indices.

Weather. Clim. Extremes

29

,

100269, (2020).

28 Zscheischler, J.

et al.

Future climate risk from compound events.

Nat. Clim. Change

8

, 469-

477, (2018).

29 Bardgett, R. D. & Caruso, T. Soil microbial community responses to climate extremes:

resistance, resilience and transitions to alternative states.

Phil. Trans. R. Soc. B

375

,

13

20190112, (2020).

30 Bennett, A. C.

et al.

Resistance of African tropical forests to an extreme climate anomaly.

Proc. Natl. Acad. Sci. U.S.A.

118

, (2021).

31 Zhao, M. & Running, S. W. Drought-induced reduction in global terrestrial net primary

production from 2000 through 2009.

Science

329

, 940-943, (2010).

32 Tao, F.

et al.

Microbial carbon use efficiency promotes global soil carbon storage.

Nature

,

618

, 981–985 (2023).

33 Eckardt, N. A.

et al.

Climate change challenges, plant science solutions.

Plant Cell

35

, 24-

66, (2023).

34 Xue, K.

et al.

Tundra soil carbon is vulnerable to rapid microbial decomposition under

climate warming.

Nat. Clim. Change

6

, 595-600, (2016).

35 Plaza, C.

et al.

Direct observation of permafrost degradation and rapid soil carbon loss in

tundra.

Nat. Geosci.

12

, 627-631, (2019).

36 Liu, Q.

et al.

Drought-induced increase in tree mortality and corresponding decrease in the

carbon sink capacity of Canada's boreal forests from 1970 to 2020.

Glob Chang Biol

29

,

2274-2285, (2023).

37 Slessarev, E. W.

et al.

Initial soil organic carbon stocks govern changes in soil carbon:

Reality or artifact?

Glob Chang Biol

29

, 1239-1247, (2023).

38 Cotrufo, M. F., Ranalli, M. G., Haddix, M. L., Six, J. & Lugato, E. Soil carbon storage

informed by particulate and mineral-associated organic matter.

Nat. Geosci.

12

, 989-994,

(2019).

39 Lugato, E., Lavallee, J. M., Haddix, M. L., Panagos, P. & Cotrufo, M. F. Different climate

sensitivity of particulate and mineral-associated soil organic matter.

Nat. Geosci.

14

, 295-

300, (2021).

40 Du, E.

et al.

Global patterns of terrestrial nitrogen and phosphorus limitation.

Nat. Geosci.

13

, 221-226, (2020).

41 Wang, H., Richardson, C. J. & Ho, M. Dual controls on carbon loss during drought

in peatlands.

Nat. Clim. Change

5

, 584-587, (2015).

42 Bao, T., Jia, G. & Xu, X. Weakening greenhouse gas sink of pristine wetlands under

warming.

Nat. Clim. Change

13

, 462-469, (2023).

43 Gauthier, S., Bernier, P., Kuuluvainen, T., Shvidenko, A. Z. & Schepaschenko, D. G. Boreal

forest health and global change.

Science

349

, 819-822, (2015).

44 Wang, J., Taylor, A. R. & D'Orangeville, L. Warming-induced tree growth may help offset

increasing disturbance across the Canadian boreal forest.

Proc. Natl. Acad. Sci. U.S.A.

120

,

14

e2212780120, (2023).

45 Jarvis, P. & Linder, S. Constraints to growth of boreal forests.

Nature

405

, 904-905, (2000).

46 Lim, H.

et al.

Boreal forest biomass accumulation is not increased by two decades of soil

warming.

Nat. Clim. Change

9

, 49-52, (2018).

47 Wang, P., Huang, K. & Hu, S. Distinct fine-root responses to precipitation changes in

herbaceous and woody plants: a meta-analysis.

New Phytol

225

, 1491-1499, (2020).

48 Vereecken, H.

et al.

Soil hydrology in the Earth system.

Nat Rev Earth Environ

3

, 573-587,

(2022).

49 García-Palacios, P.

et al.

Evidence for large microbial-mediated losses of soil carbon under

anthropogenic warming.

Nat Rev Earth Environ

2

, 507-517, (2021).

50 Ridder, N. N.

et al.

Global hotspots for the occurrence of compound events.

Nat. Commun.

11

, 5956, (2020).

51 Xu, L., Wang, A., Wang, D. & Wang, H. Hot Spots of Climate Extremes in the Future.

J

Geophys Res-Atmos

124

, 3035-3049, (2019).

52 Tschumi, E. & Zscheischler, J. Countrywide climate features during recorded climate-

related disasters.

Clim Change

158

, 593-609, (2020).

Methods

Data sources

We used the same global dataset of SOC stock (Mg C ha

–1

) as ref

. The dataset includes

110,695 SOC measurements in the 0–30 cm soil layer depth compiled into the WoSIS snapshot

of 2019, which was collated and managed by the World Soil Information Service

53

. The dataset

includes the snapshot of 2019. Another 2,778 measurements for the 0–30 cm soil layer from

permafrost-affected regions

54

were also included in the dataset. In brief, the data were quality-

assessed and standardized, using consistent procedures. The dataset covers nine major biome

groups (Fig. S1) under various climate conditions across the globe. At the locations of SOC

measurements, mean annual temperature (MAT) ranges from –20.0 to 30.7 °C and mean annual

precipitation (MAP) from 0 to 6,674 mm.

MAT and MAP at each soil location were extracted from WorldClim version 2 (ref.

55

),

which is a global mapping product at the resolution of 0.0083°×0.0083° using monthly

temperature and precipitation recorded for the period 1970-2000. Soil profiles in the same

0.0083° grid (i.e., ~1 km

2

) share the same MAT and MAP. The WWF (World Wildlife Fund)

15

map of terrestrial ecoregions of the world

56

was used to extract the biome type at each soil

profile. The MODIS land cover map

57

at the same resolution of WorldClim was used to identify

that if the soil profile is located in a cultivated land (i.e., cropland or cropland/natural vegetation

mosaic). A global landform spatial layer was obtained from Global Landform Classification

58

(https://esdac.jrc.ec.europa.eu/content/global-landform-classification) to identify the landform

at each soil location. Global terrestrial lands were divided into three general landform types:

plains including lowlands, plateaus, and mountains including hills. A global spatial layer of the

12 USDA soil orders (https://www.nrcs.usda.gov/resources/guides-and-instructions/soil-

taxonomy) were used to distinguish the soil type for each soil profile.

Climate extremes

At the location of each soil profile, 25 climate extreme indices (CEIs) defined by the expert

team on Climate Change Detection and Indices

59

were extracted from a high resolution

historical dataset of climate extreme indices

(that is a gridded mapping product at the

resolution of 0.25°×0.25°, which is equivalent to ~ 25 km at the equator). In brief, the average

CEIs were calculated based on sub-daily temperature and precipitation observations during the

period 1970–2016 using the climate data managed by the Global Land Data Assimilation

System

60

, representing baseline climate extreme conditions such as the intensity and frequency

of heatwaves, drought, and heavy precipitation. Soil profiles locate in the same 0.25°×0.25°

grid (i.e., ~25 km × ~25 km) share the same CEIs. A detailed description of CEIs is listed in

Table S1.

In addition, another eight indices indicating the magnitude and frequency of heatwave,

cold waves, extreme dry and wet events (Table S1) were also obtained from this global mapping

product of CEIs (Fig. S2)

. A heatwave event is defined as consecutive three or more days

with a positive value of the excess heat factor

61

. The magnitude and frequency of heatwaves

(HWM and HWF) are the average temperature across all individual heatwave events and the

total number of days that contribute to all heatwave events every year, respectively. The cold

wave event is identified by the excess cold factor

62

with a negative value for consecutive three

or more days. The mean magnitude and frequency of cold wave (CWM and CWF) is the mean

temperature and number of days of all cold wave events, respectively. The standardized

precipitation index (SPI) at the temporal resolution of 3-month

in this global historical dataset

of CEIs

was used to identify extreme dry and wet events. An extreme dry event (i.e., drought)

16

was identified if SPI

≤

–1, and an extreme wet event was identified if SPI > 1 (ref.

64,65

). Then,

the intensities of every extreme dry and wet events were calculated as the average of the

difference between SPI in event duration and event thresholds. The average magnitude of

extreme dry and wet event (EDM and EWM) is the average value of intensity of individual

events during the period 1970–2016. Their average frequency was calculated as the percentage

of total extreme dry and wet (EDF and EWF) months in total months during the period 1970–

2016.

Space-for-time substitution-based grouping

Taking advantage of the big dataset of more than 100 thousand soil profiles, we adopted an

approach combining space-for-time substitution with meta-analysis techniques

to estimate

the role of climate extremes in regulating the response of SOC to warming. This approach

controls the effects of other confounding variables such as soil type and landform as well as

precipitation and its seasonality, and conducts assessment based on grouped data distinguished

by the intensity and/or frequency of climate extremes of interest.

The approach has been explicitly described in ref.

, and applied to estimate the response

of whole-soil profile SOC dynamics to warming without considering climate extremes

. In

this study, we quantify the influences of climate extremes on SOC dynamics under climate

warming by comparing SOC changes under two general climate scenarios (Fig. 1): warming

alone (W), and warming with changes in climate extremes represented by the 33 CEIs (W+E).

For each CEI, the estimated response under W+E was compared to that under W. In addition,

since there are no consistent estimations on the magnitude of the changes of climate extremes

1

,

an potential change gradient with five levels is assigned to each CEI (ΔCEI, Table S1). For this

reason, there is a total of 165 comparisons (i.e., 33 CEIs × 5 change levels). For each CEI, in

detail, the estimation starts from a soil profile randomly selected from the 113,013 soil profiles.

This soil profile is used as a criterion to select a group of soils sharing the same MAT, CEI,

MAP, precipitation seasonality, landform and soil type. This group is called ambient group (A

group).

Then, another two groups (i.e., W and W+E groups) are selected according to the same

selection criterion as the A group, except that MAT in both W and W+E groups must be certain

degrees higher (1.5 °C in this study, which is the global climate target

66

) than that in the A

group to reflect warming, and CEI in the W+E group must be certain magnitude higher or lower

17

(i.e., the five change levels of CEI, Table S1) than CEI in the A group to reflect shift of extreme

climate. During the selection process, MAT, CEI, and MAP of soil profiles (due to they are

continuous variables) in the same group cannot be exactly the same, and tolerance values (τ,

which defines the allowed maximum deviation of the variable of interest in the same group) of

0.5 °C and 25 mm were assigned for MAT and MAP, respectively. To make the estimation more

robust, these groups with less than 5 soil profiles is rejected. For each CEI, an empirical

tolerance value has also been assigned by considering its distribution across the globe (Table

S1). These tolerance values are different among the CEIs. By conducting this selection and

grouping, we obtained triples including A, W and W+E groups. With the replacement of the

selected profile within the triple, the same random sampling, selection and grouping procedure

have been applied until the number of triples are large enough to obtain robust estimation of

response ratio (Fig. S14). At last, we obtain a set of triples with the same warming and climate

extreme change levels under different environmental conditions (i.e., ambient MAT and CEI,

MAP, precipitation seasonality, soil type and landform).

Compound extreme events (i.e., concurrent occurrence of multiple climate extremes) may

be also more common in the future

28,67

. To estimate the response of SOC to compound climate

extremes, several CEIs were simultaneously used to conduct the selection and grouping

procedure to represent the occurrence of multiple evets. Unfortunately, the data does not allow

us to quantify the effects of all possible compound events as more criteria have to be applied

to the selection of triples for comparison, which will substantially shrink the number of triples,

resulting in unreliable estimation. Here, we only quantified the influence of shifts in two

common compound events involving two individual events: heatwave + extreme dry, and

heatwave + extreme wet. The data only allowed us to obtain enough triples for the first two

levels of changes of the two events. For this assessment, the similar sampling and selection

processes for single event have been applied except that both CEIs were considered.

Estimation of the response of SOC to warming and climate extreme changes

Meta-analysis techniques were used to estimate the response of SOC to W and W+E by

comparing SOC stock in the A group to that in the W and W+E groups, respectively. For

i

th

selected triple (A, W and W+E), the logarithmic response ratios of SOC (

lnR

i

) to W and W+E

were calculated, respectively, as:

𝑙𝑛𝑅

𝑊,𝑖

= ln (

𝑆𝑂𝐶

𝑊,𝑖

𝑆𝑂𝐶

𝐴,𝑖

), (1)

18

𝑙𝑛𝑅

𝑊+𝐸,𝑖

= ln (

𝑆𝑂𝐶

𝑊+𝐸,𝑖

𝑆𝑂𝐶

𝐴,𝑖

), (2)

where

𝑆𝑂𝐶

𝐴,𝑖

,

𝑆𝑂𝐶

𝑊,𝑖

and

𝑆𝑂𝐶

𝑊+𝐸,𝑖

are the mean SOC stock in the

i

th

A, W and W+E groups,

respectively. The two response ratios enable us to estimate the response attributed to climate

extreme changes (ln

R

E,i

) as:

𝑙𝑛𝑅

𝐸,𝑖

= ln 𝑅

𝑊+𝐸,𝑖

− ln 𝑅

𝑊,𝑖

= ln (

𝑆𝑂𝐶

𝑊+𝐸,𝑖

𝑆𝑂𝐶

𝑊,𝑖

). (3)

Then, the individual ln

R

E,i

values calculated for each triple were used to estimate a global

mean response ratio (

𝑙𝑛𝑅

𝐸

) by a weighted mixed-effects model using the

rma.av

function in

package

metafor

in R-4.1.1. The weight (

w

) was calculated as the inverse of the sum of within-

(

v

) and between-group (

τ

2

) variances in the same set of triples:

𝑙𝑛𝑅

𝐸

=

∑(𝑙𝑛𝑅

𝐸,𝑖

× 𝑤

𝑖

)

∑ 𝑤

𝑖

, (4)

where

𝑤

𝑖

=

1

𝑣

𝑖

+𝜏

2

is the weight for the

lnR

E,i

in the

i

th

triple. The within-group variance (

v

) was

calculated as:

𝑣

𝑖

=

𝑆

𝑒

2

𝑁

𝑒

𝑋

𝑒

2

+

𝑆

𝑐

2

𝑁

𝑐

𝑋

𝑐

2

, (5)

where

S

e

and

S

c

are the standard deviation of SOC in the A and W (or W+E) groups, respectively;

X

e

and

X

c

are the mean of SOC in the A and W (or W+E) groups, respectively;

N

e

and

N

c

are

the number of samples in the A and W (or W+E) groups, respectively. The between-group

variance (

τ

2

) was estimated by restricted maximum likelihood (REML). To assist interpretation,

𝑙𝑛𝑅

𝐸

were back-transformed and reported as percentage change, i.e.,

(𝑒

𝑙𝑛𝑅

𝐸

− 1) ×

100%. The

back-transformed values were also used for subsequent data analyses.

Especially, the percentage change calculated from

𝑙𝑛𝑅

𝐸

was referred as the additional

global percentage change induced by climate extremes under W+E (ΔSOC

E

), i.e.,

∆𝑆𝑂𝐶

𝐸

= (𝑒

𝑙𝑛𝑅

𝐸

− 1) × 100%

. A linear mixed-effects model treating ecosystem type as a

random factor was fitted to assess the relationship between

∆𝑆𝑂𝐶

𝐸

and change levels of climate

extremes using the

lmer

function in the lmerTest package in R-4.1.1.

Here we note that an implicit assumption underlying the space-for-time substitution

19

approach is that critical land surface processes, which substantially change the state or

succession direction of studied groups (e.g., volcano disruption-induced new soil formation

and/or soil cultivation in one group but not in another group), are independent of space and

time

68

. That is, soils in any typical triple were not and will not be differently influenced by

other critical processes significantly regulating SOC dynamics. In addition, the approach

assumes that soils in the three groups of the same triple are at the same average succession

stage (e.g., steady state or similar succession direction). Croplands are subjected to substantial

management interventions (e.g., tillage, fertilizer application, irrigation), and these

management varies widely across space and over time. As such, croplands would not meet the

assumption of steady state or similar succession direction, providing an opportunity to test the

sensitivity of these assumptions on the results. To do so, we repeated all assessment by

excluding soil profiles from croplands and compared the results with that including cropland

soils (Fig. S13). In addition, soils sampled at different time may be also at different stages, we

repeated all assessment by excluding data measured before 1970 (Fig. S15). In addition, we

conducted a sensitivity test using a leave-out-out approach that sequentially removes one triple

to calculate the grand means using the remaining data. For all CEIs, the results show that the

grand means are not sensitive to any single triple (Fig. S16).

Driver analysis and global mapping

Focusing on ΔSOC

E

, we explained its variance by fitting a machine learning-based statistical

model - random forest (RF) driven by 67 potential environmental variables (Table S2)

including mean and extreme climate conditions, soil properties, vegetation properties, and

terrain attributes. The RF model has an advantage of taking into account nonlinear and

interactive relationships among predictors, providing a robust estimation of overall and

individual importance of the assessed predictors. Before fitting the model, collinearity among

the predictors is identified by variance inflation factor

(VIF) using the

vif

function in the car

package in R 4.1.1, and predictors with VIF of greater than 10 (28 out of the 67 variables) were

removed from the model. Model hyperparameters (mtry, randomly selected predictors; splitrule,

splitting rule; min.node.size, minimal node size) were optimized by using a tenfold cross-

validated approach using the

train

function in the caret package in R-4.1.1. That is, 80% data

were used to trained RF model, and the remaining 20% data to validate the model (Fig. S6).

The best RF model was selected with the highest coefficient of determination (R

2

).

The relative importance of each predictor was estimated by permutation variable

20

importance measurements using

varImp

function in package caret. Briefly, based on the times

a predictor selected for splitting when growing a tree, the mean square error for every given

regression tree with out-of-bag estimates is calculated. The change of model residuals due to

that splitting represents the influence of the predictor. To compare the relative importance of

individual predictors among CEIs, their influences were first normalized with the influence of

the most important predictor and then multiplied by the corresponding model performance (i.e.,

R

2

). In addition, a partial dependence analysis was used to explore the partial dependence of

ΔSOC

E

on CEI change gradients using the

partial

function in package

pdp

in R-4.1.1.

The RF model with the lowest RMSE for predicting the effects of climate extremes on

ΔSOC

E

was applied to the globe at the resolution of 0.25°×0.25° using gridded data of model

predictors (HSWD dataset

69

was used for baseline SOC level) to map the regulating effect of

selected climate extreme shifts on SOC responses under 1.5 °C warming. This mapping was

only demonstrated using the eight indices representing the magnitude and frequency of extreme

events with their median change levels (Table S1). That is, the magnitude (HWM) and

frequency of heat waves (HWF) are increased by 3 °C and 6 days, respectively; the magnitude

and frequency of both extreme dry (EDM) and wet (EWM) are increased 0.15 unit and 6%,

respectively; the magnitude and frequency of cold wave (CWM and CWF) are decreased by

3 °C and 6 days, respectively. The uncertainty of the response was estimated as the standard

error of individual predictions of 500 trees in the RF model (Figs. S6, S10 and S11).

Data availability

The

33

climate

extreme

indices

can

be

accessed

at

https://doi.pangaea.de/10.1594/PANGAEA.898014. Global mapping products generated in

this

study

are

publicly

available

and

deposited

to

https://doi.org/10.6084/m9.figshare.22317202. Other data used in this study is the same to that

used in ref.

https://doi.org/10.6084/m9.figshare.22317202.

, which are publicly accessible.

Code availability

Code (R scripts) used to assess the data and generate the results are deposited to

References

53 Batjes, N. H., Ribeiro, E. & van Oostrum, A. Standardised soil profile data to support global

mapping and modelling (WoSIS snapshot 2019).

Earth Syst. Sci. Data

12

, 299-320, (2020).

21

54 Mishra, U.

et al.

Spatial heterogeneity and environmental predictors of permafrost region

soil organic carbon stocks.

Sci Adv

7

, (2021).

55 Fick, S. E. & Hijmans, R. J. WorldClim 2: new 1-km spatial resolution climate surfaces for

global land areas.

Int J Climatol

37

, 4302-4315, (2017).

56 Olson, D. M.

et al.

Terrestrial Ecoregions of the World: A New Map of Life on Earth: A

new global map of terrestrial ecoregions provides an innovative tool for conserving

biodiversity.

Bioscience

51

, 933-938, (2001).

57 Channan, S., Collins, K. & Emanuel, W. Global mosaics of the standard MODIS land cover

type data.

University of Maryland and the Pacific Northwest National Laboratory, College

Park, Maryland

, (2014).

58 Meybeck, M., Green, P. & Vörösmarty, C. A New Typology for Mountains and Other Relief

Classes: An Application to Global Continental Water Resources and Population

Distribution.

Mt Res Dev

21

, 34-45, (2001).

59 Alexander, L. V. Global observed long-term changes in temperature and precipitation

extremes: A review of progress and limitations in IPCC assessments and beyond.

Weather.

Clim. Extremes

11

, 4-16, (2016).

60 Rodell, M.

et al.

The Global Land Data Assimilation System.

bams

85

, 381-394, (2004).

61 Alexander, L. V. & Perkins, S. E. On the Measurement of Heat Waves.

J. Clim.

26

, 4500-

4517, (2013).

62 Sheridan, S. C. & Lee, C. C. Temporal Trends in Absolute and Relative Extreme

Temperature Events Across North America.

J Geophys Res-Atmos

123

, (2018).

63 Vicente-Serrano, S. M., Beguería, S. & López-Moreno, J. I. A Multiscalar Drought Index

Sensitive to Global Warming: The Standardized Precipitation Evapotranspiration Index.

J.

Clim.

23

, 1696-1718, (2010).

64 Schwalm, C. R.

et al.

Global patterns of drought recovery.

Nature

548

, 202-205, (2017).

65 Chen, H., Wang, S., Zhu, J. & Zhang, B. Projected Changes in Abrupt Shifts Between Dry

and Wet Extremes Over China Through an Ensemble of Regional Climate Model

Simulations.

J Geophys Res-Atmos

125

, (2020).

66 IPCC

.

Global warming of 1.5°C. An IPCC Special Report on the impacts of global warming

of 1.5°C above pre-industrial levels and related global greenhouse gas emission pathways,

in the context of strengthening the global response to the threat of climate change,

sustainable development, and efforts to eradicate povertyy [Masson-Delmotte, V., P. Zhai,

H.-O. Pörtner, D. Roberts, J. Skea, P.R. Shukla, A. Pirani, W. Moufouma-Okia, C. Péan, R.

Pidcock, S. Connors, J.B.R. Matthews, Y. Chen, X. Zhou, M.I. Gomis, E. Lonnoy, T.

22

Maycock, M. Tignor, and T. Waterfield (eds.)].

Cambridge University Press, Cambridge,

UK and New York, NY, USA, 616 pp.

(2018).

67 Bevacqua, E., Zappa, G., Lehner, F. & Zscheischler, J. Precipitation trends determine future

occurrences of compound hot–dry events.

Nat. Clim. Change

12

, 350-355, (2022).

68 Strayer, D.

et al.

Long-Term Ecological Studies: an Illustrated Account of Their Design,

Operation, and Importance to Ecology.

Adv. Ecol. Acad. Sinica

5

, 51-55, (1988).

69 FAO/IIASA/ISRIC/ISS-CAS/JRC. Harmonized World Soil Database (version 1.2). FAO,

Rome, Italy and IIASA, Laxenburg, Austria., (2012).

Acknowledgement

This research has been financially supported by the National Natural Science Foundation of

China (grant no. 32241036, 32171639) and the National Key Research Program of the Ministry

of Science and Technology of China (grant no. 2021YFE0114500). Contributions of U. Mishra

were supported through a U.S. Department of Energy grant to the Sandia National Laboratories,

which is a multi-mission laboratory managed and operated by National Technology and

Engineering Solutions of Sandia, LLC, a wholly owned subsidiary of Honeywell International,

Inc., for the U.S. Department of Energy's National Nuclear Security Administration under

contract DE-NA-0003525.

Author contributions

Z.L. conceived the study; M.W. and Z.L. led data assessment; Y.Y. and U.M. contributed to

permafrost data; Z.L. and M.W. interpreted the results with the contribution of all authors; Z.L.

and M.W. led the writing of the manuscript, and all authors improved the manuscript.

Competing interests

The authors declare no competing interests.

23

Fig. 1. Schematic representation of the approach used to quantify the response of soil

organic carbon (SOC) to climate extreme and warming scenarios.

Each dot represents one

SOC observation. Dots in the

i

th

plate share the same mean annual temperature (

MAT

A,

i

) and

the frequency or intensity of the

m

th

climate extreme (

CEI

A

,

m, i

), while dots in the same color

plates indicate a group sharing the same mean annual precipitation (

MAP

), precipitation

seasonality, landform, and soil type.

ΔCEI

m

(shown in Table S1) and

ΔT

(1.5 °C in this study)

are changes in

CEI

A

, m

and

MAT

A

of interest, respectively. SOC values in the warming (W),

and warming plus more extreme (W+E) plates are compared with the values in the ambient (A)

plates in the same group to calculate a weighted average effect size (i.e., the response of SOC

to

ΔT

, and

ΔT

+

ΔCEI

m

) by the inverse of the sum of within- and between-group variances,

using meta-analytic techniques. The response of SOC to W+E attributed to climate extremes is

calculated by the comparison between W and W+E plates with meta-analysis (Methods). This

figure is modified from Wang et al.

24

Fig. 2. Global responses of soil organic carbon to climate change scenarios. a

, responses to

eight integrated climate extremes. HWM and HWF, the magnitude and frequency of heat waves,

respectively; CWM and CWF, the magnitude and frequency of cold waves, respectively; EWM

and EWF, the magnitude and frequency of extreme wet, respectively; EDM and EDF, the

magnitude and frequency of extreme dry, respectively. Numbers in the parentheses show the

five change levels.

b

, responses to compound climate extremes. Two levels of climate extreme

changes with their selected combinations (due to data limitation) are assessed. Different letters

in the same group of comparisons indicate significant difference between climate change

scenarios. A symbol of “+” indicates concurrent occurrence of the relevant changes. Error bars

show 95% confidence interval. See Table S1 for the detailed definition of climate extremes.

25

Fig. 3. The dependence of additional changes in soil organic carbon (

Δ

SOC

E

) on change

levels of climate extremes under warming.

ΔSOC

E

is defined as the difference between

percentage responses of SOC to W+E (1.5 °C warming plus more extreme climate) and that to

W (1.5 °C warming), which can also be explained as SOC changes attributed to climate

extremes under W+E.

↑

,

↓

and ↕ besides each CEIs represent that the CEI is increased,

decreased, or either under warming, respectively.

a

, global average dependence. Error bars

show one standard error. **,

P

< 0.01; *,

P

< 0.05.

b

, ecosystem-specific dependence.

Grey

grids indicate that the estimated dependence is statistically insignificant (

P

> 0.05). TS forests,

tropical/subtropical forests; Med/mon shrublands, Mediterranean/montane shrublands; TS

grasslands/savannas, tropical/subtropical grasslands/savannas.

Detailed descriptions of CEIs

are shown in Table S1.

26

Fig. 4. Influences of environmental variables on soil organic carbon (SOC) changes

attributed to climate extreme shifts (i.e.,

Δ

SOC

E

) under warming. a,

the relative influences

are estimated using a machine learning-based statistic model (i.e., random forest) taking into

account non-linear relationships. Numbers above the bars show the determination coefficients

of the model for explaining the corresponding CEI.

b-d

, partial dependence of ΔSOC

E

on SOC

stock (

b

), mean annual temperature (

c

) and precipitation (

d

). Detailed descriptions of

individual CEIs and predictor variables are shown in Table S1 and S2, respectively.

27

Fig. 5. Global spatial pattern of changes in soil organic carbon attributed to climate

extreme shifts (i.e.,

Δ

SOC

E

) under a warmer and more extreme climate.

ΔSOC

E

is defined

as the difference between percentage responses of SOC to W+E and that to W, which can also

be explained as the additional changes in SOC induced by climate extremes.

28

Table S1.

Climate extreme indices used in this study. The definition of climate extremes is adopted from ref.

and recommended by World

Meteorological Organization. TX, daily maximum near-surface air temperature; TN, daily minimum near-surface air temperature; TM, daily mean

mean-surface temperature; PR, daily near-surface total precipitation; SPI, standardized precipitation index.

Climate

extreme types

CEI

Definition

Unit

Global range

(mean, min, max)

ΔCEI

CEI

tolerance

Extreme heat

TXx

Max TX: warmest daily TX

°C

(35.52, 2.15, 49.73)

1, 2, 3, 4, 5

0.5

TXn

Min TX: coldest daily TX

°C

(4.05, -52.30, 30.14)

1, 2, 3, 4, 5

0.5

TNx

Max TN: warmest daily TN

°C

(21.90, -4.39, 33.59)

1, 2, 3, 4, 5

0.5

TNn

Min TN: coldest daily TN

°C

(-8.25, -60.42, 23.28)

1, 2, 3, 4, 5

0.5

TX90p

Amount of hot days: percentage of days when

TX > 90th percentile

%

(11.83, 7.97, 23.72)

2, 4, 6, 8, 10

1

TN90p

Amount of cold days: percentage of days

when TN > 90

th

percentile

%

(13.29, 8.22, 31.08)

2, 4, 6, 8, 10

1

TR

Number of tropical nights: number of days

when TN > 20 °C

days

(68.44, 0.00, 365.19) 5, 10, 15, 20, 25

2.5

GSL

Growing season length: annual number of

days between the first occurrence of 6

consecutive days with TM > 5 °C and the

first occurrence of 6 consecutive days with

TM < 5 °C

days

(300.93, 0.00, 365.26) 5, 10, 15,20, 25

2.5

DTR

Daily temperature range: mean difference

between daily TX and daily TN

°C

(11.50, 1.06, 19.10)

-2, -1, 0, 1, 2

0.5

WSDI

Warm spell duration index: annual number of

days contributing to events where 6 or more

consecutive days experience TX > 90

th

percentile

days

(5.08, 0.00, 50.53)

1, 2, 3, 4, 5

0.5

SU

Number of summer days: number of days

when TX > 25 °C

days

(168.03, 0.00, 365.19) 5, 10, 15, 20, 25

2.5

HWM

Heatwave magnitude: mean temperature of

all heatwaves identified by excess heat factor

°C

2

(4.03, 0.15, 17.69)

0.5, 1, 1.5, 2,

2.5

0.5

29

(EHF). See ref.

61

for more details

HWF

Heatwave frequency: number of days that

contribute to heatwaves as identified by

excess heat factor (EHF). See ref.

61

for more

details.

days

(10.11, 2.28, 36.98)

2, 4, 6, 8, 10

1

Extreme cold

TX10p

Amount of cool days: percentage of days

when TX < 10th percentile

%

(11.17, 7.49, 26.16)

-2, -4, -6, -8, -

10

1

TN10p

Amount of cool nights: percentage of days

when TN < 10

th

percentile

%

(9.94, 7.02, 21.42)

-2, -4, -6, -8, -

10

1

FD

Number of frost days: number of days when

TN < 0 °C

days

(66.26, 0.00, 365.09)

-2, -4, -6, -8, -

10

1

ID

Number of icing days: number of days when

TX < 0 °C

days

(23.26, 0.00, 335.06)

-2, -4, -6, -8, -

10

1

CSDI

Cold spell duration index: annual count of

days with at least 6 consecutive days when

daily minimum temperature is less than the

10th percentile.

days

(2.29, 0.00, 51.23)

-1, -2, -3, -4, -5

0.5

CWM

Cold wave magnitude: average temperature

across all individual cold waves identified by

excess cold factor (ECF)

62

-°C

2

(-10.44, -63.31, -0.15)

0.5, 1, 1.5, 2,

2.5

0.25

CWF

Cold wave frequency: total number of days

that contribute to all individual cold waves

identified by excess cold factor (ECF)

62

days

(12.29, 2.47, 51.13)

-2, -4, -6, -8, -

10

-1

Extreme wet

R99p

Total annual PR from very heavy

precipitation days: annual sum of daily PR >

99th percentile

mm

(663.55, 0.00,

6080.61)

50, 100, 150,

200, 250

25

R95p

Total annual PR from heavy precipitation

days: annual sum of daily PR > 95

th

percentile

mm

(1993.70, 1.20,

9091.31)

100, 200, 300,

400, 500

50

Rx5day

Maximum 5-day PR: monthly 5-day PR total

precipitation

mm

(675.78, 6.28,

2383.89)

50, 100, 150,

200, 250

25

Rx1day

Maximum 1-day PR: maximum 1-day PR

mm

(367.86, 5.40,

20, 40, 60, 80,

10

30

total

1139.58)

100

CWD

Consecutive wet days: maximum number of
consecutive wet days (when PR >= 1.0 mm)

days

(21.39, 0.77, 194.34)

2, 4, 6, 8, 10

1

SDII

Simple precipitation intensity index: annual

total PR divided by the number of wet days

(when total PR >= 1.0 mm)

mm/day

(41.25, 3.75, 125.04) 5, 10, 15, 20, 25

2.5

R20mm

Number of very heavy precipitation days:

number of days when PR >= 20 mm

days

(75.63, 0.00, 306.15)

2, 4, 6, 8, 10

1

R10mm

Number of heavy precipitation days: number

of days when PR >= 10 mm

days

(100.62, 0.30, 332.79) 5, 10, 15, 20, 25

2.5

EWM

Wet magnitude: mean SPI value of all

extreme wet events. See ref.

for more

details.

-

(0.46, 0.25, 1.44)

0.05, 0.1, 0.15,

0.2, 0.25

0.02

EWF

Wet frequency: percentage of extreme wet

months on all months

%

(17.67, 8.89, 34.97)

2, 4, 6, 8, 10

1

Extreme dry

CDD

Consecutive dry days: maximum length of

dry spell: maximum number of consecutive

days with daily precipitation less than 1 mm

days

(36.05, 1.00, 250.00) 5, 10, 15, 20, 25

2.5

EDM

Extreme dry magnitude: mean SPI value of

all extreme dry events

-

(-0.38, -1.19, -0.22)

-0.05, -0.1, -

0.15, -0.2, -0.25

-0.025

EDF

Extreme dry frequency: percentage of

extreme dry months on all months

%

(17.10, 9.68, 32.07)

2, 4, 6, 8, 10

1

31

Table S2.

Environmental variables used in this study.

Classes

Predictors

Description

Unit

Soil

SOCs

Soil organic carbon stock

Mg ha

-1

Soil types

Soil order in USDA soil classification system

Terrain

Elevation

m

Landform

Terrain landform

Vegetation

NPP

Vegetation net primary productivity

g C m

-2

Biomes

Vegetation types

Mean

climate

ai_et0

Arid index

Bio1

Annual Mean Temperature

°C

Bio2

Mean Diurnal Range (Mean of monthly (max temp - min temp))

°C

Bio3

Isothermality (BIO2/BIO7) (* 100)

°C

Bio4

Temperature Seasonality (standard deviation *100)

°C

Bio5

Max Temperature of Warmest Month

°C

Bio6

Min Temperature of Coldest Month

°C

Bio7

Temperature Annual Range (BIO5-BIO6)

°C

Bio8

Mean Temperature of Wettest Quarter

°C

Bio9

Mean Temperature of Driest Quarter

°C

Bio10

Mean Temperature of Warmest Quarter

°C

Bio11

Mean Temperature of Coldest Quarter

°C

Bio12

Annual Precipitation

mm

Bio13

Precipitation of Wettest Month

mm

Bio14

Precipitation of Driest Month

mm

Bio15

Precipitation Seasonality (Coefficient of Variation)

mm

Bio16

Precipitation of Wettest Quarter

mm

Bio17

Precipitation of Driest Quarter

mm

Bio18

Precipitation of Warmest Quarter

mm

Bio19

Precipitation of Coldest Quarter

mm

PCI

Precipitation concentration index

mm

SPC

Summer precipitation concentration index

32

WPC

Winter precipitation concentration index

PS

Precipitation seasonality

Extreme

climate

ΔCEI

climate extreme change levels

CEI

See Table S1

33

Table S3. Global absolute (mean ± se) soil organic carbon changes (i.e., ΔSOC, Pg) attributed to climate extreme shifts under warming

and more extreme climate.

HWM and HWF, the magnitude and frequency of heat waves, respectively; CWM and CWF, the magnitude and

frequency of cold waves, respectively; EWM and EWF, the magnitude and frequency of extreme wet, respectively; EDM and EDF, the magnitude

and frequency of extreme dry, respectively. T/ST, tropical/subtropical; Med/Mon, Mediterranean/Montane.

Biomes

SOC

stock (Pg,

0–30 cm)

Absolute SOC changes induced by climate extremes (Pg C) under climate extreme shifts scenarios

HWF

increase 6

days

HWM

increase 3

°

C

2

CWF

decrease 6 days

CWM

increase 3

°

C

2

EDF

increase 6

%

EDM

increase 0.15

unit

EWF

increase 6

%

EWM

increase 0.15

unit

T/ST forests

140.4

-48.6 ± 2

-36.0 ± 1.4

-19.9 ± 1.5

-22.3 ± 1.8

-32.5 ± 1.4

-50.4 ± 1.4

-45.5 ± 1.4

-47.9 ± 1.2

T/ST
grasslands/savannas

72.1

-11.7 ± 1.0

-14.6 ± 0.6

3.5 ± 0.8

12.1 ± 1.0

-6.5 ± 0.7

-6.5 ± 0.8

-7.9 ± 0.9

-10.9 ± 0.6

Temperate forests

51.8

-11.8 ± 0.6

-6.3 ± 0.3

-2.6 ± 0.4

-2.5 ± 0.5

-9.9 ± 0.5

-11.9 ± 0.4

-11.5 ± 0.4

-12.5 ± 0.5

Temperate grasslands

28.0

-2.9 ± 0.3

1.6 ± 0.2

0.8 ± 0.3

1.8 ± 0.3

-2.0 ± 0.3

-3.4 ± 0.2

-2.1 ± 0.3

-5.0 ± 0.3

Med/Mon shrublands

56.5

-6.7 ± 0.6

2.0 ± 0.5

-0.9 ± 0.5

5.5 ± 0.7

0.8 ± 0.6

-6.0 ± 0.5

-7.6 ± 0.5

Boreal forests

201.6

-34.4 ± 3.0

-39.1 ± 2.3

-14.0 ± 3.4

-14.0 ± 3.2

-57.7 ± 2.7

-10.9 ± 2.5

-35.4 ± 2.9

-41.8 ± 2.7

Tundra

60.6

-2.0 ± 1.3

1.6 ± 1.0

2.2 ± 1.3

5.0 ± 1.3

-2.2 ± 1.4

5.2 ± 1.3

-3.4 ± 1.2

-5.3 ± 1.5

Deserts

42.2

-3.3 ± 0.5

-0.5 ± 0.4

5.7 ± 0.5

9.5 ± 0.6

-1.4 ± 0.5

-3.5 ± 0.4

-0.3 ± 0.5

-4.8 ± 0.4

Croplands

97.2

-17.8 ± 1.1

-7.0 ± 0.6

-2.1 ± 0.8

0.4 ± 0.7

-13.6 ± 0.8

-18.6 ± 0.7

-17.7 ± 0.9

-23.9 ± 0.7

Globe

750.4

-139.0 ± 7.1

-98.2 ± 4.5

-27.1 ± 6.1

-4.5 ± 6.3

-125.0 ± 5.8

-105.9 ± 5.2

-129.9 ±

5.8

-159.7 ± 5.8

34

Fig. S1. Distribution of soil organic carbon data points used in this study.

The data is

derived from two data sources, and the number of data points for each biome is shown in the

parentheses. The figure is adopted from ref.

35

Fig. S2. The spatial pattern of eight climate extreme indices

. These indices represent the

magnitude and frequency of heat waves, cold waves, extreme dry and wet days. Data are

adopted from ref.