Meng et al. (2026). Nonlinear temperature change responses shape soil organic carbon loss-gain transitions in global Mollisol croplands

and play a crucial role in the long-term storage and release of

carbon globally

. The high soil organic carbon (SOC) stocks in Mol-

lisol are crucial for global food security. However, SOC stocks in Mol-
lisol show signi

ﬁ

cant spatial and temporal variation, with declining

trends mainly attributed to climate change and agricultural activities

Predictions of SOC stock dynamics in Mollisol vary signi

ﬁ

cantly due to

differences in temporal and spatial scales, technological factors, and
modeling approaches

11

, resulting in substantial uncertainties and

leaving their global dynamics insuf

ﬁ

ciently constrained

–

17

. Despite

numerous studies analyzing the causes of declining SOC stocks

18

–

20

, the

responses of SOC stocks to multiple drivers in global Mollisol remain
poorly understood. This limits our ability to generalize impact patterns
across different temporal and regional conditions, hindering efforts to
restore dynamic equilibrium in SOC stocks.

In recent decades, heterogeneous changes in global Mollisol

SOC stocks have created a pressing need for continuous, high-
accuracy spatial predictions. However, most large-scale predic-
tions rely on statistical relationships with environmental variables
and are constrained by existing SOC data, which suffer from
limited spatial resolution, incomplete spatiotemporal coverage,
and consequently inconsistent estimates of baseline stocks and
their key drivers

. Satellite remote sensing, through the con-

tinuous monitoring of soil spectral re

ﬂ

ectance

–

11

, offers a pow-

erful alternative to track SOC variations. This approach directly
addresses the critical limitations of sparse data points, coarse
spatial resolution

,

–

15

, and static modeling frameworks inherent

in traditional methods

21

–

23

. The capacity of satellite remote sen-

sing to provide long-term time series with greatly improved
spatial detail fundamentally enhances our ability to characterize

Received: 21 August 2025

Accepted: 14 May 2026

Check for updates

1

State Key Laboratory of Black Soils Conservation and Utilization, Northeast Institute of Geography and Agroecology, Chinese Academy of Sciences,

Changchun, China.

2

University of Chinese Academy of Sciences, Beijing, China.

3

Key Laboratory of Geographical Processes and Ecological Security in

Changbai Mountains, Ministry of Education, School of Geographical Sciences, Northeast Normal University, Changchun, China.

4

Center for Spatial Tech-

nologies and Remote Sensing (CSTARS), Department of Land, Air, and Water Resources, University of California, Davis, USA.

5

College of Information

Technology, Jilin Agricultural University, Changchun, China.

6

V.V. Dokuchaev Soil Science Institute, Moscow, Russian Federation.

e-mail:

byl1211neau@yeah.net

;

huanjunliu@yeah.net

| (2026) 17:6875

1

12

3

4

5

6

78

12

3

4

5

6

7

8

and quantify the dynamics of SOC stocks across global Mollisol
croplands.

Several factors, including temperature increases

19

, change in

precipitation

20

, cultivation duration

, and crop type

7

, have been

repeatedly identi

ﬁ

ed as important controls on cropland SOC stocks in

Mollisol regions. Recent syntheses and meta-analyses have examined
how SOC responds to multiple global-change drivers across climates
and soil types, suggesting that warming can amplify SOC responses to
elevated CO

2

, nitrogen addition, and altered moisture regimes, that

increases in plant productivity only partly compensate warming-
induced SOC losses, and that deeper, mineral-associated SOC pools
can also be vulnerable to sustained warming

25

27

. These studies

represent important progress by explicitly considering multiple dri-
vers, soil depth, and environmental context. Nevertheless, because
these studies typically aggregate data across soil orders, include rela-
tively few Mollisol cropland observations, and do not explicitly resolve
how SOC responses vary along gradients of warming rate at high
spatial resolution, they offer limited insight into whether, and at what
warming rates, SOC losses might be reversed in the global Mollisol
region.

Against this background, Mollisol regions provide a critical setting

for examining how soils and climate interact under intensive agri-
cultural use. Existing studies have shown that the mechanisms by
which climate change affects SOC stocks in Mollisol vary across cli-
mates and spatial-temporal scales

28

29

. For example, increased tem-

peratures may enhance microbial activity and accelerate SOC losses in
temperate climates

8

30

, whereas warming-induced soil freeze-thaw

cycles may, under some conditions, promote carbon accumulation in
cold regions

31

. We therefore hypothesize that soil-climate feedbacks in

Mollisol croplands are inherently nonlinear, and that suf

ﬁ

ciently

strong warming may, in some regimes, enhance plant productivity and
organic inputs enough to offset SOC losses. Such threshold-like
responses are likely to vary across latitudinal gradients, driven by
differences in baseline climate, SOC pools and ecosystem-level ther-
mal adaptation. In addition, increased nitrogen fertilizer use and
atmospheric deposition can promote plant growth and enhance SOC
stocks in the short term

32

, whereas water stress and elevated tem-

peratures may suppress carbon accumulation by inhibiting photo-
synthesis and altering microbial communities, ultimately reducing
SOC stocks

33

. The relative importance of climatic, biological, anthro-

pogenic and soil property drivers in Mollisol has therefore been shown
to vary from broad latitudinal zones to

ﬁ

ne-scale regional

microclimates

,

34

, but the existence, location and magnitude of

warming-rate thresholds, as well as the ways in which the dominant
drivers shift across these thresholds, remain poorly quanti

ﬁ

ed in

Mollisol regions.

To address these knowledge gaps, we integrate 5535 soil samples

(Supplementary Fig. 1a and Supplementary Table 1) with 1,086,202
Landsat images to construct a global database of surface (0

–

20 cm)

SOC stocks in Mollisol croplands at 30 m resolution for the period
1990

–

2020. Using a hybrid deep-learning framework that combines

ﬁ

eld observations with multi-decadal satellite re

ﬂ

ectance, we generate

four decadal SOC stock maps (1990, 2000, 2010, and 2020) and vali-
date their spatial and temporal dynamics through independent data-
sets and meta-analysis. We then analyze how changes in SOC stocks
relate to the rate of change in annual mean temperature (

Δ

Tavg) and to

key climatic, biological, anthropogenic, soil property, and topographic
factors. Speci

ﬁ

cally, we identify two distinct temperature change

thresholds in the

Δ

SOC-

Δ

Tavg relationship, namely a segmented

response threshold and an SOC gain-loss crossover threshold. Finally,
we use these thresholds to delineate regimes in which different drivers
dominate SOC changes, thereby elucidating the mechanisms by which
warming rate and co-varying factors jointly shape SOC dynamics in
Mollisol croplands. Our results challenge simple linear assumptions in
soil-climate feedbacks and provide an observation-based framework

for assessing the sustainability of SOC stocks under ongoing soil
degradation and declining agricultural soil quality.

Results

Spatiotemporal dynamics of SOC stocks in global Mollisol
croplands

Surface SOC stocks in croplands across the global Mollisol region have
been signi

ﬁ

cantly affected by global warming and intensive agri-

cultural activities over recent decades (Fig.

a and Supplementary

Fig. 1b, c). We used Landsat image,

ﬁ

eld-collected soil samples, and a

hybrid probability deep learning model (Supplementary Fig. 2) to
generate 30 m spatial resolution global estimates of Mollisol cropland
surface SOC stocks (root mean square error (RMSE), 17.83 MgC/ha;
coef

ﬁ

cient of determination (R

2

), 0.68) across multiple time periods

(Supplementary Fig. 3a) (Veri

ﬁ

cation results for the spatial and tem-

poral patterns of SOC stock changes are detailed in the Supplementary
Note 1) and the spatial distribution of uncertainty in SOC stock esti-
mations (Supplementary Fig. 4). From 1990 to 2020, global Mollisol
cropland surface SOC stocks declined by 5.59%, with average values
decreasing from 59.37 ± 21.30 MgC/ha to 56.05 ± 21.17 MgC/ha (Fig.

c

and Supplementary Table 2). The decadal decreases were 2.98%
(1990

–

2000),

1.35%

(2000

–

2010),

and

1.36%

(2010

–

2020),

respectively.

Signi

ﬁ

cant spatial heterogeneity is observed in the spatial dis-

tribution of SOC stocks in the global Mollisol (Fig.

a, the spatial dis-

tribution of SOC stocks in the global Mollisol for other periods are
provided in Supplementary Fig. 5). SOC stocks exhibited a clear lati-
tudinal gradient across global Mollisol croplands, with higher values at
higher latitudes and a decrease of 0.26 MgC/ha for each 1° decrease in
latitude toward the equator, highlighting the strong in

ﬂ

uence of

annual average temperature gradients on SOC stocks (Supplementary
Fig. 6). Among continents, Siberia (SI) and Europe (EU) had the highest
SOC stocks, with 71.44 ± 12.32 MgC/ha and 60.15 ± 21.15 MgC/ha,
respectively (Supplementary Table 2). Asian (AS) and North American
(NA) Mollisol had the second highest SOC stock, at 59.76 ± 16.68 MgC/
ha and 45.60 ± 16.72 MgC/ha, respectively. South American (SA) Mol-
lisol showed the lowest SOC stock, at 33.51 ± 8.17 MgC/ha.

Global Mollisol SOC stocks are still declining, but the rate of

decline has slowed from 0.18 MgC/ha/yr (1990

–

2000) to 0.08 MgC/

ha/yr (2010

–

2020) (Fig.

d). The change in global Mollisol SOC stocks

varied signi

ﬁ

cantly across continents (Fig.

b). The largest declines in

SOC stocks over the past 30 years occurred in the SI (a decline of 8.47
MgC/ha) and SA Mollisol regions (a decline of 7.69 MgC/ha). The AS (a
decline of 4.23 MgC/ha), NA (a decline of 3.46 MgC/ha), and EU (a
decline of 1.84 MgC/ha) Mollisol regions followed closely with
declines. From 1990 to 2020, SOC stocks decreased in 60.84% of the
global Mollisol region (Fig.

e), with declines concentrated in eastern

NA, SA, SI, and northeastern AS (Fig.

b). In the remaining 39.16% of the

area, SOC stocks increased, particularly in western NA and eastern EU.
The distribution of SOC stock changes indicates that new areas of
decline are mainly in the southern EU, southeastern NA, and the entire
SA Mollisol region (Supplementary Fig. 7). The proportion of regions
with declining SOC stocks from 1990 to 2020 varies across continents,
in addition to the general pattern described above. The area of
declining SOC stocks in the EU Mollisol region showed no signi

ﬁ

cant

change (Fig.

m). In contrast, the area of SOC stock decline is

increasing annually in other regions: the NA Mollisol region (up by
9.02%), SA Mollisol region (up by 22.54%), SI Mollisol region (up by
10.06%), and AS Mollisol region (up by 7.42%). However, the rate of
increase in declining SOC stock areas is slowing.

Climatic and anthropogenic controls on SOC stocks spatial
heterogeneity

We investigated the spatial heterogeneity of SOC stocks across global
Mollisol croplands over four time periods and assessed the in

ﬂ

uence

Article

https://doi.org/10.1038/s41467-026-73759-w

| (2026) 17:6875

2

of multiple environmental factors. These environmental factors were
selected because they are key determinants of SOC stocks and are
among the most widely used indicators in SOC-related studies (Sup-
plementary Table 3). A boosted regression tree (BRT) algorithm was
employed to assess the effects of climatic, anthropogenic, biological,
soil, and topographic factors on SOC stocks (Supplementary Fig. 8).
The BRT model explained 82% of the spatial heterogeneity in SOC
stocks (Fig.

a). Annual average temperature (Tavg) was the most

in

ﬂ

uential factor (53.26%), followed by nitrogen (N) (19.42%). Other

factors contributed less than 10%, including farm machinery data
(FMD) (5.73%), minimum normalized difference tillage index (min-
NDTI) (5.07%), elevation (ELE) (4.41%), annual precipitation (PR)
(1.89%), gross primary productivity (GPP) (1.75%), consecutive wet days
(CWD) (1.63%), consecutive dry days (CDD) (1.27%), and livestock
manure (1.11%). When examining the individual effects of each factor,
Tavg exhibited a strong negative effect on SOC stocks below 14 °C,
which plateaued above this threshold (Fig.

b). N below 59 g/kg had a

signi

ﬁ

cant positive effect on SOC stocks, which tended to saturate

above this threshold (Fig.

c). FMD showed a signi

ﬁ

cant negative effect

