Soil Property, Rather than Climate, Controls Subsoil Carbon Turnover Time in Forest Ecosystems across China

Abstract

Subsoil (0.2–1 m) organic carbon (C) accounts for the majority of soil organic carbon (SOC), and SOC turnover time (τ, year) is an important index of soil C stability and sequestration capacity. However, the estimation of subsoil τ and the identification of its dominant environmental factors at a regional scale is lacking in regards to forest ecosystems. Therefore, we compiled a dataset with 630 observations to investigate subsoil τ and its influencing factors in forest ecosystems across China using the structural equation model (SEM). The results showed a large variability of subsoil τ from 2.3 to 896.2 years, with a mean (± standard deviation) subsoil τ of 72.4 ± 68.6 years; however, the results of one-way analysis of variance (ANOVA) showed that subsoil τ differed significantly with forest types (p = 0.01), with the slowest subsoil τ obtained in deciduous-broadleaf forests (82.9 ± 68.7 years), followed by evergreen-needleleaf forests (77.6 ± 60.8 years), deciduous-needleleaf forests (75.3 ± 78.6 years), and needleleaf and broadleaf mixed forests (71.3 ± 80.9 years), while the fastest subsoil τ appeared in evergreen-broadleaf forests (59.9 ± 40.7 years). Subsoil τ negatively correlated with the mean annul temperature, occurring about three years faster with a one degree increase in temperature, indicating a faster subsoil SOC turnover under a warming climate. Subsoil τ significantly and positively correlated with microbial activities (indicated by microbial C and nitrogen), highlighting the importance of microbial communities in regulating subsoil C dynamics. Climate, forest types, forest origins, vegetation, and soil variables explained 37% of the variations in subsoil τ, as indicated by the SEM, and the soil property was the most important factor affecting subsoil τ. This finding challenged previous perception that climate was the most important factor driving subsoil C dynamics, and that dominant drivers varied according to climate zones. Therefore, recognizing different dominant factors in predicting subsoil C dynamics across climate zones would improve our understanding and reduce the uncertainties regarding subsoil C dynamics in biogeochemical models under ongoing climate change.

1. Introduction

As one of the most widely distributed vegetation types on the global land surface [1], forests not only contribute most of the world’s gross primary product (GPP), but also store more carbon (C) in the biomass and soil of the forest than in the atmosphere [2]. Thus, forest ecosystems play an important role in the global C cycle, providing feedback into the climate system [3]. Forests are one of the most important net C sinks in the atmosphere [4]; however, whether this C sink will continue in terms of climate change remains largely uncertain. Therefore, studying forest C dynamics is of great significance for understanding the role of forests in global C cycling.

Currently, most of the studies are focused on the dynamics and decomposition of soil organic carbon (SOC) in topsoil within the range of 0–0.2 m [5,6]. However, subsoil (below 0.2 m) has a better C sequestration capacity than topsoil [7]. The dynamics of SOC in the subsoil may be as sensitive as topsoil to environmental changes [8], but this has not currently been proved [9]. SOC turnover time (τ, years) is not only an indicator of soil C stability, but it also determines the C sequestration capacity of soil [10]. SOC turnover time is the average time from the initial photosynthetic fixation to the respiratory or non-respiratory loss of C atoms in terrestrial ecosystems, which can be calculated from C stocks relative to C inflow or outflow [11]. However, τ has an important steady-state assumption; that is, under the steady-state assumption, C output is equal to input [12]. In fact, few soils maintain a stable state because of the effects of natural and anthropogenic factors and climate change [8]. Compared with topsoil, the soil organic C in the subsoil is more stable and less susceptible to disturbance [13], indicating that subsoil can be considered to be in a stable state [8].

Previous studies have already explored the C turnover time of soil or vegetation at regional or global scales [14,15]. A recent study showed that, using the modeling data at global scales [8], soil property was the main factor driving subsoil τ (0.3 to 1 m). However, the result was estimated based on the Harmonized World Soil Database (HWSD) and remote sensing data, and this study did not determine which environmental variable was more important for subsoil τ in different forest types and climate zones. Therefore, fundamental problems related to subsoil τ remain: firstly, questions regarding whether conclusions from field observations support those from global model data arise because the mismatch between the observed data and global grid data is a well-known problem in ecological studies due to differences in spatial resolution [15]. Secondly, there is a question concerning whether the subsoil τ of different forest types or climate zones is governed by the same environmental variable. Thirdly, because vegetation is one of the sources of soil C, the determination of how vegetation variables (such as NDVI, GPP, EVI) influence subsoil τ is a major concern. Therefore, it is important to evaluate subsoil τ and identify the main driving factors by using field observations for a reliable evaluation of subsoil C cycle-climate feedback.

