Spatial Control of Carbon Dynamics in Soil by Microbial Decomposer Communities

Pagel, Holger; Kriesche, Björn; Uksa, Marie; Poll, Christian; Kandeler, Ellen; Schmidt, Volker; Streck, Thilo

doi:10.3389/fenvs.2020.00002

ORIGINAL RESEARCH article

Front. Environ. Sci., 30 January 2020
Sec. Soil Processes
Volume 8 - 2020 | https://doi.org/10.3389/fenvs.2020.00002

Spatial Control of Carbon Dynamics in Soil by Microbial Decomposer Communities

Holger Pagel¹^*

Björn Kriesche²

Marie Uksa³

Christian Poll³

Ellen Kandeler³

Volker Schmidt²

Thilo Streck¹

¹Biogeophysics, Institute of Soil Science and Land Evaluation, University of Hohenheim, Stuttgart, Germany
²Institute of Stochastics, Ulm University, Ulm, Germany
³Soil Biology, Institute of Soil Science and Land Evaluation, University of Hohenheim, Stuttgart, Germany

Trait-based models have improved the understanding and prediction of soil organic matter dynamics in terrestrial ecosystems. Microscopic observations and pore scale models are now increasingly used to quantify and elucidate the effects of soil heterogeneity on microbial processes. Combining both approaches provides a promising way to accurately capture spatial microbial-physicochemical interactions and to predict overall system behavior. The present study aims to quantify controls on carbon (C) turnover in soil due to the mm-scale spatial distribution of microbial decomposer communities in soil. A new spatially explicit trait-based model (SpatC) has been developed that captures the combined dynamics of microbes and soil organic matter (SOM) by taking into account microbial life-history traits and SOM accessibility. Samples of spatial distributions of microbes at μm-scale resolution were generated using a spatial statistical model based on Log Gaussian Cox Processes which was originally used to analyze distributions of bacterial cells in soil thin sections. These μm-scale distribution patterns were then aggregated to derive distributions of microorganisms at mm-scale. We performed Monte-Carlo simulations with microbial distributions that differ in mm-scale spatial heterogeneity and functional community composition (oligotrophs, copiotrophs, and copiotrophic cheaters). Our modeling approach revealed that the spatial distribution of soil microorganisms triggers spatiotemporal patterns of C utilization and microbial succession. Only strong spatial clustering of decomposer communities induces a diffusion limitation of the substrate supply on the microhabitat scale, which significantly reduces the total decomposition of C compounds and the overall microbial growth. However, decomposer communities act as functionally redundant microbial guilds with only slight changes in C utilization. The combined statistical and process-based modeling approach derives distribution patterns of microorganisms at the mm-scale from microbial biogeography at microhabitat scale (μm) and quantifies the emergent macroscopic (cm) microbial and C dynamics. Thus, it effectively links observable process dynamics to the spatial control by microbial communities. Our study highlights a powerful approach that can provide further insights into the biological control of soil organic matter turnover.

Introduction

Microorganisms drive biochemical processes such as C cycling in soil (Falkowski et al., 2008). There is growing consensus that soil organic matter dynamics and stability are strongly controlled by microbial processing and associated bioenergetics constraints (Schmidt et al., 2011; Lehmann and Kleber, 2015; Williams and Plante, 2018). Yet, understanding how microbial community characteristics affect rates of biogeochemical processes remains a major research challenge. Further progress to quantitatively describe spatial arrangements between microorganisms in their micro-environment and their corresponding substrate is needed (Graham et al., 2016; Baveye et al., 2018).

Categorizing microbial communities based on life-history strategies (e.g., copiotrophs/oligotrophs, r-/K-stategists, autochtonous/zymogenous microorganisms, or competitors/stress tolerators/ruderals) is useful to link microbial community characteristics to biogeochemical processes (Fierer et al., 2007; Kuzyakov et al., 2009; Martiny et al., 2015; Fierer, 2017; Blankinship et al., 2018; Hall et al., 2018). These frameworks are based on the transfer of macroscale ecology concepts to microbial ecology. A recent study refined the competitor-stress tolerator-ruderal concept from plant ecology and suggested to define three microbial life history strategies: resource acquisition, stress tolerance, and high yield (Malik et al., 2020). Life-history strategies embrace combinations and trade-offs of microbial community traits related to maximum growth rate, dormancy, substrate affinity, production of specific enzymes, or stress tolerance mechanisms (Webb et al., 2010; Fierer et al., 2014; Trivedi et al., 2016; Alster et al., 2018; Rath et al., 2019; Malik et al., 2020). Mineralization of soil C could be seen as an emergent process that is regulated by functional traits of soil microorganisms and microbiological interactions (Addiscott, 2010). Therefore, decomposition of C compounds is controlled by dynamics of assemblages of somewhat functionally redundant organisms organized in microbial guilds with characteristic life-history strategies (Schimel and Schaeffer, 2012).

Including measured functional traits of plants as well as soil microorganisms and fauna in biogeochemical modeling is a promising approach to improve predictions of biogeochemical cycling in soil (Fry et al., 2019). Biogeochemical C models increasingly include metabolic and physiological traits as well as life-history strategies to account for microbial regulation of decomposition processes (Garnier et al., 2001; Ingwersen et al., 2008; Neill and Guenet, 2010; Allison, 2012; Bouskill et al., 2012; Pagel et al., 2014, 2016; Perveen et al., 2014; Wang et al., 2014; Le Roux et al., 2016). Including microbial dormancy of microbes in models has been shown to improve the prediction of soil organic C dynamics (He et al., 2015) as did accounting for copiotrophic and oligotrophic microorganisms as physiologically distinct functional groups (Wieder et al., 2015). A model-based analysis demonstrated that adaptive microbial responses to C limitation and water stress might emerge from microbial traits related to dormancy and production of extracellular polymeric substances (Brangarí et al., 2018). The importance of community-level regulation and microbial trait trade-offs was highlighted by trait-based modeling of litter decomposition (Kaiser et al., 2015; Allison and Goulden, 2017).

Further integration of trait-based and spatial explicit approaches is, however, essential to advance the quantitative description of microbial C utilization, because microbial activity is controlled by spatial characteristics. Physical accessibility of organic compounds to microorganisms strongly affects substrate supply and microbial community functioning (Brookes et al., 2017; Nunan et al., 2017; Schimel, 2018). It has been conjectured that at the pore-scale, which is relevant for microbial processes, the supply of assimilable C (low molecular weight compounds <600 Da) to microorganisms is mainly regulated by (i) physical accessibility of soil organic matter, (ii) exoenzymatic decomposition of C compounds that are not directly assimilable (high molecular weight compounds ≥ 600 Da), and (iii) diffusive transport of assimilable C in the soil solution from locations of exoenzymatic action to microbial cells (Lehmann and Kleber, 2015; Schimel et al., 2017; Blankinship et al., 2018; Sokol et al., 2019).

Quantitative measurements of microbial distribution and processes at the pore-scale are extremely challenging. Though there is limited, albeit growing, experimental data on the spatial organization and activity of microorganisms in soils, a number of mechanistic models have been applied to understand and predict the impact of spatial heterogeneity in soil on microbial and physico-chemical processes (Baveye et al., 2018). Raynaud and Nunan (2014) analyzed the spatial distribution of bacterial cells in soil thin sections and described the spatial structure of observed bacterial distributions as aggregated point patterns using a Log Gaussian Cox process as spatial statistical model. Their analysis indicated that distributions of bacterial cells in soils are clustered and non-random at the μm-scale, most probably as a result of heterogeneity in soil structure and pore network architecture. Recent experimental evidence from combined X-ray microtomography and fluorescence microscopy at different sampled spatial scales (0.2, 1, 5 mm) indicates that pore characteristics might effectively influence the distribution of bacteria in soil mainly at a spatial scale of ~5 mm (Juyal et al., 2019). Most rapid decomposition rates were associated with pores of neck diameters of 15–90 μm. This was attributed to optimal microbial habitat conditions with respect to nutrient and oxygen supply and organism motility (Strong et al., 2004; Kravchenko and Guber, 2017). There is some experimental evidence that pore characteristics and microenvironmental conditions control the relative contributions of specific functional microbial groups to decomposition of C compounds and the extent of their functional redundancy (Ruamps et al., 2013; Negassa et al., 2015; Kravchenko and Guber, 2017; Nunan et al., 2017).