on SOC stocks below 2.2 Ps/1000 ha, which saturated beyond this level
(Fig.

d). minNDTI exhibited a signi

ﬁ

cant negative effect on SOC

stocks, which saturated at a threshold of 0.5 (Fig.

e). ELE below 160 m

had a signi

ﬁ

cant positive effect on SOC stocks, whereas between

160

–

450 m, it had a signi

ﬁ

cant negative effect (Fig.

f). PR exhibited a

signi

ﬁ

cant positive effect up to 55 mm, after which it saturated

(Fig.

g). Similar partial response patterns were observed across all

time periods (Supplementary Figs. 6

–

9). These partial response pat-

terns highlight the complex nature of SOC responses to key climatic
and management-related factors across global Mollisol croplands.

We compared the environmental factors in

ﬂ

uencing the spatial

heterogeneity of SOC stocks across different continents. Tavg was the
most important factor on all four continents, except in SA, where
nitrogen (N) was the dominant factor (Supplementary Fig. 8c). In all
cases, the factor importance exceeded 45%. The response of SOC
stocks to Tavg varied across regions. SOC stocks in NA and AS were the
most sensitive, showing the greatest decline with a 1 °C increase in
Tavg (Supplementary Figs. 13, 17). In contrast, SA was the least sensitive
(Supplementary Fig. 14).

Nonlinear response of SOC stock change to temperature
change rate

We analyzed temporal changes in cropland SOC stocks across global
Mollisol and evaluated the effects of environmental factors change
over the past three decades. Climate change was the dominant factor
of change in SOC (

Δ

SOC) (Fig.

b), with change rate in Tavg (

Δ

Tavg)

exerting the strongest in

ﬂ

uence (Fig.

c). Using a generalized additive

model (GAM), we identi

ﬁ

ed and corroborated a nonlinear relationship

between

Δ

SOC and

Δ

Tavg (Supplementary Fig. 18), characterized by a

The distribution of SOC stock
in 2020 (MgC/ha)

0–20

20–40

40–60

80–100

100–120

120–140

60–80

0

1,500

3,000 4,000 6,000

Km

The change of SOC stock
during 1990–2020 (MgC/ha)

-90–-60

-60–-30

0–30

30–60

60–90

-30–0

0

1,500

3,000 4,000 6,000

Km

1990

2000

20

C st

SOC st

ock

(MgC/h

(MgC/

a

)

Y

Year

ear

Y

2010 2020

20

0

40

60

80

100

120

0

140

1990 20

000

OC st

SOC

ock

(MgC/h

a)

Y

ear

Y

2010 2020

20

0

40

60

80

100

120

140

19

990 2000

00

SOC st

ock

oc

(MgC/ha)

Yea

ar

Y

20

2010 2020

2

20

0

4

40

60

80

100

120

140

1990 2000

SOC

s

t

S

ock

SOC st

S

ock

(MgC/h

h

C

a)

a

(MgC/h

MgC

a)

)

a

Year

Y

2010 2020

20

0

40

0

60

80

100

120

1

140

14

40

1990 2000

SOC stock

(MgC/ha)

Year

2010 2020

20

0

40

60

80

100

120

140

1990

2000

SOC stock

(MgC/ha)

Year

2010

2020

20

0

40

60

80

100

120

140

GM

NA

SA

EU

SI

AS

-5

-3

-4

0

-1

-2

2

1

5

4

3

Change rate of SOC stock (MgC/ha/y)

1990-2000

45.57%

41.76%

39.16%

54.43%

58.24%

60.84%

Area proportion (%)

0

20

40

60

80

100

1990-2010

1990-2020

1990–2

000

1990–2010

Area p

roport

ion (%

)

1990–2

2

0

020

20

0

65.37%

5.37

82.87%

2.8

2.87

87

1%

7.9

9

7

1

34.63%

4.63

17.13%

%

7.13

12.09%

2.09

40

60

80

100

1990–2

000

90–2

1990

010

Area p

rea p

roport

opo

