Microscale Heterogeneity of the Spatial Distribution of Organic Matter Can Promote Bacterial Biodiversity in Soils: Insights From Computer Simulations

There is still no satisfactory understanding of the factors that enable soil microbial populations to be as highly biodiverse as they are. The present article explores in silico the hypothesis that the heterogeneous distribution of soil organic matter, in addition to the spatial connectivity of the soil moisture, might account for the observed microbial biodiversity in soils. A multi-species, individual-based, pore-scale model is developed and parameterized with data from 3 Arthrobacter sp. strains, known to be, respectively, competitive, versatile, and poorly competitive. In the simulations, bacteria of each strain are distributed in a 3D computed tomography (CT) image of a real soil and three water saturation levels (100, 50, and 25%) and spatial heterogeneity levels (high, intermediate, and low) in the distribution of the soil organic matter are considered. High and intermediate heterogeneity levels assume, respectively, an amount of particulate organic matter (POM) distributed in a single (high heterogeneity) or in four (intermediate heterogeneity) randomly placed fragments. POM is hydrolyzed at a constant rate following a first-order kinetic, and continuously delivers dissolved organic carbon (DOC) into the liquid phase, where it is then taken up by bacteria. The low heterogeneity level assumes that the food source is available from the start as DOC. Unlike the relative abundances of the 3 strains, the total bacterial biomass and respiration are similar under the high and intermediate resource heterogeneity schemes. The key result of the simulations is that spatial heterogeneity in the distribution of organic matter influences the maintenance of bacterial biodiversity. The least competing strain, which does not reach noticeable growth for the low and intermediate spatial heterogeneities of resource distribution, can grow appreciably and even become more abundant than the other strains in the absence of direct competition, if the placement of the resource is favorable. For geodesic distances exceeding 5 mm, microbial colonies cannot grow. These conclusions are conditioned by assumptions made in the model, yet they suggest that microscale factors need to be considered to better understand the root causes of the high biodiversity of soils.

Introduction

During the last decade, soils have become increasingly central to a number of crucial debates on issues of great societal concern (e.g., Baveye et al., 2018). Because soils contain a very large stock of carbon, there is a risk that, with rising ambient temperatures associated with global climate change, soils will release vast amounts of greenhouse gases and thereby accelerate change. Biodiversity losses have also emerged as a major concern in many parts of the world. In this context, it is not surprising that in recent years, there has been a significant surge of interest into the biodiversity of soils, and the effect it has on traditional soil functions (Nannipieri et al., 2003).

Many aspects of the biodiversity of soils have proven very difficult to understand. Soils are highly complex media in which a huge number of bacteria, archaea, and fungi live. In a single gram of soil, it is not exceptional to find as many as 10¹⁰ bacterial cells and 5 × 10⁴ species (Roesch et al., 2007), with commensurate numbers found for other microorganisms. To the extent that many microorganisms (an estimated 98.5% in the case of bacteria) have never been isolated or characterized, the measurement of soil biodiversity itself raises a number of fundamental questions (e.g., Nannipieri et al., 2003; Baveye et al., 2016a,b, 2018). Functionally, it is not clear at all to what extent this very large diversity of soil microbial populations is crucial and whether it needs to be preserved at all cost. Experimental results are contradictory in this respect. Whereas, e.g., Philippot et al. (2013) show that the loss of biodiversity in soils decreases denitrification activity and nitrogen cycling, the experimental results of Werts et al. (2006) suggest on the contrary that biogeochemical functions of soil such as carbon mineralization and denitrification are not impacted by a reduction of microbial diversity. Powlson et al. (2017) have recently described as an unresolved 60-year old paradox the fact that CHCl₃ fumigation, wiping out 90% of the soil microbial population and modifying drastically its diversity, does not appear to have an effect on soil organic carbon (SOC) mineralization in soils that have a pH above 5.5. SOC mineralization continues at the same rate, after fumigant removal, once the initial decomposition flush is over (Powlson et al., 2017).

A similarly high uncertainty surrounds the features of soils that allow such a large microbial diversity to exist in the first place. Some researchers consider that diversity is mainly caused by biotic interactions between cells (Hanson et al., 2012), but experimental observations increasingly suggest that a high biodiversity is associated with soil spatial heterogeneity (Rainey and Travisano, 1998) and is caused by biotic and abiotic interactions taking place in the soil architecture. Yet the exact mechanisms involved remain elusive. The often advocated explanation that the heterogeneous, disconnected distribution of moisture in unsaturated soils causes distinct groups of microorganisms to be physically isolated from each other (Treves et al., 2003; Long and Or, 2009) is appealing, but it does not apply to fungi or filamentous bacteria (Baveye et al., 2016b) and cannot account by itself for the biodiversity of soils that are periodically saturated after rainfall events. At this point, there is no real, satisfactory explanation of how the spatial heterogeneity of soils might foster the biodiversity of their microbial populations.

In this general context, the key objective of the present article is to explore the hypothesis that the heterogeneous distribution of the basic nutrient resources used by bacteria in soils can account to some extent for their diversity. The spatial distribution of organic matter in soils is known to be highly heterogeneous. Incorporation of plant residues by tillage results in patchy distribution at the centimeter scale (e.g., Elyeznasni et al., 2012) while at smaller millimeter scales, heterogeneous distribution of soil organic matter has been visualized by Peth et al. (2014) and Kravchenko et al. (2014). On the basis of these microscale observations, the effect of the heterogeneity of the spatial distribution of soil organic matter and of the connectivity of the aqueous phase on bacterial biodiversity was examined in a series of scenarios using a 3D pore-scale carbon dynamics model. Bacterial cells of three strains of the Arthrobacter sp. including highly competitive-, generalist-, and poorly competitive-strains, were randomly placed within the 3D pore space of a small soil sample (of volume size of 314 mm³) imaged at a resolution of 68 μm, suitable to visualize meso- and macro-pores.

Materials and Methods

Soil Image

Undisturbed soil cores were sampled in the surface horizon of a cultivated soil, a silty loamy (19% clay, 75% silt, 6% sand) Albeluvisol (Vogel et al., 2015). 3D images of the samples were obtained using an X-ray CT scanner (HMX 225, NIKON metrology, Tring, United Kingdom). A global threshold according to Elyeznasni et al. (2012) was used to obtain binary images in which the voxels of the gray CT image were classified either as soil or void voxels. We selected one sub-image (called G6 in Vogel et al., 2015) of 100³ voxels size out of the set of segmented CT images. The voxel-resolution of the image is 68 μm, so that the pore space explored in this study encompasses most of the structural porosity made of meso- and macro-pores. This image corresponds to a volume size of about 314 mm³, and it has a porosity of 18.82%. This number is undoubtedly smaller than the actual porosity of the soil, because of the fact that sub-resolution pores are ignored (Baveye et al., 2017).

The localization of the fluid and gas voxels corresponding to a given water saturation index, S_w (the proportion of the pore space filled with water), is calculated using a two-phase two-relaxation-times (TRT) lattice-Boltzmann model (LBM) as described by Genty and Pot (2013) and Pot et al. (2015). Three levels of water saturation of the CT-visible porosity are assumed in the present work, S_w = 1.00, S_w = 0.50, and S_w = 0.25. After a visual inspection of the 3D distribution of the gas phase in the images, we selected a few of the smallest visible pores containing gas and recalculated from the Young-Laplace law a rough estimate of the matric potential for S_w = 0.25 and S_w = 0.50. The matric potentials estimated in this manner are about -0.6 and -0.3 kPa, respectively. Therefore, even though a water saturation level of 0.25 would suggest that the soil is relatively dry, the fact that this number refers only to the CT-visible pore space means that the scenarios reported in the present contribution correspond to wet conditions.

The image offers a tradeoff between resolution and sample size (constrained by X-ray computed tomography). It is a compromise that presents the advantage of reproducing the millimeter-scale variability of the microbial activity, as reported by Vieublé-Gonod et al. (2006), and of enabling us to work with a soil volume that is large enough for carbon mineralization to be measured in practice.

Model Description

The description below of the Ib-LBioS-Comp model follows the overview, design concepts, and details (ODD) protocol (Grimm et al., 2006, 2010), which was especially developed to communicate the features of agent and individual-based models.

Purpose

Ib-LBioS-Comp, which stands for Individual-based LBM for biodegradation affected by soil structure and competition, is a pore scale modeling approach developed to study the impact of the soil structure or, more appropriately, architecture (see discussion in Baveye et al., 2018), molecular diffusion, spatial heterogeneity of resource distribution, and competition among bacterial species on the biodegradation dynamics of organic matter in the soil.

Entities, State Variables, and Scales