A few recent models linked mechanistic descriptions of a soil's pore structure with trait-based microbial dynamics. Experimental work using artificial micrometric pore networks etched in glass combined with modeling has demonstrated that oxygen-carbon counter-gradients (as commonly found in microbial hotspots like the rhizosphere or detritusphere) induce the spatial organization of aerobic and anaerobic bacteria and promote their stable coexistence (Borer et al., 2018). Scenario simulations using a multi-species 3D pore-scale soil C model have indicated microscale (μm) control of bacterial diversity driven by the degree of heterogeneity in the spatial distribution of organic matter (Portell et al., 2018). In these simulations, the spatial heterogeneity of organic matter affected the succession of functional bacterial types differing in growth rates and substrate affinities. Irrespective of the spatial SOM distribution, however, the small-scale (mm) C turnover was similar. This indicates functional redundancy with respect to C cycling. While there are some first successful attempts to derive mechanistic effective rate laws for specific biogeochemical processes at pedon to landscape scale from pore-scale modeling (e.g., Ebrahimi and Or, 2018; Schmidt et al., 2018), the upscaling of microbial processes and their control from pore scale to macroscopic scales (pedon to landscape), which are practically relevant and accessible to direct observation, remains a largely unresolved research challenge (Baveye et al., 2018).

This theoretical study aims to elucidate the control of emerging C dynamics in soil at the macroscale (cm) by the pore-scale (μm) distribution of decomposer communities consisting of microorganisms with differing life-history traits. A new trait-based soil C model was utilized in combination with a spatial statistical model of microbial biogeography (Raynaud and Nunan, 2014) to test two hypotheses: (i) increasing spatial heterogeneity in the distribution of microbial decomposers results in an increase in diffusion-limited C availability and lower C turnover and (ii) with increasing spatial heterogeneity, the composition of decomposer communities shifts to a higher proportion of oligotrophic organisms that can outcompete copiotrophs at low C availability.

Materials and Methods

Model Rationale and Main Assumptions

The 2D spatially explicit trait-based soil C model (SpatC) has been developed to study the effects of mm-scale heterogeneous distribution of functionally diverse microbial communities on C cycling in soil. Following the conceptual soil continuum model of soil organic matter cycling (Lehmann and Kleber, 2015), SpatC distinguishes three conceptual carbon pools with respect to their assimilability by microorganisms (Figure 1). Microbial communities are grouped into three functional types that distinguish different life-history strategies according to ecological categorizations, a technique used similarly in other models (e.g., Allison, 2012; Kaiser et al., 2015). This structure reflects fundamentally different life-history strategies according to functional-ecological frameworks such as the copiotrophy–oligotrophy continuum or Grime's competitor–stress tolerator–ruderal concept (Fierer et al., 2007; Krause et al., 2014; Fierer, 2017; Ho et al., 2017; Huang et al., 2018; Fry et al., 2019; Maynard et al., 2019). The biomass of all microbial groups is regulated by growth of predators that utilize microbial pools as C and energy sources. SpatC thereby explicitly considers exploitative competition (interception of a common resource), interference competition (direct interactions between microorganisms), and predator-mediated competition (top-down control of microorganisms by selective predation) between the three functional microbial groups (see Buchkowski et al., 2017).

FIGURE 1

Figure 1. Conceptual scheme of coupled carbon turnover and biochemical interactions implemented in the 2D spatially explicit trait-based soil C model (SpatC). Solid arrows indicate carbon fluxes. Dashed green arrows depict the controls on extracellular depolymerisation reactions. C_M, C_S, and C_L stands for monomers, small polymers, and large polymers, respectively. Superscript “S” indicates sorbed phase concentration of C_M and C_S. Monomers and small polymers may be transported by 2D diffusion (not shown). Microbial communities consist of active (superscript “a”) and dormant (superscript “d”) oligotrophs (B_O), copiotrophs (B_C), and copiotophic cheaters (B_CC). P stands for predators.

Governing Equations and Fluxes

SpatC is formulated as a set of coupled partial and ordinary differential equations. All C pools are based on the C mass balance in soil and expressed in mg g⁻¹. We assumed $\frac{\partial C_{S}}{\partial n} = 0$ and $\frac{\partial C_{M}}{\partial n} = 0$ at all boundaries (with n denoting the outward facing normal vector), i.e., there was no flux of C_S and C_M out of the considered domain. Model parameters are given in Tables 1, 2. Fluxes and functions are specified in Section Fluxes and Functions. A concise description of all model equations is given in the Supplementary Material.

TABLE 1

Table 1. Parameterization of functional microbial traits.

TABLE 2

Table 2. Inhibition, maintenance, and non-microbial parameters and initial values of SpatC model simulations.

Non-microbial Carbon

Large biopolymers (C_L, Equation 1) are not directly assimilable by microorganisms, but need to be first depolymerized by extracellular enzymes to dissolved small biopolymers (C_S, Equation 2) and monomers (C_M, Equation 3). Small biopolymers are similarly prone to extracellular depolymerisation. This enzymatic process is simulated using Michaelis-Menten kinetics without explicitly considering enzyme dynamics. The depolymerization rate of large and small polymers is instead directly controlled by microbial biomass (Equation 21). Small polymers and monomers are directly consumed by microorganisms. The decay of microorganisms and predators leads to C input of non-microbial C to C_L, C_S, and C_M. While large biopolymers are not transported, SpatC accounts for transport of small polymers and monomers by diffusion. Using the approach of Streck et al. (1995), the bioavailability of C_S and C_M is further constrained by rate-limited, two-stage, non-linear sorption (Equations 4–6).

\begin{array}{l} \frac{\partial C_{L}}{\partial t} & = & - \underset{depolymerization}{\underset{︸}{r_{L}}} \\ + & \underset{decay of microorganisms and predators}{\underset{︸}{f_{P, L} \cdot \vec{1} \cdot \vec{r_{P}} + f_{m, L} \cdot \vec{1} \cdot (\vec{r_{m}^{B}} - \vec{r_{m}^{M}}) + r_{d, P}}} & (1) \end{array}

\begin{array}{l} \frac{\partial C_{S}}{\partial t} = \underset{\begin{matrix} retardation \\ factor due to \\ equillibrium \\ sorption \end{matrix}}{\underset{︸}{\frac{1}{R_{S}}}} \cdot (\begin{matrix} \overset{\begin{matrix} input from C_{L} \\ depolymerization \end{matrix}}{\overset{︷}{f_{S} \cdot r_{L}}} - \overset{depolymerization}{\overset{︷}{r_{S}}} - \overset{microbial consumption for growth}{\overset{︷}{\frac{1}{Y_{S, O}} \cdot r_{μ, O}^{S} - \frac{1}{Y_{S, C}} \cdot (r_{μ, C}^{S} + r_{μ, C C}^{S})}} + \\ \underset{microbial decay}{\underset{︸}{f_{P, S} \cdot \vec{1} \cdot \vec{r_{P}} + f_{m, S} \cdot \vec{1} \cdot (\vec{r_{m}^{B}} - \vec{r_{m}^{M}})}} - \underset{kinetic sorption}{\underset{︸}{α_{S} \cdot (C_{S, S 1} - C_{S, S 2})}} + \underset{diffusion}{\underset{︸}{D_{e, S} \cdot \nabla^{2} C_{S}}} \end{matrix}) & (2) \end{array}