n (%

ion (

)

1990–2020

20

2

0

55.21%

5.2

62.12%

2.1

64.23%

4.2

4.2 %

44.79%

4.7

4.79

37.88%

7.8

35.77%

5.7

5.7 %

40

4

60

8

80

100

0

1990–2

90–2

000

–2

1990–2

010

A

rea proportion

(%)

1990–2

19

020

20

0

40

0

60

0

80

0

100

47.23%

7.23

60.05%

0.05

47.54%

7.54

52

7%

7

2.7

77

7

2

7

39.95%

9.95

52.46%

6%

2.46

1990–2000 1990–2

010

Area p

p

roport

ropo

(%

ion (%

)

990

1990–2

020

20

0

40

0

60

80

8

100

1

69 55%

69.55%

9 5

9.5

9 5

9.5

74.17%

4.1

%

79.61%

9.6

9.6 %

30.45%

0.4

25.83%

5.8

20.39%

0.3

0.3 %

1990–2

199

000

1990–2

010

Area p

A

roport

ion (%)

1990–2

99

020

20

0

40

60

80

100

62.30%

%

2.30

74.67%

4.67

4.6

4.67

69.72%

9.72

37.70%

7.70

25.33%

5.33

33

5.33

3

30.28%

0.2

0.28

0

0.28

Positive change

Negative change

1990–2010

1990–2000

1990–2020

a

b

c

d

e

f

g

h

i

j

k

l

m

n

o

Fig. 1 | Spatial distribution of spatiotemporal changes in soil organic carbon
(SOC) stock in cropland topsoil layer in the global Mollisol (GM). a

Spatial

patterns of SOC stock in 2020.

b

Spatial patterns of changes in SOC stock during

1990

–

2020.

c

Boxplot of SOC stock in the GM in different years (the number of

pixels,

n

= 4,629,359,240).

d

Boxplot of SOC stock in different regions in different

years.

e

Area proportion with positive and negative changes SOC stock in different

periods.

f

–

j

Boxplots of SOC stock in the North America (NA,

n

= 1,436,314,361),

South America (SA,

n

= 355,062,490), Europe (EU,

n

= 1,745,314,562), Siberia (SI,

n

= 853,234,890), and Asian (AS,

n

= 239,432,937) Mollisol region in different years.

k

–

o

Area proportion with positive and negative changes SOC stock in the NA, SA,

EU, SI, and AS Mollisol region in different periods.

Article

https://doi.org/10.1038/s41467-026-73759-w

| (2026) 17:6875

3

segmented response (Fig.

a). When

Δ

Tavg was below 0.0086 °C/yr

(95% CI: 0.0077

–

0.0102 °C/yr),

Δ

SOC was negatively correlated with

Δ

Tavg (

r

=

−

0.48,

P

< 0.001) (Fig.

b), indicating that faster cooling was

associated with larger SOC gains. When

Δ

Tavg exceeded 0.0086 °C/yr,

the relationship became positive (

r

= 0.42,

P

< 0.001) (Fig.

b), imply-

ing that faster warming was associated with weaker SOC losses.
Notably, at

Δ

Tavg > 0.0309 °C/yr (95% CI: 0.0302

–

0.0316 °C/yr)

(Fig.

b),

Δ

SOC crossed from negative to positive, suggesting that

higher warming rates can shift SOC from loss to accumulation. We
term this threshold the SOC gain-loss crossover threshold. This
response was most pronounced in Mollisol regions of Southwest AS,
Northern SA, and Eastern EU (Fig.

a). We also detected a clear latitu-

dinal gradient in the gain-loss crossover, increasing from high to low
latitudes: 0.0261 °C/yr (95% CI: 0.0170

–

0.0352 °C/yr) for > 55°N,

0.0307 °C/yr (95% CI: 0.0299

–

0.0314 °C/yr) for 40°

–

55°N, and

0.0893 °C/yr (95% CI: 0.0113

–

0.1595 °C/yr) for 40°S

–

40°N (Supple-

mentary Table 4). To further assess this gradient, we conducted sup-
plementary analyses by 5° latitude bands, which corroborated these
trends (Supplementary Table 4). These results indicate that stronger
warming is required to trigger positive SOC responses in low-latitude
croplands. In addition,

Δ

SOC showed a stronger positive association

with changes in CDD (

Δ

CDD) (Supplementary Fig. 19a). Under low N,

Δ

SOC was generally negative, whereas at higher N levels,

Δ

SOC

increased with N (Supplementary Fig. 19b).

Based on the segmented response threshold (0.0086 °C/yr) and

the SOC gain-loss crossover threshold (0.0309 °C/yr), we further
examined the relationship between

Δ

SOC and change in environ-

mental factors across different stages.

During the stage where

Δ

Tavg was below 0.0086 °C/yr,

Δ

SOC

correlated with multiple environmental factors, including

Δ

PR,

Δ

Tavg,

Δ

FMD,

Δ

GPP, N, and others (Fig.

b). Structural equation modeling

(SEM) indicated that this stage was dominated by soil property

constraints (Fig.

d): N had the strongest direct negative effect on

Δ

SOC (path coef

ﬁ

cient,

PC

=

−

0.41) and also negatively affected

Δ

GPP

(

PC

=

−

0.39). Because

Δ

GPP had a positive effect on

Δ

SOC (

PC

= 0.18),

N indirectly inhibited SOC accumulation by suppressing

Δ

GPP, thereby

weakening carbon inputs and amplifying its net negative effect on

Δ

SOC. ELE indirectly effected

Δ

SOC (indirect effect = 0.16; Supple-

mentary Fig. 20) via its negative effect on N (

PC

=

−

0.43).

Δ

Tavg and

Δ

PR had negative effects on

Δ

SOC (

PC

=

−

0.27 and

−

0.15,

respectively).

During the stage with

Δ

Tavg between 0.0086 and 0.0309 °C/yr,

Δ

SOC was correlated with

Δ

Tavg,

Δ

Yield,

Δ

GPP, clay, N, and SLO

(Fig.

c).

Δ

SOC was positively affected by

Δ

Tavg, N, clay, and SLO

(

PC

s = 0.28, 0.16, 0.13, and 0.13, respectively). The negative effect of N

on

Δ

GPP persisted (

PC

=

−

0.29; Fig.

e) but was signi

ﬁ

cantly weaker

than in the < 0.0086 °C/yr stage. SLO had a direct positive effect on

Δ

SOC (

PC

= 0.13) and indirect negative effects via

Δ

Yield and N

(

PC

s =

−

0.27 and

−

0.09, respectively; Fig.

e). The negative effects of

Δ

Yield and

Δ

GPP on

Δ

SOC were weaker (

PC

s =

−

0.08 and

−

0.11,

respectively). Compared with the previous stage,

Δ

Tavg had a stronger

direct effect on

Δ

SOC (

PC

= 0.28; Supplementary Fig. 20), whereas the

effect of N weakened (

PC

= 0.16).

In the stage with

Δ

Tavg > 0.0309 °C/yr, SEM indicated that

Δ

PR,

Δ

Tavg,

Δ

GPP, and N were the primary drivers (Fig.

f).

PC

s from

Δ

PR

and

Δ

Tavg to

Δ

SOC were 0.24 and 0.16, respectively. Concurrently,

Δ

Tavg positively affected

Δ

GPP (

PC

= 0.31) and N (

PC

= 0.27), sug-

gesting that strong warming activated both biological and nutrient
processes. However, N had a stronger negative effect on

Δ

GPP (

PC

=

−

0.61) than in earlier stages. Furthermore, N mediated the relation-

ship between

Δ

GPP and

Δ

SOC (

PC

=

−

0.19), partly offsetting the direct

positive climate effect. Overall, this stage featured stronger direct
positive climate effects on

Δ

SOC, partially counteracted by negative

effects from N and

Δ

GPP.

Manure

CDD

CWD

GPP

PR

ELE

min NDTI

FMD

N

Tavg

0

10

20

30

40

50

60

Relative importance (%)

1990 (R

2

=0.85)

Average

P

redictor

v

ariable

2000 (R

2

=0.81)

2010 (R

2

=0.82)

2020 (R

2

=0.81)

40

50

60

70

0

10

20

30

Tavg (°C)

SOC stocks (MgC/ha)

40

50

60

70

0

50

100

150

N (g/kg)

SOC stocks (MgC/ha)

40

50

60

70

0

2

4

6

8

FMD (Ps/1000ha)

SOC stocks (MgC/ha)

40

50

60

70

-1.0

-0.5

0.0

0.5

1.0

min NDTI

SOC stocks (MgC/ha)

40

0

500 1000 1500 2000 2500

50

60

70

ELE (m)

SOC stocks (MgC/ha)

40

50

60

70

10

50

90

130

170

PR (mm)

SOC stocks (MgC/ha)

a

b

c

d

e

f

g

Fig. 2 | Environmental factors in

ﬂ

uencing the spatial variability of soil organic

carbon (SOC) stocks in global Mollisol. a

The relative importance of environ-

mental factors in

ﬂ

uencing the spatial variability of SOC stocks in global Mollisol in

1990, 2000, 2010, and 2020, as determined by the boosted regression tree (BRT)
model. The model was run 10 times to minimize random errors. This study
explored 19 environmental factors, with the top 10 shown in this

ﬁ

gure. The

importance of the remaining factors is presented in Supplementary Fig. 8. Tavg,
annual average temperature; N, nitrogen; FMD, farm machinery data; minNDTI,

minimum normalized difference tillage index; ELE, elevation; PR, annual total
precipitation; GPP, gross primary productivity; CWD, consecutive wet days; CDD,
consecutive dry days; Manure, livestock manure.

b

–

g

Partial dependence plots for

the

ﬁ

rst six environmental factors explaining spatial variation in SOC stocks in

2020. Partial dependence plots for other years and environmental factors are
shown in Supplementary Figs. 9

–

12. The blue line represents a smoothed response

with

ﬁ

tted values (model predictions based on original data), and the shaded areas

indicate 95% con

ﬁ

dence intervals.

Article

https://doi.org/10.1038/s41467-026-73759-w

| (2026) 17:6875

4

$Figures and page graphics, page 5$

Using the SOC gain-loss crossover threshold estimated for each 5°

latitudinal band (Fig.

5

b and Supplementary Table 4), we project that,

under the SSP2-4.5 scenario, 64.76% of Mollisol cropland areas will
experience annual warming rates exceeding this threshold by 2050
(Fig.

5

a and Supplementary Table 5). This projected spatial pattern is

particularly pronounced in southeastern NA, southern SA, and exten-
sive parts of the EU, SI, and AS Mollisol regions. In these regions, SOC
trends may shift from net losses to net gains. We hereafter refer to
these regions as warming-enabled carbon reversal zones.

Discussion

In this study, we developed a data-driven framework integrating
remote sensing, a hybrid probability deep learning model, and

ﬁ

eld-

collected soil samples to generate high-resolution (30 × 30 m) SOC
stock maps for croplands in the global Mollisol region across four
periods (1990, 2000, 2010, and 2020) (Supplementary Fig. 2). Our
maps reliably capture

ﬁ

ne-scale spatial and temporal variability of SOC

stocks (Supplementary Fig. 21), as con

ﬁ

rmed by multiple independent

validations (Supplementary Note 1). The success of high-accuracy SOC
stock prediction can be attributed to three key factors: (1) leveraging
spectral heterogeneity in visible, near-infrared, and shortwave infrared
bands to detect SOC-related re

ﬂ

ectance signatures

35

–

38

(Supplementary

Fig. 22). (2) using a bare soil pixel selecting strategy to generate bare
soil composite image (Supplementary Fig. 23), eliminating the in

ﬂ

u-

ence of surface coverings such as vegetation, snow cover, and crop
residues

,

. (3) implementing a hybrid probability deep learning

model (Supplementary Fig. 24) that partitions the study area into
homogeneous subregions (Supplementary Fig. 25), effectively redu-
cing spatial heterogeneity and boundary artifacts

,

40

(Supplementary

Fig. 26). In addition, the deep learning component of the model
effectively learns the nonlinear relationships between SOC stocks and
remote sensing variables using large-scale

ﬁ

eld-collected data

23

,

41

resulting in improved prediction performance (Supplementary Fig. 3),
with an R

2

increase of 0.16 compared to a conventional global model.

Building on these robust spatial and temporal predictions, we

further investigated how

Δ

SOC responds to climate change, with a

particular focus on different stages of temperature change. During the
stage with

Δ

Tavg < 0.0086 °C/yr (cooling or weak warming), the

association between

Δ

SOC and

Δ

Tavg likely re

ﬂ

ects decomposition

constrained by temperature, while

Δ

GPP remains at a moderate level.

This is consistent with the negative effect of

Δ

Tavg on

Δ

SOC and the

positive effect on

Δ

GPP in Fig.

c, as well as the high temperature

sensitivity of soil decomposition

8

42

. It also accords with evidence that

lower decomposition rates promote increases in SOC stocks

43

,

44

.

Cooling reduces enzyme activity and diffusivity, and

Δ

PR can weaken

aggregate-mineral associations via erosion and leaching, exposing
protected substrates and fostering decomposition and respiration
under moist conditions, thereby reducing stable SOC retention

45

Concurrently, soil pathways did not translate into proportional input
gains: higher N tended to intensify mineralization by accelerating
turnover of soluble substrates and inducing positive feedbacks,
thereby suppressing net carbon sequestration. This aligns with the

Low: -0.06

0

1,500 3,000 4,500 6,000

> 30

20–30

10–20

0–10

< -30

-30–-20

> -20–-10

> -10–0

SOC stock change (MgC/ha)

Annual average temperature
change rate in 1990–2000 (°C/yr)

High: 0.13

GM

NA

SA

EU

SI

AS

0

20

40

60

80

100

Proportion of weight (%)

GM

ΔCDD ΔCWD ΔCWF ΔHWF

ΔPR

ΔT

avg

ΔAR

ΔFMD

ΔManureΔminNDTI

ΔNrate ΔY

ield

ΔGPP ΔNDVI

Clay

N

ASP

ELE

SLO

Climatic factors

Anthropogenic factors

Biological factors Soil

Topographic factors

Correlation (

r

)

0.6

0.4

0.2

0

-0.2

-0.4

-0.6

NA

SA

EU

SI

AS

SLO

ELE

ASP

N

Clay

Ƹ

NDVI

Ƹ

GPP

Ƹ

Yield

Ƹ

Nrate

Ƹ

minNDTI

Ƹ

FMD

Ƹ

Manure

Ƹ

AR

Ƹ

Tavg

Ƹ

PR

Ƹ

HWF

Ƹ

CWF

Ƹ

CWD

Ƹ

CDD

a

c

b

Fig. 3 | Temporal changes in soil organic carbon (SOC) stocks in relation to
changes in environmental factors. a

Spatial distribution of changes in SOC stocks

(

Δ

SOC) associated with the change rate in annual average temperature (

Δ

Tavg).

b

Multiple linear regression analysis method is used to calculate the proportion of

weights that changes in environmental factors in different regions affect

Δ

SOC. GM,

NA, SA, EU, SI, and AS represent the global, North America, South America, Europe,
Siberia, and Asian Mollisol region, respectively.

c

Correlation (

r

) between changes

in environmental variables and

Δ

SOC in different regions.

Δ

AR,

Δ

FMD,

Δ

Manure,

Δ

minNDTI,

Δ

Nrate,

Δ

Yield,

Δ

GPP,

Δ

NDVI,

Δ

CDD,

Δ

CWD,

Δ

CWF,

Δ

HWF, and

Δ

PR

represent change in reclamation age, farm machinery data, livestock manure,
minimum normalized difference tillage index, N fertilizer application intensity,
yield, gross primary production, normalized difference vegetation index, con-
secutive dry days, consecutive wet days, coldwave frequency, heatwave frequency,
and annual total precipitation, respectively. Clay, N, ASP, ELE, and SLO represent
clay content, nitrogen, aspect, elevation, and slope, respectively. Black dots indi-
cate that the factor is extreme signi

ﬁ

cantly (

P

< 0.001) correlated with

Δ

SOC.

Article

https://doi.org/10.1038/s41467-026-73759-w

| (2026) 17:6875

5

negative path from N to

Δ

SOC identi

ﬁ

ed in this stage

39

46

. Collectively,

these

ﬁ

ndings support a regulatory sequence at this stage, in which

cooling leads to decomposition limitation and ultimately enhances
SOC accumulation.

During the stage where

Δ

Tavg ranged from 0.0086 to 0.0309 °C/

yr, the mean

Δ

SOC remained negative (Fig.

a). However, the overall

in

ﬂ

uence of

Δ

Tavg on

Δ

SOC shifted from negative to positive com-

pared to the previous stage, though this effect was partially offset

-40

-20

0

20

40

y = -507.83x - 5.12

Segmented response

threshold = 0.0086

SOC gain-loss

crossover

threshold = 0.0309

r

= -0.48

P

< 0.001

y = 410.35x - 12.66

r

= 0.42

P

< 0.001

-0.04

-0.02

0

0.02

0.04

ΔTavg (°C/yr)

ΔSOC stocks (MgC/ha)

0.0

0.2

0.4

0.6

0.8

Correlation (|

r

|)

ΔTavg <

0.0086

Climatic

Topographic

Soil

Biological

Anthropogenic

-0.64

**

0

.47

**

0

.18

**

0.40

**

0.3

2

**

-0.48

**

-0.16

*

-0.2

6

**

-0.26

**

0.20

**

0.46

**

0.09

0.12

ΔCDD

ΔCW

D

ΔC

W

F

ΔH

WF

ΔPR

Δ

T

avg

ΔA

R

ΔF

MD

ΔManure

ΔNra

te

ΔYi

eld

ΔGPP

ΔNDV

I

Cl

ay

N

ASP

EL

E

SLO

-0.14

0.20

**

-0.22

**

-0.20

**

-0.31

**

0.15

*

ΔTavg (°C/yr)

0.0086

0.0309

a

b

c

ΔSOC

(

R

2

=0.52)

ELE

ΔTavg

ΔGPP

ΔSOC

(

R

2

=0.58)

N

ΔFMD

ΔPR

-0.39

N

Positive

Negative

ΔGPP

Clay

ΔTavg

SLO

ΔYield

-0.29

-0.09

0.05

0.16

0.10

0.13

0.28

0.15

0.13

0.25

-0.27

-0.06

-0.07

-0.11

-0.08

N

ΔGPP

ΔPR

ΔTavg

ΔSOC

(

R

2

=0.45)

-0.61

0.16

0.24

0.16

0.31

0.27

-0.19

-0.12

d

e

f

-0.43

0.19

0.13

0.18

0.25

-0.43

-0.29

-0.13

-0.27

-0.41

-0.01

-0.08

-0.15

0.0

0.2

0.4

0.6

0.8

Correlation (|

r

|)

ΔTavg >

0.0086

Climatic

Topog

raphic

Soil

Biological

An

thropogenic

0.30

**

-0.

21

**

0.

13

**

0.42

**

0.09

*

-0

.2

4

**

-0

.21

**

0.0

6

-0.17

**

ΔCDD

ΔCWD

ΔCWF

ΔH

W

F

ΔP

R

ΔT

a

v

g

Δ

AR

ΔF

MD

ΔManure

ΔNra

te

ΔYi

eld

Δ

G

P

ΔNDVI

C

lay

N

AS

P

EL

E

SLO

-0.09

0.

20

**

0.13

**

0.0

6

0.28

**

0.12

**

0.07

-0.08

*

-0.14

**

0.18

**

ΔminNDTI

Fig. 4 | Segmented response of change in soil organic carbon (

Δ

SOC) stocks to

change rate in annual average temperature (

Δ

Tavg) and its driving mechan-

isms. a

Response of

Δ

SOC to

Δ

Tavg based on segmented linear regression; brown

points indicate observations in the low-warming regime (

Δ

Tavg <0.0086 °C/yr),

whereas green points indicate observations in the higher-warming regime
(

Δ

Tavg > 0.0086 °C/yr).

b

Correlation between changes in environmental factors

and

Δ

SOC at the

Δ

Tavg < 0.0086 °C/yr; brown bars denote correlations for this low-

warming regime.

c

Correlation between changes in environmental factors and

Δ

SOC at the

Δ

Tavg > 0.0086 °C/yr; green bars denote correlations for this higher-

warming regime.

d

Structural equation modeling (SEM) was conducted for the

region with

Δ

Tavg < 0.0086 °C/yr (R

2

= 0.58,

X

2

= 1.81,

df

= 1.00,

P

= 0.18, RMEA =

0.02).

e

SEM was conducted for region with 0.0086 <

Δ

Tavg < 0.0309 °C/yr

(

R

2

= 0.52,

X

2

= 1.33,

df

= 1.00,

P

= 0.25, RMEA = 0.01).

f

SEM was conducted for

region with

Δ

Tavg > 0.0309 °C/yr (

R

2

= 0.45,

X

2

= 1.70,

df

= 1.00,

P

= 0.07, RMEA =

0.01).

Δ

AR,

Δ

FMD,

Δ

Manure,

Δ

minNDTI,

Δ

Nrate,

Δ

Yield,

Δ

GPP,

Δ

NDVI,

Δ

CDD,

Δ

CWD,

Δ

CWF,

Δ

HWF, and

Δ

PR represent change in reclamation age, farm

machinery data, livestock manure, minimum normalized difference tillage index, N
fertilizer application intensity, yield, gross primary production, normalized differ-
ence vegetation index, consecutive dry days, consecutive wet days, coldwave fre-
quency, heatwave frequency, and annual total precipitation, respectively. Clay, N,
ASP, ELE, and SLO represent clay content, nitrogen, aspect, elevation, and slope,
respectively. In panel (

d

–

f

), green dashed arrows indicate positive path coef

ﬁ

cients,

whereas brown solid arrows indicate negative path coef

ﬁ

cients. In panel a, the red

dashed line represents the segmented linear threshold estimated through 1000
bootstrap iterations; the black dashed line represents the SOC gain-loss crossover
threshold; the shaded area denotes the 95% con

ﬁ

dence interval.

r

denotes the

Pearson correlation coef

ﬁ

cient.

P

< 0.001 indicates that the statistical association

reached the level of extreme signi

ﬁ

cance. In (

b

,

c

), black font indicates extreme

signi

ﬁ

cant correlation with

Δ

SOC and is used to establish the SEM.

*

and

**

denote

signi

ﬁ

cance at the

P

< 0.01 and

P

< 0.001 levels, respectively.

Fig. 5 | Soil organic carbon (SOC) stock reversal zones by 2050 under the SSP2-4.5 scenario. a

Spatial distribution of SOC reversal zones by 2050 under the SSP2-

4.5 scenario.

b

Latitudinal gradient of the SOC loss-gain crossover threshold.

Article

https://doi.org/10.1038/s41467-026-73759-w

| (2026) 17:6875

6

through the

Δ

GPP pathway. Slow warming without concurrent

hydrological relief increases vapor-pressure de

ﬁ

cit and thermal stress,

suppressing photosynthesis and stomatal conductance and thereby
reducing

Δ

GPP. This is consistent with reports of warming-induced

productivity declines

47

49

. The positive path from clay to

Δ

SOC indi-

cates that mineral phases began to play a larger storage role in this
stage. Higher clay content and active mineral surfaces promote con-
version of new inputs to mineral-associated organic carbon and reduce
respiratory losses

45

50

. The negative effect of N on

Δ

GPP, together with

its weak positive effect on

Δ

SOC (Fig.

e), indicates a dual role:

accelerating decomposition while offering limited support for stabili-
zation. Under warming, additional N can amplify turnover via
enhanced decomposition, reduced microbial carbon-use ef

ﬁ

ciency,

and priming, without necessarily increasing effective inputs

Conversely, higher N may also re

ﬂ

ect legacy organic matter and fer-

tilization in regions with ample clay and active mineral surfaces, pro-
viding some capacity for carbon sequestration. This helps explain why
N becomes a net negative driver of

Δ

SOC in the

Δ

Tavg < 0.0086 °C/yr

stage. Although SLO had a weak direct positive effect on

Δ

SOC, it

indirectly signaled erosion-driven losses of long-term inputs on slop-
ing land by attenuating the pathway from N to

Δ

yield

55

. Overall, this

stage re

ﬂ

ects a trade-off between input limitation from temperature-

inhibited

Δ

GPP and stabilization provided by clay and aggregates:

when the former dominates, the system tends toward an input-con-
strained, turnover-dominated de

ﬁ

cit; when mineral stabilization partly

compensates for input de

ﬁ

cits, SOC approaches relative stability.

When

Δ

Tavg exceeds 0.0309 °C/yr, the mean

Δ

SOC shifted to

positive values, and the positive correlation between

Δ

SOC and

Δ

Tavg

persisted (Fig.

a). This occurred because the direct positive effects of

Δ

PR and

Δ

Tavg on

Δ

SOC became dominant. Concurrently, climate

established a critical indirect pathway by in

ﬂ

uencing

Δ

GPP, enabling

both hydrotemperate and biogeochemical processes to jointly govern

Δ

SOC dynamics

(Fig.

f). Warming and humidi

ﬁ

cation generally

enhance photosynthesis and root turnover, thereby increasing the
potential for carbon inputs. However, the negative effect of

Δ

GPP on

Δ

SOC and the strong negative effect of N on

Δ

GPP indicate that the

pathway linking N,

Δ

GPP and

Δ

SOC acts as an inhibitory loop that

accelerates carbon turnover at this stage. This pattern is consistent
with mechanisms in which warming and moistening increase substrate
availability and mineralization, reduce microbial carbon use ef

ﬁ

ciency,

and amplify priming effects

42

,

58

,

59

. Consequently, additional carbon

inputs are more readily lost via respiration rather than stabilized in
longer-lived pools

60

,

61

. Thus, this stage essentially represents a com-

petition between hydrothermal forcing

–

driven increases in carbon

inputs and nitrogen-regulated acceleration of turnover, on the one
hand, and the capacity of soils to stabilize additional carbon, on the
other. When increased inputs are balanced by stabilization capacity,
SOC is more likely to show net accumulation; conversely, when
nitrogen-driven turnover dominates, even enhanced inputs are unli-
kely to be retained over the long term

46

. This net gain state could be

partly transient, because under continued rapid warming, the accel-
eration of turnover may eventually outpace input gains, especially if
stabilization capacity becomes limiting.

Overall, the above phased mechanism provides a conceptual basis

for explaining the latitudinal dependence of the SOC gain-loss cross-
over threshold. In cold, energy-limited high-latitude regions (Supple-
mentary Fig. 27), weak warming (

Δ

Tavg < 0.0086 °C/yr) mainly

manifests as a partial release of temperature constraints on decom-
position, accompanied by relatively moderate changes in

Δ

GPP. As

plant inputs and mineral protection gradually begin to compensate for
enhanced decomposition,

Δ

SOC can shift from net loss to net accu-

mulation even at relatively low warming rates. By contrast, in warmer,
more water- and nutrient-limited low-latitude regions (Supplementary
Fig. 27), observed

Δ

Tavg more frequently and persistently falls within

the

moderate

and

strong

warming

stages

identi

ﬁ

ed

here

(0.0086

–

0.0309 and > 0.0309 °C/yr; Fig.

a). Under these conditions,

in the absence of concurrent improvements in hydrological status,
warming often maintains an input-constrained, turnover-dominated
carbon de

ﬁ

cit state. Additional nitrogen is then more likely to rein-

force the inhibitory pathway involving N,

Δ

GPP and

Δ

SOC than to

translate into effective additional carbon inputs. Only when warming is
suf

ﬁ

ciently strong and is accompanied by adequate moisture and clay-

mediated stabilization can these constraints be relaxed, leading to net
SOC gains. This is consistent with the observed spatial pattern in which
SOC loss-gain crossover thresholds increase monotonically toward
lower latitudes.

A key advancement of this study is the empirical identi

ﬁ

cation of a

temperature change threshold that governs changes in SOC stocks
(segmented response threshold, and SOC gain-loss crossover thresh-
old). To our knowledge, this critical threshold and its potential drivers
have not been explicitly examined in global, observation-based studies
of SOC-climate feedbacks. Existing global soil carbon assessments
primarily rely on two approaches. First, spatial statistical products
such as SoilGrids 2.0 provide high-resolution, static maps of SOC
stocks but cannot directly reveal their dynamic relationships with cli-
mate change

10

62

. Second, process-based models, including Earth sys-

tem models and the RothC model, can simulate SOC stock responses
under future scenarios. However, their predictions remain sensitive to
assumptions about model structure and parameterization and are not
constrained by empirically derived thresholds in temperature
change

63

64

. Building on these advances, our study provides an

observation-based estimate of the warming-rate threshold at which
SOC transitions from continuous loss with warming to a reversal from
loss to gain. This reveals that soil-climate feedbacks are not linear but
instead exhibit critical in

ﬂ

ection points. It moves beyond traditional

empirical correlations toward a testable threshold framework, pro-
viding new observational constraints for improving SOC-climate
models and informing climate-adaptation and carbon-sequestration
strategies. Our analytical approach complements controlled incuba-
tion experiments, which offer irreplaceable advantages for elucidating
speci

ﬁ

c mechanisms

8

(e.g., temperature sensitivity) but struggle to

capture the synergistic effects of multiple drivers (climatic, biogiocal,
anthropogenic, etc.) in large-scale agricultural ecosystems. By recon-
structing multi-decadal SOC stock time series across global Mollisol
croplands, we quantitatively characterized the nonlinear relationship
between

Δ

SOC and

Δ

Tavg and precisely identi

ﬁ

ed the SOC gain-loss

crossover threshold. More importantly, using SEM, we revealed a
systematic succession in the dominant controlling factors across dif-
ferent temperature change stages. This integrated empirical frame-
work, which links thresholds, stages and mechanisms, moves SOC-
climate feedback research beyond traditional statistical correlations
and unconstrained model extrapolation. It enables systematic identi-

ﬁ

cation of real-world driving mechanisms and provides critical theo-

retical foundations and data support for developing soil carbon-
climate models with stronger observational constraints.

However, several limitations and sources of uncertainty should be

acknowledged. First, our prediction of SOC reversal zones under the
SSP2-4.5 scenario was derived from statistical inference of the
empirical relationship between

Δ

SOC and

Δ

Tavg, implicitly assuming

that other drivers (land use, crops, and management practices) remain
within their historical ranges. Therefore, these projections should be
regarded as

ﬁ

rst-order estimates of potential net SOC gains driven by

warming alone, rather than as fully coupled scenario simulations.
Second, our analysis is con

ﬁ

ned to the 0

–

20 cm surface soil layer and

does not capture SOC stock dynamics in deeper horizons, where
slower turnover and stronger physical and mineralogical protection
may dampen or delay responses to warming. In addition, the geo-
graphic coverage of independent validation remains imperfect in a few
key Mollisol cropland regions (e.g., South America/Pampas), which
may leave some residual uncertainty in transferability. Third, our

Article

https://doi.org/10.1038/s41467-026-73759-w

https://websoilsurvey.nrcs.usda.

| (2026) 17:6875

7

analysis is based on observational data and statistical models and
therefore primarily reveals correlations between

Δ

SOC and its poten-

tial drivers rather than establishing explicit causal relationships.
Finally, the AR variable, which approximates tillage duration, was cal-
culated using land-use data from 1990 to 2020 and therefore does not
fully capture pre-1990 tillage history. This limitation may partially
confound the estimated in

ﬂ

uence of the AR variable on

Δ

SOC. To

address these limitations, future research should aim to reconstruct
pre-1990 land use and soil management histories by integrating his-
torical literature and long-term statistics. Furthermore, integrating
remote sensing data with process-based models such as CENTURY or
RothC will help predict SOC stocks across multiple soil depths and
provide a means to test the causal mechanisms implied by the iden-
ti

ﬁ

ed temperature change rate thresholds.

Methods

Data sources

Global Mollisol region

. Mollisol are one of the 12 soil orders in the

USDA Soil Taxonomy and are characterized by a black, organic-rich
diagnostic topsoil horizon with high base saturation

65

. These data were

obtained from the Web Soil Survey (

gov/

). The soil map was resampled to a regular 30 m grid, and isolated

patches smaller than a single grid cell were removed.

Cropland dataset

. Land-cover datasets for 1990, 2000, 2010, and

2020 were acquired from the Data Sharing and Service Portal

66

(

https://data.casearth.cn/

). The data have a spatial resolution of 30 m.

For this study, only rainfed cropland (ID = 10) and irrigated cropland
(ID = 20) were selected to delineate the cropland extent.

We used the spatial extent of Mollisol to clip cropland areas for

each period, thereby de

ﬁ

ning the Mollisol cropland regions. All sub-

sequent analyses were restricted to these Mollisol cropland regions.

SOC stock dataset

. In this study, we collected a global SOC stock

database comprising 5535 soil samples from the topsoil layer of
cropland of Mollisol from different periods (Supplementary Fig. 2a and
Supplementary Table 1). To ensure depth consistency, we retained
only records explicitly reported as topsoil within the 0

–

20 cm interval

(e.g., 0

–

10, 0

–

15, 0

–

20 cm) and excluded samples outside this range.

When multiple sublayers within 0

–

20 cm were available for the same

site, we aggregated them to a single 0

–

20 cm value using thickness-

weighted averaging for SOC concentration, thereby ensuring that all
samples correspond to a consistent 0

–

20 cm topsoil de

ﬁ

nition. This

includes a comprehensive dataset of 4345 soil samples from the World
Soil Information Service (WoSIS), which were collected between 1986
and 2015. Another set of 113 soil samples was taken from the Chinese
Academy of Sciences Discipline Data Center for Ecosystems

67

(EDDC)

during 2000

–

2016. These samples were predominantly concentrated

in the eastern part of Northeast China, with some being continuously
monitored for more than

ﬁ

ve years. An additional 1077 soil samples

were collected from croplands within Northeast China between 2016
and 2024. These samples included speci

ﬁ

c details such as the latitude

and longitude coordinates of each sampling site, SOC content, and the
date of soil sample collection. Detailed information on each dataset,
including SOC content data standardization and comparability pro-
cedures, is provided in the Supplementary Note 2. The SOC stock was
calculated using Eq. (

), and the bulk density was computed using Eq.

(

), which was derived from 1630 soil pro

ﬁ

les and 358,772 soil samples

from cropland soils

68

69

.

SOC stock = SOC

=

10 × BD × D

ð

1

Þ

BD = 1

:

377 ×

e

0

:

0048 × SOC

R

2

= 0

:

79,

p

< 001

ð

2

Þ

There was signi

ﬁ

cant spatial heterogeneity in global Mollisol SOC

stock (Supplementary Table 6), with the following average SOC con-
tent across various Mollisol regions within croplands: Europe (73.64
MgC/ha) > Asian (52.95 MgC/ha) > North American (52.06 MgC/ha) >
South American (40.63 MgC/ha).

Satellite-based remote sensing image datasets

. We harnessed the

power of the Landsat Thematic Mapper/Operational Land Imager (TM/
OLI) image to cover the Mollisol region from January 1986 to Decem-
ber 2024. This entailed a vast dataset of 1,086,202 images accessible
via the Google Earth Engine (GEE) platform (Supplementary Fig. 28).
To ensure the quality and accuracy of our analysis, we employed the

‘

pixel_qa

’

bit mask band for each Landsat TM/OLI image, allowing us to

effectively mask out low-quality observations caused by cloud and
snow cover. Six optical bands from the Landsat TM/OLI plus the
thermal band were selected for analysis (Supplementary Table 7).
Variations in the parameters of Landsat TM and OLI sensors result in
differing re

ﬂ

ectance of the same feature, complicating the temporal

analysis of soil properties

70

. To resolve this issue, we obtained Landsat

Enhanced Thematic Mapper Plus (ETM + ) images to calibrate Landsat
TM to the OLI level (Speci

ﬁ

c step is provided in the Supplementary

Note 2). Considering the relatively stable nature of SOC stock, which
tends to exhibit minimal changes over ten-year intervals

11

, we divided

the timeframe spanning from 1986 to 2024 into four distinct periods:
1986

–

1995, 1996

–

2005, 2006

–

2015, and 2016

–

2024. This segmenta-

tion allowed us to create composite images of bare soil for each of
these four periods. The process for generating these bare soil com-
posite images involved the following steps: we established speci

ﬁ

c

criteria for each pixel, ensuring to meet certain conditions, namely, a
normalized difference vegetation index (NDVI) of less than 0.35, a
normalized burn ratio 2 (NBR2) of less than 0.125, and a bare soil index
(BSI) exceeding 0.021 (refs.

38

,

39

.) (Supplementary Note 2 and Sup-

plementary Table 8). In instances where a single observation location
contained multiple image elements, we performed an average calcu-
lation of these elements. This methodology was implemented to
mitigate the potential in

ﬂ

uence of soil moisture on our prediction

results of SOC stock, aligning with established practices

71

. The bare soil

composite images for the four periods are shown in Supplemen-
tary Fig. 23.

Climatic datasets

. For climatic conditions, we collected total pre-

cipitation, minimum and maximum temperature data from the Eur-
opean Center for Medium-Range Weather Forecasts (ECMWF)
atmospheric reanalysis of the global climate (ERA5,

https://www.

ecmwf.int/

) dataset for the period from 1986 to 2024. In this study, we

calculated the annual total precipitation (PR) and annual average
temperature (Tavg). The annual mean temperature data were derived
by averaging the annual minimum and maximum temperatures. For
example, for the period 1986

–

1995, the annual mean temperature is

the average of the annual mean temperatures for those ten years.
Consecutive dry days (CDD), consecutive wet days (CWD), coldwave
frequency (CWF), and heatwave frequency (HWF) were downloaded
from the historical dataset of climate extreme indices

72

at a resolution

of 0.25° × 0.25° (Supplementary Table 9). Finally, we obtained the PR,
Tavg, CDD, CWD, CWF, and HWF for four periods: 1986

–

1995,

1996

–

2005, 2006

–

2015, and 2016

–

2024.

Topographic datasets

. For topographic features, elevation, slope and

aspect were selected. The elevation data were obtained from the
Shuttle Radar Topography Mission-Digital Elevation Model (SRTM-
DEM) data with a spatial resolution of 30 m. Slope and aspect data were
derived from the elevation data with the SAGA-GIS software.

Anthropogenic activities datasets

. For the anthropogenic activity

data, N fertilizer application intensity (Nrate), reclamation age (AR),

Article

https://doi.org/10.1038/s41467-026-73759-w

| (2026) 17:6875

8

livestock manure (Manure), yield (Yield), farm machinery data (FMD),
and crop residue cover were acquired.

We obtained Nrate data for the period from 1986 to 2020 with 5-

arc-min spatial resolution

73

.

The AR was calculated based on land use data obtained for the

years 1990, 2000, 2010, and 2020. The speci

ﬁ

c calculation process for

the reclamation ages is described as follows. The age of reclamation
was categorized into three periods: 1) 10 years (2010

–

2020), 2) 20

