Shift in controlling factors of carbon stocks across biomes on the Qinghai-Tibetan Plateau

The Qinghai-TP lies between 26°0ʹ–39°47ʹN, 73°19ʹ–104°47ʹE, with an area of ∼2.5 × 10⁶km² [24]. Its climate is characterized by intense solar radiation, long durations of sunshine, low temperature, and extensive diurnal temperature range. Grassland as the dominant ecosystem on the plateau occupying over 60% of the total area [25]. Forests and shrublands summed together to ∼5.35% of the land surface and scattered in the central and western TP. Barren lands were mainly distributed in the northern TP with dry and cold conditions.

2.2.1. SOC and biomass data

We collected 4809 SOC observations and 2387 biomass data from multiple regional field sampling campaigns across the plateau since the 1980s [1, 2, 15–17], and established the most representative soil dataset to date, which covers the main vegetation types and spatial scope (see supplementary text in details). The sites distribution map is shown in figure 1.

Zoom In Zoom Out Reset image size

SOC stock (SOC_s, kg C m⁻²) was calculated by the following equation (1):

$$\begin{equation}{\text{SO}}{{\text{C}}_{\text{S}}}{{ = }}\frac{{\text{1}}}{{{\text{100}}}} \times \sum\limits_{j{\text{ = 1}}}^n {\text{SO}}{{\text{C}}_j} \times {\text{(1}} - {C_j}{\text{)}} \times {\text{B}}{{\text{D}}_j} \times {D_j}\end{equation} \tag{ 1 }$$

where j is the soil horizon (i.e. 1, 2, 3, ..., n), SOC is SOC concentration on a mass basis (g kg⁻¹), C_j is the volume percentage of gravel fraction, BD is the soil bulk density corrected (g cm⁻³), and D_j is the thickness of the jth horizon (cm). There were ∼7.5% (360/4809) SOC measurements missed their corresponding bulk density data at specific sites and layers, we hereby built an empirical (nonlinear) relationship between SOC and BD [18] to estimate the missing BD (supplementary figure S1).

Generally, each site in our dataset has multiple layers with different depth intervals. Sites with the maximum depth larger than 30 cm were used to calculate SOC stock of 0–30 cm. While sites with the maximum depth larger than 70 cm were used to estimate SOC stock of 30–100 cm, 0–100 cm, sites with depth larger than 100 cm were used to estimate SOC stock of 100–200 cm. In sum, we got 2019 sites for 0–30 cm, 1181 sites for 30–100 cm, 1302 sites for 0–100 cm, and 307 sites for 100–200 cm. These measurements are sufficient for us to assess the C stocks of the TP and further explore its controlling factors.

2.2.2. Environmental variables

Rotated PCA (details see supplementary text) was used to reduce dimension of the independent variables and further improve the interpretation of larger groups of controls [26]. According to the variable loadings of rotated PCA and correlations among the meteorological variables (supplementary figure S2, tables S2–S5), we found temperature related 11 bioclimatic factors and precipitation related eight bioclimatic factors were always classified into different components. In addition, MAT has high correlation with variables, such as solar irradiation (RD), vapor pressure deficit (VPD) and ratio of MAT to solar irradiation (MAT_RD), and mean annual precipitation (MAP) has high correlation with variables like aridity (AI), evapotranspiration (ET), soil moisture (SM) and ratio of ET to solar irradiation (ET_RD) (supplementary tables S2–S5), we thereby grouped these principal components (PCs) related to these climatic variables into two larger categories: temperature related factors, including 11 temperature derived bioclimatic factors, RD, VPD and MAT_RD, and water related factors, including eight precipitation derived bioclimatic variables and AI, ET, SM and ET_RD. The rest PCs were further grouped into three different environmental categories including AGB, which always related with diurnal temperature range, soil physico-chemical properties and topographic factors (supplementary tables S2–S5). Machine learning algorithms was then used to build the relationships between these PCs and different C stocks and further identify the importance of each PC.

Combined with the spatial environmental covariates, we trained and tuned (supplementary text 1.3.1) the models with several machine learning algorithms, including support vector machine (SVM), random forest (RF), extreme gradient boosting and feedforward neural network. To test the model performance, we used the five-fold cross-validation approach and calculated the mean square errors as the criteria [27]. RF finally outperformed other algorithms (supplementary figure S3, table S7) and was adopted to mapping C stocks (at a resolution of 5 km) on the TP.