\begin{array}{l} \frac{\partial C_{M}}{\partial t} = \underset{\begin{matrix} retardation \\ factor due to \\ equillibrium \\ sorption \end{matrix}}{\underset{︸}{\frac{1}{R_{M}}}} \cdot (\begin{matrix} \overset{\begin{matrix} input from C_{L} and C_{S} \\ depolymerization \end{matrix}}{\overset{︷}{(1 - f_{S}) \cdot r_{L} + r_{S}}} - \overset{microbial consumption for growth and maintenance}{\overset{︷}{\frac{1}{Y_{M, O}} \cdot r_{μ, O}^{M} - \frac{1}{Y_{M, C}} \cdot (r_{μ, C}^{M} + r_{μ, C C}^{M}) - \vec{1} \cdot \vec{r_{m}^{M}}}} + \\ \underset{microbial decay}{\underset{︸}{f_{P, M} \cdot \vec{1} \cdot \vec{r_{P}} + f_{m, M} \cdot \vec{1} \cdot (\vec{r_{m}^{B}} - \vec{r_{m}^{M}})}} - \underset{kinetic sorption}{\underset{︸}{α_{M} \cdot (C_{M, S 1} - C_{M, S 2})}} + \underset{diffusion}{\underset{︸}{D_{e, M} \cdot \nabla^{2} C_{M}}} \end{matrix}) & (3) \end{array}

Equilibrium sorption is considered using a Freundlich isotherm. Sorbed phase concentrations of small biopolymers and monomers at sorption sites in region 1 are accordingly expressed as:

\begin{array}{l} \begin{matrix} C_{S, S 1} & = & K_{F, S} \cdot {(\frac{ρ_{B}}{θ} \cdot C_{S})}^{m_{S}} \\ C_{M, S 1} & = & K_{F, M} \cdot {(\frac{ρ_{B}}{θ} \cdot C_{M})}^{m_{M}} \end{matrix} & (4) \end{array}

Kinetic sorption is expressed as mass transfer between sorption sites in region 1 and region 2:

\begin{array}{l} \begin{matrix} \frac{\partial C_{S, S 2}}{\partial t} & = & \frac{1}{1 - f_{S, S}} \cdot α_{S} \cdot (C_{S, S 1} - C_{S, S 2}) \\ \frac{\partial C_{M, S 2}}{\partial t} & = & \frac{α_{M}}{1 - f_{M, S}} \cdot (C_{M, S 1} - C_{M, S 2}) \end{matrix} & (5) \end{array}

Total sorbed phase concentrations are given by the sum of sorbed phase concentrations in region 1 and 2, each weighted by the fraction of sorption sites in both regions:.

\begin{array}{l} \begin{matrix} C_{S}^{S} & = & \overset{\begin{matrix} fraction of \\ region 1 sites \\ for small \\ biopolymers \end{matrix}}{\overset{︷}{f_{S, S}}} \cdot C_{S, S 1} + \overset{\begin{matrix} fraction of \\ region 2 sites \\ for small \\ biopolymers \end{matrix}}{\overset{︷}{(1 - f_{S, S})}} \cdot C_{S, S 2} \\ C_{M}^{S} & = & \underset{\begin{matrix} fraction of \\ region 1 sites \\ for monomers \end{matrix}}{\underset{︸}{f_{M, S}}} \cdot C_{M, S 1} + \underset{\begin{matrix} fraction of \\ region 2 sites \\ for monomers \end{matrix}}{\underset{︸}{(1 - f_{M, S})}} \cdot C_{M, S 2} \end{matrix} & (6) \end{array}

Functional Microbial Groups

SpatC accounts for three functional microbial types: oligotrophs (B_O), copiotrophs (B_C) and copiotrophic cheaters (B_CC) (Equations 7–11). All microbial groups are considered to be able to switch from an active to a dormant physiological state (Lennon and Jones, 2011; Blagodatskaya and Kuzyakov, 2013; Joergensen and Wichern, 2018) with different parameterizations for different functional types (Table 1, Figure 2). Active microorganisms use dissolved small biopolymers and monomers for growth, while dormant microorganisms do not grow. Maintenance energy requirements of microorganisms are assumed to be fulfilled through the uptake of monomers at sufficient substrate supply and are met from biomass when monomers become limiting (Wang and Post, 2012). That is, microorganisms switch from exogenous to endogenous maintenance (see Equation 18) leading to microbial decay at low substrate availability. We consider that endogenous maintenance proportionally results in the formation of dead microbial biomass and CO₂ (Equations 1–3 and 14). Additionally, microbial biomass decays due to predation. Thereby, microbial C is used for growth of predators (Equation 13), reallocated to non-microbial C pools in soil (Equations 1–3) and lost to CO₂ (Equation 14).

FIGURE 2

Figure 2. Schematic illustration of trade-off in functional microbial traits as implemented in SpatC (see also Table 1).

Dynamics of active microorganisms are expressed as follows:

\begin{array}{l} \frac{\partial B_{O}^{a}}{\partial t} = \underset{growth}{\underset{︸}{r_{μ, O}^{S} + r_{μ, O}^{M}}} - \underset{deactivation + reactivation}{\underset{︸}{r_{d, O} + r_{r, O}}} \\ - \underset{\begin{matrix} microbial decay \\ due to maintenance \end{matrix}}{\underset{︸}{\frac{1}{Y_{m}} \cdot (r_{m, O}^{a, B} - r_{m, O}^{a, M})}} - \underset{microbial decay by predation}{\underset{︸}{\frac{1}{Y_{P}} \cdot (1 - f_{P}) \cdot r_{P, O}^{a} - f_{P} \cdot r_{P, O}^{a}}} & (7) \end{array}

\begin{array}{l} \frac{\partial B_{C}^{a}}{\partial t} = \underset{growth}{\underset{︸}{r_{μ, C}^{S} + r_{μ, C}^{M}}} - \underset{deactivation + reactivation}{\underset{︸}{r_{d, C} + r_{r, C}}} \\ - \underset{\begin{matrix} microbial decay \\ due to maintenance \end{matrix}}{\underset{︸}{\frac{1}{Y_{m}} \cdot (r_{m, C}^{a, B} - r_{m, C}^{a, M})}} - \underset{microbial decay by predation}{\underset{︸}{\frac{1}{Y_{P}} \cdot (1 - f_{P}) \cdot r_{P, C}^{a} - f_{P} \cdot r_{P, C}^{a}}} & (8) \end{array}

\begin{array}{l} \frac{\partial B_{C C}^{a}}{\partial t} = \underset{growth}{\underset{︸}{r_{μ, C C}^{S} + r_{μ, C C}^{M}}} - \underset{deactivation + reactivation}{\underset{︸}{r_{d, C C} + r_{r, C C}}} \\ - \underset{\begin{matrix} microbial decay \\ due to maintenance \end{matrix}}{\underset{︸}{\frac{1}{Y_{m}} \cdot (r_{m, C C}^{a, B} - r_{m, C C}^{a, M})}} - \underset{microbial decay by predation}{\underset{︸}{\frac{1}{Y_{P}} \cdot (1 - f_{P}) \cdot r_{P, C C}^{a} - f_{P} \cdot r_{P, C C}^{a}}} & (9) \end{array}

Dynamics of dormant microorganisms are given by:

\begin{array}{l} \frac{\partial B_{O}^{d}}{\partial t} = \underset{deactivation - reactivation}{\underset{︸}{r_{d, O} - r_{r, O}}} - \underset{maintenance}{\underset{︸}{\frac{1}{Y_{m}} \cdot (r_{m, O}^{d, B} - r_{m, O}^{d, M})}} \\ - \underset{microbial decay by predation}{\underset{︸}{\frac{1}{Y_{P}} \cdot (1 - f_{P}) \cdot r_{P, O}^{d} - f_{P} \cdot r_{P, O}^{d}}} & (10) \end{array}

