Microscale Heterogeneity Explains Experimental Variability and Non-Linearity in Soil Organic Matter Mineralisation

Soil respiration represents the second largest CO2 flux from terrestrial ecosystems to the atmosphere, and a small rise could significantly contribute to further increase in atmospheric CO2. Unfortunately, the extent of this effect cannot be quantified reliably, and the outcomes of experiments designed to study soil respiration remain notoriously unpredictable. In this context, the mathematical simulations described in this article suggest that assumptions of linearity and presumed irrelevance of micro-scale heterogeneity, commonly made in quantitative models of microbial growth in subsurface environments and used in carbon stock models, do not appear warranted. Results indicate that microbial growth is non-linear and, at given average nutrient concentrations, strongly dependent on the microscale distribution of both nutrients and microbes. These observations have far-reaching consequences, in terms of both experiments and theory. They indicate that traditional, macroscopic soil measurements are inadequate to predict microbial responses, in particular to rising temperature conditions, and that an explicit account is required of microscale heterogeneity. Furthermore, models should evolve beyond traditional, but overly simplistic, assumptions of linearity of microbial responses to bulk nutrient concentrations. The development of a new generation of models along these lines, and in particular incorporating upscaled information about microscale processes, will undoubtedly be challenging, but appears to be key to understanding the extent to which soil carbon mineralization could further accelerate climate change.


Introduction
Soil respiration has been acknowledged as the second largest CO 2 flux from ecosystems to the atmosphere, and researchers have recognized that a small change in soil respiration could contribute significantly to a further rise in atmospheric CO 2 , increasing the threat of a 'runaway climate' scenario [1]. In spite of that, understanding the environmental controls of soil respiration, and especially of whether soil respiration is likely to be enhanced by rising soil temperatures and changes in the frequency of rainfalls, remains contradictory.
Convergent evidence suggests that this limited understanding, and the ensuing reality that the outcomes of experiments in this area remain largely unpredictable, stem from the fact that we are missing something essential, which occurs in soils at the spatial scale at which microorganisms operate, and that existing models are based on assumptions, which need to be alleviated. Experiments typically involve samples at the core-scale (5 cm height and diameter), characterized in terms of bulk properties, without taking into account heterogeneities that are known to exist at micro-scale with regard to the distribution of microorganisms, nutrients, or the geometry of the pore space in which both microorganisms and nutrients co-evolve [2][3][4][5][6][7]. Models of the fate of soil organic matter likewise ignore details of the heterogeneity of microbial populations and their physiology, overwhelmingly describing microbial processes by 1 st order kinetics relations, which often conflict with experimental results that exhibit high variability [1], a clear influence of microbial physiology [2], and a marked non-linear character [3]. For example, experiments with fungi have shown their mode of spread and growth switches from invasive to non-invasive, which can be interpreted as autocatalytic and conservative behaviours. However, further work has shown in spatially explicit soil, plant and nutrient environments that this behaviour is affected not just by fungal physiology but also by environmental cues including the density, quantity and heterogeneity of spatially distributed resources [4] or the pore network of a soil environment through which fungi spread [5]. We argue that accurate predictions can be obtained only when considering the physiology of microorganisms and their interplay with carbon distribution within a physically heterogeneous environment. Fortunately, tremendous technological advances, in terms of imaging, simulation [8][9][10][11] and experimental techniques, now allow unprecedented quantification of the microscale heterogeneity of the physical, chemical and biological characteristics of soils. This information can be used to develop a new generation of models to better describe the fate of soil organic matter, which in turn help predict the macroscopic outcomes associated with various scenarios of microscopic heterogeneity.
In this general context, a first objective of this article is to develop a spatially-explicit model of the fate of soil Particulate Organic Matter (POM), a fraction of the soil organic matter that typically has residence times ranging from 1 to 20 years in temperate areas [6]. The model relies on information about the 3-dimensional pore space of soils obtained via X-ray computed tomography. It accounts for physical and chemical processes at the micro-scale, as well as for the growth and metabolism of fungi, whose significance in POM decomposition is well acknowledged, especially in forest ecosystems [13][14][15]. Using this model, a second objective described in the following is to assess the extent to which the amount and spatial distribution of soil POM influences the fungal colony growth and spread, and the amount of CO 2 that evolves over time.

Materials and Methods
The soil samples on which this study was based were obtained from land belonging to French National Institute for Agricultural Research (INRA) and permission was granted. Integration of soil physical, chemical and biological components requires three components outlined in turn below: (1) characterisation of the pore volume; (2) location of POM at the pore/solid interface; (3) formulations of SOM-microbial feedback model.