Ib-LBioS-Comp combines a lattice-Boltzmann solute diffusion model and an individual-based biological module describing bacterial activity. The model involves single bacterial cells (biotic agents) of up to three different species or strains, dissolved organic carbon or DOC (abiotic lattice-Boltzmann populations), and particulate organic matter (POM, abiotic agents). Bacterial cells and DOC are distributed in the 3D volume of the soil pore volume and POM is distributed in the 3D volume of the soil solid matrix. Enzymatic hydrolysis of POM is assumed to continuously supply DOC to the liquid phase. Then DOC diffuses through the liquid phase toward the microenvironments where it is taken up by bacteria and used as nutrient source. The diffusion and uptake of oxygen are not accounted for explicitly in the model at the moment and are assumed not to limit microbial activity.

The simulated space is divided into a regular 3D grid made of cubic voxels that can be either solid or filled with a fluid (air or water).

A bacterial cell or individual (I_i) is defined by the variables Π_i, identifying its position in the domain, its species j_i, its mass m_i (mg C), its specific uptake rate υ_DOC,i (tu^-1), and its mass at reproduction m_R,i (mg C). All masses in Ib-LBioS-Comp are expressed in terms of mass of carbon. Letting P = I(t) denote the number of bacteria at time t, one obtains for the state of the population at time t:

$\mathit{P}\mathit{n}\mathrm{=}{\mathrm{\{}{\mathit{I}}_{\mathit{i}}\mathit{[}{\mathit{\Pi }}_{\mathit{i}}\mathit{,}{\mathit{j}}_{\mathit{i}}\mathit{,}{\mathit{m}}_{\mathit{i}},{\mathrm{\upsilon }}_{\mathit{D}\mathit{O}\mathit{C}\mathit{,}\mathit{i},}{\mathit{m}}_{\mathit{R}\mathit{,}\mathit{i}}\mathit{]}\mathit{\}}}_{\mathit{i}\mathrm{=}\mathrm{1}\mathrm{,}\mathrm{2}\mathrm{,}\mathrm{...}\mathrm{,}\mathit{n}\mathrm{(}\mathit{t}\mathrm{)}}\mathrm{\left(}\mathrm{1}\mathrm{\right)}$

where i is the index of individual bacteria and n(t) the total number of cells at time t.

A POM agent (A_l) is defined by a variable indicating its position in the domain, Π_l, and its carbon mass, m_POM,l (mg C). Letting A(t) denote the number of POM agents, the POM population state at time t is:

${\mathit{P}\mathit{O}\mathit{M}}_{\mathit{m}}\mathrm{=}{\mathrm{\{}{\mathit{A}}_{\mathit{l}}\mathrm{[}{\mathrm{\Pi }}_{\mathit{l}}\mathrm{,}{\mathit{m}}_{\mathit{P}\mathit{O}\mathit{M}\mathit{,}\mathit{l}}\mathrm{]}\mathrm{\}}}_{\mathit{l}\mathrm{=}\mathrm{1}\mathrm{,}\mathrm{2}\mathrm{,}\mathrm{...}\mathrm{,}\mathit{m}\mathrm{(}\mathit{t}\mathrm{)}}\mathrm{\left(}\mathrm{2}\mathrm{\right)}$

where l is the index of the POM agents and m(t) the total number of POM agents at time t.

The DOC solute is simulated by microscopic lattice-Boltzmann populations, f_q, that are microscopic solute entities defined in the Q microscopic directions at each liquid node of the 3D grid. The Q directions are defined by the unit microscopic velocity vectors, $\overrightarrow{{\mathit{C}}_{\mathit{q}}}\mathrm{=}\mathrm{\left\{}{\mathit{C}}_{\mathit{q}\mathrm{\alpha }}\mathrm{\right\}}\mathit{q}\mathrm{=}\mathrm{0}\mathrm{,}\mathrm{...}\mathrm{,}\mathit{Q}\mathrm{-}\mathrm{1}\mathrm{;}\mathrm{\alpha }\mathrm{=}\mathrm{1}\mathrm{,}\mathrm{...}\mathrm{,}\mathit{d}$ where d is the dimension of the grid or lattice. We used the model D3Q7.

The DOC concentration in the liquid phase, C_DOC (t) (mg C lu^-3) can be calculated at each liquid node of the grid as the sum of the f_q populations:

${\mathit{C}}_{\mathit{D}\mathit{O}\mathit{C}}\mathrm{=}\frac{\mathrm{1}}{\mathrm{\Delta }{\mathit{V}}_{\mathit{x}\mathit{y}\mathit{z}}}\underset{\mathit{q}\mathrm{=}\mathrm{0}}{\overset{\mathit{Q}\mathrm{-}\mathrm{1}}{\mathrm{\Sigma }}}{\mathit{f}}_{\mathit{q}}\mathrm{}\mathrm{\left(}\mathrm{3}\mathrm{\right)}$

with ΔV_xyz the volume of one voxel expressed in lu³ where lu is the spatial unit of the lattice, in our case determined by the scanning resolution so that 1 lu = 68 μm.

The temporal evolution of the system is divided into equal intervals associated with time steps or units (tu) of a time step length dictated by the lattice-Boltzmann submodel. The temporal extent of the simulations was set to 10 days according to previous simulations made by Vogel et al. (2015) in which exponential growth and decline of biomass were observed within this duration. The time step length is 3.44 s (see Sections “Abiotic Processes” and “Model Parameterization”).

Process Overview and Scheduling

Global simulation comprises three sections (Figure 1): (i) the initialization of the simulated system, (ii) the time step loop, which is repeated until the end of the defined time steps, and (iii) the model output section, where the system-level (aggregated) and individual-level (non-aggregated) data are saved in files. Initialization of the system includes: reading of model parameters, initialization of the bacterial agents, and initialization of the LBM parameters and DOC populations. The time step loop includes, chronologically: (ii.i) storage of the simulation state in temporary data structures, (ii.ii) the POM agents actions loop, (ii.iii) the bacterial actions loop (Figure 2), and (ii.iv) the lattice-Boltzmann actions. Output files of aggregated and state variables are created from the temporary data structures saved previously.

FIGURE 1

FIGURE 1. General workflow of the Ib-LBioS-Comp model.

FIGURE 2

FIGURE 2. General workflow of the bacterial agents loop of the Ib-LBioS-Comp model.

At each time step, bacterial cells perform the following set of actions: uptake, metabolism, reproduction, and mortality. The order in which bacteria act is changed randomly every time step to avoid privileging always the same first-acting bacteria. At each time step, the existing POM agents undergo hydrolysis to produce DOC. The DOC lattice-Boltzmann populations, f_q, are then updated through the following set of actions: collision, propagation, and bounce-back when they encounter a solid or a gas neighbor. This last action is motivated by the premise that DOC occurs only in the water phase.

Design Concepts

Emergence

Bacterial processes (uptake, metabolism, reproduction, mortality) are defined at a single-cell level and the population or system level behavior emerges from the interactions among individuals and between the individuals and the media. The main population-level emerging characteristics are the population density (bacteria present in the media), the population biomass, the DOC taken up by the population, the CO₂ produced by the population and the bacterial size distribution. Solute diffusion processes (collision, propagation) are defined at the microscopic level of the Q directions of the lattice, and the solute behavior (diffusion) emerges from the interactions among the solute lattice-Boltzmann populations.

Interaction

Bacteria are considered to be immobile but they compete directly for space through the maximum carrying capacity of a voxel. Therefore, the presence of other bacterial cells in the local space can directly affect the fate of new-born bacteria. Indirect interactions among individual bacteria arises from competition for the DOC available in the local environment.

Stochasticity

Randomness is considered when the rules are applied to individuals by using probabilistic distributions to deal with or manage individual events. Random processes or events include the assignment of the position of new-born bacteria in the physical domain near the mother when the grid element of the mother reached maximum occupancy, and the occurrence for cell death. The sequence of actions of the simulated bacteria changes randomly at each time step to avoid privileging one over the others. The model can also introduce further stochasticity when setting the mass of the initial individuals, the individual specific uptake rate, and the mass at reproduction of the individual bacteria using folded normal distributions, but this is not considered in the current study.

Observation

Global and single-cell outputs are recorded at the beginning of the simulation, at regular intervals during the simulation, as well as at the end of the simulation. Global variables calculated at the scale of the entire domain include: mass of DOC, POM, and CO₂ produced in the media, and, for each of the three bacterial species, the number of bacteria, and total bacterial biomass. The state of all the individual bacteria is saved at sampling times. These single-cell data include the position within the domain, species, mass, uptake rate, and mass at reproduction of the individuals. The carbon mass and position of all the POM spots (abiotic agents) are also saved at sampling times. The DOC concentration of the liquid voxels holding at least one bacterium is also recorded at sampling times. The final DOC concentration of all the liquid grid cells of the domain is recorded at the end of the simulation.

Initialization