\begin{array}{l} \frac{\partial B_{C}^{d}}{\partial t} = \underset{deactivation - reactivation}{\underset{︸}{r_{d, C} - r_{r, C}}} - \underset{maintenance}{\underset{︸}{\frac{1}{Y_{m}} \cdot (r_{m, C}^{d, B} - r_{m, C}^{d, M})}} \\ - \underset{microbial decay by predation}{\underset{︸}{\frac{1}{Y_{P}} \cdot (1 - f_{P}) \cdot r_{P, C}^{d} - f_{P} \cdot r_{P, C}^{d}}} & (11) \end{array}

\begin{array}{l} \frac{\partial B_{C C}^{d}}{\partial t} = \underset{deactivation - reactivation}{\underset{︸}{r_{d, C C} - r_{r, C C}}} - \underset{maintenance}{\underset{︸}{\frac{1}{Y_{m}} \cdot (r_{m, C C}^{d, B} - r_{m, C C}^{d, M})}} \\ - \underset{microbial decay by predation}{\underset{︸}{\frac{1}{Y_{P}} \cdot (1 - f_{P}) \cdot r_{P, C C}^{d} - f_{P} \cdot r_{P, C C}^{d}}} & (12) \end{array}

Dynamics of predators are modeled using first-order growth and decay. It is considered that only part of the killed microbial biomass is actually taken up by predators. A fraction of C from killed microorganisms (f_P, Equation 20) is directly released to the soil solution and reallocated to non-microbial soil pools:

\begin{array}{l} \frac{\partial P}{\partial t} = \underset{growth}{\underset{︸}{(1 - f_{P}) \cdot \vec{1} \cdot \vec{r_{P}}}} - \underset{decay}{\underset{︸}{r_{d, P}}} & (13) \end{array}

Formation of carbon dioxide (CO₂) results from energy metabolism by aerobic respiration during microbial growth and maintenance as well as growth of predators:

\begin{array}{l} \begin{matrix} \frac{\partial C O_{2}}{\partial t} = \overset{microbial growth respiration}{\overset{︷}{\frac{1 - Y_{S, O}}{Y_{S, O}} \cdot r_{μ, O}^{S} + \frac{1 - Y_{M, O}}{Y_{M, O}} \cdot r_{μ, O}^{M} + \frac{1 - Y_{S, C}}{Y_{S, C}} \cdot (r_{μ, C}^{S} + r_{μ, C C}^{S}) + \frac{1 - Y_{M, C}}{Y_{M, C}} \cdot (r_{μ, C}^{M} + r_{μ, C C}^{M})}} + \\ \underset{microbial maintenance respiration}{\underset{︸}{\frac{1 - Y_{m}}{Y_{m}} \cdot \vec{1} \cdot (\vec{r_{m}^{B}} - \vec{r_{m}^{M}}) + \vec{1} \cdot \vec{r_{m}^{M}}}} + \underset{growth respiration of predators}{\underset{︸}{\frac{1 - Y_{P}}{Y_{P}} \cdot (1 - f_{P}) \cdot \vec{1} \cdot \vec{r_{P}}}} \end{matrix} & (14) \end{array}

Fluxes and Functions

The following flux equations define the C flows between soil organic matter pools and soil biota. All fluxes are expressed in mg g⁻¹d⁻¹.

Predation and maintenance fluxes were combined into column vectors. These were then used in the governing equations (Equations 1–14) as a scalar product with a row vector of ones for an effective description of the model:

\begin{array}{l} \vec{r_{P}} = (\begin{matrix} r_{P, O}^{a} \\ r_{P, O}^{d} \\ r_{P, C}^{a} \\ r_{P, C}^{d} \\ r_{P, C C}^{a} \\ r_{P, C C}^{d} \end{matrix}) \vec{r_{m}^{M}} = (\begin{matrix} r_{m, O}^{a, M} \\ r_{m, O}^{d, M} \\ r_{m, C}^{a, M} \\ r_{m, C}^{d, M} \\ r_{m, C C}^{a, M} \\ r_{m, C C}^{d, M} \end{matrix}) \vec{r_{m}^{B}} = (\begin{matrix} r_{m, O}^{a, B} \\ r_{m, O}^{d, B} \\ r_{m, C}^{a, B} \\ r_{m, C}^{d, B} \\ r_{m, C C}^{a, B} \\ r_{m, C C}^{d, B} \end{matrix}) \\ \vec{1} = (1, 1, 1, 1, 1, 1) \end{array}

A multi-substrate Monod kinetic (Lendenmann and Egli, 1998) is used to simulate grow of functional microbial types on small polymers and monomers (Equation 15). Following the proposed application of Grime's competitor–stress tolerator–ruderal concept to soil bacterial heterotrophs (Fierer, 2017), copiotrophs are parameterized as competitors. They are assumed to be most competitive by inhibiting the growth of oligotrophs and copitrophic cheaters. This is implemented using a first-order inhibition term (Buchkowski et al., 2017):

\begin{array}{l} \begin{matrix} r_{μ, O}^{S} & = & \underset{growth on small polymers}{\underset{︸}{\frac{μ_{max, O} \cdot C_{S} \cdot k_{O, S}}{μ_{max, O} + C_{S} \cdot k_{O, S} + C_{M} \cdot k_{O, M}} \cdot B_{O}^{a}}} - \underset{\begin{matrix} inhibition by \\ copiotrophs \end{matrix}}{\underset{︸}{k_{I} \cdot B_{C}^{a} \cdot B_{O}^{a}}} \\ r_{μ, O}^{M} & = & \underset{growth on monomers}{\underset{︸}{\frac{μ_{max, O} \cdot C_{M} \cdot k_{O, M}}{μ_{max, O} + C_{S} \cdot k_{O, S} + C_{M} \cdot k_{O, M}} \cdot B_{O}^{a}}} - \underset{\begin{matrix} inhibition by \\ copiotrophs \end{matrix}}{\underset{︸}{k_{I} \cdot B_{C}^{a} \cdot B_{O}^{a}}} \\ r_{μ, C}^{S} & = & \frac{μ_{max, C} \cdot C_{S} \cdot k_{C, S}}{μ_{max, C} + C_{S} \cdot k_{C, S} + C_{M} \cdot k_{C, M}} \cdot B_{C}^{a} \\ r_{μ, C}^{M} & = & \frac{μ_{max, C} \cdot C_{M} \cdot k_{C, M}}{μ_{max, C} + C_{S} \cdot k_{C, S} + C_{M} \cdot k_{C, M}} \cdot B_{C}^{a} \\ r_{μ, C C}^{S} & = & \frac{μ_{max, C C} \cdot C_{S} \cdot k_{C C, S}}{μ_{max, C C} + C_{S} \cdot k_{C C, S} + C_{M} \cdot k_{C C, M}} \cdot B_{C C}^{a} \\ - k_{I} \cdot B_{C}^{a} \cdot B_{C C}^{a} \\ r_{μ, C C}^{M} & = & \frac{μ_{max, C C} \cdot C_{M} \cdot k_{C C, M}}{μ_{max, C C} + C_{S} \cdot k_{C C, S} + C_{M} \cdot k_{C C, M}} \cdot B_{C C}^{a} \\ - k_{I} \cdot B_{C}^{a} \cdot B_{C C}^{a} \end{matrix} & (15) \end{array}

Switching between dormant and active state was modeled as first-order process (Equation 16) based on the approach of Mellage et al. (2015). Deactivation and reactivation rates are triggered by the concentration of dissolved monomers using a switching function (Equation 17). This function approaches zero if the monomer concentration is below a trait-specific threshold value and takes a value of one above the threshold. The shape parameter α controls the sharpness of the transition. It was fixed to a value of 0.1 to reflect a relatively sharp switching from and to dormancy.

\begin{array}{l} \begin{matrix} r_{d, O} = \underset{Deactivation}{\underset{︸}{(1 - ϕ_{O}) \cdot k_{d, O} \cdot B_{O}^{a}}} \\ r_{r, O} = \underset{Re a c t i v a t i o n}{\underset{︸}{ϕ_{O} \cdot k_{r, O} \cdot B_{O}^{d}}} \\ r_{d, C} = (1 - ϕ_{C}) \cdot k_{d, C} \cdot B_{C}^{a} \\ r_{r, C} = ϕ_{C} \cdot k_{r, C} \cdot B_{C}^{d} \\ r_{d, C C} = (1 - ϕ_{C C}) \cdot k_{d, C C} \cdot B_{C C}^{a} \\ r_{r, C C} = ϕ_{C C} \cdot k_{r, C C} \cdot B_{C C}^{d} \end{matrix} & (16) \end{array}