We used two methods to quantify the uncertainty of our estimate, absolute uncertainty and relative uncertainty (details see supplementary text). Briefly, for the former one, we used a quantile regression forest approach to calculated the interquartile of each pixel (difference between the 75th and 25th percentiles), further the summed and averaged value of the 25th and 75th percentiles were used to represent the uncertainty of OC stocks on the TP [2]. To get the relative uncertainty, the sensitivity and model structure uncertainty were then summed up to represent the total uncertainty, which then was divided by the median values of each pixel to obtain the relative uncertainty [2, 28, 29].

If the P-value of Kolmogorov Smirnov [30] test was smaller than 0.05, we used log transformation for the data for further model establishing and analysis. Biomass C, including AGB and BGB, in this study were all transformed into base-10 log scale for further analysis. They were transformed back after model prediction. All statistical analysis were conducted in R platform [31]. The packages we mainly used were listed in supplementary text.

The relationship between C stocks (bin-averaged) and climatic factors varied with vegetation types (figure 2). Specifically, both grassland AGB and BGB showed a significant positive tendency with MAP (p < 0.01, figure 2(a)), but their correlations with MAT were sharply weakened (figure 2(d)). While for shrubland, a decreasing relationship with MAT was found for both its AGB and BGB (p < 0.01, figure 2(e)). In addition, only AGB in forest showed a positive relationship with MAP (p < 0.01, figures 2(c) and (f)).

Zoom In Zoom Out Reset image size

Moreover, the correlations between SOC and MAP gradually shifted from positive in grassland (p < 0.01) to negative in forest (p < 0.01) for both topsoil and subsoil (figures 2(g) and (i)). SOC stock in shrubland showed an opposite trend to MAP at the two depth intervals, increasing in topsoil and decreasing in subsoil. However, the response of SOC stock to MAT was almost synchronized at both soil layers (figures 2(i)–(l)), which all showed a negative relationship with MAT for the three vegetation types.

The primary governing environmental category differs for vegetation and soil C stocks on the TP. For AGB, temperature and water related factors contributed the most to its spatial distribution, taking half of its variation across the TP (figure 3(a)), indicating the dominant role of climatic variables in shaping amount of AGB. While for the individual land covers, the main driving environmental category changed from precipitations related factors in grassland to soil properties and temperature related factors in shrubland and forest (figure 3(a)), respectively.

Zoom In Zoom Out Reset image size

For SOC stock, AGB played a primary role for both topsoil (0–30 cm) and subsoil (30–100 cm) (figures 3(b) and (c)), demonstrating the importance of plant-derived C input to SOC stocks on the TP. For topsoil, both AGB and water related factors (MAP, etc) showed positive effect on the SOC stocks across the TP, while temperature related factors like temperature fluctuation ratio of monthly to yearly (isothermality, supplementary figure S5(e)) almost showed a linear negative effect on SOC stocks. However, the dominant driving environmental category varied among different vegetation types (figure 3(b); supplementary figures S4(e)–(h)), which changed from water related factors in grassland to topographic factors, temperature related factors in shrubland and forest (figure 3(b)), suggesting divergent controlling factors among land covers on the TP.

Unlike topsoil, besides AGB, there were no exceptional environmental categories contributing to subsoil SOC stocks. Categories such as temperature related and soil physico-chemical properties had almost the same contribution to SOC in sublayers throughout the TP (figure 3(c)). Moreover, the dominant regulatory categories were relatively consistent for the individual land cover (figure 3(c)), which were all referring to water and temperature related factors, revealing a relative complex mechanism in regulating the SOC in deeper soils. Specifically, the influence of individual PCs on different SOC stocks with depth exhibited the same trend, which all increased with precipitation and decreased with temperature related factors (supplementary figures S4(i)–(l)).

The spatial distribution of both vegetation and soil C densities presented a similar pattern. Namely, the C density decreased from the east and southeast parts to the central and western parts and then increased again at the northwest edge within Xinjiang Uygur autonomous region (figures 4(a)–(c) and (g)–(i)). While, their uncertainties showed an opposite tendency between plant and soil carbon pools (insets of figures 4(a)–(c) and (g)–(i)). The uncertainties of vegetation C density exhibited a similar pattern with its density distribution, suggesting proportional uncertainties to their C density. While the highest uncertainty for soil C densities within 1 m depth was concentrated in the northwest part and then decreased to the eastern part (insets of figures 4(g) and (h)).