The specific uptake rate (v_DOC,i), and the reproduction mass (m_R,i) of the initial individuals (i = 4,...,N_B0) are assumed to be specific for each jth species. For the individuals of each species j, these properties are set using the model parameters ${\mathit{v}}_{\mathit{D}\mathit{O}\mathit{C}}^{\mathit{j}}$ and ${\mathit{m}}_{\mathit{t}\mathrm{0}}^{\mathit{j}}$, respectively, for the uptake rate and the reproduction mass. No intraspecific variability is considered in the present study. Similarly, the initial mass of the individuals starting the simulations (m_i) is initialized according to the model parameters ${\mathit{m}}_{\mathit{t}\mathrm{0}}^{\mathit{j}}$.

Since it is generally assumed that bacteria in soil microenvironments tend to be sorbed to, or be at least very near, solid surfaces, the model assumes that bacteria can be located only in liquid voxels having at least one solid neighbor. What defines a neighbor here is the particular lattice-Boltzmann connectivity that is adopted in the calculations (D3Q7 in this case). The initial N_B0 bacterial cells are randomly distributed among a number of bacterial spots (N_SP0T) that, in turn, are randomly chosen, with replacement, from the liquid voxels having one or more solid neighbors.

The POM agents are situated in the solid matrix of the soil. The initial POM agents are randomly distributed among solid voxels that have at least one liquid neighbor.

An initial amount of dissolved organic carbon, DOC₀, is distributed homogeneously among the liquid voxels of the image.

Input Data

The model uses soil structural data as described in see Section Soil Image.

Biological Processes

Several separate submodels describe quantitatively the uptake, metabolism, reproduction, and mortality, respectively, of individual bacteria.

In the Uptake submodel, the uptake (U_i) of carbon substrate by bacterium i, belonging to the species j = j_i, is given by the equation depending on the mass m_i (t) of the bacterium

${\mathit{U}}_{\mathit{i}}\mathrm{\left(}\mathit{t}\mathrm{\right)}\mathrm{=}\frac{{\mathit{v}}_{\mathit{D}\mathit{O}\mathit{C}\mathit{,}\mathit{i}}{\mathit{C}}_{\mathit{D}\mathit{O}\mathit{C}}^{\mathrm{\left(}\mathit{t}\mathrm{\right)}}}{{\mathit{C}}_{\mathit{D}\mathit{O}\mathit{C}}^{\mathrm{\left(}\mathit{t}\mathrm{\right)}}\mathrm{+}{\mathit{k}}_{\mathit{D}\mathit{O}\mathit{C}}^{\mathit{j}}}{\mathit{m}}_{\mathit{i}}\mathrm{\left(}\mathit{t}\mathrm{\right)}\mathrm{\left(}\mathrm{4}\mathrm{\right)}$

where v_DOC,i is the specific uptake rate of the ith bacterium (tu^-1), and the parameter ${\mathit{k}}_{\mathit{D}\mathit{O}\mathit{C}}^{\mathit{j}}$ is the half saturation constant (mg C lu^-3) of the jth species.

In the Metabolism submodel, the mass of carbon taken up is used by the cell to create new biomass. Since catabolic reactions need energy, due to the respiration process, a fraction of the carbon that is taken up is released again to the media in the form of CO₂ emissions. Accordingly, the growth of the bacterial cell is modeled using the following equation:

${\mathit{m}}_{\mathit{i}}\mathrm{(}\mathit{t}\mathrm{+}\mathrm{1}\mathrm{)}\mathrm{=}{\mathit{m}}_{\mathit{i}}\mathrm{\left(}\mathit{t}\mathrm{\right)}\mathrm{+}{\mathit{U}}_{\mathit{i}}\mathrm{\left(}\mathit{t}\mathrm{\right)}\mathrm{-}{\mathit{k}}_{\mathit{r}}^{\mathit{j}}{\mathit{m}}_{\mathit{i}}\mathrm{\left(}\mathit{t}\mathrm{\right)}\mathrm{\left(}\mathrm{5}\mathrm{\right)}$

where ${\mathit{k}}_{\mathit{r}}^{\mathit{j}}$ is the respiration rate (tu^-1) of the species j. Equations (4) and (5) assume that growth respiration (traditionally calculated as a function of substrate during uptake) and maintenance respiration (proportional to the biomass) are not distinguished.

The Reproduction submodel adopts a simple bipartition condition (Banitz et al., 2015) according to which the bacterial cell has to attain a specific individual reproduction mass before dividing. Every time step, the mass of a bacterium is compared to m_R,i (mg C), the reproduction mass of the ith individual. If m_i(t + 1) > m_R,i, the mass of the bacterium is halved, and simultaneously a daughter cell of the same mass is created. If the number of bacteria occupying the voxel of the mother cell is less than N_VOX, the maximum carrying capacity of a voxel, the cell is created in the current voxel. Otherwise, the simulator chooses randomly a voxel that has not attained the maximum carrying capacity. The searching algorithm looks progressively to neighbors situated at increasing distances (in voxels) from the mother’s voxel of origin, until all the tridimensional space is inspected. If all the voxels reach the maximum carrying capacity, the simulation stops. The daughter cell remains active but does not act until the next time step is reached. The specific uptake rate of the new-born individual (v_DOC,i) and its mass at reproduction (m_R,i) are inherited from the mother.

Finally, the Mortality submodel accounts for bacterial cell death derived from internal and external events (e.g., predation by other organisms). The cell cannot survive anymore due to internal events when the cell size decreases below the minimal cell size characteristic of its species j,${\mathit{m}}_{\mathit{M}\mathit{I}\mathit{N}}^{\mathit{j}}$ which can be attained due to a starvation process. Cell death due to external events is accounted by a probability, ${\mathit{p}}_{\mathit{M}}^{\mathit{j}}$(dimensionless), independent of the cell state. At every time step, the submodel compares r_p with ${\mathit{p}}_{\mathit{M}}^{\mathit{j}}$, where r_p is a random realization coming from a uniform distribution between 0 and 1. If ${\mathit{r}}_{\mathit{p}}\mathrm{>}{\mathit{p}}_{\mathit{M}}^{\mathit{j}}$ the individual bacterium dies. The cell carbon lyses and creates new DOC in the current voxel.

Abiotic Processes

The main abiotic processes simulated by the model are the hydrolysis of POM and the diffusion of DOC in the 3D pore space. No convective movement of DOC was considered in the model. POM agents release DOC, decreasing as a result the mass of the POM agent. The model assumes only one homogeneous fraction of POM with a unique hydrolysis rate. The hydrolysis process is modeled assuming a first order kinetics of constant rate, k_POM. The underlying hypothesis is that exo-enzymes are ubiquitous in soil (Folse and Allison, 2012). Nannipieri et al. (2003) reported that newly produced exo-enzymes by bacterial cells are short-lived molecules because, for instance, proteases can degrade them. Thus, the ubiquitous enzymes are probably those physically protected though their adsorption to clay particles or humic molecules (Burns, 1982). Burns (1982) considers it a likely scenario that these enzymes become active when a POM fragment comes into contact with them. Then, the DOC released by the agent is distributed equally among the liquid voxels neighboring the solid voxels containing the POM agent.

We implemented the TRT lattice-Boltzmann approach of Ginzburg (2005). The evolution equation of the DOC microscopic entities at the liquid nodes (grid cell with k = lq), $\overrightarrow{{\mathit{V}}_{\mathit{x}\mathit{y}\mathit{z}}}$ from time t to t+ 1 is given by:

${\mathit{f}}_{\mathit{q}}\mathrm{(}\overrightarrow{{\mathit{V}}_{\mathit{x}\mathit{y}\mathit{z}}}\mathrm{+}\overrightarrow{{\mathit{c}}_{\mathit{q}}}\mathrm{,}\mathit{t}\mathrm{+}\mathrm{1}\mathrm{)}\mathrm{-}{\mathit{f}}_{\mathit{q}}\mathrm{(}\overrightarrow{{\mathit{V}}_{\mathit{x}\mathit{y}\mathit{z}}}\mathrm{,}\mathit{t}\mathrm{)}\mathrm{=}{\mathrm{\lambda }}_{\mathit{e}}\mathrm{\left[}{\mathit{f}}_{\mathit{q}}^{\mathrm{+}}\mathrm{\right(}\overrightarrow{{\mathit{V}}_{\mathit{x}\mathit{y}\mathit{z}}}\mathrm{,}\mathit{t}\mathrm{)}\mathrm{-}{\mathit{e}}_{\mathit{q}}^{\mathrm{+}}\mathrm{(}\overrightarrow{{\mathit{V}}_{\mathit{x}\mathit{y}\mathit{z}}}\mathrm{,}\mathit{t}\mathrm{\left)}\mathrm{\right]}\mathrm{+}{\mathrm{\lambda }}_{\mathit{o}}\mathrm{\left[}{\mathit{f}}_{\mathit{q}}^{\mathrm{-}}\mathrm{\right(}\overrightarrow{{\mathit{V}}_{\mathit{x}\mathit{y}\mathit{z}}}\mathrm{,}\mathit{t}\mathrm{)}\mathrm{-}{\mathit{e}}_{\mathit{q}}^{\mathrm{-}}\mathrm{(}\overrightarrow{{\mathit{V}}_{\mathit{x}\mathit{y}\mathit{z}}}\mathrm{,}\mathit{t}\mathrm{\left)}\mathrm{\right]}\mathrm{+}{\mathit{S}}_{\mathit{q}}\mathrm{}\mathrm{\left(}\mathrm{6}\mathrm{\right)}$