Soil sampling and X-ray CT scanning
Soil was sampled from the surface 15 cm of a Eutric Cambisol (FAO; 17% clay, 56% silt and 27% sand) in the "La Cage" long-term trial at the INRA research centre of Versailles (France). Soil was sieved to aggregates <5 mm and packed at a bulk density of 1.4 g cm -3 . The soil samples were scanned with a Nikon Metrology X-Tek HMX CT scanner at energies of 100 kV 131 mA, with 3010 projections. A Varian Paxscan 2520 V detector, a 225 kV X-ray source (Nikon Metrology X-Tek Systems Ltd, Tring, UK) and a focal spot of 5 mm were used. A molybdenum target was used with a 0.25 mm thick aluminium filter to reduce beam hardening. NIKON software CT Pro v2.0 was used for reconstruction using a filtered back-projection algorithm. The voxel resolution was 50 microns. 2D cross sections (top view) from the 3D volumetric images of soil samples were saved as image sequences (slice thickness = voxel size) using VGStudio MAX v2.1 software. These images were cropped to a size of 200 × 200 × 200 voxels. The cropped images were converted to binary images using Indicator Kriging to separate the volume into a pore and a solid phase, and to produce the 3D structures within which POM was placed and simulations were conducted [7].

Placement of POM in the 3D soil structures
Since the true distribution of POM is unattainable using benchtop XRay μCT, a prescribed amount of Particulate Organic Matter (POM), defined as being larger than 50 μm, was distributed along the pore-solid interface, assuming homogenous or heterogeneous coverage controlled by magnitude and degree of aggregation (clustering).
A method was developed to produce random spatial distributions of POM within the 3D soil structure, adhering to global constraints, such as the POM (g) inputs per cm 3 of soil, and the particle size class distribution of POM. This permits detailed testing of microbial response to contrasting POM distributions that can be described by general characteristics obtained from soil analysis, such as quantity of POM and size fraction of POM. POM is distributed randomly (from a uniform random distribution) in the solid matrix at the pore-solid interface, according to the following set of simple rules and recognizing that POM is part of the solid phase. POM is generated from a central voxel randomly chosen in the set of solid voxels that have at least one porous neighbour. The POM aggregate is defined as the set of solid voxels located not further than a given distance from this centre. In the case of a flat surface POM aggregates would have a semi-spherical shape. Each voxel in the 3D volume can therefore have one of three states: (i) pore volume, (ii) inert solid phase, or (iii) POM. The number of particles of POM (N aggreg ) contained in the sample is set equal to the ratio between the volume of POM and the volume of a sphere that has the specified diameter d pom : where m pom is the specified amount of POM, and ρ pom is the mass density of POM. In the same manner, the number voxels of POM voxels (N pom_voxel ) in the sample required is given by: where h is the voxel edge length (50μm).

Scenario Modelling
Having developed a flexible way to place POM within a given soil architecture, one can then investigate the interplay between soil carbon dynamics and colony dynamics. The scenarios will investigate the effect of POM inputs and size class and distribution on fungal colonisation and activity. The amounts of POM and certain POM size class distributions are obtained from experimental data [8]. The POM amounts were varied between 0-5% sample weights (g/g) [9]. The allocation of POM into size classes has been described in numerous papers [10,11]. POM is defined as being greater than 50 microns in size and is typically grouped into 51-200 microns and 201-2000 micron classes. The same size classes were adopted in this study with the inclusion of a fine POM size being equal to 50 microns [10,11]. Following [12], Minkowski functionals were used to characterise the distribution of POM in space, including measures for POM volume fraction and % accessible POM (Table 1). POM volume fraction calculates the fraction of POM voxels with respect to total number of voxels in the sample. % accessible POM calculates the fraction of POM voxels that are located on a solid-pore interface-hence are accessible-, with respect to the total number of POM voxels.

Modelling Scenario 1-Assessing the effect of POM inputs
In the first scenario, we tested the impact of randomly distributed POM within the range of 0-5% C (g/g soil). We assume all POM to be of the smallest size (= 50 μm), therefore all POM is contained within a single voxel layer between the solid and pore phase (as the voxel size of the data was 50 μm). The conditions are given in scenario 1-7 in Table 1. Each simulation was repeated 5 times with different realisations in POM distributions, resulting from the random placement of POM on the pore/solid interface (except for when POM = 0 and 5% when there is no and complete coverage of solid/pore interface with POM).