years (2000

–

2020), and 3) 30 years (1990

–

2020). The land cover

classi

ﬁ

cation for regions that were cropland in 2020 but not cropland

in 2010 was designated as the reclamation period for years 10 (ref.

https://doi.org/10.1038/s41597-

).

Other reclamation periods were identi

ﬁ

ed by this approach. It needs to

be emphasized that these land use data were utilized to delineate the
cultivated land regions for four distinct periods.

Livestock manure data were obtained from

https://www.fao.org/

home/en/

, where data from 1986 to 2020 for all countries within the

Mollisol region were downloaded and processed into raster data for
analysis.

Yield data were sourced from

020-0433-7

, with a grid resolution of 0.5-degree grid-cell.

FMD data were obtained from

https://ourworldindata.org/

grapher/machinery-per-agricultural-land?tab=table

, where data from

1986 to 2019 for all countries within the Mollisol region were down-
loaded and processed into raster data for analysis.

Crop residue cover was represented using the minimum normal-

ized difference tillage index (minNDTI), an index designed to minimize
the confounding effects of green vegetation independent of
moisture

62

. The index was calculated based on bare soil composite

images.

Biological datasets

. For the biological data, we obtained the gross

primary production (GPP) from 1986 to 2024 using the Landsat GPP
CONUS dataset on the GEE platform, with a spatial resolution of
30 × 30 m. In addition, we obtained the annual average NDVI from 1986
to 2024 on the GEE platform.

