Introduction

Globally, soils hold a large amount of carbon (C); the size of the SOC pool is twice that of the atmosphere and greater than the atmospheric and terrestrial vegetation C pools combined1,2. Small changes in the balance between inputs to and outputs from the SOC pool (especially the stable C pool) could have a significant impact on atmospheric CO2 concentration3,4,5, which may either reduce or exacerbate the consequences of burning of fossil fuels. For example, global climate change due to rising atmospheric CO2 concentration will likely change aboveground vegetation dynamics. Such changes can then impact the quantity of litter C inputs to soil by altering plant net primary productivity (NPP) and/or affect the quality of litter C inputs by altering plant community structure (e.g., population sizes and species composition) and/or chemical composition (e.g., nitrogen and phosphorus concentrations) of plant litter6,7,8,9.

Impacts on plant-C inputs to soil are fairly well represented in current climate and Earth system models (CESMs) that are driven largely by plant productivity responses10,11. However, microbial roles in soil C cycling are still poorly reflected in these models12, despite that dynamics of SOC is driven ultimately by microbial catabolic and anabolic activities2,13,14. In practice, microbial models can be both plausible and straightforward (but not easy) to parameterize; importantly they also show promise for improving overall knowledge and our ability to predict the effects of global changes12,15.

Calls for the explicit consideration of microorganisms and their activities in models are increasing15,16,17, the necessary theoretical supporting research are growing18,19,20,21 and results are beginning to emerge14,22,23,24,25,26. For example, a recent study showed that performance of the Community Land Model (a global land surface model) was substantially improved by addition of microbial processes24. While microbial models are not absent in current studies, unfortunately most of them only target microbial biomass. Inevitably, microbial controls over SOC formation, transformation and stabilization are engaged by numerous functional species that constitute the microbial community. Yet so far, virtually few published models have incorporated dynamics of microbial community and examine its relevance for SOC cycling27. Therefore, accounting for the responses of the microbial community and its physiology in CESMs may be necessary to reliably predict SOC dynamics28.

Emerging opinion suggests that soil microbes act as important agents of SOC formation13,29,30, in part, because of growing evidence that microbial-derived C forms are primary constituents of the stable SOC pool14,31,32,33,34,35. In this view, microbial activities simultaneously lead to (1) significant CO2 emissions via decomposition of plant residues and soil organic matter and (2) deposition of microbial-derived C into the SOC reservoir via turnover of microbial biomass and necromass. Hence, we believe that the balance of C response between microbial decomposition and deposition channels will be a key factor determining changes in the magnitude, composition (labile vs recalcitrant) and source configuration (microbial-derived vs plant-derived) of SOC.

Accordingly, we develop a mechanistic microbial model (Fig. 1) that represents and explores the coupled mechanisms behind SOC dynamics resulting from variations in plant litter inputs over decadal time scales by taking into account (1) the origins of SOC (i.e., microbial-derived and plant-derived), (2) the responses of various SOC and microbial biomass C pools to plant litter inputs and (3) the feedbacks between SOC pools and microbial functional groups (fungi and bacteria). A global optimization strategy36 and several strictly defined constraints (e.g., biomass, maintenance respiration) were used to ensure the validity and reliability of the calibrated model parameters and enable reliable model simulations. After model parameters and the corresponding equilibrium state were determined, we systematically explored the responses of various SOC pools to different litter C input scenarios by using our microbial model (Table 1). Our model is a biologically realistic model based on established theory and microbial parameters. We expect our work serve as the first step towards a new generation of models that include key physical and chemical mechanisms in the SOC cycling.

Table 1 The summary of simulations for the set of modeling exercises (EqualFungi) with the assumption that the ratio between fungal and bacterial biomass is 1:1. The other two sets of modeling exercises with the assumption that the ratio between fungal and bacterial biomass is 1:2 (LowFungi) and 2:1 (HighFungi) have the same simulation setups as 1:1 (EqualFungi)
Full size table
Figure 1
figure 1

Schematic of the model framework that was used to simultaneously simulate the interactions among dynamics of litter C input (i.e., quantity and quality), dynamics of biomass and necromass of microbial functional groups (i.e., fungi and bacteria) and dynamics of soil non-living C pools (i.e., labile and recalcitrant C).

Full size image

Results

