Introduction

Forests have been proposed as a natural-solution to enhance carbon (C) sequestration in terrestrial ecosystems1,2 because they cover ~4200 million hectares of the world’s land surface3, and 50–90% of the total annual C flux of terrestrial ecosystems occurs at the forest-atmosphere interface4. However, primary productivity and other biological processes in forest ecosystems are widely constrained by soil nutrient availabilities especially by soil nitrogen (N) and phosphorus (P)5,6,7,8. Microbial metabolism serves as the base of detrital food-webs, mediates soil nutrient cycles and releases soil carbon (C) to atmosphere9. Biogeochemical constraints on catabolic processes, relative to primary production, determine the magnitude of soil respiration and the retention of detrital organic C in soils9,10. Therefore, nutrient limitation to microbial metabolism is a potential control in soil C turnover and release in forest ecosystems9,11. It is especially relevant to quantify the potential of forests to capture C in their soils as restoration efforts are expected to plant millions of trees over the next few decades.

Previous studies suggest that microbial metabolism is commonly limited by P in tropical forests12,13, while limited by N in temperate forests14,15. However, recent studies suggest that soil microbial metabolism can be limited by P in non-tropical zones16,17, and by N in subtropical soils6. In general, forest N and P availability vary with distance from the equator, plant productivity, climate, and soil depth7,18,19. Accurately predicting the role of soil nutrients in controlling microbial activity must account for all these environmental factors simultaneously. For example, C inputs and nutrient availability are generally lower in subsoil than topsoil so that limiting factors for microbial metabolism often vary with depth within a single profile17,19. Therefore, exploring the patterns and drivers of microbial metabolic limitation in soil profiles across large geographical scales could help decipher the effects of both local and geographical controls.

The forested areas of China account for 5% of the world’s forests and currently store over 4.75 Pg of C20. They are experiencing rapid revegetation in recent decades and their soils are considered significant C sinks with large C sequestration potential21,22. Although several studies have estimated pervasive soil P limitation to aboveground plant production in this region7,8, comprehensive estimates of belowground microbial nutrient limitations are unavailable. Ecoenzymatic stoichiometry based on the activities of extracellular enzymes (ecoenzymes), the proximate agents of microbial metabolism, connects the metabolic theory of ecology and ecological stoichiometry theory to predict microbial nutrient assimilation and growth9, and thus is an appropriate method to predict microbial metabolic limitation at a community level in soils17,23,24,25,26.

Here, we used ecoenzymatic stoichiometry methods to investigate patterns of microbial metabolic limitations in contrasting forests across a wide latitudinal gradient, (from tropical to cold regions) and at different soil depths (organic (O), eluvial (A), and parent material (C) layers) in China. Our results revealed widespread soil P limitation to microbial metabolism across Chinese forest ecosystems with the lowest limitation in warm-temperate regions. Significant correlations between microbial and plant P limitations over a wide geographic scale indicate synchronous P constraints in both microbial metabolism and primary productivity. These results imply that soil P availability is a long-term and inherent constraint on the turnover of soil organic C and the potential of ecosystem C sink in Chinese forests.

Results and discussion

Patterns of microbial metabolic limitation in soil profiles and latitudes across Chinese forests

Soil samples from 181 sites of Chinese forests showed microbial P limitation (microbial N/P limitation value > 0) at most locations (Fig. 1a), from tropical to boreal forests, with significant variability across soil horizons, vegetation types and climate zones (P < 0.001; Fig. 1, Supplementary Tables 4 and 5). Although 31% of all sampling sites from warm-temperate to cold-temperate zones were N limited (microbial N/P limitation value <0) (Supplementary Table 6), P limitation was pervasive across these Chinese forests and had much higher magnitude than N limitation. We thus focused on microbial P rather than N limitations in subsequent discussions.

Fig. 1: Spatial patterns of relative C limitation and N/P limitation of soil microbial metabolism in Chinese forests.
figure 1