\begin{array}{l} \begin{matrix} ϕ_{O} = \frac{1}{e^{\frac{C_{t h r e s, O} - C_{M}}{α \cdot C_{t h r e s, O}}} + 1} \\ ϕ_{C} = \frac{1}{e^{\frac{C_{t h r e s, C} - C_{M}}{α \cdot C_{t h r e s, C}}} + 1} \\ ϕ_{C C} = \frac{1}{e^{\frac{C_{t h r e s, C C} - C_{M}}{α \cdot C_{t h r e s, C C}}} + 1} \\ with α = 0.1 \end{matrix} & (17) \end{array}

Total required maintenance uptake is given by the product of the trait-specific maximum maintenance rate coefficient and microbial biomass. Reduced maintenance needs of dormant microorganisms are considered using a reduction factor (β) of maximum maintenance rate coefficients. The relative C flux needed for maintenance that can be fulfilled from dissolved monomers (exogenous maintenance) is calculated using a Michealis-Menten type rate law (Lendenmann and Egli, 1998).

\begin{array}{l} \begin{matrix} r_{m, O}^{a, B} = \overset{\begin{matrix} Total C-flux required \\ to fulfill maintenance needs \end{matrix}}{\overset{︷}{m_{max, O} \cdot B_{O}^{a}}} \\ r_{m, O}^{a, M} = \overset{\begin{matrix} Proportional C-flux \\ of exogenous maintenance \end{matrix}}{\overset{︷}{(\frac{m_{max, O} \cdot C_{M} \cdot k_{O, M}}{m_{max, O} + C_{M} \cdot k_{O, M}}) \cdot B_{O}^{a}}} \end{matrix} \\ \begin{matrix} r_{m, O}^{d, B} = m_{max, O} \cdot β_{O} \cdot B_{O}^{d} r_{m, O}^{d, M} \\ = (\frac{m_{max, O} \cdot C_{M} \cdot k_{O, M}}{m_{max, O} + C_{M} \cdot k_{O, M}}) \cdot β_{O} \cdot B_{O}^{d} \\ r_{m, C}^{a, B} = m_{max, C} \cdot B_{C}^{a} r_{m, C}^{a, M} = (\frac{m_{max, C} \cdot C_{M} \cdot k_{C}}{m_{max, C} + C_{M} \cdot k_{C}}) \cdot B_{C}^{a} \\ r_{m, C}^{d, B} = m_{max, C} \cdot β_{C} \cdot B_{C}^{d} \end{matrix} \\ \begin{matrix} | r_{m, C}^{d, M} = (\frac{m_{max, C} \cdot C_{M} \cdot k_{C}}{m_{max, C} + C_{M} \cdot k_{C}}) \cdot β_{C} \cdot B_{C}^{d} \\ r_{m, C C}^{a, B} = m_{max, C C} \cdot B_{C C}^{a} r_{m, C C}^{a, M} \\ = (\frac{m_{max, C C} \cdot C_{M} \cdot k_{C C}}{m_{max, C C} + C_{M} \cdot k_{C C}}) \cdot B_{C C}^{a} \\ r_{m, C C}^{d, B} = m_{max, C C} \cdot β_{C C} \cdot B_{C C}^{d} r_{m, C C}^{d, M} \\ = (\frac{m_{max, C C} \cdot C_{M} \cdot k_{C C}}{m_{max, C C} + C_{M} \cdot k_{C C}}) \cdot β_{C C} \cdot B_{C C}^{d} \end{matrix} & (18) \end{array}

Predation of microorganisms and the associated growth of predators as well as the decay of predators is reflected by first-order expressions. Decreased predation of dormant microorganisms is considered by reduction factors (γ) of predation rate coefficients:

\begin{array}{l} \begin{matrix} r_{P, O}^{a} & = & k_{P, O} \cdot P \cdot B_{O}^{a} \\ r_{P, O}^{d} & = & k_{P, O} \cdot γ_{O} \cdot P \cdot B_{O}^{d} \\ r_{P, C}^{a} & = & k_{P, C} \cdot P \cdot B_{C}^{a} \\ r_{P, C}^{d} & = & k_{P, C} \cdot γ_{C} \cdot P \cdot B_{C}^{d} \\ r_{P, C C}^{a} & = & k_{P, C C} \cdot P \cdot B_{C C}^{a} \\ r_{P, C C}^{d} & = & k_{P, C C} \cdot γ_{C C} \cdot P \cdot B_{C C}^{d} \\ r_{d, P} & = & k_{P} \cdot P \end{matrix} & (19) \end{array}

The proportion of C lost to non-microbial C pools by predation is given by:

\begin{array}{l} f_{P} = f_{P, L} + f_{P, S} + f_{P, M} & (20) \end{array}

Enzymatic breakdown of large and small biopolymers is modeled using Michaelis-Menten kinetics. Oligotrophs control the depolymerisation of large and small polymers, while copiotrophs only affect the depolymerisation of small polymers (Equation 21). This was done to implicitly reflect a higher metabolic versatility of oligotrops than copiotrophs. Copiotrophic cheaters fully rely on the direct uptake of small polymers and monomers and do not affect extracellular depolymerization of polymers:

\begin{array}{l} \begin{matrix} r_{L} = v_{max, L} \cdot \frac{C_{L}}{K_{L} + C_{L}} \cdot B_{O}^{a} \\ r_{S} = v_{max, S} \cdot \frac{C_{S}}{K_{S} + C_{S}} \cdot (B_{O}^{a} + B_{C}^{a}) \end{matrix} & (21) \end{array}

Retardation factors of dissolved small polymers and monomers to consider non-linear equilibrium sorption are calculated as follows (see Jury and Horton, 2004):

\begin{array}{l} \begin{matrix} R_{S} = 1 + f_{S, S} \cdot K_{F, S} \cdot m_{S} \cdot {(\frac{ρ_{B}}{θ})}^{m_{S}} \cdot {C_{S}}^{(m_{S} - 1)} \\ R_{M} = 1 + f_{M, S} \cdot K_{F, M} \cdot m_{M} \cdot {(\frac{ρ_{B}}{θ})}^{m_{M}} \cdot {C_{M}}^{(m_{M} - 1)} \end{matrix} & (22) \end{array}

Effective diffusion coefficients of small polymers and monomers in soil are derived from corresponding aqueous diffusion coefficients by accounting for unsaturated porous media permeability (after Millington and Quirk, 1961):

\begin{array}{l} \begin{matrix} D_{e, S} = \overset{\begin{matrix} Correction term to account \\ for unsaturated permeability \end{matrix}}{\overset{︷}{\frac{θ^{\frac{7}{3}}}{{(1 - \frac{ρ_{B}}{ρ_{S}})}^{2}}}} \cdot \overset{\begin{matrix} aquaeus \\ diffusion \\ coefficient \end{matrix}}{\overset{︷}{D_{S}}} \\ D_{e, M} = \frac{θ^{\frac{7}{3}}}{{(1 - \frac{ρ_{B}}{ρ_{S}})}^{2}} \cdot D_{M} \end{matrix} & (23) \end{array}

Parameterization

Parameters of SpatC and default values used in all simulations are given in Tables 1, 2. We used uniformly distributed initial concentrations of SOM pools and predators. Parameter values were derived from available data if possible and based on logical consideration elsewhere. All microorganisms were assumed to be initially in a dormant state, i.e., initial values of active microorganisms were set to zero. We set a low initial abundance of dormant microbial biomass in the order of 10⁻⁴ mg g⁻¹ (C soil⁻¹) to assure the detection of emerging behavior of microbial groups due to growth in the simulation. Uniform initial SOM pools and a homogeneous medium with isotropic transport and sorption properties were assumed in order to clearly derive effects of spatial distribution of functional microbial groups on C dynamics. Spatial heterogeneity was restricted to microbial distributions. Diffusion coefficients (Table 2) were set to values which reflect a higher molecular weight of C_S than of C_M (Worch, 1993; Hendry et al., 2003).