Soil properties dataset

. For the soil properties data, we downloaded

the clay content data and nitrogen (N) data with a spatial resolution of
250 m via

https://data.isric.org/

.

To maintain the consistency of the spatial analysis, all raster data

were resampled in ArcGIS 10.8 software employing the nearest-
neighbor method to 30 m spatial resolution.

Determination of input variables for SOC stock prediction

The input variables speci

ﬁ

cally include (Supplementary Table 10):

the six raw bands of Landsat TM/OLI (B, G, R, NIR, SWIR1, SWIR2), the
normalized difference index (NDI), the ratio index (RI), and the
difference index (DI), which are the most sensitive to SOC stock
(Supplementary Fig. 29), and the soil-related soil adjusted vegeta-
tion index (SAVI) and modi

ﬁ

ed soil adjusted vegetation index

(MSAVI) indices

74

. According to the information provided on the

of

ﬁ

cial website of the USGS, the BT band of Landsat TM/OLI was

utilized to calculate the soil temperature (land surface temperature,
LST), since temperature is typically a key driver for change in SOC
content, with elevated temperatures expedite the decomposition of
SOC

17

75

. In addition, land surface wetness (LSW) and the normalized

difference moisture index (NDMI) re

ﬂ

ect soil moisture status, which

can indirectly affect soil carbon stocks by reducing total plant pri-
mary productivity under water stress and increasing vegetation
mortality

76

. Land surface brightness (LSB) and brightness index (BI)

can re

ﬂ

ect the brightness of the soil, which, in general, is usually

lower in areas with high SOC content, whereas the re

ﬂ

ectance is

usually higher in areas with low SOC content (Supplementary
Fig. 22). Lower SOC content causes a brighter soil color and this
indicates lower light absorption and higher re

ﬂ

ectance

77

. The clay

index (CI) re

ﬂ

ects the soil texture

78

, which is considered to be one of

the most promising indicators for SOC content

79

. Consequently, we

selected a total of 17 input variables for SOC stock predicting
(Supplementary Table 10).

Determination of the factor affecting the spatial distribution of
SOC stock and its changes

The factors affecting the spatial distribution of SOC stock include:
climatic factors (CDD, CWD, CWF, HWF, PR, and Tavg), anthropogenic
factors (AR, FMD, Manure, minNDTI, Nrate, and Yield), biological fac-
tors (GPP, and NDVI), soil properties (Clay, and N), and topographic
factors (ASP, ELE, and SLO) (Supplementary Table 9).

The factors affecting the changes in SOC stock include: change in

climatic factors (

Δ

CDD,

Δ

CWD,

Δ

CWF,

Δ

HWF,

Δ

PR, and

Δ

Tavg),

change in anthropogenic factors (

Δ

AR,

Δ

FMD,

Δ

Manure,

Δ

minNDTI,

Δ

Nrate, and

Δ

Yield), change in biological factors (

Δ

GPP, and

Δ

NDVI),

soil properties (Clay, and N), and topographic factors (ASP, ELE, and
SLO). In this study, changes in SOC stock, climatic factors, anthro-
pogenic factors, and biological factors are calculated by comparing the
periods 1990 and 2020.

To assess the validity and reliability of the environmental fac-

tors used to interpret the spatial distribution and temporal change
of SOC stock, we applied three complementary multicollinearity
diagnostics: (i) computed the Pearson correlation coef

ﬁ

cient matrix

among environmental factors (Supplementary Figs. 30 and 31) to
preliminarily identify highly correlated pairs, using |

r

|

≥

0.7 as the

reference threshold

80

; (ii) Calculated the variance in

ﬂ

ation factors

(VIF) (Supplementary Tables 11 and 12), using a threshold of VIF

≤

5

to indicate the absence of signi

ﬁ

cant multicollinearity, and reported

the maximum VIF value

81

; (iii) Applied the Belsley diagnostic

method, reporting the condition index (CI) values (Supplementary
Tables 13 and 14). Moderate and severe multicollinearity were
de

ﬁ

ned as CI > 15 and CI > 30, respectively

82

. Based on the diagnostic

results, no signi

ﬁ

cant multicollinearity was detected among envir-

onmental factors explaining either the spatial distribution or the
temporal change of SOC stock.

Modeling SOC stock dynamics

To link SOC stock observations with time-varying input variables, each
soil sample was

ﬁ

rst assigned to a period according to its sampling year

(for example, 1986

–

1995, 1996

–

2005, 2006

–

2015, and 2016

–

2024).

For every soil sample, we then extracted the values of all raster cov-
ariates from the input variable datasets corresponding to the same
period (Supplementary Fig. 32). Although some positional uncertainty
is inevitable, we used accurate geographic coordinates and a 30 m
pixel size to ensure a close match between sampling locations and
selected pixels.

After calibrating the SOC stock prediction model on all samples

from all regions and periods jointly, we sequentially applied the

ﬁ

tted

model to the raster covariates for each target period. This procedure
yields spatially explicit SOC stock maps through time, including for
regions and periods where no soil samples were available, by relying on
the learned relationships between SOC stocks and input variables
within the observed range of covariate conditions.

Hybrid probability deep learning model

The hybrid probability deep learning model employed in this study
integrates the C-means model with the convolutional neural network
model (CNN). In the hybrid probability deep learning model, the soil
samples were grouped into several clusters using the C-means model
(Supplementary Fig. 24). Next, the CNN model was employed to train
in different clusters separately. Then, the probability that each image
pixel is grouped into different clusters was calculated using RStudio
software. Finally, the predicted SOC stock for a pixel is calculated by
summing the product of the probability of that pixel in each cluster