Zoom In Zoom Out Reset image size

There was a total of ${\text{2}}{\text{.4}}_{{\text{0}}{\text{.95}}}^{{\text{5}}{\text{.16}}}$ Pg C (confidence interval of upper and lower estimate indicated 75th and 25th percentiles to quantify the prediction uncertainty) stored in vegetation on the TP. Of which ${\text{1}}{\text{.03}}_{{\text{0}}{\text{.2}}}^{{\text{2}}{\text{.7}}}$ Pg was allocated to AGB and ${\text{1}}{\text{.37}}_{{\text{0}}{\text{.75}}}^{{\text{2}}{\text{.45}}}$ Pg was assigned to BGB (figures 4(a)–(c)). The average carbon density of AGB was increased from 0.10 kg m⁻² in grassland to 0.52 and 2.86 kg m⁻² in shrubland and forest, respectively. Correspondingly, their carbon stock size increased from ${\text{0}}{\text{.16}}_{{\text{0}}{\text{.04}}}^{{\text{0}}{\text{.5}}}$ Pg in grassland to ${\text{0}}{\text{.27}}_{{\text{0}}{\text{.06}}}^{{\text{0}}{\text{.8}}}$ and ${\text{0}}{\text{.56}}_{{\text{0}}{\text{.09}}}^{{\text{1}}{\text{.22}}}$ Pg in shrubland and forest. In contrast, the carbon stored in BGB was decreased from ${\text{0}}{\text{.7}}_{{\text{0}}{\text{.36}}}^{{\text{1}}{\text{.31}}}$ Pg in grassland to ${\text{0}}{\text{.29}}_{{\text{0}}{\text{.18}}}^{{\text{0}}{\text{.49}}}$ and ${\text{0}}{\text{.26}}_{{\text{0}}{\text{.16}}}^{{\text{0}}{\text{.41}}}$ Pg in forest (figures 4(d)–(f)).

The estimated SOC stock on the entire TP was ${\text{32}}{\text{.01}}_{{\text{19}}{\text{.69}}}^{{\text{47}}{\text{.9}}}$ Pg (${\text{11}}{\text{.72}}_{{\text{7}}{\text{.2}}}^{{\text{17}}{\text{.53}}}$ kg m⁻² on average) for 0–100 cm soil and ${\text{40}}{\text{.11}}_{{\text{24}}{\text{.55}}}^{{\text{65}}{\text{.06}}}$ Pg for 0–200 cm depth (${\text{14}}{\text{.68}}_{{\text{8}}{\text{.98}}}^{{\text{23}}{\text{.81}}}$ kg m⁻² on average). Though the average SOC stock (0–100 cm) increased from grassland (9.95 kg m⁻²) to forest (19.12 kg m⁻²), the opposite trend of their total storage was found, the largest SOC pool was grassland soils (${\text{15}}{\text{.52}}_{{\text{9}}{\text{.91}}}^{{\text{23}}{\text{.5}}}$ Pg, 48.5%), of which ${\text{10}}{\text{.2}}_{{\text{6}}{\text{.62}}}^{{\text{15}}{\text{.3}}}$ Pg stored in alpine meadow and ${\text{3}}{\text{.87}}_{{\text{2}}{\text{.36}}}^{{\text{5}}{\text{.88}}}$ Pg stored in alpine steppe, followed by shrubland (${\text{7}}{\text{.52}}_{{\text{4}}{\text{.8}}}^{11.6}$ Pg, 23.5%) and forest (${\text{3}}{\text{.72}}_{{\text{2}}{\text{.5}}}^{{\text{5}}{\text{.36}}}$ Pg, 11.6%) (figures 4(f)–(h)). There was ${\text{14}}{\text{.27}}_{{\text{9}}{\text{.77}}}^{{\text{20}}{\text{.15}}}$ Pg SOC stored in 0–30 cm depth, which accounted for 44.57% of 0–100 cm SOC storage, and of which ${\text{7}}{\text{.81}}_{{\text{5}}{\text{.4}}}^{{\text{11}}{\text{.0}}}$ Pg stored in grassland (${\text{5}}{\text{.47}}_{{\text{3}}{\text{.89}}}^{7.4}$ Pg stored in alpine meadow and ${\text{1}}{\text{.81}}_{{\text{1}}{\text{.18}}}^{{\text{2}}{\text{.82}}}$ Pg stored in alpine steppe), ${\text{3}}{\text{.13}}_{{\text{2}}{\text{.2}}}^{{\text{4}}{\text{.34}}}$ Pg in shrubland and ${\text{1}}{\text{.6}}_{{\text{1}}{\text{.11}}}^{{\text{2}}{\text{.23}}}$ Pg in forest, respectively. All indicated that grassland, especially its alpine meadow system, played as the largest C reservoir to the TP.