in which the collision and propagation steps are described, respectively, by the two first terms of right hand side and left hand side of (7), respectively. The sink/source term of DOC, S_q, is calculated from the hydrolysis and bacterial processes. In the TRT scheme, the microscopic entities f_q are decomposed into symmetric, ${\mathit{f}}_{\mathit{q}}^{\mathrm{+}}$ and antisymmetric components, ${\mathit{f}}_{\mathit{q}}^{\mathrm{-}}$ along their opposite velocities $\overrightarrow{{\mathit{c}}_{\mathit{q}}}\mathrm{=}\mathrm{-}{\overrightarrow{\mathit{c}}}_{\overline{\mathit{q}}}$ (Ginzburg, 2005). During the collision step, the relaxation of moments resulting from the entities’ distribution at time t toward an equilibrium state ${\mathit{e}}_{\mathit{q}}\mathrm{=}{\mathit{e}}_{\mathit{q}}^{\mathrm{+}}\mathrm{+}{\mathit{e}}_{\mathit{q}}^{\mathrm{-}}$ governs the reorganization of the entities. The relaxation parameter λ_e is a free parameter and the relaxation parameter λ_o is related to the molecular diffusion coefficient, ${\mathit{D}}_{\mathit{m}}^{\mathit{L}\mathit{B}\mathit{M}}$ (lu²tu^-1) along:

${\mathit{D}}_{\mathit{m}}^{\mathit{L}\mathit{B}\mathit{M}}\mathrm{=}\frac{\mathrm{-}\mathrm{1}}{\mathrm{3}}\mathrm{(}\frac{\mathrm{1}}{\mathrm{2}}\mathrm{+}\frac{\mathrm{1}}{\lambda \mathit{o}}\mathrm{)}\mathrm{}\mathrm{\left(}\mathrm{7}\mathrm{\right)}$

Both parameters must be comprised between -2 and 0 for stability. The rescaling of lattice Boltzmann time units (tu) in real time units [T] is made through the relation:

${\mathit{T}}_{\mathit{R}}\mathrm{=}\frac{{\mathit{D}}_{\mathit{m}}^{\mathit{L}\mathit{B}\mathit{M}}{\mathit{L}}_{{\mathit{R}}^{\mathrm{2}}}}{{\mathit{D}}_{\mathit{m}}^{\mathit{R}}{\mathit{L}}_{\mathit{L}\mathit{B}\mathit{M}}^{\mathrm{2}}}{\mathit{T}}_{\mathit{L}\mathit{B}\mathit{M}}\mathrm{}\mathrm{\left(}\mathrm{8}\mathrm{\right)}$

where T_LBM and L_LBM are the space and time in lattice units (respectively, tu and lu) and T_R and L_R are their corresponding time and space units in real units [respectively (T) and (L)].

Model Parameterization

Typical parameters for the three different strains 3R, 7R, and 9R of Arthrobacter sp. (see Table 1) were taken from the literature. The specific uptake rate, ${\mathit{v}}_{\mathit{D}\mathit{O}\mathit{C}}^{\mathit{j}}$, the half saturation constant of the uptake rate, ${\mathit{k}}_{\mathit{D}\mathit{O}\mathit{C}}^{\mathit{j}}$, and the respiration rate, ${\mathit{k}}_{\mathit{r}}^{\mathit{j}}$, of the three species were directly taken from Monga et al. (2014), while the population mortality rate, provided by the same authors, was reinterpreted as probability, ${\mathit{P}}_{\mathit{m}}^{\mathit{j}}$, dependent on the duration of the time step. These strains have different growth patterns, with 3R expected to be the fastest growing strain and strain 9R the least competitive one due to its highest value of the half saturation constant of the uptake rate. Strain 7R is supposed to be representative of a more generalist strain.

TABLE 1

TABLE 1. Model parameters values used in the simulations. See text for an explanation of the computations and the bibliographic references used.

The maximum carrying capacity of a lattice-Boltzmann node, N_VOX, was calculated from the volume of a single image voxel (683 μm³) and the mean cell volume of a bacterial cell. The mean volume of an Arthrobacter cell was calculated to be 0.418 μmł according to data from Erlebach et al. (2000), assuming a spherical shape for the bacterial cell. We used this mean volume for the three strains of Arthrobacter sp. The mean reproduction diameter (i.e., the diameter attained by the cell before division) was also estimated from Arthrobacter cell size distributions (Erlebach et al., 2000) and assumed to be in the range 1.25 ± 0.15 μm, which includes the biggest diameters of the size distribution measured with a Coulter counter. The central value of 1.25 μm has been used in the present study.

Then, the mean carbon content of a single cell (${\mathit{m}}_{\mathit{t}\mathrm{0}}^{\mathit{j}}$, Table 1) and the value for the reproduction mass, ${\mathit{m}}_{\mathit{R}}^{\mathit{j}}$, were calculated from the mean cell volume and the mean reproduction volume calculated earlier, assuming a density of 1.1 g/cm³, a ratio of dry to wet cell weight of 0.25, and a carbon content of 0.47 g C per gram of dry cells (Gras et al., 2011).

The value of the decay rate of the POM agents, k_POM, was set to 0.25 day^-1 as reported by Iqbal et al. (2014) for the decomposition rate of maize (Zea mays) stem residues. These authors found “optimal decomposition conditions similar to those obtained with ground material” for fragments of POM of 0.02 cm length. Assuming a density of POM of 0.12 g cm^-3 (Iqbal et al., 2013) and a volume of POM residue of about 0.02 cm × 0.02 cm × 0.01 cm, we calculated an initial POM mass of carbon, POM₀, of 1.92 10^-4 mg C. We translated this fragment of POM into four fragments of parallelepiped shape (1 × 1 × 3 voxels), located at the solid/liquid interface (Figure 3). When the POM hydrolyses, the DOC produced is included in the neighboring fluid site.

FIGURE 3

FIGURE 3. 2D view of the POM initialization scheme used in the scenarios. In the figure, the POM (red), the DOC (yellow), the soil solid matrix (black), and the pore space (white) are depicted. In (a), POM disaggregated, four parallelepiped POM fragments of 1 × 1 × 3 image voxels, were assumed to be present and connected to the pore space through a single voxel. In (b), POM aggregated, a single 3D POM fragment of 4 × 4 × 3 was assumed to be present and connected to the pore pace through 4 of its voxels.

The molecular diffusion coefficient in lattice-Boltzmann units, ${\mathit{D}}_{\mathit{m}}^{\mathit{L}\mathit{B}\mathit{M}}$, is fixed to 0.5 lu²tu^-1 (Vogel et al., 2015, 2018). Since the DOC molecular diffusion coefficient, ${\mathit{D}}_{\mathit{m}}^{\mathit{R}}$, is 6.73 10^-6 cm²s^-1 (Weast et al., 1986) and 1 lu = 68 μm, the rescaling equation (8) gives 1 tu = 3.44 s.

Simulation Scenarios

The scenarios were designed to investigate the effect of the local micro-environments of bacteria on the strains abundance. For all scenarios we randomly placed 230 bacteria of each strain in the medium which resulted in 690 bacterial spots. In one particular scenario (S4) the three strains systematically co-exist in the same bacterial spot, so that only 230 bacterial spots were generated. The local micro-environments of these bacterial spots were modified by introducing heterogeneity in the placement of the resource (POM and DOC) in the medium. The water content level was also modified. Table 2 summarizes the characteristics of the different scenarios. Scenarios S1, S2, and S3 test the effect of the spatial distribution of organic matter on the global organic matter degradation and strain abundance. Scenario S4 tests the effect of the direct interspecific bacterial competition on strain abundance, whereas scenarios S5 and S6 are designed to test the effect of the water saturation level on organic matter degradation and strain abundance. Details of these scenarios are given below.

TABLE 2

TABLE 2. Simulation scenarios overview.

Scenario S1

The initial 690 bacteria were distributed among 690 randomly selected liquid voxels neighboring the soil solid matrix. Four POM fragments were randomly distributed in the medium (Figure 3a). Ten replicated simulations were performed with the position of POM fragments controlled by a random seed. The positions of the 690 bacterial spots were left unchanged for the ten replicated simulations. The water saturation level, S_w, was 0.50.

Scenario S2