Article

https://doi.org/10.1038/s41467-026-73759-w

| (2026) 17:6875

9

and the predicted SOC stock of that cluster.

SOC = C

1

% × CNN

C

1

+ C

2

% × CNN

C

2

+

+ C

n

% × CNN

C

n

ð

3

Þ

where SOC represents the SOC stock, CNN is the convolutional neural
network model,

n

represents the number of clusters; and

C

1 represents

the probability that a soil sample or image pixel is classi

ﬁ

ed as cluster 1.

The C-means clustering associates each data point to the nearest
cluster center while maintaining the greatest possible distance
between cluster centers

83

. The algorithm groups soil samples based

on the similarity of the input quantities and performs iterative
calculations to minimize the value of the objective function (J) (Eq.
(

)). The advantage of this algorithm is the introduction of the fuzzy

performance index (FPI), modi

ﬁ

ed packet entropy (MPE), and cluster

independence index (S) to determine the optimal number of clusters.
The FPI index was employed to evaluate the misclassi

ﬁ

cation of soil

samples, while the MPE index assesses the differences among samples
within the same cluster, with

S

denotes the differences among

different clusters. The optimal number of clusters is obtained when
the values of all three indexes are minimized (Supplementary Fig. 33).
In this study, the number of clusters ranged from a minimum of 2 to a
maximum of 12. In addition, the number of soil samples in each cluster
had to be more than 30 to meet the minimum number of sample
points required for modeling. Finally, the optimum number of clusters
was determined to be

ﬁ

ve (Supplementary Fig. 33 and Supplementary

Table 15).

J

=

X

c

i

= 1

X

n

j

= 1

u

m

ij

ð

x

j

c

i

Þ

2

ð

4

Þ

Here,

c

represents the number of principal components,

n

denotes the number of clusters, and

m

signi

ﬁ

es the number of soil

samples,

u

ij

(1

≤

i

≤

c

, 1

≤

j

≤

n

) indicates the membership of soil sample

c

i

to cluster center

x

j

, while

ð

x

j

c

i

Þ

2

is the squared distance between

x

j

and

c

i

in the eigenvector space.

The CNN algorithm consists of multiple layers, each of which

processes the input differently

84

. The convolutional layer serves as the

initial layer of the CNN, which extracts local features of the input using
a

ﬁ

lter (Supplementary Fig. 34), and then passes these features

through a

ReLU

function that handles the nonlinear relationship

between these features and the SOC stock

85

. In the subsequent layers

of the CNN, convolutional layers are generally alternated with pooling
layers. This structure reduces the spatial dimensionality of local fea-
tures, enhancing computational ef

ﬁ

ciency

86

. Finally, the fully con-

nected layer performs SOC stock prediction based on the features
extracted in the previous layer. It is worth noting that different
neighborhood sizes play different roles in SOC stock prediction results
when the input is processed in the input layer

87

. Therefore, to validate

the effect of spatial neighborhood, we tested different neighborhood
sizes. The results of SOC stock prediction for various neighborhood
sizes are displayed in Supplementary Fig. 35, with a

ﬁ

nal neighborhood

size of 9*9 chosen for SOC stock prediction.

In this study, we developed a CNN model using the TensorFlow

framework, comprising three convolutional layers for feature extrac-
tion, three max pooling layers for spatial dimensionality reduction, and
three fully connected layers for

ﬁ

nal prediction (Supplementary

Table 16). A 3D convolutional layer with 3*3*3 pixel

ﬁ

lters and a max

pooling layer with 2*2 pixel

ﬁ

lters were employed. To prevent over-

ﬁ

tting and enhance network generalization, we utilized a dropout

layer

85

and monitored model performance on the validation dataset

during training. In practice, we

ﬁ

rst carried out preliminary experi-

ments with different combinations of learning rate, number of training
epochs and dropout rate and selected the con

ﬁ

guration that mini-

mized the validation RMSE. An

Adam

optimizer was employed to

determine the optimal parameters of the model. The learning rate,

stochastic gradient descent momentum, number of training sessions,
and dropout rate were

ﬁ

nally established at 1e-5, 0.5, 50 and 0.5,

respectively. All hyperparameter choices were made based only on the
calibration and validation dataset, while the independent test dataset
was reserved for

ﬁ

nal evaluation.

To validate the bene

ﬁ

ts of the proposed hybrid probability model,

a global CNN and a random forest (RF) model constructed directly
using the entire soil samples were established for comparison with the
hybrid probability model (Supplementary Fig. 3). The process of
clustering and calculating the probability of dividing into different
clusters using the C-means model is implemented in GEE platform, and
the SOC stock prediction process using the RF and CNN model were
conducted in Python software.

Model evaluation and uncertainty analysis

The soil data, comprising 5535 samples in total, were divided into
calibration (2767 samples), validation (1384 samples), and test
(1384 samples) datasets at a 2:1:1 ratio. The validation dataset is used to
optimize model parameters. The calibration model

’

s accuracy in pre-

dicting SOC stock was assessed using RMSE (Eq. (

5

)), R

2

(Eq. (

6

)), and

RPIQ (Eq. (

7

)) using the test dataset

88

. In this study, we used the stan-

dard deviation (SD) of the 10-fold cross-validation results to assess the
uncertainty of the predictive model in SOC stock mapping. The SD
metric is critical for evaluating the robustness of our model and the
con

ﬁ

dence in our predictions

89

. Firstly, we generated spatial distribu-

tion maps of SOC stock for each period using cross-validation, we then
averaged over the 10 cross-validation iterations to produce the

ﬁ

nal

spatiotemporal distribution maps of SOC stock. Finally, the prediction
uncertainty was quanti

ﬁ

ed by calculating the SD of the 10 SOC stock

spatial distribution maps for each period. A lower uncertainty indicates
higher con

ﬁ

dence in the predictions

90

.

RMSE =

ﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ

1

n

X

m

i = 1

ð

p

i

m

i

Þ

2

s

ð

5

Þ

R

2

= 1

P

n

i = 1

ð

p

i

m

i

Þ

2

P

n

i = 1

ð

p

i

m

i

Þ

2

ð

6

Þ

RPIQ =

IQ

RMSE

ð

7

Þ

where

n

represents the number of samples in

D

T

;

p

i

and

m

i

denote the

predicted and measured SOC contents of the

ith

sample point in

D

T

,

respectively; and

m

is the mean value of the SOC stock of all soil

samples in

D

T

.

IQ

is the interquartile range, speci

ﬁ

cally the difference

between the 3rd and 1st quartiles (

IQ

=

Q

3

Q

1).

SOC stock dataset accuracy validation

To verify the spatial accuracy of the SOC stock dataset, we obtained
SoilGrids250m data

, combined it with lab-measured sample points,

and compared them for accuracy of

ﬁ

t and spatial resolution. In

addition, we compiled and recalculated 63 soil samples to indepen-
dently validate the accuracy of the predicted spatial distribution of
SOC stocks (Supplementary Fig. 36a). These samples were obtained
from the Russian Soil Register and represent SOC stocks in the top
20 cm of soil across Russia

’

s Mollisol regions.

To verify the SOC stock dataset trend accuracy from 1990

–

2020,

we obtained literature from Web of Science, Science direct, Google
Scholar, and Scopus using keywords Organic carbon, Organic matter,
Mollisol, and Chernozem. We examined the primary texts and sup-
plementary materials of these literatures to locate SOC stock data-
bases. Raw data were obtained directly from tables when available, or
alternatively extracted by digitizing graphs using Origin (version
2021b; OriginLab Corporation)

91

. After collecting data related to SOC

Article

https://doi.org/10.1038/s41467-026-73759-w

| (2026) 17:6875

10

stock, a meta-analysis was performed on these data, a process that
mainly involves data standardization. Due to variations in SOC metrics
reported across the literature (e.g., depth intervals and measurement
units), we employed a mass-preserving pro

ﬁ

le function for depth

standardization. This method uniformly converts all records to SOC
stocks within the 0

–

20 cm depth range (MgC/ha)

,

92

, ensuring com-

parability between existing literature SOC stocks and those predicted
in this study. In addition, we simulated annual SOC stocks from 1990
onward using the RothC model, which represents ecosystem carbon
cycling through input, decomposition, transformation, and output
processes. The rates and

ﬂ

uxes of these processes are in

ﬂ

uenced by

environmental conditions, soil properties, vegetation types, land use,
and management practices, making RothC well suited for character-
izing SOC stock changes

93

.

Statistical analysis

Drivers of spatial and temporal variation in SOC stocks

. To eliminate

the effects of spatial autocorrelation, we calculated the minimum
distances to avoid spatial autocorrelation based on the global Moran

’

s

I indices among different environmental factors (Supplementary
Table 17). Ultimately, we re-aggregated the data into 50 km × 50 km
grids

94

, generating 5253 cells for analyzing the drivers of spatial and

temporal variation in SOC stocks.

We used BRT analysis to quantify the relative contributions of

environmental factors, including climatic, anthropogenic, biological,
soil property, and topographic factors, to the spatial variation in SOC
stocks. This approach has been widely applied in ecological research to
explore complex relationships between predictors and response
variables

95

96

. BRT is a machine learning ensemble method that itera-

tively

ﬁ

ts and combines multiple regression trees to enhance pre-

dictive performance

97

. To minimize random variation, we performed

10 iterations of the BRT model. Partial dependence plots (PDP) derived
from BRT analysis were used to illustrate the marginal effects of indi-
vidual environmental factors on SOC stocks, isolating the in

ﬂ

uence of

each variable from others. To identify, for each environmental factor,
the value at which SOC stock is most sensitive, we located the point on
its PDP where the slope was maximal and de

ﬁ

ned this as the value

exerting the greatest in

ﬂ

uence on SOC stocks.

We quanti

ﬁ

ed

Δ

SOC in croplands across the global Mollisol region

by calculating the difference in SOC stocks between 1990 and 2020 at
the pixel level. Multiple linear regression models were constructed to
examine the response of

Δ

SOC to changes in environmental factors,

with

Δ

SOC as the dependent variable and topography, soil properties,

changes in climatic, changes in anthropogenic, and changes in biolo-
gical factors as predictors. All predictor variables were standardized
prior to modeling to eliminate the in

ﬂ

uence of differing units and

magnitudes. Statistical signi

ﬁ

cance was evaluated at the 95% con-

ﬁ

dence level. In addition, GAM was used to test whether a nonlinear

relationship exists between

Δ

SOC and

Δ

Tavg. Using the

ﬁ

rst- and

second-order derivatives of the smoothing function, together with a
residual minimization criterion, we identi

ﬁ

ed segmented response

thresholds. A 95% con

ﬁ

dence interval for these thresholds was esti-

mated using 1000 bootstrap iterations (Supplementary Note 3).

SEM was used to analyze the direct and indirect relationships

between

Δ

SOC and changes in environmental factors. Building on

previous studies

98

https://esgf.github.io/nodes.html

100

, we speci

ﬁ

ed an a priori SEM that hypothesized

causal relationships among these factors (Supplementary Fig. 37).
Prior to modeling, we assessed correlations between

Δ

SOC and

changes in environmental factors to identify key in

ﬂ

uencing factors.

Model

ﬁ

t was evaluated using the root mean square error of approx-

imation (RMEA; acceptable

ﬁ

t: 0.00

≤

RMSEA

≤

0.05), chi-squared test

(X

2

), and probability level (

P

; acceptable

ﬁ

t:

P

> 0.05)

101

.

Data analysis and visualization were conducted using the follow-

ing R packages: caret

102

, gbm

103

, piecewiseSEM, lme4 (ref.

104

) and

ggplot2 (ref.

105

.).

Potential areas for reversal of

Δ

SOC stocks in 2050

. To identify

areas with a potential reversal in

Δ

SOC stock by 2050, we

ﬁ

rst obtained

projections of Tavg for 2020

–

2050 under the SSP2-4.5 scenario from

seven CMIP6 models (AWI-CM-1-1-MR, BCC-CSM2-MR, CMCC-ESM2,
MPI-ESM1-2-HR, MRI-ESM2-0, NorESM2-MM and TaiESM1) via the ESGF
data portal (

). To reduce inter-model

uncertainty, we calculated the multi-model

Δ

Tavg over 2020

–

2050 for

each grid cell. We then applied the empirically derived thresholds
linking