Parameter values of functional microbial groups were chosen to reflect ecological trade-offs between growth, dormancy and maintenance traits (Figure 2, Table 1). Oligotrophs were parameterized as slowest growers with most efficient substrate uptake and usage. In contrast, copiotrophic cheaters can grow fastest, but are characterized by least efficient substrate uptake and usage. Copiotrophs grow slower than cheaters and have higher maintenance requirements, but are more competitive due to their more efficient substrate uptake in combination with their ability to depolymerize small polymers and inhibit other microorganisms. Oligotrophs were considered to stay active at low substrate supply with lowest maintenance requirements in active state, but highest in dormant state. Copiotrophic cheaters can switch fastest from and to dormancy and switching is triggered already at a low monomer threshold, i.e., they respond fastest to monomer supply. Copiotrophs reactivate and deactivate at a relatively high monomer threshold concentration, but respond much more slowly to insufficient substrate supply than copiotrophic cheaters.

Initialization and Scenario Simulations

Initial pool sizes of dormant functional microbial types were set up in two steps based on a spatial statistical model of microbial biogeography. A Log Gaussian Cox process (LGCP) (Moller et al., 1998) was used as a spatial stochastic model to generate point patterns of microbial cells in a 100 × 100 mm² soil domain. The LGCP model is characterized by three parameters; the mean (μ), the variance (σ²), and the scale (β) of the Gaussian random measure. Following Raynaud and Nunan (2014), an isotropic exponential covariance function C(r) = σ²e^−r/β with distance variable r was used to model the Gaussian process. All parameters were related to the μm-scale. The mean initial density of microbial cells was set to 20 cells mm⁻² (close to the lower limit observed by Raynaud and Nunan, 2014; and Juyal et al., 2019). This is equivalent to an average intensity of the LGCP $λ = e^{μ + \frac{σ^{2}}{2}} = 20 \times 10^{- 6}$ points μm⁻². The spatial heterogeneity of microbial cell distributions was determined by σ² values. Point patterns of increasing spatial heterogeneity and clustering were simulated using four different σ² values: 0.1, 0.5, 2, 6. Corresponding μ values were calculated as $μ = ln (λ) - \frac{σ^{2}}{2}$ to: −10.82, −10.94, −12.82, −28.82. The scale parameter β was fixed to 25 μm in all simulations corresponding to average estimates of Raynaud and Nunan (2014).

The generated point patterns of total microbial cells were then aggregated to 1 mm² resolution by discretizing the 100 × 100 mm² soil domain into 10,000 squares of 1 mm². The total number of cells at 1 mm² resolution was then randomly split into three subsets to derive average cell densities (cells mm⁻²) for the three functional microbial groups (B_O, B_C, B_CC). Initial pool sizes of dormant functional microbial types in mg g⁻¹ (C soil⁻¹) were calculated from these cell densities by assuming a soil bulk density (ρ_S) of 1.2 g cm⁻³, a bacterial cell mass of 10⁻¹¹ mg (McMahon and Parnell, 2014), and a representative layer thickness of 10⁻³ mm (see also Raynaud and Nunan, 2014). Thus, the average total initial microbial biomass was 1.67 × 10⁻⁴ mg g⁻¹ (C soil⁻¹).

In total, 400 simulations comprising 100 realizations per σ² value were performed. All simulations were run for 100 days. This simulation time was chosen as an adequate trade-off between computational effort and process insight.

Technical Implementation

Simulations of the described LGCPs and the aggregation of the generated point patterns were performed using the package spatstat (Baddeley, 2015) and the statistical computing environment R (R Core Team, 2018). The coupled system of partial and ordinary differential equations was implemented and solved using the multipurpose finite element code COMSOL Multiphysics^® in combination with the COMSOL^® module LiveLink™ for MATLAB^®.

Continuous spatial distributions of all state variables were discretized using finite elements. The computational mesh was constructed by converting and refining a regular quadrilateral mesh with 10,000 elements of 1 mm edge length such that every 1 mm square is further discretized by 16 tetrahedral elements (Supplementary Figure 1). As a result the 100 × 100 mm² domain was represented by 160,000 tetrahedral finite elements with an area of 62.5 μm² each. The initial distribution of microorganisms at 1 mm² resolution was mapped on the finite element mesh using nearest-neighbor interpolation as implemented in COMSOL. Test simulations using finer and coarser meshes showed that the chosen mesh resolution provided accurate results at a reasonable computation time.

The equations were solved numerically using an adaptive implicit time-stepping scheme with a backward differentiation formula of varying order from 1 to 5. Newton's method was used to linearize the system of equations. A flexible generalized minimum residual iterative method (Saad, 1993) was used in combination with a geometric multigrid solver (Hackbusch, 1985) to solve the final system of linear equations. The multigrid solver utilized successive over-relaxation for pre- and post-smoothing and a parallel sparse direct method as coarse solver. MATLAB^® was used to set the initial distribution patterns of dormant functional microbial pools, to control the model runs and for post-processing of simulation results.

The derived discrete initial pool sizes of functional microbial groups at 1 mm² resolution could not be directly used for initializing the simulation, because strong differences between individual 1 mm² squares would have required a highly resolved computational mesh for numerical accuracy. Therefore, the initial discrete spatial bacterial distributions were slightly smoothed by running a reduced version of the full model that only simulated slight diffusion of bacterial cells. By this procedure, sharp fronts were removed by an initializing COMSOL model run. The resulting smooth bacterial density fields were then used to initialize the functional microbial types for running the actual SpatC COMSOL^® model.

We explored the effect of biokinetic parameterization by varying some key biokinetic parameters within reasonable bounds by running SpatC with one stochastic realization in a 1 × 10 mm² soil domain (Supplementary Figures 2, 25–27).

Results

Spatiotemporal Dynamics

Spatial clustering of initial microbial communities resulted in the emergence of coupled spatial patterns of C pools and microbial succession (Figure 3, see also Supplementary Figures 4–24 for spatial distributions of soil pools at all degrees of heterogeneity). The spatial distribution of large polymers (Supplementary Figures 6, 7), however, was largely unaffected by microbial distribution. Largely homogenously distributed initial microbial communities (σ² = 0.1) led to a uniform decline of monomers and small polymers. Strong spatial clustering (σ² = 6) induced local depletion zones of small polymers and zones of increased monomer concentrations after 20 and 40 days at spots of high abundance of microbial biomass (Figure 3, Supplementary Figures 25, 26). Higher diffusive transport of monomers compared to small polymers resulted in sharper spatial concentration gradients of small polymers at certain local spots.

FIGURE 3

Figure 3. Microbial biogeography triggers the emergence of spatial patterns of carbon utilization and microbial succession. Each square exemplifies the spatial distribution of C pools (left) and the fraction of microbial functional groups of the total microbial biomass (right) for low (σ² = 0.1) and strong (σ² = 6) initial spatial clustering of microbial communities within a 100 × 100 mm² soil domain for one stochastic realization.

Spatial clustering of initial microbial communities (σ² = 6) led to distinct spots of high microbial abundance. At these spots, also predators became highly abundant (Supplementary Figure 24). The distribution of oligotrophs was characterized by relatively large and more uniformly distributed spots in comparison to the other microbial functional groups (Figure 3). Spots of high abundant copiotrophs were most segregated (Figure 3) and associated with low abundances of the other two functional groups (Supplementary Figure 27). Thus, copiotrophs outcompeted oligotrophs and copiotrophic cheaters and their biomass increased until the end of the simulation period. This pattern emerged as a direct consequence of the simulated interference competition of copiotrophs' inhibition of microbial growth.