The four fragments of POM were gathered in one fragment made of 2 × 2 × 3 voxels with the base 2 × 2 voxels being contiguous solid sites neighboring 2 × 2 fluid sites (Figure 3b). Bacteria were located in the same 690 positions as for scenario S1 and the same water saturation level, S_w = 0.50, was adopted. Ten replicated simulations were again performed with the position of POM fragments controlled by a random seed Comparison of scenario 1 and 2 made it possible to assess the effect of the degree of spatial heterogeneity of POM.

Scenario S3

We assumed that all the carbon that can, potentially, be hydrolyzed and released to the liquid phase of the soil in the scenarios S1 and S2, was already homogeneously distributed in the liquid phase of the medium at the beginning of the simulations as a DOC. A water saturation index of S_w = 0.50 was also adopted. In ten replicated simulations, a random seed controlled the position of 690 bacterial spots. This resulted in 9 extra configurations of the 690 bacterial positions tested in S1 and S2. Comparison of S3 with S2 and S1 allowed us to evaluate the impact of homogenously- vs. heterogeneously distributed C within the soil.

Scenario S4

This scenario was aimed at the effect of direct interspecific competition on strain abundance. Two simulations (identified as S4a and S4b) involved 3 bacteria, one of each strain, in selected 230 bacterial spots. In S4a, we took one particular repetition of scenario S2 in which the three strains had a noticeable growth (repetition S2r2) and we classified the 690 spots according to the amount of biomass growth of the colony in a descending order. The first 230 bacterial spots were selected and used to place the initial bacteria. In S4b, we took the repetition S2r2 but we just selected the 230 bacterial positions occupied by the less competitive strain (9R). In each of these spots, we placed initially 3 bacteria, one of each strain.

Scenarios S5 and S6

These scenarios are identical to scenario S2 except for the water saturation level that was fixed to 1.00 (S5) and 0.25 (S6). Because the positions of the 690 bacterial spots are left unchanged, in scenario S6 there were 165 bacterial spots that were found to be in the gas phase, and thus did not grow. In the 525 spots still placed in water filled grid cells, the three strains were found to be equally distributed as 178 cells of strain 3R, 172 cells of strain 9R and 175 cells of strain 7R.

Results and Discussion

Effect of Spatial Distribution of Organic Matter on Global Organic Matter Degradation and Strain Abundance

Simulation scenarios S1 and S2 assume that the soil organic matter is found in a number of POM fragments that are distributed (S1) or aggregated (S2). In the proposed scenarios, the hydrolysis rate of POM is set to a constant value, so that for k_POM = 0.25 d^-1, about 92% of the initial mass of carbon of the POM fragments is hydrolyzed at the end of the simulations. The predicted time evolutions of the POM, DOC, total biomass (B), and CO₂ are shown in Figure 4. At the beginning of the simulations, the carbon hydrolyzed from the POM particles cannot be totally consumed and accumulates in the liquid phase. After one day, when exponential growth of bacteria starts, the major part of the dissolved organic carbon is metabolized by the cells. At the end of the simulations, the CO₂ emitted has not yet reached a plateau although the slope of its cumulated increase begins to decrease. Except for one repetition of scenario S2, no differences are observed in the time evolutions of the POM, DOC, total biomass, and CO₂ between the two POM distributions.

FIGURE 4

FIGURE 4. Global biodegradation kinetics (left) and strains abundance (right) of scenarios with dispersed POM fragments (S1) and aggregated POM fragments (S2). M/M₀ is the ratio of the carbon mass (mg C) of each output over the total initial carbon mass M₀.

Although the global model outputs are similar for both spatial distributions of POM, differences in the strains abundance are observed (Figure 4, right). When the POM is fragmented into 4 pieces, the less competitive strain 9R cannot grow significantly (Figure 4, scenario S1). The fastest growing strain 3R experiences an exponential growth up to two days before starting to decline. The more generalist strain 7R has a much smoother exponential increase compared to 3R, and even surpasses the biomass of 3R after five days, before starting to decline after about the day seven. When the POM fragments are gathered into a single piece, 3R has the same overall dynamics, with comparable mean abundances but with higher dispersion between replicates. The coefficient of variation (CV) of the abundance peak is 0.03 for S1 and 0.10 for S2. The strain 7R has also the same dynamics as in scenario S1 but with a lower mean abundance peak that does not exceed 78% of the 7R mean abundance peak observed for the scenario S1. Again, the CV is highest in scenario S2 compared to S1 with the values of 0.24 and 0.10, respectively. On the contrary, the strain 9R shows a much higher growth variation among replicates when the POM is aggregated in a single piece of POM. In particular, two replicates (S2r2 and S2r5) show growth kinetics similar to those of strain 7R. In three other replicates, 9R grows without exceeding the abundance of the other species. In the last five replicates, 9R presents a similar growth as in the case of POM fragmented. The S2 replicate with lower global total biomass and emitted CO₂ (S2r9) corresponds to the lowest abundances of the strains 3R and 9R, and to the highest abundance of 7R. In that case, the highest abundance of 7R does not compensate the low abundances of strains 3R and 9R, while compensation is observed for the other replicates at the scale of the whole soil sample.

A very different biodegradation kinetics is observed when the DOC is homogeneously distributed in the pore space (Figure 5, scenario S3). The global growth of the total biomass is much faster (Figure 5, left) with a maximum peak reached at day one. The value of the peak is 4.5 times higher than the peaks reached by the total biomass in scenarios with POM fragments. The higher amount of DOC available initially permits a higher bacterial uptake that translates into a faster growth of the bacteria. As a result, the amount of DOC quickly decreases within the same time interval (one day) and the cumulated CO₂ emitted is higher, with values of about 72% (Figure 4) and 91% (Figure 5), respectively. The different picture shown by the biodegradation kinetics has relevant implications for the species abundance (Figure 5, right). The most competitive 3R strain dominates from the very beginning. The high value of its maximal uptake rate (Table 1) makes 3R cells benefit more from the larger initial DOC concentration, and experience a quick and large exponential growth. This strain shows a more pronounced decline phase than in scenarios with POM fragments. Only when DOC becomes scarcer, the more generalist strain 7R with the lowest k_DOC values succeeds to grow at about the same extent as observed in previous scenarios with POM fragments. The strain 9R with intermediate growth rate but with a higher k_DOC value is not competitive enough to grow. A probable explanation is that, after two days, the DOC concentration at the local microbial habitat of strain 9R remains always lower than in the case of the POM fragments. In the latter case, some microbial habitats of strain 9R that are located close to the POM fragments can still benefit from a sufficient DOC concentration for them to grow. On the contrary, the more generalist strain 7R is not impacted by the lower DOC concentrations. Figure 5 also shows that the position of the bacterial spots does not affect the predicted model outputs, and the ten replicated simulations almost overlap. Simulation scenario 3 is consistent with the nutritional state of the soil after a sudden flush of nutrients, which is typical of the anthropic addition of fertilizers to agricultural land, or as observed after rainfall following a period of drought.

FIGURE 5

FIGURE 5. Global biodegradation kinetics (left) and strains abundance (right) of the scenario with the C available as DOC and a water saturation level of 0.5 (S3). M/M₀ is the ratio of the carbon mass (mg C) of each output over the total initial carbon mass M₀.

Effect of Spatial Distribution of Organic Matter on Spatial Distribution of Strain Abundance

The spatial distribution of the abundances of the different strains in the scenarios S1 and S2 reveals a number of interesting patterns. In terms of the biomass in each of the 690 spots, we observe that the maximum peak abundance reached in the spots is 2.43 ± 1.01 higher in the case of scenario S2 compared to scenario S1 (Figure 6). In both scenarios, a few spots containing cells of the less competitive strain 9R can surpass spots of strain 3R and 7R and even be the most active spots (simulations S2r5 Figure 6, and simulation S2r2, not shown). When POM is fragmented, some spots containing 9R cells also end up with an amount of biomass that is similar to what is found with the strains 3R and 7R (simulations S1r2 and, to a lesser extent, simulation S1r9, Figure 6) showing that the global, per strain representation displayed in the graphic at right in Figure 4 hides the very large dispersion of kinetics at the local scale (that of the microbial habitat).

FIGURE 6

FIGURE 6. Biomass growth kinetics of the three strains (3R, 9R, and 7R) of each of the 690 spots in scenarios of dispersed POM fragments (S1, Top figures) and aggregated POM fragments (S2, Bottom figures). M/M₀ is the ratio of the carbon mass (mg C) of the cells over the total initial carbon mass M₀. For each scenario, the replicate having the minimal, mean and maximum biomass growth is displayed from left to right.