Δ

SOC to

Δ

Tavg, including latitude-speci

ﬁ

c SOC gain-loss

crossover thresholds estimated in 5° latitudinal bands, to statistically
extrapolate the spatial distribution of potential SOC stock reversals in
global Mollisol croplands by 2050 under SSP2-4.5.

To evaluate the robustness of the identi

ﬁ

ed segmented response

threshold, SOC gain-loss crossover threshold, and potential

Δ

SOC

reversal areas, we excluded regions with high prediction uncertainty
and recalculated these thresholds and reversal areas. We then com-
pared the two sets of results. Detailed

ﬁ

ndings are presented in Sup-

plementary Note 4.

Data availability

All data used to support the

ﬁ

ndings of this study are publicly available.

The Landsat surface re

ﬂ

ectance data used in this study are freely

available and can be obtained from

https://earthengine.google.com/

https://isric.org/explore/

The SOC sample dataset are available from

wosis

and

https://www.resdc.cn/

. Global Mollisol region datasets are

available from

https://websoilsurvey.nrcs.usda.gov/

. Cropland region

datasets are available from

https://data.casearth.cn/

. Climatic datasets

are available from

https://www.ecmwf.int/

. Livestock manure data are

available from

https://www.fao.org/home/en/

. Yield data are available

from

https://doi.org/10.1038/s41597-020-0433-7

(ref.

106

). Farm

machinery data are available from

https://ourworldindata.org/

grapher/machinery-per-agricultural-land?tab=table

.

Soil

nitrogen

data are available from

https://data.isric.org/

. Topographic datasets

are available from

http://glovis.usgs.gov/

. The dataset used to predict

and analyze soil organic carbon stocks in global Mollisols is publicly
available

on

Figshare

at

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

(ref.

). Source data are provided in this paper.

Code availability

All code used in the analyses is available via GitHub at

https://github.

com/mxt07150123-hue/Soil-Organic-Carbon-Density-Prediction-and-
Analysis.git

.

References

1.

Scharlemann, J. P. W. et al. Global soil carbon: understanding and
managing the largest terrestrial carbon pool.

Carbon Manag.

5

,

81

–

91 (2014).

2.

Lal, R. Soil carbon sequestration impacts on global climate
change and food security.

Science

304

, 1623

–

1627 (2004).

3.

Padarian, J. et al. Soil carbon sequestration potential in global
croplands.

PeerJ

10

, e13740 (2022).

4.

Lal, R. World cropland soils as a source or sink for atmospheric
carbon.

Adv. Agron.

71

, 145

–

191 (2001).

5.

Rosa, D. D. et al. Soil organic carbon stocks in European croplands
and grasslands: How much have we lost in the past decade?

Glob.

Change Biol.

30

, e16992 (2023).

6.

Sanderman, J., Hengl, T. & Fiske, G. J. Soil carbon debt of 12,000
years of human land use.

PNAS

114

, 9575

–

9580 (2017).

7.

Dietz, C. L. et al. Soil carbon maintained by perennial grasslands
but lost in

ﬁ

eld crop systems over 30 years in a temperate Mollisol

according to longitudinal, compaction-corrected, full-soil pro

ﬁ

le

analysis.

Commun. Earth Environ.

5

, 1

–

8 (2024).

8.

Davidson, E. A. & Janssens, I. A. Temperature sensitivity of soil
carbon decomposition and feedbacks to climate change.

Nature

440

, 165

–

173 (2006).

Article

https://doi.org/10.1038/s41467-026-73759-w

https://doi.org/10.4060/cb9015en

| (2026) 17:6875

11

9.

Hengl, T. et al. SoilGrids250m: global gridded soil information
based on machine learning.

PLoS ONE

12

, e0169748 (2017).

10.

Poggio, L. et al. SoilGrids 2.0: producing soil information for the
globe with quanti

ﬁ

ed spatial uncertainty.

Soil

7

, 217

–

240 (2021).

11.

Meng, X. T. et al. SOC content of global Mollisols at a 30 m spatial
resolution from 1984 to 2021 generated by the novel ML-CNN
prediction model.

Remote Sens. Environ.

300

, 113911 (2024).

12.

Bao, Y. L. et al. Dynamic modeling of topsoil organic carbon and its
scenarios forecast in global Mollisols regions.

J. Clean. Prod.

421

,

138544 (2023).

13.

Shangguan, W. et al. A global soil data set for earth system
modeling.

J. Adv. Model. Earth Sy.

6

, 249

–

263 (2014).

14.

Batjes, N. H. Harmonized soil property values for broad-scale
modelling (WISE30sec) with estimates of global soil carbon
stocks.

Geoderma

269

, 61

–

68 (2016).

15.

FAO and ITPS: Global Soil Organic Carbon Map

–

GSOCmap v.1.6.

(2018).

16.

Georgiou, K. et al. Global stocks and capacity of mineral-
associated soil organic carbon.

Nat. Commun.

13

, 3797 (2022).

17.

Guo, L. B. & Gifford, R. M. Soil carbon stocks and land use change:
a meta-analysis.

Glob. Change Biol.

8

, 345

–

360 (2002).

18.

Bravo, O., Blanco, M. D. C. & Amiotti, N. Control factors in the
segregation of Mollisols and Aridisols of the semiarid

–

arid transi-

tion of Argentina.

CATENA

70

, 220

–

228 (2007).

19.

David, M. B. et al. Long-term changes in mollisol organic carbon
and nitrogen.

J. Environ. Qual.

38

, 200

–

211 (2008).

20.

Liu, H. W. et al. The driving mechanism of soil organic carbon
biodegradability in the black soil region of Northeast China.

Sci.

Total Environ.

884

, 163835 (2023).

21.

Luo, C. et al. Spatial prediction of soil organic matter content
using multiyear synthetic images and partitioning algorithms.

CATENA

211

, 106023 (2022).

22.

He, X. L. et al. Soil organic carbon prediction using phenological
parameters and remote sensing variables generated from
Sentinel-2 images.

CATENA

205

, 105442 (2021).

23.

Meng, X. T. et al. An advanced soil organic carbon content pre-
diction model via fused temporal-spatial-spectral (TSS) informa-
tion based on machine learning and deep learning algorithms.

Remote Sens. Environ.

280

, 113166 (2022).

24.

Wang, X. et al. Effects of cropland reclamation on soil organic
carbon in China

’

s black soil region over the past 35 years.

Glob.

Change Biol.

29

, 5460

–

5477 (2023).

25.

Liu, X. J. et al. Long-term moderate warming shifts soil carbon
cycling but maintains carbon sinks in a subtropical forest.

One

Earth

9

, 101465(2025).

26.

Sun, Y. et al. Warming ampli

ﬁ

es responses of soil organic carbon

to multiple global change drivers.

Glob. Change Biol.

31

,

e70612 (2025).

27.

Wang, M. M. et al. Increased plant productivity exacerbates sub-
soil carbon losses under warming via nitrogen mining.

Nat.

Geosci.

18

, 510

–

517 (2025).

28.

Gao, Y. et al. Precipitation increment reinforced warming-induced
increases in soil mineral-associated and particulate organic
matter under agricultural ecosystem.

Appl. Soil Ecol.

196

,

105301 (2024).

29.

Fang, R. et al. Warming offsets the bene

ﬁ

cial effect of elevated

CO

2

on maize plant-carbon accumulation in particulate organic

carbon pools in a Mollisol.

CATENA

213

, 106219 (2022).

30.

Zosso, C. U. et al. Rapid loss of complex polymers and pyrogenic
carbon in subsoils under whole-soil warming.

Nat. Geosci.

16

,

344

–

348 (2023).

31.

Wang, J. Y. et al. Effects of freezing

–

thawing cycle on peatland

active organic carbon fractions and enzyme activities in the Da
Xing

’

anling Mountains Northeast China.

Environ. Earth Sci.

72

,

1853

–

1860 (2014).

32.

Li, Y. & Li, Y. Nitrogen addition enhances soil carbon and nutrient
dynamics in Chinese croplands: a machine learning and nation-
wide synthesis.

Carbon Bal. Manage.

20

, 15 (2025).

33.

Zhang, L. M. et al. Combined effects of temperature and pre-
cipitation on soil organic carbon changes in the uplands of eastern
China.

Geoderma

337

, 1105

–

1115 (2019).

34.

Bennie, J. J. et al. Characterising inter-annual variation in the
spatial pattern of thermal microclimate in a UK upland using a
combined empirical-physical model.

Agr. Forest Meteorol.

150

,

12

–

19 (2010).

35.

Rossel, R. A. V. & Chen, C. Digitally mapping the information
content of visible

–

near infrared spectra of sur

ﬁ

cial Australian soils.

Remote Sens. Environ.

115

, 1443

–

1455 (2011).

36.

Meng, X. et al. Regional soil organic carbon prediction model
based on a discrete wavelet analysis of hyperspectral satellite
data.

Int. J. Appl. Earth Obs.

89

, 102111 (2020).

37.

Akbari, M. et al. Predicting soil organic carbon by integrating
Landsat 8 OLI, GIS and data mining techniques in semi-arid region.

Earth Sci. Inform.

14

, 2113

–

2122 (2021).

38.

Zepp, S. et al. Optimized bare soil compositing for soil organic
carbon prediction of topsoil croplands in Bavaria using Landsat.

ISPRS J. Photogramm. Remote Sens.

202

, 287

–

302 (2023).

39.

Demattê, J. A. M. et al. Geospatial soil sensing system (GEOS3): a
powerful data mining procedure to retrieve soil spectral re

ﬂ

ec-

tance from satellite images.

Remote Sens. Environ.

212

,

161

–

175 (2018).

40.

Meng, X. T. et al. A long-term global Mollisols SOC content pre-
diction framework: Integrating prior knowledge, geographical
partitioning, and deep learning models with spatio-temporal
validation.

Remote Sens. Environ.

318

, 114592 (2025).

41.

Odebiri, O., Mutanga, O. & Odindi, J. Deep learning-based national
scale soil organic carbon mapping with Sentinel-3 data.

Geo-

derma

411

, 115695 (2022).

42.

Allison, S. D., Wallenstein, M. D. & Bradford, M. A. Soil carbon
response to warming dependent on microbial physiology.

Nat.

Geosci.

3

, 336

–

340 (2010).

43.

Kirschbaum, M. U. F. The temperature dependence of soil organic
matter decomposition, and the effect of global warming on soil
organic C storage.

Soil Biol. Biochem.

27

, 753

–

760 (1995).

44.

Conant, R. T. et al. Temperature and soil organic matter decom-
position rates

–

synthesis of current knowledge and a way forward.

Glob. Change Biol.

17

, 3392

–

3404 (2011).

45.

Lehmann, J. & Kleber, M. The contentious nature of soil organic
matter.

Nature

528

, 60

–

68 (2015).

46.

Janssens, I. A. et al. Reduction of forest soil respiration in response
to nitrogen deposition.

Nat. Geosci.

3

, 315

–

322 (2010).

47.

Novick, K. A. et al. The increasing importance of atmospheric
demand for ecosystem water and carbon

ﬂ

uxes.

Nat. Clim.

Change.

6

, 1023

–

1027 (2016).

48.

Yuan, W. P. et al. Increased atmospheric vapor pressure de

ﬁ

cit

reduces global vegetation growth.

Sci. Adv.

5

, eaax1396 (2019).

49.

Humphrey, V. et al. Soil moisture-atmosphere feedback dom-
inates land carbon uptake variability.

Nature

592

, 65

–

69 (2021).

50.

Crowther, T. W. et al. Quantifying global soil carbon losses in
response to warming.

Nature

540

, 104

–

108 (2016).

51.

Frey, S. D. et al. The temperature response of soil microbial ef

ﬁ

-

ciency and its feedback to climate.

Nat. Clim. Change.

3

,

395

–

398 (2013).

52.

Ye, J. S. et al. Increasing microbial carbon use ef

ﬁ

ciency with

warming predicts soil heterotrophic respiration globally.

Glob.

Change Biol.

25

, 3354

–

3364 (2019).

53.

Tao, F. et al. Microbial carbon use ef

ﬁ

ciency promotes global soil

carbon storage.

Nature

618

, 981

–

985 (2023).

54.

Kuzyakov, Y. Priming effects: Interactions between living and dead
organic matter.

Soil Biol. Biochem.

42

, 1363

–

1371 (2010).

Article

https://doi.org/10.1038/s41467-026-73759-w

http://download.nust.na/pub3/cran/web/