Initially, we simulated the impacts of a change in the quantity of litter C inputs (doubled amount; i.e., LF) on the dynamics of SOC pools (Table 1). The simulation results showed that increasing the amount of litter C inputs alone eventually increased the size of living C pools (fungal and bacterial biomass C; Fig. 2, panels d and e), but had no impacts on the ratio of fungal to bacterial biomass C, the size of non-living C pools (labile and recalcitrant C; Fig. 2, panels a, b and q) and the source configuration of non-living C pools (i.e., the relative contributions of C derived from fungi, bacteria and plant sources in labile and recalcitrant C pools; Fig. 2, panels g – l). These results were consistent among three sets of modeling exercises in which we assumed that living fungal biomass C was the same as (EqualFungi), greater than (LowFungi), or less than (HighFungi) living bacterial biomass C (Fig. 2, Supplementary S1–S5). Our simulations with LF suggest that increases in litter C inputs due to increased NPP may change the size of soil non-living C pools in the short term (assuming plant community structure maintains stable), but that this response is transient. After ca. 3 years of increased litter C inputs alone, the non-living SOC pools returned to steady states equivalent to the equilibrium levels before NPP was increased. Basically, the extra C from doubled litter inputs increased microbial growth and eventually sustained a larger pool of living microbial biomass, without causing significant changes in non-living SOC pools.

Figure 2
figure 2

Modelled responses of SOC dynamics to different plant litter input scenarios.

The simulations were conducted with the assumption that 1) fungal biomass C equals to bacterial biomass C and 2) fungi and bacteria are equally sensitive to the litter input. B – Bacteria; F – Fungi; M – Microbial; Rec – Recalcitrant; P – Plant; LF – doubled quantity of litter input; LQ – lower quality of litter input; HQ – higher quality of litter input.

Full size image

Next we simulated how a change in the quality of litter C (i.e., HQ or LQ; hereafter referred to as ChangeQl) affected SOC dynamics (Table 1). We conducted these model runs under the assumption that fungal and bacterial communities have synchronous responses to a change in litter quality (i.e., fungi and bacteria had the same positive or negative percentage changes in biosynthesis rates), that is to say that fungal and bacterial groups were assumed to respond with the same MBA. Under this assumption, ChangeQl had large impacts on the size and composition of the living C pools and also altered the composition and source configuration of non-living SOC pools (Fig. 2, panels a, b, d, e, g – l; Supplementary S2 and S4). However, changes in the size of total non-living SOC were less than <5% (Fig. 2, panel c; Fig. 3) and this limited response was unaffected by the ratio of fungal-to-bacterial biomass (i.e., the three sets of modeling exercises: EqualFungi, LowFungi and HighFungi). Thus, when plant community structure is altered by changes in land-use or climate, these results suggest that a corresponding change in litter quality (whether increased or decreased) is likely to have limited impacts on the size of total non-living SOC, if microbial functional groups (fungi and bacteria) exhibit the same MBA.

Figure 3
figure 3

A stylized illustration integrating modelled results under various scenarios.

It depicts the model simulation results with the assumption that the ratio between fungal and bacterial biomass C is 1.0 at the initial equilibrium state. The filled square (biomass C pool) and circle (non-living C pool) shapes are drawn to scale with the length or diameter derived from the simulated pool sizes shown in Figure 2 and Appendix; however, different scale number was used for the biomass (square shape) and non-living (circle shape) C pools. The ratios (for example, 1.0:1.7) indicate the ratios between plant-derived and microbial-derived non-living C. S0 denotes the control simulation with no changes in quality and quantity of litter C input; St denotes the simulation with change in quantity of litter C input; Sl denotes the simulation with change in quality of litter C input. The red and blue line color denotes the simulation results with positive and negative microbial biosynthesis acceleration (MBA) to litter C input, respectively and fungi and bacteria have same magnitude (indicated as same filled color) or different magnitude (indicated as different color with darker color for greater magnitude) of such MBA.

Full size image

Finally, we conducted a series of simulations in which we assumed that fungal and bacterial communities react asynchronously to a change in litter quality (i.e., fungal and bacterial groups can respond with different percentage changes in their biosynthesis rates, that is to say, different MBA) (Table 1). Under this assumption, our simulations for all sets of modeling exercises (EqualFungi, LowFungi and HighFungi) suggested that ChangeQl significantly impacted the size, composition and source configuration of non-living C pools over the long term (Fig. 3, Supplementary S1 and S3).