SOC and AGB were closely related on the TP, especially for grassland and shrubland (figure 5). SOC increased significantly with the increase of AGB in grassland for the topsoil (p < 0.01) (figure 5(a)). Though the strength of the correlation (bin averaged) decreased in deeper soil (30–100 cm), the coefficient still at a high level (0.78) (figure 5(b)). There was also a strong positive correlation between AGB and SOC in shrubland in topsoil (0–30 cm), of which the coefficient was 0.70 (p < 0.01). For the forest, though there was a negative trend between AGB and SOC, they are not statistically significant (p > 0.05) (figure 5).

Zoom In Zoom Out Reset image size

Based on the unique compiled dataset, widely collected environmental variables, and the advances of machine learning (RF) approach, this study provided a comprehensive assessment of the TP for both its biomass and soil carbon pools. The total estimated biomass carbon stock was ${\text{2}}{\text{.4}}_{{\text{0}}{\text{.95}}}^{{\text{5}}{\text{.16}}}$ Pg, which accounted for 23.4% of China's total vegetation carbon stock (figure 6(a)). The estimated biomass C stock size of grassland was 0.86 Pg, which accounted for 63.7% of China's total grassland biomass C size (1.35 Pg from Tang et al [34]), this proportion was similar with the findings by Fan et al (56.4%) [35], but the amount was much lower than theirs (1.87 Pg). Our results of biomass C densities were higher in shrubland and lower in forest compared with China's average value by Tang et al [34] (1.07 vs 0.64 kg m⁻² in shrubland, and 4.2 vs 5.32 kg m⁻² in forest).

Zoom In Zoom Out Reset image size

The estimated SOC stock of the entire TP was ${\text{32}}{\text{.01}}_{{\text{19}}{\text{.69}}}^{{\text{47}}{\text{.9}}}$ Pg C (1 m depth), with an average density of ${\text{11}}{\text{.72}}_{{\text{7}}{\text{.2}}}^{{\text{17}}{\text{.53}}}$ kg C m⁻², correspondingly accounting for ${\text{37}}{\text{.2}}_{{\text{22}}{\text{.9}}}^{{\text{55}}{\text{.6}}}$% (25%–75% quantile) of China's SOC stocks (median = 86.09 Pg) (figure 6(a)), and ∼2.2% of global SOC stock (synthesized ∼1448.9 Pg) [36, 37]. This size was much similar to Jiang's et al [14] report (32.79 Pg with 12.68 kg m⁻² on average) but lower than the estimates of Fang [13] (38.4 Pg for 0–0.72 m depth) and Tian et al [12] (49 Pg for 0–0.75 m depth) (table 1). To compare with the results from earth system models (ESMs), we chose ten ESMs derived from Coupled Model Intercomparison Project v.5 over the TP, and found our estimates (0–100 cm) was similar with the averaged prediction (∼29.32 Pg C) of the ten ESMs (supplementary figure S6). However, almost all of these ten ESMs failed to capture the spatial pattern of SOCs distribution [38], our outputs could be provided as a benchmark for future carbon model prediction on the TP.

Though the larger proportion of SOC on the TP compare to its area (figure 6(a)), the sensitivity of which to warming such as Q₁₀ was found to be ∼20.7% larger than the average level of China (2.32–1.93) (figure 6(b)), suggesting its vulnerability in response to climate warming on the TP. Process-based model (ORCHIDEE) has indicated that a 2 °C warming could invoking ∼10% of SOC stock lost (1.2 Pg) in grassland, which would make the TP a weak C source [39]. While Chen et al [40] synthesized 65 publications and found the total amount of SOC on the TP would not change under warming, the lost C induced by enhanced respiration could be compensated by increased plant biomass and dissolved OC. Concerning the uneven distribution of current warming experimental sites on the TP, more attention should be further focused on its central and west part and besides grassland to obtain a more comprehensive assessment of its SOC feedbacks to global warming.