A closer look at the microbial habitats suggests that 10 ± 1% of the spots do not experience any bacterial growth in scenarios S1. Among those spots, about 45% contain one initial bacterial cell of strain 3R, while 21% contain one initial 9R cell, and 34% contain one initial 7R cell. A higher proportion of non-active spots is found in scenario S2, 18.5 ± 5%. A similar distribution of strains among those spots is observed for the scenario S1. When it comes to the active spots, when POM is fragmented in four pieces, only at 9.5 ± 4.0% of the spots does biomass exceed 10% of the maximum biomass registered among all the spots. Among these spots, about 60% contain 3R cells, 38% contain 7R cells and only 2% contain 9R cells. When POM is present in a single piece, the proportion of active spots exceeding 10% of the maximum simulated biomass, drops to 4.2 ± 2.4%. Among those spots, about 58% contain 3R cells, 33% contain 7R cells, and 8% contain 9R cells. The higher proportion of 9R cells spots among the most active spots in scenario S2 compared to scenario S1 explains the observed higher abundance of this strain.

The general trend suggested by these results is that when POM is gathered into one piece, significantly fewer microbial habitats (almost half) are prone to grow, and a larger dispersion of the abundance is found compared to scenarios in which POM is fragmented. Therefore, under these conditions, it is not surprising that similar biomass growth is observed at the scale of a soil sample (Figure 4, left). In our scenarios, since the POM hydrolysis is constant and independent of the spatial position of bacteria, a lower amount of dissolved organic carbon is produced locally compared to the aggregated POM fragments when there are 4 fragments. Even if the probability of having more spots closer to these local sources of DOC is higher with dispersed POM, the available DOC concentration remains lower. Gaillard et al. (1999) show that an aggregated distribution of organic matter mineralized a lower amount of carbon than a dispersed distribution. They explain their results on the basis of a higher exchange surface with soil of the dispersed distribution. The simulation scenarios performed in this contribution were designated to create a similar contact surface area between the soil solid and liquid phases and, therefore, cannot account for the reported outputs. Nevertheless, in general, differences in the contact area still constitute a plausible explanation for differences in mineralization rate in soils.

Using an algorithm developed by Dijkstra (1971) and based on the 6-connexity of the lattice-Boltzmann grid, we further calculated the geodesic distance between each of the 690 spots and the POM fragments in order to relate biomass growth of the microbial habitats to their spatial remoteness of POM. The geodesic distance is the shortest pathway included in the liquid phase that connects two points in the pore space (Gommes et al., 2009). Divided by the Euclidian distance, the direct pathway between the two points, it gives the geometrical tortuosity as defined by Clenell (1997). Geometrical tortuosity and constriction are often considered to be good descriptors of the diffusive transport of solutes in complex pore spaces (Berg and Held, 2016). The mean geodesic distances of the bacterial spots to the POM fragments are about 1.6 times longer when POM fragments are gathered with 5,041 ± 3,190 μm and 8,144 ± 3,983 μm for scenarios S1 and S2, respectively. The mean smallest geodesic distances are 252 ± 111 μm and 490 ± 262 μm for S1 and S2, respectively. Within one replicate of either scenario, the most active spots are those having the shortest geodesic distance (Figure 7), however, when comparing the repetitions within a scenario it appears that spots where growth is maximum (S1r2) has a longer or similar minimal geodesic distance than spots having lower growth (S1r9 and S1r4). A large number of spots that do not have significant growth can have a very large range of geodesic distance values (Figure 7), which is an unexpected result. In some cases, the geodesic distance can even be close to the shortest values. However, when geodesic distances are longer than about 5000 μm in the case of S1 and about 7500 μm in the case of S2, the spots do not noticeably grow. When POM fragments are dispersed, the geodesic distances are shorter but also the DOC concentration that is emitted from the POM spots is lower.

FIGURE 7

FIGURE 7. Total biomass growth of each of the 690 spots against the geodesic distances between the spots and the POM fragments of scenarios S1 (Top figures) and S2 (Bottom figures). M/M₀ is the ratio of the carbon mass (mg C) of the cells over the total initial carbon mass M₀. For each scenario, the replicate having the minimal, mean and maximum biomass growth is displayed from left to right.

To our knowledge, it is the first time that geodesic distances between the nutrient resource and microbial habitats are calculated in 3D modeling scenarios of soil carbon dynamics. They reveal that beyond a distance of 5 mm to the POM fragments, the microbial colonies cannot grow. Interestingly, Euclidian distances of around 4 mm from straw labeled with ¹³C have been reported to hold the sites of higher microbial assimilation and referred as “residusphere” (Gaillard et al., 1999). However, although the most active microbial spots are correlated with the lowest geodesic distances, a low geodesic distance is not a sufficient condition for the microbial colony to grow. We suggest that the size of the pores also matters (e.g., discussion in Baveye et al., 2018). Large cavities can dilute the concentration of DOC that reaches the bacteria, impacting ultimately the growth in the microbial habitat. Therefore, a bacterial spot can experience a microscale environment promoting more the bacterial growth than a spot placed at a shorter geodesic distance to the POM. Calculation of constriction in addition to diffusion length would be more appropriate. More effort is thus needed to calculate other metrics of importance for diffusive transport such as the constriction factor and the diffusion length in order to characterize and fully understand species abundance and functioning at pore scale.

Effect of Direct Interspecific Competition on Strain Abundance

Simulation scenarios S1, S2, and S3 suggest that when the strain 9R is closely competing with 3R or 7R, growth of 9R can occur only when it is located in an advantageous point with respect to a sufficient source of DOC (for instance in the case of gathered POM fragments, Figure 4, scenario S2). To test this hypothesis, we chose one replicate of the scenario S2 (S2r2) in which a few 9R microbial habitats experienced growth, and we initially placed three cells in now 230 spots, one of each strains. The global outputs of the model are similar (Figure 8, left), but the strain abundances are very different (Figure 8, right). When the three strains co-exist in the same microbial habitat, the strain 9R cannot grow, even when it is located in the spots close to the POM fragments. The most competitive strain 3R grows to a much higher extent and the strain 7R has a delayed growth. There are no prominent differences observed between scenario S4a and S4b.

FIGURE 8

FIGURE 8. Global biodegradation kinetics (left) and strains abundance (right) of replicate 2 of the scenario with aggregated POM fragments (S4) in which the three strains are initially alone in the spots (S2r2, solid lines) or gathered in the spots (S4a dotted lines and S4b dashed lines). M/M₀ is the ratio of the carbon mass (mg C) of each output over the total initial carbon mass M₀. Simulation outputs obtained in the second repetition of the scenario S2 (solid lines) are compared to the two simulations of the scenario S4 (dashed lines).

Effect of Water Saturation Level on Organic Matter Degradation and Strain Abundance

The role of the water saturation level of the pore space on the carbon dynamics was investigated in the scenarios S5 and S6. The global model outputs are similar when the pore space is fully saturated with water (S_w = 1.00, Figure 9, scenario S5). As mentioned in Section Entities, State Variables, and Scales, oxygen limitations are not yet considered in Ib-LBioS-Comp, so that the different water saturation levels do not impact the role of oxygen in the bacterial activity in these scenarios. Furthermore matric potentials considered in the scenarios S2 and S6 are very close (about -0.3 and -0.6 kPa, respectively) so that oxygen limitations are not expected to happen under the tested conditions. Since the local positions of the aggregated POM fragments and the 690 microbial habitats are not changed, the only effect of the higher water saturation is to decrease the DOC concentrations. However, the diffusive transport of DOC is also accelerated because all the pore space is now connected. The strain abundances are rather similar although this environment is comparatively more favorable for the generalist strain 7R, which becomes the most abundant strain after day 6. The establishment on the system for the simulations showing a noticeable growth of the less competing strain appear mostly driven by the spatial heterogeneity of POM distribution.

FIGURE 9

FIGURE 9. Global biodegradation kinetics (left) and strains abundance (right) of the scenario with aggregated POM fragments and a water saturation level of 1.0 (S5). M/M₀ is the ratio of the carbon mass (mg C) of each output over the total initial carbon mass M₀.

When the water saturation level decreases, the pore space filled with air increases and more disconnected aqueous regions appear. Consequently, some POM spots may not hydrolyze resource to the liquid phase, and, some bacterial microhabitats may not have access to the DOC hydrolyzed by the connected POM agents. In particular, when the water saturation is divided by two (from S_w = 0.50 to S_w = 0.25), 525 of the original bacterial spots are still in the aqueous phase. In spite of that, five of the replicates (S6r2, S6r4, S6r5, S6r6, and S6r10) produced global model outputs similar to the ones observed for higher water saturation levels (Figure 10, upper panel). In one replicate (S6r1), the aggregated POM fragments are located in solid voxels whose neighbors are disconnected from the aqueous phase, preventing the release of DOC, and thus bacterial growth (Figure 10, lower panel). In three other replicates (S6r3, S6r7, and S6r9) one of the four gathered solid voxels containing the aggregated POM fragments also has a dry neighbor voxel. Furthermore, for two of them (S6r7 and S6r9), and for the replicate S6r8, no bacterial growth is recorded, resulting in the accumulation of DOC in the liquid phase (Figure 10, bottom left). In these three last situations, spatial disconnections between the microbial spots and the POM fragments are made possible by the aqueous phase fragmentation. Spatial disconnections are also observed in the repetition S6r3 where total biomass growth is very small (Figure 10, bottom). The strain abundances in the active spots are very similar to those observed when the water saturation is 0.50. In one repetition, the growth of strain 9R of all microbial spots surpasses the growth of strain 7R (Figure 10, upper right).