a Geographical distribution of 181 Chinese forest study sites, the effects of both climate types and soil layers on relative C limitation and N/P limitation. b, c Identified relationships between C-acquiring enzyme activities and N-acquiring or P-acquiring enzyme activities in Chinese forests (n = 1520) via standardized major axis (type II) regression. The 1:1 dotted line is gray. *** represents significant relationships at P < 0.001. d, f The effects of soil layers within forest types on relative C limitation and N/P limitation. e, g The overall effects of soil layers on relative C limitation and N/P limitation in Chinese forests. The effects of climate, forest type and soil layers on microbial metabolic limitations were analyzed via linear mixed-effects models. *** represent significant effects at P < 0.001. The asterisk above the bar and box charts represents the significance of the effects of soil layers within climate (a) and forest (d, f) types on relative C limitation and N/P limitation. In the box charts: the top of box is 75th percentile, the center of box is median, the bottom of box is 25th percentile, and the upper and lower ends of the vertical line are the maximum and minimum data values, respectively. Soil layers: O, organic horizon; A, eluvial horizon; C, parent material horizon. Vegetation types: CF coniferous forest, CBF coniferous-broadleaf mixed forest, BF broadleaf forest. C-acquiring enzymes: β−1,4-glucosidase (BG) and cellobiohydrolase (CBH)), N-acquiring enzymes: β−1,4-N-acetylglucosaminidase (NAG) and L-leucine aminopeptidase (LAP), and P-acquiring enzyme: alkaline or acid phosphatase (AP). The full names of all 31 mountain forests sampled (DX, XX, ME, etc.) are listed in Supplementary Data.

Within soil profiles, microbial P limitation significantly increased with soil depth across climate zones and vegetation types, with the O-layer having the lowest microbial P limitation (P < 0.001; Fig. 1a, f, g), consistent with a recent study by Jing et al.17. In contrast, relative C limitation was significantly higher in the surface soils than in the subsoils for all vegetation types and climate zones (P < 0.001; Fig. 1a, d, e). Because relative C limitation indicates the hydrolysis potential of organic C rather than the metric of absolute microbial C limitation25, this result suggests that microorganisms in surface soils have stronger potential to drive organic C decomposition than in subsoils. This profile pattern of relative C limitation can be attributed mainly to higher organic C content in surface than subsoil layers providing more organic C to support microbial metabolism (Supplementary Fig. 4).

Geographically, piecewise linear regression analyses identified significant latitudinal breakpoints at about latitude 27°N and 41°N for relative C limitation and N/P limitation, respectively (P < 0.001; Fig. 2a, b), although these breakpoints have obvious variations among soil layers (Fig. 2c). Relative C limitation significantly increased from 19°N to 27°N and then decreased from 27°N to 54°N, suggesting the highest potential for organic C decomposition in the subtropical zone (latitude 24–30°N) (Fig. 1a). An opposite latitude pattern in N/P limitation indicated higher microbial P limitation in tropical and cold-temperature zones, with the lowest microbial P limitation in the warm-temperate zone (latitude 41–42°N) (Fig. 1b). These results suggested that soil P limitation to microbial metabolism remains unexpectedly high at relatively high latitudes in Chinese forests. It provides novel evidence that P limitation to soil microorganisms might be far more important than expected for boreal forests, compared to forests in other regions of the world usually thought to be limited by N27,28.

Fig. 2: The latitudinal patterns and breakpoint estimation for soil microbial metabolic limitations across Chinese forests.
figure 2

a, b The latitudinal breakpoints for microbial metabolic limitations were estimated using piecewise regression analyses, and relationships of relative C limitation (a) or N/P limitation (b) with latitude in Chinese forests were identified using generalized linear models. The shaded circles indicate breakpoint latitudes, while shaded areas are the 97.5% confidence intervals of the breakpoints. Solid black lines indicate model fits between relative C limitation or N/P limitation and latitude. c Estimated latitudinal breakpoints for soil microbial metabolic limitations among different soil layers across Chinese forests.

Drivers of microbial P limitation in Chinese forests