Discussion

Our study demonstrates that the equilibrium response of total non-living SOC to plant litter inputs might not be controlled by litter quantity and microbial biomass and limitedly controlled by change in litter quality, but most significantly controlled by how MBA of microbial functional groups (fungi and bacteria) differs from each other. Considering the small fraction (<5%) of microbial living biomass in the total organic C in soils20, our results suggest that the total SOC at an equilibrium state may be an “intrinsic” property of a given ecosystem and such intrinsic property is ultimately controlled by asynchronous change in MBA of different microbial functional groups.

So far, microbial biomass is rarely considered as an additional C pool37,38,39 and even when such a microbial C pool exists it is simply represented as a C inventory rather than a driving factor of C decomposition (similar to soil temperature and moisture) in the existing CESMs13. Our model structure includes the under-explored microbial biomass as a fast turnover C pool and a major driver of SOC decomposition and more importantly our model also captures microbial functional groups. Furthermore, our model has capability of differentiating the microbial-derived C from the C of other origins, from which we could explore the extent of how turnover of microbial community biomass changes allocations of SOC between plant left-over and microbial recycled products. The clarification and quantification of the percentage of total SOC that is microbially derived have significant scientific implications as well as practical value, e.g., it is a prerequisite to understanding soil processes on the Earth and particularly useful for broad-scale calculations of global C model parameters29. This is fundamental yet under-explored knowledge, as few of ecologists and soil scientists have even attempted to measure it in a meaningful way. Our model simulations show that the contribution of microbial-derived C to total non-living C stock at equilibrium state ranges from 47% to 80%, which is consistent with previous studies indicating that microbial-derived C is important in SOC stabilization14,30,32. Although this range might be different in reality, these results do suggest that shifts in the source configuration (microbial-derived vs. plant-derived) of soil non-living C may be significantly changed due to environmental changes and such changes are controlled by plant-C quality and asynchronous MBA.

We argue that our model framework provides a baseline to investigate microbial-centered soil processes in SOC cycling and is ready to be coupled with models that simulate other abiotic factors such as oxidation and hydrolysis of SOC, soil habitat, physicochemical protection and spatial inaccessibility of SOC to microbes40. In particular, physical protections of SOC in micro-aggregates and chemically-bound to mineral particles have long been known to affect soil C storage41,42,43 and are drawing increasing attentions in recent years. New studies have shown that interactions between SOC and microbes and minerals play a more important role in SOC sequestration than chemical properties of SOC17,44. Critical evaluation of these mechanisms (i.e., physical and chemical protections) needs to be addressed in future model development. Also, it should be noted that the quality of fungal-, bacterial- and plant-derived labile and recalcitrant C was assumed to be the same in the model, while many studies have shown that microbial- and plant-derived C likely have different vulnerability to decomposition5,45. Nonetheless, our study highlights the important role of microbial community and that shifts in community biomass turnover change the size, composition and source configuration of non-living C owning to “acceleration effect” (here is asynchronous MBA to C input between microbial functional groups) and further may have great impacts on the global C cycling.

More broadly, the driving force (plant litter inputs) in our microbial model could be extended to investigate SOC dynamics due to other events such as root exudate inputs and C amendment (e.g., fertilization) or disturbances (e.g., land-use change, warming, fire and species invasion). It is likely that, qualitatively similarly to the changes of plant litter inputs, those events and disturbances will have great impacts on total SOC if MBA of fungal physiology differs from that of bacterial physiology. Overall, our model simulations suggest that future experimental studies may improve our understanding of plant-microbe-environment interactions by concentrating on the metabolic response of different microbial functional groups to global change and, most critically, the relative differences among those responses.

Methods

We developed a mechanistic model to simulate and investigate the microbial control over SOC cycling. In our modeling framework, the soil C pools are separated into two categories: living and non-living C pools. The living C pools (i.e., biomass) are further divided into living fungal and bacterial C pools; the non-living C pools are further divided into labile and recalcitrant C pools. The labile C pool depends on C loss due to microbial anabolism and catabolism (defined as the sum of biosynthesis and loss to respiration) and transformation of labile C to recalcitrant C and C gain due to both microbial necromass and external litter C inputs. The microbial biomass C pool depends on C loss (due to microbial death and maintenance respiration) and C gain (due to biosynthesis). The recalcitrant C pool depends on C loss due to the microbial anabolism and catabolism and C gain due to microbial necromass input and transformation of labile C to recalcitrant C. Below we describe the processes that affect the dynamics of each C pool and the corresponding mathematical equations. Table 2 provides a list of model parameters and their definitions and units.