FIGURE 10

FIGURE 10. Global biodegradation kinetics (left) and strains abundance (right) of the scenario with aggregated POM fragments and a water saturation level of 0.25 (S6). M/M₀ is the ratio of the carbon mass (mg C) of each output over the total initial carbon mass M₀. The repetition number generating the trends depicted are detailed in the figure.

It is known that lack of competition resulting from spatial disconnection of soil microhabitats may promote biodiversity (Kim et al., 2008; Vos et al., 2013). Complete or partial spatial separation due to disconnection of liquid soil volumes in unsaturated soils is a common hypothesis used to explain biodiversity. For instance, one study using experimental setups with two bacterial strains competing for a dissolved resource in sand showed dominance of the more competitive strain under water-saturated conditions while drier conditions allowed the less competitive bacteria to establish (Treves et al., 2003). Regardless of the high variability on the model outputs found for the drier conditions (S_w = 0.25), our scenarios do not show a clearly improved establishment of the less competing strain. Another study (Zhou et al., 2002), based on an rRNA-based cloning approach, reported differing biodiversity distributions in the microbial communities living in four geographically distinct sites at different soil depths. A uniform biodiversity distribution, which is thought to arise from a lack of microbial competition, was obtained for the saturated subsurface of both high and low carbon soils. Since the hypothesis of resource disconnection is difficult to hold under water-saturated conditions, the authors explain their observed pattern by a lack of competition due to specialization for different substrates. Although this hypothesis remains plausible, our simulations suggest that the spatial heterogeneity of the resource placement could also explain part of this biodiversity.

Conclusion

We have coupled a multi-species individual-based model describing bacterial growth to a 3D lattice-Boltzmann diffusion model to simulate organic matter dynamics in soil pore space. The resulting model, Ib-LBioS-Comp, has been used to study the influence of the spatial heterogeneity of a nutrient source on the organic matter degradation and species abundance of a competitive-, a generalist-, and a poorly-competitive bacterial strain. The scenarios used three resource placements showing a gradient in the spatial heterogeneity of its distribution: organic matter dissolved in the aqueous phase (low heterogeneity), four fragments of organic matter (intermediate heterogeneity), and a single fragment of organic matter (high heterogeneity).

A number of results can be highlighted from the modeling scenarios performed: (i) In general terms, the greater the spatial heterogeneity of the location of the resource, the greater the variability in the output at the level of the soil volume imaged; (ii) When the resource is found as particulate organic matter, the fastest growing strain tends to dominate at first but then, when the resources becomes scarcer, it is overtaken by the generalist strain, showing that the spatial distribution of organic matter affects bacterial succession; (iii) the global bacterial growth is faster when the nutrient resource is available in the liquid phase at the beginning of the simulation. Under these circumstances, the fastest growing strain is able to reach much higher relative abundances, having a negative effect on biodiversity; (iv) When the resource is present as particulate organic matter, the total biomass created does not differ noticeably between the intermediate and high spatial heterogeneity schemes but in contrast the species abundance is impacted; (v) The least competing strain, which does not reach noticeable growth for the low and intermediate resource spatial heterogeneity schemes, is able to grow appreciably in the absence of direct competition, if the position of the nutrient resource is favorable. According to this observation, heterogeneity of the spatial distribution of the organic matter in soil would promote microbial diversity; (vi) the geodesic distance among the nutrient resources and the bacteria alone is not sufficient to explain this phenomenon; and, (vii) In the scenarios tested, the water saturation level does not seem to change much the observed biodiversity.

As a cautionary note, one needs to remember that the various predictions made in this article emanate from a relatively simple model, which ignores many different aspects of soils. In actual soils, microbial diversity and growth are potentially affected by a myriad of factors (e.g., resource quality, different nutritional requirements of individual bacteria, mutation, predation, transport by soil fauna, toxin production). In particular, bacterial motility and the diffusion of oxygen in the pore space are likely to have a significant influence on microbial activity. Subresolution pores, too small to be visible in X-ray CT images, are most probably playing a role in this context as well, especially in the dry range of the hydrological regime of soils. Future research will determine if the relationship between bacterial diversity and the spatial heterogeneity in resource distribution, highlighted in the present article, still holds when some of the assumptions of the model are lifted. It is possible that the same relationship will be observed, or that other, more complex behaviors will unfold.

Data Availability

Data underlying this paper can be accessed at https://doi.org/10.17862/cranfield.rd.6744158.

Author Contributions

XP implemented the models. XP and VP ran simulations and analyzed the data. All authors conceived, designed the study, and contributed to the writing of the manuscript.

Funding

The research reported in this article was made possible by the financial support of the ANR Project Soil-μ3D (Number ANR-15-CE01-0006-01), and the Natural Environment Research Council (NE/P014208/1).

Conflict of Interest Statement

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.

References

Banitz, T., Gras, A., and Ginovart, M. (2015). Individual-based modeling of soil organic matter in NetLogo: transparent, user-friendly, and open. Environ. Model. Softw. 71, 39–45. doi: 10.1016/j.envsoft.2015.05.007

CrossRef Full Text | Google Scholar

Baveye, P. C., Baveye, J., and Gowdy, J. (2016a). Soil “ecosystem” services and natural capital: critical appraisal of research on uncertain ground. Front. Environ. Sci. 4:41. doi: 10.3389/fenvs.2016.00041

CrossRef Full Text | Google Scholar

Baveye, P. C., Berthelin, J., and Munch, J.-C. (2016b). Too much or not enough: reflection on two contrasting perspectives on soil biodiversity. Soil Biol. Biochem. 103, 320–326. doi: 10.1016/j.soilbio.2016.09.008

CrossRef Full Text | Google Scholar

Baveye, P. C., Pot, V., and Garnier, P. (2017). Accounting for sub-resolution pores in models of water and solute transport in soils based on computed tomography images: are we there yet? J. Hydrol. 555, 253–256. doi: 10.1016/j.jhydrol.2017.10.021

CrossRef Full Text | Google Scholar

Baveye, P. C., Otten, W., Kravchenko, A., Balseiro Romero, M., Beckers,É., Chalhoub, M., et al. (2018). Emergent properties of microbial activity in heterogeneous soil microenvironments: different research approaches are slowly converging, yet major challenges remain. Front. Microbiol. 8:1364. doi: 10.3389/fmicb.2017.01364

CrossRef Full Text

Berg, C. F., and Held, R. (2016). Fundamental transport property relations in porous media incorporating detailed pore structure description. Transp. Porous Media 112, 467–487. doi: 10.1007/s11242-016-0661-7

CrossRef Full Text | Google Scholar

Burns, R. G. (1982). Enzyme activity in soil: location and a possible role in microbial ecology. Soil Biol. Biochem. 14, 423–427. doi: 10.1016/0038-0717(82)90099-2

CrossRef Full Text | Google Scholar

Clenell, M. B. (1997). Tortuosity: a guide through the maze. Geol. Soc. Lond. Spec. Publ. 122, 299–344. doi: 10.1144/GSL.SP.1997.122.01.18

CrossRef Full Text | Google Scholar

Dijkstra, E. W. (1971). Ewd316: A Short Introduction to the Art of Programming. Eindhoven: T.U. Eindhoven.

Google Scholar

Elyeznasni, N., Sellami, F., Pot, V., Benoit, P., Vieublé-Gonod, L., Young, I., et al. (2012). Exploration of soil micromorphology to identify coarse-sized OM assemblages in X-ray CT images of undisturbed cultivated soil sores. Geoderma 179–180, 38–45. doi: 10.1016/j.geoderma.2012.02.023

CrossRef Full Text

Erlebach, C. E., Illmer, P., and Schinner, F. (2000). Changes of cell size distribution during the batch culture of Arthrobacter strain PI/1-95. Antonie Van Leeuwenhoek 77, 329–335. doi: 10.1023/A:1002632128921

PubMed Abstract | CrossRef Full Text | Google Scholar

Folse, H. J. III., and Allison, S. D. (2012). Cooperation, competition, and coalitions in enzyme-producing microbes: social evolution and nutrient depolymerization rates. Front. Microbiol. 3:338. doi: 10.3389/fmicb02012.00338

PubMed Abstract | CrossRef Full Text | Google Scholar

Gaillard, V., Chenu, C., Recous, S., and Richard, G. (1999). Carbon, nitrogen and microbial gradients induced by plant residues decomposing in soil. Eur. J. Soil Sci. 50, 567–578. doi: 10.1046/j.1365-2389.1999.00266.x

CrossRef Full Text | Google Scholar