Modelling Scenario 2-Assessing the effect of POM size distribution
To determine the effect of spatial distribution and POM size on evolved CO 2 , the quantity of POM was held constant while the POM sizes and their distributions varied. A POM content of 3% was selected as this was where a switch appeared to occur from non-invasive (limited spatial extent) to invasive growth (filled the accessible pore volume available). To represent POM  Table 1), (ii) a scenario with all POM of a large size (scenario 8 where POM diameter is 1 mm (20 voxels) ( Table 1)), (iii) and (iv) assume distributions of POM with scenario iii skewed towards the smaller size and with a distribution consistent with experimentally obtained POM size class distributions, and scenario (iv) with a similar distribution but now skewed to the larger POM sizes (Table 1).
In both scenarios the spatial distribution of POM is unknown, and therefore was randomly distributed within the solid phase as described in section Placement of POM.

Microbial SOM-Microbial Feedback Model
The key physiological processes of the fungal growth model are uptake, biomass recycling (inter-conversions amongst three fungal biomass types required for biomass recycling), and respiration associated with uptake of Dissolved Organic Carbon (DOC) (Fig 1A). Internal resource is translocated through the colony and is converted into biomass-hyphal tips (NIB) or an insulated biomass (IB), the latter possessing a reduced uptake capacity of DOC. NIB is converted into IB over time which is analogous to rigidification of hyphal structure [13]. NIB and IR spread via diffusion, being constrained by pore space and biomass respectively. Extensions to the model included the formulation of SOM-fungal feedback embedded in a soil context. The SOM-microbial feedback can be described as how the type, size and distribution of POM affect the partitioning of biomass between hyphal tips (NIB) or biomass for conservation (IB), which affects enzyme production and therefore the subsequent rate of carbon decomposition. The feedback is mediated via enzyme production which is assumed to be in proportion to the density of hyphal tips (NIB). The enzymes transform SOM (POM ! DOC), and the DOC can subsequently be utilised by hyphal tips (NIB) and converted into internal resources (IR). Owing to the large density of carbon in POM compared to fungal biomass, POM degradation is assumed to be driven by a high affinity process. The degradation rate is therefore proportional to the enzyme concentration. Enzymes are released from the hyphal tips (NIB) and are assumed to have a limited range of diffusion. POM degradation follows a first order kinetic process with respect to amount of NIB. As the simulation proceeds (simulation run equivalent to 3 weeks real time), POM amounts remain nearly constant, which is consistent with persistence time of POM in soils e.g. months for particulate organic matter such as maize [6]. Once DOC is produced it can diffuse freely in the porous media but it is also available for fungal uptake. Uptake is described by a Monod equation. The metabolic cost for NIB is represented by a loss of IR proportional to NIB amount that results in the production of CO 2 . Again the sustainability cost for IB is considered to be negligible.
The system of equations can be written as: DOC ðK doc þ DOCÞ À x þ NIB:ðb n :p À a n :p y Þþ IB:ðb i :p À a i :p y Þ @NIB=@t ¼ r:ðD nib :rNIBÞ À NIB:ðb n :p À a n :p y Þ À z:NIB @IB=@t ¼ ÀIB:ðb i :p À a i :p y Þ þ z:NIB Where α n , α i , β n , β i and θ are the recycling parameters, π stands for the internal resource concentration in the mycelium (π = IR/(IB+NIB)) and z is the insulation parameter (see [13]). D doc , D ir and D nib are the diffusion parameters respectively for DOC, IR and NIB which characterise the spread of the fungus through the pore volume as it finds new C resources. V pom , K pom, V doc and K doc are the Michaelis-Menton (MM) degradation and uptake parameters. ξ and δ u are metabolic costs. They are commonly referred to as maintenance and growth respiration, respectively. ϕ b is a function of ∑ B = NIB + IB + IR bounded between 0 and 1. It limits the total amount of fungal biomass in each single voxel to a maximal value B max being such as ϕ b = 0 when ∑ B = B max and ϕ b = 1 when ∑ B = 0. This function can be given by the empirical expression ϕ b = 1 -∑ B/B max .