SOC stored in the grassland was ${\text{15}}{\text{.52}}_{{\text{9}}{\text{.91}}}^{{\text{23}}{\text{.52}}}$ Pg, which accounted for 64.5% of China's grassland carbon size (24.03 Pg), indicating the large contribution (∼64%) of TP's grassland to China's grassland C stock on ∼55% of its area [34]. The large size of SOC reservoir in grassland primarily due to the slow OC decomposition rate by low temperature in most alpine grassland regions [2, 10, 19]. To make a further comparison, SOC storage of alpine grassland was estimated (14.07 Pg, 9.92 kg m⁻² on average), and the amount was almost two times the size of Yang et al [1] results for 0–1 m (7.36 Pg, 6.52 kg m⁻² on average). Additionally, we overlay a digital permafrost distribution map on our predicted maps, and found that there were ${\text{4}}{\text{.71}}_{{\text{3}}{\text{.07}}}^{{\text{6}}{\text{.87}}}$, ${\text{10}}{\text{.74}}_{{\text{5}}{\text{.83}}}^{{\text{16}}{\text{.06}}}$, and ${\text{3}}{\text{.41}}_{{\text{1}}{\text{.94}}}^{{\text{8}}{\text{.73}}}$ Pg C stored in 0–0.3 m, 0–1 m and 1–2 m in permafrost regions, respectively, which were all fell into the range of previous estimations, but mainly at a high level for each depth interval (table 1).

Differences between our estimates to previous studies could be mainly induced by three possible reasons. Firstly, the dataset was unprecedently improved for both its data volume and geographical coverage, making the data more representative and minimizing the uncertainties [10, 15]. In particular, we have appended a large number of measurements from regional field works conducted in shrub and forest ecosystems, which were rarely concerned by former estimates. Moreover, number of observations from permafrost regions were complemented recent years, which has largely increased the average value of SOC density (SOCD) in grassland than before (increased to 6.96 and 9.95 kg C m⁻² for 0–0.3 m and 0–1 m, respectively). Secondly, algorithms have a non-negligible effect on regional SOC prediction. The machine learning approach used in this study was expected to have better performance than traditional up-scaling methods [2]. A recent work by Liu et al [41] used ensemble machine learning algorithms mapping China's soil information, and achieved more accurate predictions compared with previous soil maps like the harmonized world soil database (HWSD, v1.2), which produced by linking raster to harmonized soil property data. Thirdly, differences in scale may also result in estimation deviations among vegetation types. For instance, area of alpine grassland used in this study was ∼1.42 × 10⁶ km² which was larger than 1.13 × 10⁶ km² that used in Yang et al [1].

Though machine learning has proven to have great advantages in OC mapping [2, 17], this study still has uncertainties. For biomass C stocks, the large uncertainty from forest AGB and grassland BGB suggested that additional observations in above ground forest and belowground grassland are needed in the future (figures 4(a)–(c)). For SOC stocks, the uncertainty decreased for 0–100 cm compared with 0–30 cm, which was consistent with the tendency from Ding et al [2] and Jiang et al [42]. While the uncertainty increased sharply when it concerned 100–200 cm depth, which was in accordance with the findings by Mishra et al [10], who claimed the uncertainty was much higher at 100–200 cm depth interval on the TP's permafrost region. These above uncertainties could mainly be resulted from the lack of sampling points, which again suggested the necessity of supplemental samplings for the non-sampling areas and deeper soils.

Ecological correlation does not necessarily imply causation. Temperature related factors, specifically the mean diurnal temperature range (Tdr), was shown as the primarily factor negatively controlling AGB across the TP (figure 3(a), supplementary figures S4 and S5). Usually, vegetation through physiological processes like evapotranspiration and biomass energy absorption modifying ambient air temperature during daytime, and releasing the absorbed energy to atmosphere during nighttime affecting ambient air temperature, and thus influence the Tdr. Therefore, Tdr, as an indicator, is regulated by interaction among plant and surrounding factors like solar irradiation, VPD etc [43]. On the TP, the spatial pattern of AGB was asynchronous with the Tdr, which increased from forest in southwest to barren desert steppe in west central part (figure 4(a)). From this point of view, the dominant driving factor for AGB should be water related factors on the TP (supplementary figure S4(a)). However, our result implied that Tdr can be used as an effective covariate for AGB prediction on the TP.