Increased microbial P limitation with soil depth in our study seems to contradict the prevailing idea of soil parent material being the main source of P29,30. However, two mechanisms could generate this vertical pattern of microbial P limitation in our study. First, soil organic matter is a potential pool of available P9,31. Higher organic matter decomposition reflected by higher relative C limitation in the O-layer than in A- and C-layers (Fig. 1a, d, e) suggests higher P release from organic matter in the O-layer compared to A- and C-layers. Thus, organic matter decomposition could help alleviate microbial P limitation in the O-layer in a P-deficient environment. Second, differential uptake of P by root systems with depth could be another mechanism driving vertical patterns of microbial P limitation. Vegetation index showed positive and increasing effects on microbial P limitation from O- to C-layers (Fig. 3d–j), suggesting an increased competition in P-acquisition between plants and microbes with increasing roots densities from O- to C-layers. A consistent pattern of P limitation between plants and microbes across this large geographical scale and soil layers (Fig. 4a–c) supports this argument and suggests stronger nutrient competition of plants with microbes in subsoils with abundant root systems than in topsoils32,33. Our previous study on rhizospheres also found that plants and microbes show strong competition for P when microbes are P constrained34. Therefore, P competition between roots and microbes could contribute to shaping these vertical patterns of microbial P limitation.

Fig. 3: Driving factors of soil microbial N/P limitation across Chinese forests.
figure 3

ac Variation partitioning analysis showing the effects of environmental drivers (climate, vegetation and soil) on microbial N/P limitation in different soil layers. Soil layers: O, organic horizon; A, eluvial horizon; C, parent material horizon. Climate includes aridity index (AI), mean annual precipitation (MAP), mean annual temperature (MAT), precipitation seasonality (PSEA), temperature seasonality (TSEA), acid rain (pH) (Acid rain), total N deposition (N deposition). Vegetation includes vegetation type and vegetation index (NDVI). Soil includes soil bulk density (BD), soil pH (pH), soil organic C (SOC), soil total N (TN), soil total P (TP). df Random-forest models identifying the contributions of each predictor variable in explaining microbial N/P limitation in different soil layers. MSE, mean squared error. * represents significant effects at P < 0.05, ** represents significant effects at P < 0.01, *** represents significant effects at P < 0.001. gj Partial least squares path modeling disentangling major pathways of the influences of climate, vegetation, and soil properties on microbial N/P limitation in different soil layers (gi), and the total effects of these variables on microbial N/P limitations (j). Blue and yellow arrows indicate positive and negative flows of causality, respectively. Numbers on the arrow indicate significant standardized path coefficients. R2 indicates the variance of dependent variable explained by the models. * represents significant effects at P < 0.05; ** represents significant effects at P < 0.01; *** represents significant effects at P < 0.001.

Fig. 4: Relationships between microbial N/P and plant N/P limitations in different soil layers across Chinese forests.
figure 4

a O-layer, soil organic horizon; b A layer, soil eluvial horizon; c C layer, soil parent material horizon. Solid line indicates model fit between microbial N/P limitation and plant N/P limitation (n = 55), while gray areas are the 95% confidence intervals of the model. *** represents significant relationships at P < 0.001.

Overall, we found that the combined effects of climate, vegetation and soil factors played the most important roles in latitudinal patterns of microbial P limitation across soil profiles (Fig. 3a–c). Random-forest models and partial least squares path modeling further found that vegetation index driven by climate conditions (including aridity index and mean annual precipitation) has the largest positive effects on microbial P limitation (Fig. 3d–g). These results suggested that primary productivity has a central role in contributing to the geographic pattern of microbial P limitation. Many processes relevant to plant production can limit the availability of soil P in forests including biomass accumulation and its effect on soil pH6,30. The most direct process is that the fixation of P in biomass by plant growth reduces the availability of soil P, which is supported by the consistent pattern we found between P limitation in plants and soil microbial communities (Fig. 4a–c).

Highly productive ecosystems often develop on acidic soils, which can further decrease P availability for plants and microbes. This is especially important in warm, wet ecosystems in the tropics12. However, high P limitation at high latitudes can result from cold and dry conditions that suppress microbial activities and decomposition19,35 (Fig. 2a), subsequently slowing P release from organic matter9. In addition, P limitation should be exacerbated by anthropogenic N deposition (Fig. 3 and Supplementary Figs. 13), because it can create an imbalance in soil N:P stoichiometry and increase soil acidification36. Indeed, some low latitude forests around the world, especially in China, are receiving high N deposition due to geographic distributions of industry and human populations37,38,39, which implies an even greater increase in microbial P limitation at lower latitudes in the future.