Table 2 Descriptions of Model parameters
Full size table

Model development

Labile C pool

The labile C pool is quantitatively expressed as.

where CL is the labile C concentration (mg g−1); t is the time (month); Cinput is the external litter C input (mg g−1 month−1) that was set to 0.72 mg g−1 month−1; fd is the fraction of bacterial and fungal necromass that transfers to the labile C pools during turnover (unitless), respectively; kd,b and kd,f are the death rates of bacteria and fungi (month−1), respectively; B and F are the C concentration of bacterial and fungal biomass (mg g−1), respectively; δb and δf are the factors describing the bacterial and fungal activities, respectively and ranges from 0 to 1 (unitless); kb,l and kf,l are the labile C biosynthesis rates by bacteria and fungi (month−1), respectively. Equation (1) and the following equations indicate that bacterial and fungal biomass not only contribute to soil C pools but also are drivers of decomposition of labile and recalcitrant C pools46. The factors controlling the dynamics of bacterial and fungal activities (i.e., δb and δf) depend on the amount of substrate (e.g., CL) and the dynamical change of microbial activity follows a Michaelis-Menten type response function46,47:

where kI is the microbial inhibition constant (mg g−1).

Microbial biomass C pool

The microbial biomass C pool is quantitatively expressed as:

where CMBC is the microbial biomass C concentration (mg g−1); fu,l and fu,r are the fractions of assimilated labile and recalcitrant C that transfers to the microbial biomass C pool (unitless), respectively; kb,r and kf,r are the biosynthesis rates of the recalcitrant C by bacteria and fungi (month−1), respectively; CR is the recalcitrant C concentration (mg g−1); km is the microbial maintenance respiration rate (month−1).

Recalcitrant C pool

The recalcitrant C pool is quantitatively expressed as:

where ft is the fraction of labile C directly transformed to recalcitrant C. Other parameters in equation (7) have been defined earlier.

CO2 efflux

The total CO2 efflux (CCO2) from soil is composed of three components: CO2 efflux due to biosynthesis of both labile and recalcitrant C pools and CO2 efflux due to microbial maintenance respiration and is quantitatively expressed as:

where CCO2,l, CCO2,r and CCO2,m are the CO2 efflux derived from the biosynthesis of both labile and recalcitrant C pools and from the microbial maintenance respiration (mg g−1), respectively.

Model parameterization and modeling excercises

There are seven unknown parameters in the developed model: the death rates of fungi and bacteria (i.e., kd,f and kd,b), the biosynthesis rates of the labile C by fungi and bacteria (i.e., kf,l and kb,l), the biosynthesis rates of the recalcitrant C by fungi and bacteria (i.e., kf,r and kb,r) and microbial maintenance respiration rates (i.e., km). Five objective functions and three constraints were used to determine those unknown model parameters (see Supplementary for details). Briefly, these five objective functions were the ratio between labile and recalcitrant C pools, the ratio of microbial biomass C to total soil C, the ratio of microbial maintenance respiration to total soil heterotrophic respiration, the ratio between fungal and bacterial biomass C and the total soil non-living C content. The three constraints were that 1) turnover rate of bacterial biomass is constrained to be greater than that of fungal biomass, 2) biosynthesis rate of labile C by bacteria is constrained to be greater than by fungi and 3) biosynthesis rate of recalcitrant C by fungi is constrained to be greater than by bacteria. It is important to have these objective functions and constraints in order to determine a set of parameters that produce reasonable simulations on all of pools and processes included in our model.

After the objective functions and constraints were defined, a global optimization method (stochastic ranking evolutionary strategy36) was used to determine the model parameters. This was done by minimizing the difference between the objective functions and the simulations after model spinups to reach equilibrium state. After model parameters and the corresponding equilibrium state were determined, we systematically explored the temporal responses of various SOC pools to different litter C input scenarios by using our microbial model (Table 1). Alternative reasonable objective functions (e.g., different objective function on the ratio between fungal and bacterial biomass) might change the estimated parameters and quantitative results, but would not change the qualitative results. Please see the supplementary materials for the detailed descriptions of model parameterization and modeling exercises.