| (2026) 17:6875

12

55.

Shi, W. H., Huang, M. B. & Wu, L. H. Prediction of storm-based
nutrient loss incorporating the estimated runoff and soil loss at a
slope scale on the Loess Plateau.

Land Degrad. Dev.

29

,

2899

–

2910 (2018).

56.

Keenan, T. F. et al. Net carbon uptake has increased through
warming-induced changes in temperate forest phenology.

Nat.

Clim. Change.

4

, 598

–

604 (2014).

57.

Liu, Z. H. et al. Precipitation thresholds regulate continental-scale
net carbon exchange of the terrestrial biosphere.

Nat. Commun.

9

,

3596 (2018).

58.

Ren, C. et al. Thermal sensitivity of soil microbial carbon use ef

ﬁ

-

ciency across forest biomes.

Nat. Commun.

15

, 6269 (2024).

59.

Wang, G. et al. Enhanced response of soil respiration to experi-
mental warming upon thermokarst formation.

Nat. Geosci.

17

,

532

–

538 (2024).

60.

Tao, X. Y. et al. Experimental warming accelerates positive soil
priming in a temperate grassland ecosystem.

Nat. Commun.

15

,

45277 (2024).

61.

Bond-Lamberty, B. & Thomson, A. Temperature-associated
increases in the global soil respiration record.

Nature

464

,

579

–

582 (2010).

62.

Tang, Q. et al. Global soil organic and inorganic carbon vulner-
ability in response to irrigation.

Commun. Earth Environ.

6

,

599 (2025).

63.

Todd-Brown, K. E. O. et al. Causes of variation in soil carbon
simulations from CMIP5 Earth system models and comparison
with observations.

Biogeosciences

10

, 1717

–

1736 (2013).

64.

Coleman, K., & Jenkinson, D. S. RothC-26.3 - A Model for the
turnover of carbon in soil. in

Evaluation of Soil Organic Matter

Model

(1996).

65.

Soil Survey Staff.

Soil Taxonomy: A Basic System of Soil Classi

ﬁ

-

cation for Making and Interpreting Soil Surveys.

(1999).

66.

Liu, L. Y., & Zhang, X. Global 30-meter

ﬁ

ne surface coverage

products in 2020. in

Aerospace Information Research Institute,

Chinese Academy of Sciences

, (2020).

67.

Xu, L. et al. Carbon storage in China

’

s terrestrial ecosystems: A

synthesis.

Sci. Rep-UK.

8

, 2806 (2018).

68.

Song, G. H. et al. Topsoil organic carbon storage of China and its
loss by cultivation.

Biogeochemistry

74

, 47

–

62 (2005).

69.

Rambaut, L. A. E. et al. 15-year fertilization increased soil
organic carbon stock even in systems reputed to be saturated
like permanent grassland on andosols.

Geoderma

425

,

116025 (2022).

70.

Mishra, N. et al. Continuous calibration improvement in solar
re

ﬂ

ective bands: Landsat 5 through Landsat 8.

Remote Sens.

Environ.

185

, 7

–

15 (2016).

71.

Gasmi, A. et al. Mean spectral re

ﬂ

ectance from bare soil pixels

along a Landsat-TM time series to increase both the prediction
accuracy of soil clay content and mapping coverage.

Geoderma

388

, 114864 (2021).

72.

Mistry, M. N. A high-resolution global gridded historical dataset of
climate extreme indices.

Data

4

, 41 (2019).

73.

Adalibieke, W. et al. Global crop-speci

ﬁ

c nitrogen fertilization

dataset in 1961

–

2020.

Sci. Data.

10

, 1

–

17 (2023).

74.

Qi, J. et al. A modi

ﬁ

ed soil adjusted vegetation index.

Remote

Sens. Environ.

48

, 119

–

126 (1994).

75.

van Gestel, N. et al. Predicting soil carbon loss with warming.

Nature

554

, E4

–

E5 (2018).

76.

Green, J. K. et al. Large in

ﬂ

uence of soil moisture on long-term

terrestrial carbon uptake.

Nature

565

, 476

–

479 (2019).

77.

Kuang, B. & Mouazen, A. M. Non-biased prediction of soil organic
carbon and total nitrogen with vis-NIR spectroscopy, as affected
by soil moisture content and texture.

Biosyst. Eng.

114

,

249

–

258 (2013).

78.

Sabetizade, M. et al. Combination of MIR spectroscopy and
environmental covariates to predict soil organic carbon in a semi-
arid region.

CATENA

196

, 104844 (2020).

79.

Meersmans, J. et al. A multiple regression approach to assess the
spatial distribution of soil organic carbon (SOC) at the regional
scale (

ﬂ

anders, belgium).

Geoderma

143

, 1

–

13 (2008).

80.

Li, H. W. et al. Decipher soil organic carbon dynamics and driving
forces across China using machine learning.

Glob. Change Biol.

28

, 3394

–

3410 (2022).

81.

O

’

brien, R. M. A Caution Regarding Rules of Thumb for Variance

In

ﬂ

ation Factors.

Qual. Quant.

41

, 673

–

690 (2007).

82.

Belsley, D. A., Kuh, E. & Welsch, R. E.

Regression Diagnostics;

Identifying In

ﬂ

uence Data and Source of Collinearity.

(1980).

83.

Orhan, U., Hekim, M. & Ozer, M. EEG signals classi

ﬁ

cation using

the K-means clustering and a multilayer perceptron neural net-
work model.

Expert Syst. Appl.

38

, 13475

–

13481 (2011).

84.

LeCun, Y., Bengio, Y. & Hinton, G. Deep learning.

Nature

521

,

436

–

444 (2015).

85.

Krizhevsky, A., Sutskever, I. & Hinton, G. E. Imagenet classi-

ﬁ

cation with deep convolutional neural networks. In

Pro-

ceedings Advances in Neural Information Processing Systems
25

(2012).

86.

Nogueira, K., Penatti, O. A. & dos Santos, J. A. Towards better
exploiting convolutional neural networks for remote sensing
scene classi

ﬁ

cation.

Pattern Recogn.

61

, 539

–

556 (2017).

87.

Yang, L. et al. A deep learning method to predict soil organic
carbon content at a regional scale using satellite-based phenol-
ogy variables.

Int. J. Appl. Earth Obs.

102

, 102428 (2021).

88.

Williams, P. C. Variables affecting near-infrared re

ﬂ

ectance

spectroscopic analysis. In

Near-Infrared Technology in the Agri-

cultural and Food Industries

, 143

–

167 (1987).

89.

Goidts, E., van Wesemael, B. & Cruci

ﬁ

x, M. Magnitude and sources

of uncertainties in soil organic carbon (SOC) stock assessments at
various scales.

Eur. J. Soil Sci.

60

, 723

–

739 (2009).

90.

Nauman, T. W. & Duniway, M. C. Relative prediction intervals
reveal larger uncertainty in 3D approaches to predictive digital soil
mapping of soil properties with legacy data.

Geoderma

347

,

170

–

184 (2019).

91.

He, X. et al. Global patterns and drivers of soil total phosphorus
concentration.

Earth Syst. Sci. Data.

13

, 5831

–

5846 (2021).

92.

Bishop, T. F. A., McBratney, A. B. & Laslett, G. M. Modelling soil
attribute depth functions with equal-area quadratic smoothing
splines.

Geoderma

91

, 27

–

45 (1999).

93.

Jenkinson, D. S. & Rayner, J. H. Turnover of soil organic matter in
some of the Rothamsted classical experiments.

Soil Sci.

123

,

298

–

305 (1977).

94.

Xi, Y. B. et al. Boreal tree species diversity increases with global
warming but is reversed by extremes.

Nat.Plants.

10

,

1473

–

1483 (2024).

95.

Elith, J., Leathwick, J. R. & Hastie, T. A working guide to boosted
regression trees.

J. Anim. Ecol.

77

, 802

–

813 (2008).

96.

Elith, J. & Leathwick, J. Boosted Regression Trees for ecological
modeling. CRAN

packages/dismo/vignettes/brt.pdf

(2017).

97.

Zhang, W. et al. Socio-economic and climatic changes lead to
contrasting global urban vegetation trends.

Global Environ.

Change.

71

, 102385 (2021).

98.

Brahim, N. et al. Application of structural equation modeling for
assessing relationships between organic carbon and soil proper-
ties in semiarid Mediterranean region.

Int. J. Environ. Sci. Te.

8

,

305

–

320 (2011).

99.

Rabbi, S. M. et al. Climate and soil properties limit the positive
effects of land use reversion on carbon storage in Eastern Aus-
tralia.

Sci. Rep.

5

, 17866 (2015).

Article

https://doi.org/10.1038/s41467-026-73759-w

https://github.com/tidyverse/ggplot2/releases/tag/v2.0.
0

| (2026) 17:6875

13

100. Dai, L. et al. In

ﬂ

uence of soil properties, topography, and land

cover on soil organic carbon and total nitrogen concentration: A
case study in Qinghai-Tibet plateau based on random forest
regression and structural equation modeling.

Sci. Total Environ.

821

, 153440 (2022).

101.

Delgado-Baquerizo, M. et al. Multiple elements of soil biodiversity
drive ecosystem functions across biomes.

Nat. Ecol. Evol.

4

,

210

–

220 (2020).

102. Kuhn, M. Building predictive models in R using the caret package.

J. Stat. Softw.

28

, 1

–

26 (2008).

103. Ridgeway, G.

Generalized Boosted Models: A Guide to the gbm

Package

(CRAN, 2007).

104. Bates, D. et al. Package

‘

lme4

’

.

Convergence

12

, 2 (2015).

105. Wickham, H. et al. ggplot2: Create elegant data visualisations

using the grammar of graphics. R package version 2.0.0. GitHub

(2015).

106. Iizumi, T. & Sakai, T. The global dataset of historical yields for major

crops 1981

–

2016.

Sci. Data.

7

, 97 (2020).

107. Meng, X. A dataset used to predict and analyze soil organic carbon

stocks in global Mollisol.

ﬁ

gshare

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

ﬁ

gshare.31980039

(2026).

Acknowledgements

The authors thank all colleagues and collaborators who contributed to
sample collection, data processing, and manuscript revision.

Author contributions

All authors contributed intellectual input and assistance to this study
and manuscript preparation. X.T.M. and H.J.L. designed the original
concept and experiment strategy. Y.L.B. and C.L. collected the data.
X.L.Z. analyzed the data. X.T.M., Y.L.B., S.L.U., C.L., X.L.Z., M.S., and H.J.L.
drafted the manuscript, and all authors involved in critical revisions and
gave

ﬁ

nal approval for publication.

Funding

This work was supported by the National Key Research and
Development Program of China (grant no. 2021YFD1500100, H.J.L.),
the Young Scientists Innovation Fund of the State Key Laboratory
of Black Soils Conservation and Utilization (grant no. 2023HTDGZ-
QN-01, H.J.L.), and the Postdoctoral Fellowship Program of

the China Postdoctoral Science Foundation (grant no.
GZB20240737, X.T.M.).

Competing interests

The authors declare no competing interests.

Additional information

Supplementary information

The online version contains

supplementary material available at

https://doi.org/10.1038/s41467-026-73759-w

http://www.nature.com/reprints

Correspondence

and requests for materials should be addressed to

Yilin Bao or Huanjun Liu.

Peer review information

Nature Communications

thanks Zhou Shi and

the other anonymous reviewer(s) for their contribution to the peer
review of this work. A peer review

ﬁ

le is available.

Reprints and permissions information

is available at

Publisher

’

s note

Springer Nature remains neutral with regard to jur-

isdictional claims in published maps and institutional af

ﬁ

liations.

Open Access

This article is licensed under a Creative Commons

Attribution-NonCommercial-NoDerivatives 4.0 International License,
which permits any non-commercial use, sharing, distribution and
reproduction in any medium or format, as long as you give appropriate
credit to the original author(s) and the source, provide a link to the
Creative Commons licence, and indicate if you modi

ﬁ

ed the licensed

material. You do not have permission under this licence to share adapted
material derived from this article or parts of it. The images or other third
party material in this article are included in the article

’

s Creative

Commons licence, unless indicated otherwise in a credit line to the
material. If material is not included in the article

’

s Creative Commons

licence and your intended use is not permitted by statutory regulation or
exceeds the permitted use, you will need to obtain permission directly
from the copyright holder. To view a copy of this licence, visit

http://

creativecommons.org/licenses/by-nc-nd/4.0/

Article

https://doi.org/10.1038/s41467-026-73759-w