Implications and uncertainties

Our recent incubation experiments found that P addition significantly alleviated microbial P limitation, and resulted in an 11.5–29.1% increase in C emissions from organic C decomposition40. This suggested that constraining the decomposition of organic matter under P limitations can potentially decrease soil C loss. In plant-microbe systems, however, soil P limitation constrains not only microbial metabolism but also plant growth18,41. In contrast, soils with sufficient P supply could maintain both high primary productivity and high soil C emissions via microbial decomposition. This argument is also supported in part by the counterbalance hypothesis of plant-microbe interaction that purports a homeostatic mechanism stabilizing microbial metabolism in relation to gross primary production wherein the greatest release of C by soil microorganisms occurs in ecosystems having the highest levels of plant production41. Therefore, soil P supply may be a key driver of C turnover in forest ecosystems, i.e., low turnover rate under P limitation and higher rate with sufficient P, whereas the potential soil C sink depends on relative availabilities of soil P to plants versus microbes.

The consistent effects of soil P limitation on both microbial metabolism and primary productivity in our study provide cues to explaining the “soil C saturation” hypothesis (i.e., soils have a limited capacity to store C as organic inputs increase)42 (Fig. 4a–c). Recent studies found that tree plantings in organic soils did not increase net C sequestration after 12 to 39 years43, and CO2 enrichment did not increase C sequestration at the ecosystem level in mature forests44. An inherent constraint suggested by our study is that soil P availability simultaneously limits both plant growth and microbial catabolism, thus constraining both ecosystem C balance and soil C sink. Consequently, forest ecosystems may not sequester as much C from vegetation restoration or afforestation as expected due to inherent P constraints.

Although our study reveals geographical and vertical profile patterns of microbial metabolic limitations in the forest soils of China, the present study has uncertainties. First, the ecoenzymatic stoichiometry method we used is based on the activities of relatively few enzymes to represent microbial nutrient acquisition, and is bound to include uncertainties as suggested by several recent studies45,46. Second, we used plant nutrient limitation data based on geographic coordinates that corresponded to the microbial nutrient limitations derived in our study to explore their spatial couplings, but we did not account for likely temporal and spatial differences in both sets of data. The lack of precise space-time pairing of sites limits our knowledge of how nutrient limitation affects plant-microbial interaction on micro-scales such as in the rhizosphere32. Third, global change and human activity are increasing uncertainties in the coupling of soil nutrient constraints to ecosystem C cycles43,47. For example, different adaptabilities of plants and microbes to environmental changes41,47 means that N deposition may affect plant growth and microbial catabolism in different ways and magnitudes through changes soil nutrient limitations. This may result in a dynamic non-equilibrium state of the C cycle in microbial-plant systems, resulting in a large C source or sink effect within any terrestrial ecosystem. Regardless of these uncertainties, our study is of direct relevance to current prediction of forest C sinks, and has significant implications for future study and policy decisions in forest restoration.

Materials and methods

Study area and soil sampling

We sampled almost all accessible mountains having Chinese forest ecosystems (Fig. 1a and Supplementary Data), which span a wide geographic range (18.9°–53.5°N, 101.0°–129.7°E, 210-3830 m above sea level). These forests are found in national parks or nature reserves with limited human interference, and include broadleaf forests (BF), coniferous-broadleaf mixed forests (CBF) and coniferous forests (CF). Climates include tropical, subtropical, warm-temperate, temperate, and cold-temperate zones. The geographical coordinates of each sample site were recorded using a global positioning system (GPS) device (eTrex Venture, USA). Mean annual precipitation (MAP) and mean annual temperature (MAT) ranged from −5.91 °C to 23.18 °C and 242 mm to 2667 mm, respectively.