Genty, A., and Pot, V. (2013). Numerical simulation of 3D liquid-gas distribution in porous media by a two-phase TRT lattice boltzmann method. Transp. Porous Media 96, 271–294. doi: 10.1007/s11242-012-0087-89

CrossRef Full Text | Google Scholar

Ginzburg, I. (2005). Equilibrium-type and link-type lattice Boltzmann models for generic advection and anisotropic-dispersion equation. Adv. Water Resour. 28, 1171–1195. doi: 10.1016/j.advwatres.2005.03.004

CrossRef Full Text | Google Scholar

Gommes, C. J., Bons, H. J., Blacher, S., Dunsmuir, J. H., and Tsou, A. H. (2009). Practical methods for measuring the tortuosity of porous materials from binary or gray-tone tomographic reconstructions. AIChE J. 55, 2000–2012. doi: 10.1002/aic11812

CrossRef Full Text | Google Scholar

Gras, A., Ginovart, M., Valls, J., and Baveye, P. C. (2011). Individual-based modelling of carbon and nitrogen dynamics in soils: parameterization and sensitivity analysis of microbial components. Ecol. Modell. 222, 1998–2010. doi: 10.1016/j.ecolmodel.2011.03.009

CrossRef Full Text | Google Scholar

Grimm, V., Berger, U., Bastiansen, F., Eliassen, S., Ginot, V., Giske, J., et al. (2006). A standard protocol for describing individual-based and agent-based models. Ecol. Modell. 198, 115–126. doi: 10.1016/j.ecolmodel.2006.04.023

CrossRef Full Text | Google Scholar

Grimm, V., Berger, U., DeAngelis, D. L., Polhill, J. G., Giske, J., and Railsback, S. F. (2010). The ODD protocol: a review and first update. Ecol. Modell. 221, 2760–2768. doi: 10.1016/j.ecolmodel.2010.08.019

PubMed Abstract | CrossRef Full Text | Google Scholar

Hanson, C. A., Furhman, J. A., Horner-Devine, M. C., and Martiny, J. B. (2012). Beyond biogeographic patterns: processes shaping the microbial landscape. Nat. Rev. Microbiol. 10, 497–506. doi: 10.1038/nrmicro2795

PubMed Abstract | CrossRef Full Text | Google Scholar

Iqbal, A., Beaugrand, J., Garnier, P., and Recous, S. (2013). Tissue density determines the water storage characteristics of crop residues. Plant Soil 367, 285–299. doi: 10.1007/s11104-012-1460-8

CrossRef Full Text | Google Scholar

Iqbal, A., Garnier, P., Lashermes, G., and Recous, S. (2014). A new equation to simulate the contact between soil and maize residues of different sizes during their decomposition. Biol. Fertil. Soils 50, 645–655, doi: 10.1007/s00374-013-0876-875

CrossRef Full Text | Google Scholar

Kim, H. J., Boedicker, J. Q., Choi, J. W., and Ismagilov, R. F. (2008). Defined spatial structure stabilizes a synthetic multispecies bacterial community. Proc. Natl. Acad. Sci. U.S.A. 105, 18188–18193. doi: 10.1073/pnas.0807935105

PubMed Abstract | CrossRef Full Text | Google Scholar

Kravchenko, A., Negassa, W., Guber, A. K., and Schmidt, S. (2014). New approach to measure soil particulate organic matter in intact samples using x-ray computed microtomography. Soil Sci. Soc. Am. J. 78, 1177–1185. doi: 10.2136/sssaj2014.01.0039

CrossRef Full Text | Google Scholar

Long, T., and Or, D. (2009). Dynamics of microbial growth and coexistence on variably saturated rough surfaces. Microb. Ecol. 58, 262–275. doi: 10.1007/s00248-009-9510-3

PubMed Abstract | CrossRef Full Text | Google Scholar

Monga, O., Garnier, P., Pot, V., Coucheney, E., Nunan, N., Otten, W., et al. (2014). Simulating microbial degradation of organic matter in a simple porous system using the 3-D diffusion-based model MOSAIC. Biogeosciences 11, 2201–2209. doi: 10.5194/bg-11-2201-2014

CrossRef Full Text | Google Scholar

Nannipieri, P., Ascher, J., Ceccherini, M. T., Landi, G., Pietramellara, G., and Renella, G. (2003). Microbial diversity and soil functions. Eur. J. Soil Sci. 54, 655–670. doi: 10.1046/j.1351-0754.2003.0556.x

CrossRef Full Text | Google Scholar

Peth, S., Chenu, C., Leblond, N., Mordhost, A., Garnier, P., Nunan, N., et al. (2014). Localization of soil organic matter in soil aggregates using synchrotron-based X-ray microtomography. Soil Biol. Biochem. 78, 189–194. doi: 10.1016/j.soilbio.2014.07.024

CrossRef Full Text | Google Scholar

Philippot, L., Spor, A., Hénault, C., Bru, D., Bizouard, F., Jones, C. M., et al. (2013). Loss in microbial diversity affects nitrogen cycling in soil. ISME J. 7, 1609–1619. doi: 10.1038/ismej.2013.34

PubMed Abstract | CrossRef Full Text | Google Scholar

Pot, V., Peth, S., Monga, O., Vogel, L. E., Genty, A., Garnier, P., et al. (2015). Three-dimensional distribution of water and air in soil pores: comparison of two-phase two-relaxation-times lattice-Boltzmann and morphological model outputs with synchrotron X-ray computed tomography data. Adv. Water Res. 84, 87–102. doi: 10.1016/j.advwatres.2015.08.006

CrossRef Full Text | Google Scholar

Powlson, D., Xu, J., and Brookes, P. (2017). “Throught the eye of the needle: the story of the soil microbial biomass,” in Microbial Biomass: A Paradigm Shift in Terrestrial Biogeochemistry, ed. K. R. Tate (London: World Scientific), 1–40.

Rainey, P. B., and Travisano, M. (1998). Adaptive radiation in a heterogeneous environment. Nature 394, 69–72. doi: 10.1038/27900

PubMed Abstract | CrossRef Full Text | Google Scholar

Roesch, L. F. W., Fulthorpe, R. R., Riva, A., Casella, G., Hadwin, A. K. M., Kent, A. D., et al. (2007). Pyrosequencing enumerates and contrasts soil microbial diversity. ISME J. 1, 283–290. doi: 10.1038/ismej.2007.53

PubMed Abstract | CrossRef Full Text | Google Scholar

Treves, D. S., Xia, B., Zhou, J., and Tiedje, J. M. (2003). A two-species test of the hypothesis that spatial isolation influences microbial diversity in soil. Microb. Ecol. 45, 20–28. doi: 10.1007/s00248-002-1044-x

PubMed Abstract | CrossRef Full Text | Google Scholar

Vieublé-Gonod, L., Chadoeuf, J., and Chenu, C. (2006). Spatial distribution of microbial 2, 4-dichlorophenoxyacetic acid (2,4-D) mineralization potential at a millimetre scale in soil. Soil Biol. Biochem. 35, 373–382. doi: 10.1016/S0038-0717(02)00287-0

CrossRef Full Text

Vogel, L., Makowski, D., Garnier, P., Vieublé, L., Coquet, Y., Raynaud, X., et al. (2015). Modeling the effect of soil meso- and macropores topology on the biodegradation of a soluble carbon substrate. Adv. Water Res. 83, 123–136. doi: 10.1016/j.advwatres.2015.05.020

CrossRef Full Text | Google Scholar

Vogel, L. E., Pot, V., Makowski, D., Garnier, P., and Baveye, P. C. (2018). To what extent do uncertainty and sensitivity analyses help unravel the influence of microscale physical and biological drivers in soil carbon dynamics models? Ecol. Modell. 383, 10–22. doi: 10.1016/j.ecolmodel.2018.05.007

CrossRef Full Text | Google Scholar

Vos, M., Wolf, A. B., Jennings, S. J., Kowalchuk, G. A. (2013). Micro-scale determinants of bacterial diversity in soil. FEMS Microbiol. Rev. 37, 936–954. doi: 10.1111/1574-6976.12023

PubMed Abstract | CrossRef Full Text | Google Scholar

Weast, R. C., Astle, M. J., Beyer, W. H. (eds). (1986). CRC Handbook of Chemistry and Physics, 66th Edn. Boca Raton, FL: CRC Press.

Google Scholar

Werts, S., Degrange, V., Prosser, J. I., Poly, F., Commeaux, C., Freitag, T., et al. (2006). Maintenance of soil functioning following erosion of microbial diversity. Environ. Microbiol. 8, 2162–2169. doi: 10.1111/j.1462-2920.2006.01098.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhou, J., Xia, B., Treves, D. S., Wu, L.-Y., Marsh, T. L., O’Neill, R. V., et al. (2002). Spatial and resource factors influencing high microbial diversity in soil. Appl. Environ. Microbiol. 68, 326–334. doi: 10.1128/AEM.68.1.326-334.2002

CrossRef Full Text | Google Scholar