Our results suggested that AGB and water related variables presented as the key categories for the topsoil (0–30 cm) SOC stock, which was consistent with former studies [1, 44, 45]. SOC of topsoil on the TP was susceptible to precipitation because it could increase organic C input into the soil by enhancing plant productivity. For subsoil (30–100 cm), we found AGB and climatic variables still played as a dominant role in driving its SOC, though the accumulative contribution of soil physico-chemical properties increased (figure 3(c)). This result appears to be inconsistent with the early findings on a global scale (with ∼2700 soil profile) [19] which claimed a switched importance of SOC driving factors with depth, and emphasized the importance of edaphic factors for deep soils. However, Luo et al [46] recently used an updated global soil dataset (141 581 soil profile) detecting the controlling factors and found the importance of climatic variables did not show a decreasing trend with depth. Additionally, Mishra et al [10] used the regression kriging method (AGB was not included) also suggested the importance of climate to deep SOC which found precipitation could explain the largest portion of the SOC variation for both 0–0.3 m (47%) and 0–1 m (24%) SOC on the TP. It is important to note that the relative high SOCD the TP was mainly attributed to its cold weather. As the C output pathway was inhibited by low temperature, the mass C balance should be dominated by C input from plant which was primarily controlled by climatic variables (figure 3(a), supplementary figure S4(a)). In addition, from the data aspect, the topsoil (0–30 cm) and subsoil (30–100 cm) showed a high correlation (figure S7) and also a similar spatial pattern (figures 4(g)–(h)), therefore our findings, these two layers were synchronized on the main driving factor, should be reasonable. However, C stability in subsoils might still be controlled by edaphic related factors rather than climatic variables. A recent study used radiocarbon (Δ¹⁴C) measurements detecting controlling factors on C persistence on the TP's grassland, and found the dominant factor for SOC persistence changed to mineral properties in subsoil due to a decreased root C input [45]. While evidence from temperature control experiments have also demonstrated that deep SOC had higher sensitivity (Q₁₀) to warming than topsoils [47, 48], and was vulnerable to climate warming [49]. Given the large uncertainty and vulnerability of SOC to climate, more field observations and manipulative experiments of whole-soil profile should be paid to clarify the feedbacks of deep soils to the ongoing climate change.

For vegetation types, the dominant driving factors of topsoil SOC in grassland were coming from water related categories (supplementary table S5), which was in accordance with former studies [1, 10]. While for shrubland, the key driving factors transformed into temperature related variables. The significant negative relationship between temperature related variables (supplementary table S5) and SOC stock in shrubland in this study was consistent with the findings from shrubland on a global scale [19, 50], which indicated an increasing risk for C loss in shrubland due to warming stimulated biological degradation. The reason might be that alpine shrubs are mainly distributed on shady slopes and foothills with lower solar radiation and thus low temperature and water consumption. Therefore, water becomes no longer a limit factor and the temperature becomes the main constraint of the shrub ecosystems.

The tendency between plant productivity and SOC varies with vegetation types. The correlation coefficients decreased from grassland to shrubland and disappeared in forest (figures 5(a) and (b)). Plant productivity (i.e. AGB) in grassland was mainly controlled by water related variables (figure 3(a)), and the low temperature further constrains the microbe activities for SOM decomposition. Therefore, plant productivity was closely related to its SOC stocks in grassland, which has also been demonstrated by remote sensing observations [1, 18]. Whereas, for shrubland, temperature related variables was the main driver on its topsoil SOC which mainly affect the C loss path [50]. Therefore, the relationship between AGB, which determine C input, and SOC was weakened. For forest, only a small portion of its AGB, mainly leaves, fall on the ground and has been involved in the C cycle every year, resulting in no significant relationship between them.

Besides, the strength of the correlations between AGB and SOC weakened with depth. For subsoil (30–100 cm), only grassland exhibited a significant (p < 0.01) relationship between AGB and SOC, with a decreased R² (figure 5(b)). Which confirmed by previous findings that the importance of C input from plant decrease with increased soil depth [19, 45]. However, AGB still played a dominant role in shaping subsoil SOC distribution, indicating that the importance of plant C input for deep SOC (to 1 m) was non-negligible on the TP. The similar contribution of other environmental categories (i.e. temperature and soil properties) suggested a relative complex regulatory environment for the deep soils [51]. In sum, our findings implicated that the responses of SOC to future climate change may differs depending on interactions between vegetation types and regional climate conditions.