Field sampling was conducted from July 2012 to March 2013. At each selected forest, one to fifteen sampling sites were located along an altitudinal transect based on the range of community types available (Supplementary Data). At each site, three 10 m × 10 m plots located 50 m away from each other were established. Five samples of each O, A and C-horizon were collected from each plot and pooled as one composite sample. In total, 504 soil samples were collected from 181 sites (some sites had only 1 or 2 soil horizons). After removing roots, litter, debris, and stones, composite soil samples were divided into two subsamples. One subsample was stored at 4 °C for later analysis of enzymatic activity. Because extracellular enzymes (i.e., ecoenzymes) are the products of long-term rather than instantaneous metabolism of microorganisms and they are relatively stable in soils9, their seasonal variability could be much smaller than the differences between environmental gradients. Thus, the enzymatic activities from our study would be representative of the large-scale sampling across Chinese forests. The other subsample was passed through a 2 mm sieve and air-dried for physicochemical analysis. Soil bulk density (BD) was measured using a soil bulk sampler and a stainless-steel cutting ring (5 cm in diameter by 5 cm long), adjacent to the soil sampling points in each plot.

Environmental factors

Our study included 14 environmental variables (climate, vegetation, and soil), obtained either from field investigations or from satellites and databases. MAT, MAP and both temperature and precipitation seasonality (TSEA and PSEA) were obtained from the Worldclim database (https://www.worldclim.org; ~1 km resolution). The aridity index of each site was estimated by the ratio of MAP to potential evapotranspiration (PET) (aridity index = MAP/PET). The MAT, MAP and aridity index values ranged from 5.91 °C to 23.18 °C, 242 mm to 2667 mm, and 0.24 to 2.42, respectively. We obtained bulk N deposition rates and acid rain data for each region in China from 1980 to 2015 from Yu et al.39 and the National Earth System Science Data Center (http://www.geodata.cn), respectively. Vegetation index (NDVI, 2001-2015) was obtained from the Advanced Very High Resolution Radiometer (AVHRR) developed by the Global Inventory Modeling and Mapping Studies (GIMMS) group (https://ecocast.arc.nasa.gov/data/pub/gimms/3g.v1/) at a 1/12° resolution.

Soil properties of samples (soil pH, soil organic C (SOC), total N and total P) were measured for each horizon (O, A, C) using standardized protocols according to methods described in a previous study34. Each sample was measured in triplicate.

Measurements of microbial metabolic parameters

The activities of two extracellular C-acquiring enzymes (β−1,4-glucosidase (BG) and β-d-cellobiosidase (CBH)), two extracellular N-acquiring enzymes (β−1,4-N-acetylglucosaminidase (NAG) and L-leucine aminopeptidase (LAP)), and one extracellular organic P-acquiring enzyme (alkaline or acid phosphatase (AP)) (Supplementary Table 1) were determined using the method of Saiya-Cork et al.48 and German et al.49. Enzyme activities were expressed as nanomoles of substrate released per hour per gram of dry soil (nmol g−1 h−1). The experimental procedures are described in our previous study16.

Calculating microbial metabolic limitations

Microbial metabolic limitations were quantified by calculating the vector lengths and angles of enzymatic activity based on untransformed proportional activities (e.g. [BG + CBH]/[BG + CBH + NAG + LAP]). Vector length, representing relative C limitation of microbial metabolism, was calculated as the square root of the sum of x2 and y2 (Eq. 1), where x is the relative activity of C vs. P-acquiring enzymes, and y is the relative activity of C vs. N-acquiring enzymes25. Vector angle, representing microbial N/P limitation, was calculated as the arctangent of the line extending from the plot origin to point (x, y). The relative C limitation to microbes increases with vector length whereas a vector angle of 45° roughly represents the boundary between microbial P limitation (> 45°) and N limitation (< 45°). This definition is based on the global ecoenzymatic stoichiometry law for C:N:P activities near 1:1:1 at a global scale, i.e., the slopes of their regression relationships are close to 45°23,24. In the present study, the slopes of the relationships between enzyme activities for C vs. P and C vs. N acquisition were very close to 1.0 for all sites (Fig. 1b, c), indicating that the ecoenzymatic stoichiometry method and the 45° boundary of N/P limitation using in Chinese forests is reasonable. The relative strength of microbial P limitation increases with vector angle and strength of microbial N limitation increases as the angle decreases.

$${{{{{\rm{Relative}}}}}}\; {{{{{\rm{C}}}}}}\; {{{{{\rm{limitation}}}}}}={Lengt}h=\sqrt{({x}^{2}+{y}^{2})}$$
(1)
$${{{{{\rm{N}}}}}}/{{{{{\rm{P}}}}}}\; {{{{{\rm{limitation}}}}}}={Angle}\,({\,\!}^\circ)={{{{{\rm{atan}}}}}}2(x,y)\;\times 180/{{{{{\rm{\pi}}}}}}-45{{\,\!}^\circ}$$
(2)

Considering potential competition of microbes with plants for soil N and P33, we further linked our predictions of microbial metabolic limitation with the plant nutrient limitations from Du et al.8 to evaluate the relationships of microbial N/P limitation with plant N/P limitation. We extracted the ln-transformed NREdom/PREdom (ln NREdom/PREdom) from the database of Du et al.8 on the same or adjacent sites to our sampling sites. In total, values of ln NREdom/PREdom from 55 matched sites were obtained with a wide latitudinal range across Chinese Forests (21.9°–51.6°N). Du et al.8 defines ln NREdom/PREdom > 0 indicating plant N limitation and the value < 0 indicating plant P limitation. To facilitate interpretations in the present study, we redefined the plant nutrient limitation as follows: plant N/P limitation = −(ln NREdom/PREdom). Thus, in our study, the value of plant N/P limitation > 0 indicates plant P limitation, while the value < 0 indicates plant N limitation. As a result, consistent with soil microbial P and N limitations, the strength of plant P limitation increases with value of plant N/P limitation and strength of plant N limitation increases as the value decreases.

Statistical analysis

We first tested the variable normality and homogeneity of variances for microbial metabolic limitations, and found that neither relative C limitation nor N/P limitation conformed to normal distributions or homogeneity of variances under different vegetation types and climate zones (Supplementary Tables 2 and 3). Thus, we adopted the linear mixed-effects models to identify the effects of soil layer, forest type and climate type, as well as their interaction effects, on microbial metabolic limitations. The model was constructed using the “nlme” package in R with “soil layer”, “forest type” and “climate type” as the fixed factors and “sampling site” as a random factor50. The standardized major axis (type II) regression via the “lmodel2” package was used to identify relationships between extracellular enzymatic activities51 (Supplementary Note 1).

For examining latitudinal patterns of microbial metabolic limitations, we first identified the latitudinal breakpoints from a piecewise linear regression analysis. The regression relationships of the relative C and N/P limitations with latitude were fitted and tested with linear models using the “segmented” package52. The confidence intervals of the breakpoints were calculated via 1000 bootstrap samples using the “SiZer” package53. We then used general linear models to fit the relationships of the relative C and N/P limitations with upper and lower latitudes (before and after breakpoints), respectively. Generalized linear models also examined relationships of microbial N/P limitation with plant N/P limitation. Pearson correlations examined relationships of microbial metabolic limitations with potential predictor variables.

Variation partitioning analysis was conducted using the “vegan” package54 to identify the relative importance of climate, vegetation and soil properties to variations in microbial N/P limitation in different soil layers. We further determined the relative influence of each variable using the random-forest algorithm using the “RandomForest” package55. In these random-forest models, fourteen variables served as predictors for the microbial N/P limitation (Supplementary Table 7). To estimate the importance of these variables, we used percentage increases in the MSE (mean squared error) of variables, and higher MSE values indicate more important variables. Significance of the models and cross-validated R2 values were assessed with 1500 permutations of the response variable, by using the “A3” package. Similarly, the significance of each predictor to the response variables was assessed with the “rfPermute” package (Supplementary Note 2).

Furthermore, to explore cascading relationships between microbial N/P limitation and these key drivers, all significant variables (drivers) were divided into seven categories (hydrothermal condition (aridity index and MAT), N deposition, climate seasonality (PSEA and TSEA), vegetation, soil pH, SOC and TP). Partial least squares path models were used to identify possible pathways whereby variables control microbial N/P limitation in different soil layers. The VIF of all variables in the model requires <10 for removing the internal factor collinearity of the block (module). Meanwhile, loading in Outer Models requires >0.7 for removing those independent variables with small contributions. The models were constructed using the “plspm” package56. All statistical analyses were performed using the R software (v.3.3.2)57.

Reporting summary

Further information on research design is available in the Nature Research Reporting Summary linked to this article.