Although temperature and precipitation are the most important climatic factors affecting soil τ [16], soil and vegetation variables may also influence τ. For example, the leaf area index (LAI) is the main vegetation factor affecting C turnover [17,18], and soil bulk density (BD), soil silt, and clay indirectly affect subsoil τ by affecting fine roots and C storage, respectively [19,20]. In addition to climatic and soil factors, biological traits also affect subsoil τ (such as vegetation longevity and types, as well as biomass) [12]. Different vegetation types may be accompanied by different climatic conditions, which greatly affect vegetation productivity, impacting τ [20]. Therefore, climate, vegetation and soil variables, which are the most important factors affecting τ, are still poorly understood.

In this study, we used 630 observations across China that included most of the forest types in the forest ecosystem. The main research objectives of this paper were to: (1) estimate and compare the subsoil τ of different climate zones, forest origins, forest types, and forest ages across China; (2) analyze the influence of climate, soil, and vegetation factors on subsoil τ; (3) explore the dominant environmental factors influencing subsoil τ.

2. Materials and Methods

2.1. Data Sources

The dataset used in this study came from Xu et al. (2018) [21], which was based on field measurements obtained between 2004 to 2014 from published studies of the China National Knowledge Infrastructure (CNKI, http://www.cnki.net, accessed on 10 March 2020) and the Institute for Scientific Information (ISI, http://apps.webofknowledge.com, accessed on 10 March 2020). Our search keywords when using the database were: forest, SOC (0–0.2 m), and SOC (0–1 m). Three criteria were used for data selection: (1) SOC content/concentration data should be obtained through field investigations, excluding observations predicted by other models; (2) field investigations should have been performed after 2000; and (3) methods for soil organic C determination should be comparable. More details can be found in the literature of Xu et al. (2018) [21]. Since we focused on SOC turnover time in forests, we maintained field observations of SOC in the 0–0.2 m and 0–1 m ranges, and subsoil SOC stock (0.2–1 m) was calculated as the difference between the two soil layers. Finally, the compiled dataset contained 630 observations across China, and a general data description is found in

Table A1. General descriptions of subsoil τ (year) used in this study.

Forest Types	df	Mean	Standard Deviation	Minimum	Maximum	Coefficient of Variance (%)
DBF	154	82.9	68.7	10.1	411.0	82.9%
DNF	55	75.3	78.6	2.85	312.0	104%
EBF	148	59.9	40.7	4.59	221.0	68%
ENF	53	77.6	60.8	13.0	311.0	78.4%
NBF	220	71.3	80.9	2.32	896.2	113%
Climate zones
Boreal	79	89.5	68.6	16.4	318.0	76.7%
Subtropical	285	60.1	41.0	2.32	253.0	68.2%
Temperate	192	92.2	96.8	9.25	896.2	105%
Tropical	74	49.7	45.6	4.59	221.0	91.8%
Forest origin
Natural	337	78.2	73.2	2.32	896.2	93.6%
Plantation	293	65.6	62.6	2.85	481.0	95.3%
Forest age (486 observations)
Young	233	67.5	60.3	2.32	411.0	89.4%
Middle-aged	152	71.1	57.7	12.2	315.0	81.1%
Mature	101	92.3	111	14.8	896.2	121%
All	630	72.4	68.6	2.32	896.2	94.8%

. We divided the observed data into five forest types (Figure 1): deciduous-broadleaf forests (DBF), deciduous-needleleaf forests (DNF), evergreen-broadleaf forests (EBF), evergreen-needleleaf forests (ENF), and needleleaf and broadleaf mixed forests (NBF).

2.2. Global Environmental Variables

In this study, we investigate the correlation of subsoil (0.2–1 m) C turnover with climate, vegetation, and soil variables. Climate variables included mean annual temperature (MAT) and mean annual precipitation (MAP) from the Climatic Research Unit [22], and relative humidity (RH) from the published literature [23]. The vegetation variables included the normalized difference vegetation index (NDVI), the GPP and enhanced vegetation index (EVI) from 2004–2014 from the Moderate Resolution Imaging Spectroradiometer (MODIS). The soil variables, including soil microbial biomass C (SMC) and soil microbial biomass nitrogen (SMN) data, were obtained from the published literature [24], as was the soil nitrogen content [25]. Due to the lack of SMC and SMN in the 0–0.2 m range, the difference between 0–0.3 m and 0–1 m was used to calculate the soil microbial biomass C and N.

2.3. Subsoil τ Calculation

At a steady state, the C input was equal to the output, and the subsoil (0.2–1 m) τ was calculated as follows [8]:

$$\mathsf{\tau }=\frac{\mathrm{SOC}}{{\mathrm{C}}_{\mathrm{input}}{ \mathrm{or} \mathrm{C}}_{\mathrm{output}}}$$

(1)

At a steady state hypothesis, C_input and C_output can be considered as the net primary productivity (NPP) allocated to the soil. Therefore, the calculation formula of subsoil (0.2–1 m) τ was transformed into [26]:

$$\mathsf{\tau }=\frac{{\mathrm{SOC}}_{0.2\text{–}1}}{{\mathrm{BNPP}}_{0.2\text{–}1}}=\frac{{\mathrm{SOC}}_{0.2\text{–}1}}{\mathrm{NPP}·{\mathrm{f}}_{\mathrm{BNPP}}·{\mathrm{fr}}_{0.2\text{–}1}}$$

(2)

SOC_0.2–1 was SOC 0.2–1 m, and BNPP_0.2–1 was underground NPP 0.2–1 m; f_BNPP indicates the proportion of underground NPP to total NPP. fr_0.2–1 is the root ratio within 0.2–1 m to the total roots. In our study, the mean MODIS NPP from 2001 to 2014 represented the total NPP and the extracted site specific NPP from the site coordinates.

Based on Equation (2), the NPP distribution of the subsoil was needed to calculate subsoil τ. Therefore, according to the distribution of root biomass in the soil, the total NPP was divided into 0.2–1 m subsoil [27]. The root distribution was calculated as:

$${\mathrm{r}}_{\mathrm{D}}=\frac{{\mathrm{R}}_{\max }}{1+{(\frac{\mathrm{D}}{{\mathrm{D}}_{50}})}^{\mathrm{c}}}$$

(3)

${\mathrm{r}}_{\mathrm{D}}$ is the total number of roots across soil depth D (m); ${\mathrm{R}}_{\max }$ is the estimated value of total root C in the soil profile, and ${\mathrm{D}}_{50}$ is the depth (m) where ${\mathrm{r}}_{\mathrm{D}}=0.5·{\mathrm{R}}_{\max }$. $\mathrm{c}$ calculated as [8]:

$$\mathrm{c}=\frac{-1.27875}{{\log }_{10}{\mathrm{D}}_{95}-{\log }_{10}{\mathrm{D}}_{50}}$$

(4)

${\mathrm{D}}_{95}$ is the depth (m) at ${\mathrm{r}}_{\mathrm{D}}=0.95·{\mathrm{R}}_{\max }$. Based on Equation (3), the ratio of roots in the subsoil (fr_0.2–1) was reckoned as:

$${\mathrm{f}}_{\mathrm{r}0.2\text{–}1}=\frac{{\mathrm{r}}_1}{{\mathrm{R}}_{\max }}-\frac{{\mathrm{r}}_{0.2}}{{\mathrm{R}}_{\max }}=\frac{1}{1+{(\frac{1}{{\mathrm{D}}_{50}})}^{\mathrm{c}}}-\frac{1}{1+{(\frac{0.2}{{\mathrm{D}}_{50}})}^{\mathrm{c}}}$$

(5)

fr_0.2–1 was 0.526 for DBF, 0.596 for DNF, 0.556 for EBF, 0.618 for ENF, and 0.649 for NBF [8].

2.4. Data Analysis

Before statistical analysis, the Shapiro–Wilk test was used to assess the normality of the data [28]. Since the subsoil τ we calculated was not a normal distribution, subsoil τ was naturally log-transformed to approximate normality.

Firstly, one-way analysis of variance (ANOVA) was used to analyze the subsoil τ of forest types, forest origins, climate zones and forest ages at p = 0.05. The Tukey-HSD (honestly significant difference) test was used for multiple comparisons. A two-way analysis of variance was performed for evaluate the interaction between forest types and climate zones and their effect on subsoil τ.

Secondly, a generalized linear mixed model (GLMM) was used to analyze the correlation between climate, vegetation, and soil variables and the variation of τ, because GLMM is a more flexible method for analyzing non-normal data with random effects [29]. Since the subsoil τ of stand types may produce random effects, stand types were used as random factors. Moreover, study sites and the interactions between stand types and sites may also produce potential random effects, therefore, we included sites and stand types as random effects as well (

Table A2. Statistics of GLMM relating environmental variables with subsoil τ, treating the sites as a random effect.

Variable	Slope	95% CI	R2 (m/c)
MAT	−0.014 ***	(−0.018, −0.010)	0.09/0.10
Log (MAP)	−0.367 ***	(−0.469, −0.266)	0.07/0.09
Log (RH)	−0.558 ***	(−0.684, −0.432)	0.11/0.13
Log (GPP)	−0.370 ***	(−0.474, −0.265)	0.07/0.43
Log (EVI)	−0.769 ***	(−1.048, −0.443)	0.04/0.43
Log (NDVI)	−0.816 ***	(−1.130, −0.502)	0.04/0.44
Log (Soil C:N)	0.672 ***	(0.465, 0.665)	0.44/0.56
Log (SMC)	0.485 ***	(0.362, 0.609)	0.09/0.34
Log (SMN)	0.454 ***	(0.301, 0.607)	0.05/0.32

Note: R² (m/c): variance of fixed effects explanations (R² marginal)/variance of fixed and random effects explanations (R² conditional); 95% CI; 95% confidence interval. *** p < 0.001.

and

Table A3. Statistics of GLMM relating environmental variables with subsoil τ, treating the stand types by sites as a random effect.

Variable	Slope	95% CI	R2 (m/c)
MAT	−0.016 ***	(−0.020, −0.012)	0.10/0.23
Log (MAP)	−0.367 ***	(−0.469, −0.266)	0.07/0.07
Log (RH)	−0.601 ***	(−0.739, −0.461)	0.12/0.12
Log (GPP)	−0.370 ***	(−0.475, −0.266)	0.07/0.07
Log (EVI)	−0.769 ***	(−1.076, −0.463)	0.04/0.04
Log (NDVI)	−0.816 ***	(−1.130, −0.503)	0.04/0.04
Log (Soil C:N)	0.675 ***	(0.628, 0.732)	0.44/0.69
Log (SMC)	0.512 ***	(0.389, 0.635)	0.10/0.12
Log (SMN)	0.521 ***	(0.365, 0.676)	0.06/0.06

Note: R² (m/c): variance of fixed effects explanations (R² marginal)/variance of fixed and random effects explanations (R² conditional); 95% CI; 95% confidence interval. *** p < 0.001.

, Figure A1 and Figure A2). Pseudo-R-squared (R²) was used to evaluate the variance of the fixed effects explanations (R² marginal), and the variance of fixed and random effects explanations (R² conditional). ${\mathrm{R}}_{\mathrm{m}}^2$ and ${\mathrm{R}}_{\mathrm{c}}^2$ were calculated using the following formula [30]:

$${\mathrm{R}}_{\mathrm{m}}^2=\frac{{\mathsf{\sigma }}_{\mathrm{f}}^2}{{\mathsf{\sigma }}_{\mathrm{f}}^2+{{\displaystyle \sum }}_{\mathrm{l}=1}^{\mathrm{u}}{\mathsf{\sigma }}_{\mathrm{l}}^2+{\mathsf{\sigma }}_{\mathrm{e}}^2+{\mathsf{\sigma }}_{\mathrm{d}}^2}$$

(6)

$${\mathrm{R}}_{\mathrm{c}}^2=\frac{{\mathsf{\sigma }}_{\mathrm{f}}^2+{{\displaystyle \sum }}_{\mathrm{l}=1}^{\mathrm{u}}{\mathsf{\sigma }}_{\mathrm{l}}^2}{{\mathsf{\sigma }}_{\mathrm{f}}^2+{{\displaystyle \sum }}_{\mathrm{l}=1}^{\mathrm{u}}{\mathsf{\sigma }}_{\mathrm{l}}^2+{\mathsf{\sigma }}_{\mathrm{e}}^2+{\mathsf{\sigma }}_{\mathrm{d}}^2}$$

(7)

${\mathsf{\sigma }}_{\mathrm{f}}^2$ is the calculation variance of the fixed effect components of GLMM; u means the number of random factors; ${\mathsf{\sigma }}_{\mathrm{u}}^2$ is the variance component of the lth random factor; ${\mathsf{\sigma }}_{\mathrm{e}}^2$ is the additive discrete variance in GLMM; and ${\mathsf{\sigma }}_{\mathrm{d}}^2$ is the distribution-specific variance in GLMM.

Finally, the effects of climate, vegetation, forest types, forest origins, and soil variables on subsoil τ were investigated using the structural equation model (SEM). As a causal reasoning tool, SEM is often used in ecological research [31]. Different from regression or ANOVA, SEM can separate multiple influencing pathways, which is an effective method for studying complex relationships in ecology [32]. At the same time, SEM can test the rationality of a hypothetical model according to priori information about the relationships between specific variables. Based on prior knowledge and theoretical knowledge, the conceptual model and hypothesis mechanism were built (Figure A3 and Figure A4).

Before SEM analysis, MAT, MAP, and RH were used to characterize the influence of climate on τ; EVI, GPP, and NDVI were used to characterize the influence of vegetation on τ, while soil variables included SMC, the ratio of SOC to soil nitrogen, and SMN. Principal component analysis (PCA) was performed to construct a multivariate indicator representing each group [33]. The first principal component (PC1) explained 63–92% of the variations in the three potential groups and was used in the structural equation model (

Table A4. Results of PCA of environmental variables.

Variable	Loading Factor
Climate
MAT (°C)	0.96
Log (MAP) (mm)	0.93
Log (RH) (%)	0.99
Cumulative variance explained (%)	92%
Vegetation
Log (NDVI)	0.97
Log (GPP) (g C m−2 year−1)	0.90
Log (EVI)	0.98
Cumulative variance explained (%)	91%
Soil variables
Log (SMC) (g C m−2)	0.82
Log (C:N)	0.50
Log (SMN) (g N m−2)	0.79
Cumulative variance explained (%)	63%

). In addition, a regular numerical covariable (e.g., 1, 2, 3, 4, 5) was used to represent forest types and forest origins in the SEM analysis, based on the research of others [34]. SEM used the maximum likelihood method to fit the model, and the Chi-square (χ²) test and the root mean square errors of approximation (RMSEA) were used to evaluate the performance of the model. When the non-significant χ² test (p > 0.05), and the RMSEA < 0.08, the result of the SEM model was acceptable [35]. All analysis was performed in R 3.6.0 [36].

3. Results

3.1. Subsoil τ Varied with Forest Types and Climate Zones

The subsoil τ varied greatly, ranging from 2.3 to 896.2 years with an average (± standard deviation) of 72.4 ± 68.6 years (Table A1). However, there were significant differences in subsoil τ among different forest types (p = 0.01, Figure 2A), the slowest subsoil τ was found in DBF (82.9 ± 68.7 year), followed by ENF (77.6 ± 60.8 year), DNF (75.3 ± 78.6 year), and NBF (71.3 ± 80.9 year), while the fastest subsoil τ appears in EBF (59.9 ± 40.7 year, Figure 2A and Table A1). Similarly, the subsoil τ of different climatic zones was significantly different (p < 0.001, B), with the slowest obtained in temperate areas (92.2 ± 96.8 year) and the fastest in tropical areas (49.7 ± 45.6 year), while boreal and subtropical zones exhibited intermediate τ (Table A1). Two-way ANOVA showed that the interactions between forest types and climatic zones had a significant impact on subsoil τ (p < 0.001,

Table 1. Two-way ANOVA of forest types, climatic zones, and their interaction effect on subsoil τ.

	df	SS	MS	F Value	p
Forest types	4	15.9	3.96	3.67	<0.01
Climatic zones	3	55.4	18.45	17.09	<0.001
Forest types: Climatic zones	10	33.1	3.31	3.07	<0.001

Note: SS—sum of squares; MS—mean square; df—degree of freedom.

). In regards to forest origins and stand ages, the subsoil τ in natural forests was slower than that in plantations (p = 0.002, Figure 2C), and the subsoil τ in mature forests was slower than that in young forests (p = 0.01, Figure 2D).

3.2. Influencing Factors on Subsoil τ

GLMM analysis showed that there was a significant negative correlation between climate variables (MAT, MAP, and RH) and subsoil τ (p < 0.001, Figure 3A–C and

Table 2. Statistics of GLMM relating climatic, vegetation, and soil variables with τ.

Variable	Slope	95% CI	R2 (m/c)
MAT	−0.016 ***	(−0.020, −0.012)	0.10/0.12
Log (MAP)	−0.366 ***	(−0.468, −0.265)	0.07/0.07
Log (RH)	−0.560 ***	(−0.740, −0.460)	0.12/0.13
Log (GPP)	−0.370 ***	(−0.475, −0.265)	0.07/0.07
Log (EVI)	−0.770 ***	(−1.076, −0.463)	0.04/0.04
Log (NDVI)	−0.817 ***	(−1.130, −0.503)	0.04/0.04
Log (Soil C:N)	0.672 ***	(0.615, 0.730)	0.44/0.56
Log (SMC)	0.509 ***	(0.386, 0.632)	0.09/0.12
Log (SMN)	0.516 ***	(0.361, 0.672)	0.07/0.10

Note: R² (m/c): variance of fixed effects explanations (R² marginal)/variance of fixed and random effects explanations (R² conditional); 95% CI; 95% confidence interval. *** p < 0.001.

). Vegetation variables (GPP, NDVI, and EVI) were also significantly correlated with subsoil τ, leading to a negative impact on subsoil τ (p < 0.001, Figure 3D–F). Subsoil τ increased with the increase in C:N, SMC, and SMN (p < 0.001, Figure 3G–I).

In China’s terrestrial ecosystem, climate, forest types, forest origins, vegetation, and soil variables explained 37% of the variations in subsoil τ, as indicated by the SEM (Figure 4). Regardless of forest types, the SEM showed that soil property was the most important variable, exerting a directly positive impact on subsoil τ, followed by vegetation, climate, and forest origins, which all have negative effects on subsoil τ. However, the relationship between forest types and subsoil C turnover time was not significant. Climate has a significantly positive effect on vegetation and forest types, and vegetation variables caused exerted a positive influence on soil variables and forest types, but exhibited a negative effect on forest origins.

SEM models were also established for different forest types and climatic zones (Figure A5, Figure A6, Figure A7, Figure A8, Figure A9, Figure A10, Figure A11, Figure A12 and Figure A13). Although the SEM results of different forest types showed that the main controlling factors were soil properties (Figure A5, Figure A6, Figure A7, Figure A8 and Figure A9), the subsoil τ of different climatic zones was driven by different environmental factors (Figure A10, Figure A11, Figure A12 and Figure A13). Climate was the main driving factor for temperate forests, and soil property was the primary influencing factor for boreal, tropical, and subtropical forests.

4. Discussion

4.1. Subsoil τ in China’s Forests

The mean subsoil τ in China’s forests was 72.4 ± 68.6 years, which was slower than the subsoil τ (0–1 m) of tree plantations (30 years) and China’s forests (17.7 years) [16,37], but faster than the subsoil τ (0.3–1 m) of global forests (168 years) [8]. The main reasons for this difference include the following three points. Firstly, Wang et al. (2017) [16] estimated soil τ within 0–1 m, while Luo et al. (2019) [8] defined subsoil ranging from 0.3 to 1 m from global forests. We defined the subsoil as 0.2–1 m, and observational data from Chinese forests were used. The differences in these results indirectly indicated that SOC in deeper soils had a slower turnover time compared to that in topsoil. Secondly, the distribution of NPP in different soil depths was not clearly estimated [16]. Thirdly, the sources of the dataset and the estimated methods of subsoil τ also yielded different results. For example, soil τ was defined as the ratio of SOC and the CO₂ flux from SOC decomposition [37], while we defined subsoil τ as the ratio of SOC and NPP. For these reasons, our results are not directly comparable to those of others. Therefore, this study may be the first using field observation data to estimate forest subsoil τ across China, which can further increase our understanding of subsoil C dynamics.

4.2. Influence of Environmental Factors on Subsoil τ

4.2.1. Climatic Effects

Our results demonstrated that subsoil τ decreased with the increase in MAT (Figure 3A, p < 0.001), which was inconsistent with previous results [38], which hypothesized that the subsoil τ had no significant correlation with temperature in global forests. In addition, our study further quantified the effects of temperature on subsoil τ and found that the mean of subsoil τ decreased by three years with a one degree increase in temperature (Figure A14), which indicated that subsoil τ will potentially become faster under a warming climate [39]. This negative correlation can be partly attributed to the fact that firstly, increasing MAT will accelerate the rate of SOC mineralization by stimulating soil microbial or enzyme activity, thus making the subsoil τ faster [40]. Secondly, NPP allocated to the subsoil will increase with the increase in temperature (Figure A15), leading to a faster subsoil τ. We also found that MAP had a significantly negative impact on subsoil C turnover in forest ecosystems across China (Figure 3B, p < 0.001), which may be due to the accelerated decomposition of soil C by increasing precipitation or soil water content (SWC) [41], in turn, making the subsoil τ faster.

4.2.2. Vegetation Effects

In our study, subsoil τ decreased with the increase in vegetation variables (GPP, NDVI, EVI) (Figure 3D–F, p < 0.001), which was similar to previous research results that vegetation could affect subsoil C turnover time [42,43]. The influence of vegetation on subsoil τ mainly included three aspects. (1). Vegetation can indirectly affect the subsoil τ by changing photosynthate [42]; Street et al. (2020) [42] also found that photosynthate drives SOC mineralization in mature forests (>50 years old), which may have an effect on subsoil τ. (2). Vegetation type affects the subsoil C turnover time by changing the litter quality; the better the litter quality, the faster the subsoil τ [44]. (3). Vegetation could also affect subsoil τ by changing microclimate, soil water content, and forest structure [45]. However, in the results of our study, subsoil τ of DBF with better litter quality was slower than EBF, which was related to the distribution of the sample points. In the terrestrial ecosystem of China, EBF was mainly distributed in subtropical areas, while DBF was widely distributed (Figure 1). However, in a fixed study area (such as subtropical areas), the subsoil τ of DBF was 55.9 years, which was faster than that of EBF (59.9 years).

We also found that forest age was an important biotic factor determining subsoil τ (Figure 2D), with the fastest subsoil τ observed in young forests and the slowest in mature forests (Table A1). This result may be associated with C allocation patterns at different forest ages [46]; the C storage in vegetation increases with the development of forest age, which affects the C turnover in the subsoil [47]. In addition, the decrease in stand NPP with the increase in tree age may be another reason [48], which is caused by the decrease in soil nutrients and photosynthesis during vegetation development [49].

4.2.3. Soil Effects

Our study also found that soil variables have significant effects on subsoil τ (Figure 3G–I, p < 0.001), which is consistent with the results of previous studies. Soil properties (including physical and chemical properties) directly or indirectly affect the subsoil τ by influencing soil structure or C input and output. For example, soil C:N was significantly positively correlated with subsoil C turnover [50]. Soil clay content was an important variable affecting subsoil τ by influencing SOC storage [51]. It is generally believed that the changes in soil nutrients and water have significant effects on vegetation growth [52], but soil nutrients and water are also affected by soil chemical and physical properties, especially in arid areas [53]. Therefore, soil properties have an indirect effect on vegetation growth, and thus affect C input and output [8]. We also found a significant increasing trend between subsoil τ and SMC and SMN (Figure 3H,I), because soil microbial properties were also important factors affecting C turnover in the subsoil [16]. In addition, microbial activity can predict the change in soil properties and thus influence the C cycle [54,55].

According to SEM, soil variables (SMN, SMC, Soil C:N) were the dominant influencing factor on subsoil τ (Figure 4), which correlated with recent results showing that soil property was the most important environmental variable of subsoil τ [8]. The GLMM results also showed that soil properties (SMN, SMC, and soil C:N) had the greatest influence on the subsoil τ. Many studies found significant differences in C turnover time among different forest types [5,43]. However, forest types were not the main factor determining SOC turnover time, as indicated by the SEM model (Figure 4), because forest types and their distributions were coupled with climate, soil, and their interactions. On the other hand, our SEM results showed that although soil properties were the main controlling factors for different stand types, the main driving factors were different in different climatic zones (Figure A5, Figure A6, Figure A7, Figure A8, Figure A9, Figure A10, Figure A11, Figure A12 and Figure A13), indicating that different dominant environmental drivers for subsoil τ among different climate zones should be considered to estimate soil C dynamics, particularly for subsoil. In our study, soil property was the main driving factor of subsoil τ in subtropical and tropical regions. However, climate was the dominant driving factor for temperate areas, which was consistent with previous findings, where soil C turnover times were driven by climate in forest ecosystems [16].

4.3. Implications for Biogeochemical Modeling

Our study provided evidence for the causes of subsoil τ in forest ecosystems across China, which has significant implications for biogeochemical models to study subsoil τ at regional or global scales. Firstly, subsoil τ played an important role in evaluating soil C dynamics. The size of subsoil C storage can be determined by soil C input and output. However, subsoil τ was not well characterized in the existing biogeochemical models, leading to great uncertainties regarding the estimation of the SOC in biogeochemical models [14]. Therefore, the results of this study may provide a reference for model parameterization.

Secondly, soil C dynamics were expressed in terms of plant functional types, rather than forest types, in biogeochemical models, which might hinder our study of the relationships between soil C and ecosystem processes in different forest types, especially subsoil. Although the forest type is different, they may exhibit the same plant functional type. Thus, the C dynamics of each forest type are unclear. Our results showed that the subsoil τ of forest types was different (Figure 2A). Therefore, an accurate representation of vegetation types in biogeochemical models was a necessary condition for predicting soil C dynamics under continuous climate change, and climate change will specifically affect the structure of vegetation [45,56], thus ultimately affecting C turnover.

Thirdly, soil variables that influencing subsoil τ have not been adequately represented in the existing biogeochemical models. Considering that soil τ was usually affected by climate [57,58], many models focused on studying and analyzing the effects of soil temperature on soil C dynamics [59]. Our results showed that soil variables were the dominant driving factor of subsoil τ (Figure 4). Therefore, we suggest that more consideration should be given to the influence of soil variables on the dynamic changes of C, such as soil microbial biomass C and soil microbial biomass nitrogen, when establishing the biogeochemical model.

4.4. Uncertainties and Limitations

Although we obtained some results by analyzing the effects of climate, vegetation, and soil variables on the subsoil τ of forest ecosystems across China using field observations, we also acknowledge that there were still some uncertainties and limitations. On the one hand, since subsoil NPP is rarely observed, we had to estimate the distribution of NPP in the range of 0.2–1 m indirectly, based on the distribution of root biomass in the soil profile by using MODIS NPP. However, the distribution of root biomass may be different considering the soil profile [8]. Thus, the accuracy of the root biomass and productivity measurements and its distribution in the subsoil would be an important step in understanding the SOC turnover time and subsoil C dynamics. However, the scale mismatch in the spatial resolution between the measured SOC and the MODIS NPP was another uncertainty. This kind of scale mismatch presents an ongoing challenge in the study of large-scale ecology.

On the other hand, we calculated subsoil τ as the ratio of soil total C storage to NPP and did not consider SOC turnover time for different factions of SOC. However, it has been well documented that subsoil SOC is still a mixture of decadal cycling SOC and stable, slower cycling SOC, with the turnover times more than several thousand years [60]. Moreover, DOC is an important source of subsoil SOC [9,61], which may affect the stability and dynamics of subsoil SOC [62,63]. Therefore, estimating C turnover times for different soil C fractions, e.g., low-density, high-density, and non-oxidizable C [60], is of great significance for improving our understanding of subsoil dynamics.

5. Conclusions

SOC turnover time and its influencing factors were estimated based on field observations in forest ecosystems across China. The robust conclusions include: firstly, the subsoil τ (0.2–1 m) varied from 2.3 to 896.2 years, with an average of 72.4 ± 68.6 years, and it differed significantly among different forest types (p = 0.01), 82.9 ± 68.7 years for DBF, 75.3 ± 78.6 years for DNF, 59.9 ± 40.7 years for EBF, 77.6 ± 60.8 years for ENF, 71.3 ± 80.9 years for NBF, as well as climate zones (p < 0.001), 89.5 ± 68.6 years for boreal, 60.1 ± 41.0 years for subtropical, 92.2 ± 96.8 years for temperate, and 49.7 ± 45.6 years for tropical. The subsoil τ in natural and mature forests was slower than that in plantations and young forests. Secondly, the GLMM showed that there was a significant negative correlation between subsoil τ and MAT and MAP, indicating the sensitivity to climate change. Thirdly, the results of the SEM model indicated that soil variables (SMN, SMC, and C:N) were the most important controlling environmental factors of forest subsoil τ change on a national scale, which was different from the previous perception that climate was the main dominant factor for soil τ. Finally, the dominant driver differed with climate zones (boreal, temperate, subtropical, tropical), which could improve our understanding of subsoil τ. This study further emphasizes the importance of soil variables, particularly for soil microbial properties, on subsoil τ.