Aggregated C Turnover

Heterogeneity in the initial distribution of microbial communities affected aggregated C turnover in soil, but microbial distribution triggered only slight changes in C utilization (Figure 4). For all initial spatial distributions of microorganisms (σ² = 0.1, 0.5, 2, 6), decomposition of small polymers coincided with microbial growth. The concentration of large polymers remained close to the initial value of 10 mg g⁻¹. As a result of microbial death, it showed only a slight increase of <0.015 mg g⁻¹. Monomers showed a concentration peak after about 50 days as a result of enzymatic breakdown of small polymers triggered by the activity of copiotrophs and oligotrophs. While the maximum monomer concentration decreased from homogenous (σ² = 0.1, 0.5) to heterogeneous (σ² = 2, 6) microbial distributions, the monomer concentration peak became broader with increasing spatial clustering. The variability of all C pools increased with increasing spatial heterogeneity of decomposer communities.

FIGURE 4

Figure 4. Microbial biogeography triggers only small changes in carbon utilization. Plots show C turnover dynamics (left) and final values (right) of dissolved monomers and small polymers as well as total microbial biomass C in response to spatial heterogeneity of the initial distribution of microorganisms (σ² = 0.1, 0.5, 2, 6). Values are aggregated over the 100 × 100 mm² soil domain. Lines indicate the medians of 100 realizations and shaded areas show minimum and maximum values (left). Violin plots (right) are scaled to the same width and show the relative distribution of final values. In the inserted box plots, horizontal lines indicate median values, boxes show interquartile ranges (IQR) and whiskers reflect values within maximum 1.5 × IQR.

Moderate spatial clustering (σ² = 2) led to fastest monomer production, degradation of small polymers, and microbial growth. Strong spatial clustering (σ² = 6) resulted in slowest decomposition of small polymers and monomers in combination with the slowest increase in total microbial biomass. As a consequence, final aggregated concentrations of monomers and small polymers were higher and final microbial biomass was lower at σ² = 6 compared to the other scenarios.

Microbial Succession

The aggregated SpatC simulation results revealed a characteristic succession of microbial functional groups in response to available substrates (Figure 5). Copiotrophic cheaters reacted first, directly followed by copiotrophs. These two functional groups simultaneously grew most rapidly on the available monomers and small polymers. Copiotropic cheaters were outcompeted by copiotrophs and oligotrophs as monomers and small polymers became limiting. Copiotrophs switched from active to dormant and maintained the largest portion of their biomass in a dormant state at the end of the simulation. In contrast, active oligotrophs and copiotrophic cheaters showed net growth until the end of the simulation.

FIGURE 5

Figure 5. Microbial biogeography most strongly affects dynamics of fast-growing copiothrophs. Plots show the succession (left) and final values (right) of microbial functional groups (total biomass) in response to spatial heterogeneity of the initial distribution of microorganisms (σ² = 0.1, 0.5, 2, 6). Values are aggregated over the 100 × 100 mm² soil domain. Lines indicate the medians of 100 realizations and shaded areas show minimum and maximum values (left). Violin plots (right) are scaled to the same width and show the relative distribution of final values. In the inserted box plots, horizontal lines indicate median values, boxes show interquartile ranges (IQR), and whiskers reflect values within maximum 1.5 × IQR.

The top-down control by predators played only a minor role. While the median abundance of predators was only slightly affected by microbial distribution, strong spatial clustering of microorganisms resulted in relatively high variability in simulated predator biomass (data not shown).

Moderate spatial clustering (σ² = 2) promoted the growth of copiotrophs and triggered the fastest growth response of copiotrophic cheaters. Strong spatial clustering (σ² = 6) delayed and reduced growth for all microbial functional groups. The variability of all microbial functional groups increased proportional to the initial degree of spatial heterogeneity. Copiotrophs showed the highest sensitivity to spatial heterogeneity of their initial localization. This was evident by the highest variability of the stochastic simulation output compared to oligotrophic and copiotrophic cheaters (Figure 5).

Spatial clustering of microbial communities only slightly affected the relative contribution of functional groups to total biomass (Figure 6, first row). Oligotrophs clearly dominated and were similarly competitive independent of spatial clustering. While copiotrophs reached maximum contribution to total biomass at moderate spatial clustering (σ² = 2), copiotrophic cheaters gained highest maximum contributions at low spatial clustering (σ² = 0.1, 0.5).

FIGURE 6

Figure 6. Spatial clustering of microbial decomposers limits activity and access to monomers by copiotrophic cheaters. Moderate clustering facilitates the access to monomers of copiotrophs and their contribution to total biomass. The first row shows the contribution of microbial functional groups (active and dormant biomass) to total microbial biomass with respect to time. The second and the third row show a phase-space plot of microbial functional group fractions against dissolved monomers (second row) and dissolved small polymers (third row). Each model output is shown in response to the spatial heterogeneity of the initial distribution of microorganisms (σ² = 0.1, 0.5, 2, 6). Lines indicate medians of 100 realizations (aggregated over the 100 × 100 mm² soil domain). Shaded areas (first row) show minimum and maximum values.

The relative contributions of microbial functional groups with respect to dissolved monomer and small polymer concentrations (Figure 6, second and third row) highlights that spatial clustering of microorganisms differently affects the access of microbial functional groups to substrate. Oligotrophs were relatively more competitive at monomer concentrations >0.1 mg g⁻¹ with decreasing spatial clustering and at concentration of small polymers <0.6 mg g⁻¹ with strong spatial clustering (σ² = 6). Copiotrophs benefited most from moderate spatial clustering (σ² = 2) with monomers >0.1 mg g⁻¹ and small polymers <0.75 mg g⁻¹. Copiotrophic cheaters performed best at low spatial clustering (σ² = 0.1, 0.5), independent of substrate concentration.

Discussion

Simulation results indicate that low and moderate initial spatial clustering of microbial decomposers exert some control over the functional composition of microbial communities, whereas the overall C turnover is only slightly affected. Oligotrophs, copiotrophs, and copiotrophic cheaters predominantly act as functionally redundant microbial guilds with respect to decomposition of C compounds. This fits well with conceptual view that C turnover is a “broad” soil process that is carried out by phylogenetically diverse but functionally redundant organisms (Schimel and Schaeffer, 2012; Nunan et al., 2017). Strong spatial clustering of microbial communities, however, induces diffusion-limited C availability at the microhabitat scale which translates into lower decomposition of C compounds and microbial growth at the cm scale. This finding corroborates previous results indicating that the spatial separation of substrates and decomposers can be compensated to a certain degree by shifts in the functional composition of the microbial community (Kaiser et al., 2015), but that if critical diffusion lengths are reached, diffusive transport strongly controls C turnover at the microhabitat scale (Folse and Allison, 2012; Manzoni et al., 2014; Portell et al., 2018).

Oligotrophs are observed to be most competitive regardless of spatial organization. Their competitive advantage results from higher substrate affinities to small polymers and monomers in combination with lower maintenance costs and predation than copiotrophs and copiotrophic cheaters. Copiotrophic cheaters successfully compete with oligotrophs for monomers and small polymers as long as substrate availability remains high enough. They can only sustain relatively low total biomass under unfavorable conditions by switching to dormancy. Interestingly, our results suggest that moderate spatial heterogeneity (σ² = 2) is beneficial for copiotrophs. Moderate spatial clustering induces the formation of large areas of high monomer concentration by extracellular decomposition of small polymers. Copiotrophs become active and grow rapidly under relatively high concentrations of monomers while inhibiting the growth of other microorganisms. Thus, relatively more micro-environments of competitive advantage to copiotrophs against oligotrophs and copiotrophic cheaters are created in comparison to lower and higher spatial clustering. In addition, copiotrophs sustain themselves under less beneficial conditions by quickly switching to a dormant state, which drastically reduces maintenance costs and biomass decay by predation.