Model parameterisation and calibration
The model was parameterised for Rhizoctinia solani Kühn, a soil pathogen and saprotrophic organism, and parameters were estimated via a literature review and sensitivity analysis, and were defined in terms of carbon transfer and flow [12]. Previous research has established that the fungal growth model, based on established theory and microbial parameters, is robust. Briefly, some parameters, governing specific physiological processes, were estimated via an extensive literature search and experimentation. Parameter values that could not be determined in this manner were identified, and the consequences of uncertainty in parameter values were assessed through sensitivity analysis. This analysis showed that predictions of biomass yield and extent are most sensitive to uncertainties in the recycling processes.
In the microbial-SOM model presented in this study, the microbial component now interacts with POM. The parameters associated with the breakdown of POM into DOC are unfortunately unknown. Values associated with POM inputs and size distributions can be obtained from the literature [8][9] or measured experimentally. The main parameter uncertainty is in the distribution of POM, as this cannot reliably be classified using existing segmentation methods applied to soils. The introduced parameters associated with decay of POM into a more accessible form (DOC) are presented in Table 2 together with their values. As an initial starting point we assumed the Michaelis-Menten (MM) parameters (V pom and K pom ) to be an order of magnitude smaller than V Doc and K Doc, and the latter values were determined and provided in [14].
The selected Michaelis-Menten (MM) parameter values (V POM and V DOC ) produced reasonable simulated outputs in terms of matching slow degradation of POM and cumulative respired CO 2 amounts, consistent with the literature [8][9]. The values for POM_QUANTITY, POM_DIAMETER, and POM size distributions were derived from the literature and experimentation as justified in "Scenario Modelling" section.
The spatial location of POM within the soil is unknown therefore it is distributed randomly adhering to the criteria stipulated in section "Placement of POM". Uncertainty in how POM is distributed can be assessed via system outputs (CO 2 respired) propagated from uncertain inputs. We will use the variability in model output (cumulative CO 2 respired) as a criterion to identify uncertainty and importance associated with the POM distributions.

Integration of soil physical, chemical and biological components
On the basis of the X-ray μCT data, the pore geometry visualised at a resolution of 50 μm is depicted in Fig 2. The visible porosity, which includes pores larger than 50 μm, was 10%, with a pore-solid surface area of 74.7 mm 2 . The pore volume was dominated by two large interconnected pore clusters and many smaller disconnected clusters, represented by different colors in Fig 2B. The large disconnected clusters imply that much of the pore volume can potentially become fully occupied by fungi if sufficient carbon is available for growth and spread.

CO 2 production relates in a non-linear way to C content in soil
Simulations assessed the impact of POM, distributed at the pore-solid interface, on fungal growth and evolved CO 2 . With increasing mass of POM there is an increase in the number of POM voxels on the pore/solid interface, ranging from zero to all interface voxels being  Table 3 POM spatial distribution Location of POM in soil sample Randomly distributed at pore-solid interface The main parameters of the microbial model are obtained from experimental data and literature research as described in [14].
doi:10.1371/journal.pone.0123774.t002  3A). Depending on the POM input, two distinct responses of fungal behaviours were observed with fungi either decreasing in biomass  Soil heterogeneity affecting fungal invasion and CO 2 production. a) A 2D schematic of the soil illustrating the distribution of POM (brown pixels) solid (black pixels) and pore (white pixels) phases. Enzymes produced by fungi transform POM into DOC (yellow) which is subsequently taken up by fungi and used for growth and CO 2 production b) The complex pore geometry of the soil sample used in this study is visualised using X-ray μCT. The sample shown has a pore volume of 10% (pores > 50μm) and a pore-solid surface area of 74.7 mm 2 . Different colors represent disconnected pore volumes, and the volume is dominated by two large disconnected connected clusters (gray and green). (conservative growth for low POM inputs, Fig 3B) or increasing (autocatalytic growth for high POM inputs, Fig 3C). Above a POM content of 3% g/g, the cumulative evolved CO 2 follows an exponential trend with time (resulting from autocatalytic growth of fungi). Below a POM content of 3% a levelling off of cumulative CO 2 is observed (Fig 3D and 3E).

Micro-scale distribution of POM causes substantial variability in respiration
The threshold value of POM that promoted a transition from conservative to autocatalytic growth was further investigated. This threshold was where 60% pore/solid interface voxels were POM, which corresponds to a POM content of 3% (g/g). The investigation involved the spatial allocation and clustering of POM by manipulating the POM size classes (Table 1). In particular the POM size classes were varied ranging from fine POM (1 voxel or 50 μm) to aggregated POM (20 voxels or 1000 μm). The realised POM distribution was characterised by geometric measures such as volume and % accessible surface area (see Table 1 in MM). The micro-scale heterogeneity of POM had a large impact on mean CO 2 production as well as the variability. For samples with an identical POM content of 3%, the cumulative mean CO 2 produced ranged from 0.034 to 2.38 mg C per cm 3 . This is approaching a striking 100 fold difference in CO 2 production from samples with identical bulk POM content and porosity resulting  from differences in the spatial distribution of POM at the solid-pore interface at the microscale. In addition, and in particular for samples where POM has relatively large particles (>200 μm) substantial variability is observed between replicated samples. Due to large variability in simulation results a total of 15 replicates were run. Interestingly, and perhaps counter intuitive, higher respiration is observed for a population containing larger POM size classes despite that these will have a smaller cumulative surface area (Fig 4) accessible to fungal colonisation.