The simulated behavior of microbial functional groups supports experimental evidence of the importance of metabolic activation/deactivation strategies by microbial functional groups for regulating C turnover in soil (Placella et al., 2012; Joergensen and Wichern, 2018; Salazar et al., 2019). Our finding that interactions between microbial functional groups are controlled by the spatial localization of microorganisms is in agreement with previous results from individual-based modeling (Allison, 2005; Kaiser et al., 2015; Portell et al., 2018). SpatC model results, however, suggest a less severe impact of cheaters on microbial functioning and C turnover. In addition, our approach is able to considerably extend the total spatial dimension typically covered by individual-based modeling approaches (Allison, 2005; Folse and Allison, 2012; Kaiser et al., 2015) by several orders of magnitude, from ≤1 mm² to 100 cm². The INDISIM-SOM model (Banitz et al., 2015) is conceptually similar to SpatC but adopts an individual-based approach. INDISIM-SOM-NL simulates C and N dynamics in soil and splits the spatial domain into 30 × 30 grid cells of 310 μm, each containing two functional groups of microorganisms: heterotrophs and autotrophs. Despite assuming an individual-based approach, the spatial extend of INDISIM-SOM-NL is similar to the one covered in the present contribution. This is achieved by using the “superindividual” formalism, which assumes fewer individuals summarizing the properties and being representative of larger ensembles of microorganisms. Nonetheless, compared to INDISIM-SOM-NL, SpatC provides a higher temporal and spatial resolution, and considers three functional types of heterotrophic microorganisms.

SpatC scenario simulations provide predictions of the emergent macroscopic (cm) microbial and C dynamics resulting from small-scale (mm) distribution characteristics of microbial functional decomposer communities. Microbial biogeography at the microhabitat scale (μm) is thereby considered by using a spatial stochastic model to derive microbial distribution patterns at the μm-scale, which are aggregated to mm-scale distributions of microbial communities. The combination of statistical and process-based modeling applied with SpatC provides an upscaling approach that can consider feedbacks between microhabitats, microbial communities and soil microbial and physical processes up to the pedon scale. Hence, our study contributes to resolving the challenge of upscaling microbial regulation mechanisms from the microhabitat scale to larger scales relevant for soil management and global environmental change (Baveye et al., 2018). SpatC predictions of microbial and C dynamics are, however, dependent on the assumed biokinetic rate laws at the mm-scale, which have been shown to differ from rate laws at μm-scale (Chakrawal et al., 2019; Wang and Allison, 2019). Similarly, an exploratory analysis of SpatC revealed a very high sensitivity of model dynamics to some key biokinetic parameters that partly increased the observed mild effect of spatial heterogeneity (Supplementary Figures 28–30). The maximum growth rate of oligotrophs (μ_{max, O}) strongly controlled the coexistence of functional microbial pools independent of spatial heterogeneity (Supplementary Figure 28). Higher values of μ_{max, O} resulted in the fast depletion of small polymers and monomers and strongly reduced growth of copiotrophs and copiotrophic cheaters. In contrast, increasing depolymerisation rate coefficients of small polymers induced dominant growth of copiotrophs in combination with stronger effects of spatial distributions of microorganisms on microbial and C dynamics (Supplementary Figure 28). Parameters related to dormancy of oligotrophs and copiotrophic cheaters had a minor impact on model results, only a very high dormancy concentration threshold slightly affected dynamics (Supplementary Figure 29). As expected, less maintenance requirements of oligotrophs and less severe inhibition by copiotrophs resulted in higher biomass of oligotrophs in combination with faster and more complete degradation of small polymers and monomers (Supplementary Figure 30), but the effect of spatial heterogeneity did not significantly change.

The developed SpatC model considers the control of C turnover by spatial heterogeneity of functional microbial groups. However, SpatC currently simplifies the micro-scale distribution of organic C, which probably has a strong impact on C dynamics at larger scales. The simulated spatial patterns in organic matter decomposition are in alignment with experimentally observed patterns of extracellular enzyme activity (Kravchenko et al., 2019). Experimental evidence further suggests that C turnover is strongly determined by pore characteristics (Kravchenko and Guber, 2017; Juyal et al., 2019) and microbial activity is highest in pores between 10 and 300 μm (Kravchenko et al., 2019). Thus, an improved description of microbial C turnover could be gained by integrating realistic descriptions of soil pore structure based on X-ray computed tomography data (see e.g., Portell et al., 2018) in combination with a meaningful correlation structure of substrate and microbial group distribution using evidence-based spatial statistical modeling. In addition, the representation of biological community interactions remains limited. Crucial extensions could include the explicit representation of enzyme dynamics (Burns et al., 2013; Moyano et al., 2018; Wang and Allison, 2019) and the implementation specific fungal traits (Yang and van Elsas, 2018). Similarly, microbial dispersal and chemotactic behavior (Valdés-Parada et al., 2009; see e.g., Gharasoo et al., 2014; Locey et al., 2017; König et al., 2018) should be included in future. Other promising extensions are quorum sensing (Williams et al., 2007; Melke et al., 2010; Mund et al., 2016; McBride and Strickland, 2019; Schmidt et al., 2019) as regulator of biological interactions, as well as to improve the modeling of top-down control of microbial communities by predators and viruses (Pratama and van Elsas, 2018; Thakur and Geisen, 2019). Extensions along these lines will provide further insights into the biological controls on soil organic matter turnover by generating model-based hypotheses that can be tested against experimental evidence.

Soil organic matter formation is an emergent process. It cannot be directly predicted from community composition, but arises from non-linear feedbacks and interactions between microbial community members. To understand and predict these biogeochemical feedbacks it is crucial to combine microbial traits with the spatial arrangements between microorganisms in their micro-environment and their corresponding substrate. A key finding of our work is that the degree of spatial heterogeneity of microbial communities may control the relative contribution of functional microbial groups to biogeochemical processes and the degree of functional redundancy within microbial communities. Our simulation results suggest that metabolic activation/deactivation strategies of microbial functional groups may be a key control of C turnover in soil. These model-based implications could be tested with targeted experiments that enable spatially resolved measurements of microbial community composition and C fluxes at the microhabitat scale by extending existing approaches (e.g., Poll et al., 2006) and using novel techniques such as flow cells (Krueger et al., 2018) in combination with functional multilayered omics approaches (Jansson and Hofmockel, 2018; Sergaki et al., 2018).

Data Availability Statement

The raw data supporting the conclusions of this study will be made available by the authors, without undue reservation, to any qualified researcher.

Author Contributions

MU, EK, CP, TS, VS, and HP contributed the conception and design of the study. BK and VS performed Log Gaussian Cox Process simulations and provided aggregated point patterns. HP developed the model, performed the simulations, analyzed the data, and wrote the first draft of the manuscript. All authors contributed to manuscript revision, read, and approved the submitted version.

Funding

This study was financially supported by the Ellrichshausen Foundation and the Collaborative Research Center 1253 CAMPOS (DFG Grant Agreement SFB 1253/1 2017). The authors acknowledge support by the state of Baden-Württemberg through bwHPC.

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Acknowledgments

The authors are grateful for the excellent and stimulating feedback from three reviewers on the first manuscript version.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fenvs.2020.00002/full#supplementary-material

References

Addiscott, T. M. (2010). Soil mineralization: an emergent process? Geoderma 160, 31–35. doi: 10.1016/j.geoderma.2010.03.016

ORIGINAL RESEARCH article

Spatial Control of Carbon Dynamics in Soil by Microbial Decomposer Communities

Introduction

Materials and Methods

Model Rationale and Main Assumptions

Governing Equations and Fluxes

Non-microbial Carbon

Functional Microbial Groups

Fluxes and Functions

Parameterization

Initialization and Scenario Simulations

Technical Implementation

Results

Spatiotemporal Dynamics

Aggregated C Turnover

Microbial Succession

Discussion

Data Availability Statement

Author Contributions

Funding

Conflict of Interest

Acknowledgments

Supplementary Material

References

This article is part of the Research Topic

People also looked at