Discussion
Central to the development of an understanding of SOM decomposition by fungi is an understanding of fungal ecology, influenced by the physical environment, type and distribution of substrates. Many carbon decomposition models, which do not have an explicit consideration of microbial ecology, do not observe a nonlinear response of respiration to carbon inputs, most likely due to the assumptions of the underlying models i.e. decomposition is affected by carbon pool only (1st order) [15] as opposed to the interplay amongst microbial physiology, the structured environment in which these organisms reside and C distributions. These results therefore show that the assumptions made in soil C models critically affect the predictions of CO 2 in relation to available C and that when underlying heterogeneity is considered a different relationship can emerge. Our second main finding which shows that substantial variability in POM turnover rate occurs at identical bulk properties is far reaching; it demonstrates that under certain conditions the bulk measurement (e.g. average porosity and average C content) of soils is insufficient to predict the production of CO 2 . This variability is attributed to microscale heterogeneity and microbial dynamics. Although this is only one component of the C cycle in soils it suggests that a different approach to larger scale soil C models based upon trends derived from models with explicit consideration of microbial dynamics and heterogeneity may be required to improve the current predictions of the roles of soils in climate change models. Finally, it highlights the importance of obtaining experimental data at appropriate spatial scales. Two modes of collective fungal behaviour, conservative and autocatalytic, both sensitive to POM inputs, are observed. This behaviour is determined by a combination of the spatial distribution of POM within the pore volume as this conditions the accessibility to organic matter, and of the intrinsic physiological processes that regulate fungal growth, spread and biomass partitioning within the colony. This affects the efficiency of POM decomposition within the pore volume. When conservation occurs, the little POM that is available is used for maintenance, leaving a lesser amount available for enzyme and biomass production. Biomass becomes insulated over time and spread is insufficient to access new POM; this further limits enzyme production, uptake of DOC and hence production of CO 2 . The stationary phase is realised fast in comparison with autocatalytic behaviour. For the autocatalytic behaviour the system is yet to be limited by carbon and fungi are in an exponential growth phase. This critical switch in collective behaviour of fungi is consistent with other experimental and theoretical work examining the effect of spatial distribution of C sources in a 2D environment where a small change in density [16] or connectivity of pore network [17] made fungi switch from non-invasive to invasive spread. Non-invasive spread is when the spatial extent of fungi is limited, whereas invasive spread occurs when the colony spreads and proliferates in the available pore space.
The large differences observed in POM turnover rate for identical bulk C inputs are a consequence of including the microbial physiology. Interestingly the larger POM pieces are metabolically more efficient to decompose than many separate smaller carbon hot spots voxels, an, at first surprising result. The separated smaller POM pieces do not invoke such a rapid decomposition of POM as enzyme production and DOC uptake is lower and the decomposition rate is hampered further by colony scale insulation (colony conservation). This colony-scale insulation acts to convert active (NIB) biomass into non-active form (IB) therefore reducing the colony's ability to produce enzymes. By altering the parameters that govern the physiological processes of colony growth, reflecting different, ecological strategies other results, i.e. more biomass and activity with smaller chunks of POM, may be obtained. This would likely be achieved by faster diffusion rates and less metabolic costs thus promoting more explorative colony behaviour as opposed to the observed exploitative behaviour. Here the fungal parameters used were calibrated for Rhizoctania solani Kühn based on literature values and parameter estimation [14].

Conclusions
Predictive modelling of soil carbon dynamics is a key factor in future climate change prediction. It has been suggested that uncertainties in predictions can be associated with a lack of representation of microbial dynamics at appropriate spatial scales in current soil carbon models. Our results show how biophysical interactions predict non-linear responses to carbon inputs in soil and explain the huge variability observed for samples with similar bulk characteristics. The results highlight the need for the development of a novel class of soil carbon models and a significant departure from existing SOM decomposition models to form a comprehensive understanding of soil-carbon responses to changing abiotic factors, as this response depends on efficiency of C usage by soil microbes. Further efforts should now focus on the effect of other ecological strategies of fungi and on considering microbial diversity, chemical composition and physical characteristics to identify the key drivers of C dynamics, and to upscale the observed patterns of behaviour to determine core-scale model parameters (CO 2 efflux, biomass yields) from measures of structure and carbon distribution. These patterns of behaviours can then be used in larger scale C dynamics models. Major advancements in soil science will be made if it is possible to translate micro-scale soil processes to macro-scale soil properties.