Stochastic processes in the structure and functioning of soil biodiversity

Lindo, Zoë; Bolger, Thomas; Caruso, Tancredi

doi:10.3389/fevo.2023.1055336

SYSTEMATIC REVIEW article

Front. Ecol. Evol., 09 May 2023

Sec. Population, Community, and Ecosystem Dynamics

Volume 11 - 2023 | https://doi.org/10.3389/fevo.2023.1055336

This article is part of the Research Topic Models in Population, Community and Ecosystem Dynamics View all 11 articles

Stochastic processes in the structure and functioning of soil biodiversity

$\r\nZoë Lindo$ Zoë Lindo¹^*

Thomas Bolger^2,3

Tancredi Caruso^2,3

¹Department of Biology, The University of Western Ontario, London, ON, Canada
²School of Biology and Environmental Science, University College Dublin, Dublin, Ireland
³Earth Institute, University College Dublin, Dublin, Ireland

Ecologists are increasingly recognizing the importance of stochastic processes in generating spatial and temporal variation in biological communities. This variation is very high in soil, which hosts not <¼ of all biodiversity on Earth and is central to how terrestrial ecosystems respond to perturbations. Measurement errors, demographic stochasticity (individual variability in traits such as birth and death rates), and environmental stochasticity (fluctuations in environmental properties) are the three main sources of stochasticity in ecology. Here, we synthesize how these three sources of stochasticity are quantified and incorporated in the study of soil biodiversity, highlighting current limits, possible solutions, and future research needs. We stress the relevance of all these factors to our future understanding of terrestrial ecosystems via plant-soil and soil-climate interactions and feedbacks. In soil, measurement errors are due to the small size, high abundance, and broad distributions of soil organisms, which limit sampling in space and especially over time. We argue that positive autocorrelation is a main characteristic of soil environmental properties, which may have important consequences on the response of soil biota to perturbations. At a local scale, large populations of soil organisms also imply a minor role of demographic stochasticity. Despite demographic stochasticity being a less significant source of variability than environmental stochasticity, we show that demographic stochasticity can be sizeable, but that within soil systems, stochasticity of environmental conditions must be accounted for. Explicit consideration of stochastic processes in soil biodiversity research is essential to our future understanding of the processes that control soil biodiversity. In classical ecology, stochasticity implies probabilistic predictions in terms of population growth, extinction, species coexistence, and community diversity. In soil, stochasticity implies very variable responses to climate change and the soil-climate feedback. Future studies will have to identify the major sources of environmental stochasticity with a particular focus on the interaction between multiple global change factors.

1. Introduction

Individuals, populations, and ecological communities display very large variability over space and time. Since a fully deterministic description of all the factors that control this variability is not feasible, ecologists have developed a large number of models that explicitly incorporate stochastic processes in modeling this variability (May, 1973; Caswell, 2000; Lande et al., 2003; Vellend, 2010). Shoemaker et al. (2020) have recently clarified many of the misconceptions and confusions on stochastic processes in ecology, highlighting the importance of an integrative framework to make the study of ecological communities more robust and predictive. In fact, a major misconception is that stochastic approaches imply uncertainty and unpredictability while, instead, they simply imply that predictions are probabilistic, and that there are different sources of variance that contribute to the outcome of an ecological process. Stochastic ecological processes can thus be quantified in terms of probability distributions and expectations (Lande et al., 2003). In other words, models that incorporate stochastic processes make “the unpredictable” predictable, if we accept that most predictions in ecology will be of a probabilistic nature. There is a misconception that a stochastic process implies total ignorance of a process, which is rarely the case. In most cases, there is information available to us that allows for predictions while calculating some degree of uncertainty around the mean. More generally, modeling a process using a stochastic approach does not mean that the process is intrinsically random and has no deterministic causes. Rather that the complexity of the factors that determine the trajectory of the measured variables is so high that a probability description of the possible trajectories is tractable while a deterministic one is not (Karlin and Taylor, 1975).

Soil is one of the most variable and diverse components in terrestrial ecosystem and hosts not <¼ of all biodiversity on earth (Bardgett and Van der Putten, 2014; FAO, 2020). The variability of soil biodiversity over space and time is enormous and observed from very small (a few mm) to global scales (e.g., Ettema and Wardle, 2002; Bardgett et al., 2005; Delgado-Baquerizo et al., 2018; Phillips et al., 2019; White et al., 2020; Caruso and Bardgett, 2021). Many studies have clarified the roles of the abiotic and biotic factors that control the distribution of soil organisms from global (Delgado-Baquerizo et al., 2018; Phillips et al., 2019; Van Den Hoogen et al., 2019) to intermediate and local scales (Lindo and Winchester, 2009; Caruso et al., 2019). But the abiotic and biotic factors that control soil organism distribution, and the ecosystem functions that these organisms mediate, vary over space and time, and stochastic processes offer the quantitative framework to link fluctuations in abiotic and biotic factors to variation in population and communities (Lande et al., 2003). This is particularly important for soil communities, that are central to ecosystem functioning, especially biogeochemical cycles (Bardgett and Wardle, 2010; Crowther et al., 2019), which are connected via plants and the atmosphere to global climate change dynamics (Bardgett et al., 2008), and are subjected to a complex interaction of multiple global change factors (e.g., warming, changes in soil moisture, nitrogen deposition) (Rillig et al., 2019; Bardgett and Caruso, 2020). Interactions with multiple factors imply important temporal and spatial fluctuations in populations, biomass, microbial and faunal traits, and biological rates, including energy fluxes. Given that a fully deterministic description of all these interactions is not feasible, we argue that stochastic approaches offer a solution to conceptualize and model this complexity. We thus offer a synthesis and perspective on how stochastic processes have been considered in past soil biodiversity research, highlighting some important aspects of the state-the-of-the-art, current limitations in terms of data availability, and solutions that will boost, in our view, a deeper understanding of the response of terrestrial ecosystems to environmental variation and perturbation regimes.

2. Sources of stochasticity in soil biodiversity

Over the last 15 years, the factors and processes that structure soil biodiversity have been studied intensively and over multiple spatial scales. For example there have been some general literature reviews and synthesis papers on different biodiversity theory and methods as applied to soil biota (e.g., Vályi et al., 2016; Thakur et al., 2020; White et al., 2020), large scale studies investigating the factors that control the distribution of soil organisms (e.g., Delgado-Baquerizo et al., 2018; Crowther et al., 2019; Phillips et al., 2019; Tedersoo et al., 2022), and studies that have tested the predictions of competing community assembly models (e.g., Lekberg et al., 2007, 2011; Lindo and Winchester, 2009; Dumbrell et al., 2010; Caruso et al., 2011, 2012). Much of this work was meta-analyzed by Guerra et al. (2020). In many of these studies, and also as recently summarized in Vályi et al. (2016), Thakur et al. (2020), and White et al. (2020), there is a relatively large proportion of the variance in the data that cannot be attributed to the factors analyzed. There is typically variance observed in the abundance of organisms, but also in relative species abundance and species composition, as well as in ecosystem functions (for example, total soil respiration). Much of this variance is typically spatial because most studies have focused on the spatial dimension, and time series of soil biota are rare. Caruso et al. (2020) recently reviewed the literature for soil animals and highlighted the rarity of time series for soil biota. For example, in the BioTIME database (Dornelas et al., 2018) the only time series for soil organisms are the two we recently submitted.

Depending on the specific context of the different research, various authors have offered multiple explanations on the sources underlying this “unexplained” variance, but in general there are three main sources in ecology (Figure 1): measurement errors, environmental, and demographic stochasticity (Shoemaker et al., 2020). The term “unexplained” is thus not really appropriate in our view, because the sources of variances are known, but what it is not known are their relative contributions.

FIGURE 1

Figure 1. Stochastic processes affect population fluctuations (black lines of different types), and so species intraspecific and interspecific interaction and eventually fluctuations in any properties emerging from these interactions (for example, the simple total density as the sum of individual species density, blue solid line). The two key types of stochastic processes are demographic stochasticity and environmental stochasticity. A third source of variance is measurement errors, which propagates when variable measured with uncertainty (species fluctuation) are combined into an aggregated properties (for example, total sum of species densities).

2.1. Measurement errors

The process of measurement in science is best modeled using probability distributions, given that measurements are subjected to a myriad of sources of error that together affect the accuracy and precision of measurements (Taylor, 1997). In this sense, measurement errors can be measured and accounted for (e.g., through instrument calibration, detection limits), such that, in principle we should always be able to quantify them. In ecology, at a theoretical level, measurement errors are perhaps the least interesting source of stochastic variation, but in soil research this type of error represents an important source of variance for two reasons: firstly, soil biota are very small sized and patchily distributed from local to broad scales (Ettema and Wardle, 2002), meaning the sampling regime should be very intense at multiple scales in order to provide precise estimates of ecological variables such as population density and biomass; secondly, basically all taxonomic groups in soil are sampled destructively, meaning that repeated measurements on the “same object” are not possible, and that measurements at one particular spatial scale must be based on either a composite sample from multiple samples or technical replicates later pooled or averaged to represent the community at the chosen spatial scale (Caruso and Bardgett, 2021). While this is well-known in the practice of soil ecology, we propose that there are also other sources of measurement errors that are underestimated in soil biodiversity research, and that soil ecologists and ecologists in general, rarely apply the rules of error propagation in their estimates of measurement errors. In ecology this topic has been considered explicitly only in a few works mostly related to the estimate of the biomass of plants (Deutschman et al., 1999; Lo, 2005; Molto et al., 2013).

We propose that there are two major areas where error propagation can become an important approach to estimating the contribution of measurement error to the total variance observed in soil ecological systems. One is the estimate of population biomass, which is usually based on compounding information on body size with information on population density (Turnbull et al., 2014). The other is the related estimate of energy fluxes (e.g., Potapov et al., 2016) that is based on the allometric scaling of metabolism with body size and temperature. This latter case is particularly important, given that an energetic description of soil food webs potentially allows the estimation of matter fluxes and so, for example, rates of CO₂ emissions (Moore and de Ruiter, 2012). In energetic food webs estimated through allometric metabolic scaling (Barnes et al., 2018), the total metabolic loss of a trophic species is the sum of the losses of each individual within that species population. The metabolic loss I of one individual of biomass M is typically estimated with an equation of the type: I∝M^b. There is a measurement error σ_M in M, which is often due to estimating individual biomass M by body size, with uncertainty both in length measurements and the parameters that link body size to body dry or wet weight. But, even more fundamentally, there is also uncertainty in the scaling exponent b. The exact value of the scaling exponent, and the implication of this value, have been a subject to heated discussion in the more general debate on the metabolic scaling theory of ecology (Enquist et al., 2003; Brown et al., 2004; Makarieva et al., 2006, 2008; Glazier, 2010). Regardless of the different theories and predictions on the exact value of the scaling exponent, experimentally there is a well-known systematic variation in the scaling exponent, which varies between taxa and phylogenetic lineages (Ehnes et al., 2011). For example, Ehnes et al. (2011) estimated (mean ± SE) b as 0.68 (± 0.4) for oribatid mites and 0.69 (± 0.09) for mesostigmatid mites. Our point is that the experimental error σ_b around b as well as the error σ_M around M are both needed to assess the measurement error of I, which is a combination of the two errors.

This can be demonstrated through an equation of the type I∝M^b to estimate I, given uncertainty σ_b and σ_M in b and M, respectively. From error propagation theory (Taylor, 1997), the uncertainty in I would then be $σ_{I} \approx I^{2} [{(\frac{M}{b} σ_{M})}^{2} + {(ln (M) σ_{b})}^{2}]$ under the assumption that there is no correlation in the error structure of M and b. If there was some non-negligible correlation between the two measurements, one more term should be added to the error of I, with the general effect of further increasing measurement uncertainty. Even more, usually the general equation used to estimate I is $I = i 0 M^{b} e x p (- \frac{E}{k T})$ where k is the Boltzman constant, T is temperature, E activation energy, and i0 is a taxon-specific constant. The latter two parameters (E and i0) are estimated experimentally with some errors, which should also be incorporated in the estimate of σ_I. Given that the I_j of each individual j is estimated with an error σ_I, and that the total metabolism of a trophic species is the sum of all individual metabolisms, that is I_tot =_j, the error of I_tot will be obtained by propagating the errors from the sum of all I_j. That would simply be σ_tot = ∑σ_I assuming no correlation between the measurements of the metabolic loss of each individual (a correlation would further increase the combined error for the total).

These estimates of individual biomass based on allometric scaling become compound errors with calculations into population level biomass, which are then used in calculations of flux in food web models. The estimate of the energy fluxes between two trophic species finally requires the knowledge of other coefficients, which express the efficiency of energetic transfers between trophic levels. Arguably, these coefficients are known with some errors and the errors should thus propagate to the estimate of the energy flux together with the error of the metabolic loss. All these errors are rarely, if ever, taken into account, and are (implicitly) assumed to be negligible (Barnes et al., 2014, 2018; Potapov et al., 2019, 2021), or are not taken into account (Gauzens et al., 2019). Similar issues apply to conversions from abundance to biomass data, assignment of species into trophic groups (see Buchkowski and Lindo, 2021), extrapolation of density data across scales, and in general any conversion based on measurements and coefficients that, being based on experimental measurements, are necessarily known with a degree of uncertainty intrinsic to the measurement process. We thus recommend that future soil biodiversity studies explicitly consider the issue of error propagation in the quantification of error measurements, as this is important to estimate lower and upper bounds for key ecosystem level quantities such as respiration and fluxes of C and N.

2.2. Demographic stochasticity

Individuals within a species differ in vital rates, especially survival and fecundity, that together determine individual fitness (Caswell, 2000; Lande et al., 2003). In very small populations (typically <100), this variability in vital rates is the source of demographic stochasticity, and can be measured by following cohorts of individuals, their survival and reproduction output (Lande et al., 2003). This is relatively straightforward, while time consuming, for vertebrates such as birds (Engen and Sæther, 1998; Sæther et al., 2000; Engen et al., 2001, 2003) but much more challenging for soil organisms, apart from laboratory studies on soil animals (Siepel, 1994; Søvik and Leinaas, 2003; Stamou, 2012). Despite the fact that direct measurements of demographic stochasticity are difficult to obtain for soil biota, there is theoretical ground to expect demographic stochasticity plays a minor role in structuring soil biodiversity. For example, following Lande et al. (2003), in simple, unstructured density independent stochastic models, the total variance $σ_{λ}^{2}$ of the finite rate of population growth λ equals $σ_{e}^{2} + \frac{σ_{d}^{2}}{N}$ , where $σ_{e}^{2}$ is environmental variance (defined in the next section), $σ_{d}^{2}$ is demographic variance, and N is total population size. Density independence is appropriate because demographic stochasticity applies to small populations, which are arguably far away from their carrying capacity. In large populations, population size $N ≫ \frac{σ_{d}^{2}}{σ_{e}^{2}}$ and demographic stochasticity plays a much smaller role than environmental stochasticity. A critical population size $N_{c} = 10 \frac{σ_{d}^{2}}{σ_{e}^{2}}$ can be defined, above which demographic stochasticity can be ignored. In birds, this critical population size can range from a few tens to various hundreds of individuals (Lande et al., 2003). Clearly, the issue is that we do not know the magnitude of environmental and demographic variance in soil organisms. However, it is textbook knowledge that the density per square meter of most soil organisms is typically well above 100 or 1,000 units for animals and up to 1 billion bacterial cells and 200 m of fungal hyphae per gram of soil microbes (Coleman et al., 2004).

It is thus unlikely that demographic stochasticity plays any role at all for most species at scales larger than 1 m² (see also Section 3). At scales smaller than 1 m², however, various soil animal species, such as Collembola and mites, can have population sizes of a few tens of individuals, making demographic variance a likely source of stochasticity. If we assume a small population size at this scale, and thus a density independent model following Lande (1998), and the formulation in Lande et al. (2003), we can also formulate an unstable stochastic equilibrium for population size that establishes a critical threshold below which the probability of extinction approaches unity. This critical equilibrium value can be calculated as $N_{^{*}} = \frac{\frac{σ_{d}^{2}}{4}}{\overset{}{λ} - 1 - \frac{σ_{e}^{2}}{2}}$ . For example, if $σ_{d}^{2} = 1$ and $σ_{e}^{2} = 0.04$ , and the average λ = 1.03, the critical population size would be 25 individuals. This is a likely scenario at some spatial scales (1–100 m) for some groups of soil fauna given existing estimates of population size and finite rates of increases (e.g., Caruso et al., 2020). However, at this scale interactions between individuals in soil become very likely, as also shown by highly aggregated distribution of soil fauna and microbes at small scales (Ettema and Wardle, 2002). The implication is that other biological mechanisms will potentially contribute to population fluctuations, especially Allee effects, which makes again pure demographic stochasticity unlikely to be a large and only source of fluctuations in soil biota.

These considerations are very relevant to the large body of literature that has investigated the relative roles of stochastic and deterministic determinants of soil biodiversity. A recent synthesis has been offered by Thakur et al. (2020), who found about 100 papers explicitly testing community ecology assembly theory in soil communities, and concluded that theory such as the neutral theory, which has demographic stochasticity as the core process, are generally not well-supported suggesting niche-based explanations underly community composition, and the high levels of unexplained variance are likely due to other stochastic processes described here (see next section). This conclusion is consistent with the general consideration that for most soil species, large population sizes over relatively broad scale (>100 m) imply a small role of demographic stochasticity. It is, however, important to remark that in many of the papers reviewed in Thakur et al. (2020) there is, in our opinion, very often confusion between the assumption of neutrality in certain models and the role of stochastic processes. As made clear by Adler (2007), demographic stochasticity is just a demographic process that may be at play both in neutral and niche community dynamics (see also the simple niche stochastic model example by Tilman, 2004). This means that the rejection of a neutral model does not imply that community fluctuations are not affected by demographic stochasticity, but only that the assumption of neutrality is not sufficient to explain the tested community patterns (Vellend, 2010).

2.3. Environmental stochasticity

Besides individual-level variability in vital rates such as birth and survival rate, the average expression of rates in the population also depends on the environment, which fluctuates. Sometimes conditions are favorable, sometimes they are not. After removing structural, periodic fluctuations due to, for example, the day and night cycle, or the seasonal cycle, or trends due, for example, to global warming, many environmental fluctuations are best modeled as probability distributions (Figures 2, 3), that is, stochastic fluctuations. The implication is that the vital rates, too, will fluctuate stochastically. A simple density independent model for this type of stochastic fluctuations is $N_{t} = N_{0} \prod_{i = 1}^{t} λ_{i}$ where the population size N grows geometrically from time 0 to time t according to the finite rate of increase λ, but at each time step the value of λ changes, because the environment changes (Botsford et al., 2019). This implies fluctuations in the rate of growth and so population size over time. It will also imply fluctuations in interspecific interactions, and all the ecosystem processes that depends on these interactions (Figure 1).

FIGURE 2

Figure 2. The “observed” population stochastically fluctuate around an equilibrium. The fluctuation follows a normal distribution with autocorrelation parameter “phi”. The parameter phi and the mean and variance of the observed population can be estimated from a sample of the population. The estimate of the mean, variance and autocorrelation parameter can be used to draw a desired number of “simulated” populations. In this example, three simulated populations (Simulated R1, R2, and R3) are shown together with the originally observed population. The three simulated population and the observed one represents four realizations of the same stochastic process, with the probability distribution of the process been estimated from the observed population. This time series displays positive autocorrelation, that is “red” noise.

FIGURE 3

Figure 3. A Fourier transformation decomposes the overall pattern of the time series to reveal the underlying frequencies present in the original function. As such, the Fourier transformation can be used to provide information on seemingly stochastic processes. Environmental variables typically display “red” noise or positive autocorrelation over time. In this figure, a time series of soil moisture with red and white noise (no autocorrelation) are simulated. A Fourier transformation of the two timeseries (top panels, red and white noise) demonstrates that in the spectrum of the red noise most of the variance is displayed over short time scales, with a negative correlation between variance and periods. The spectrum of the white noise is flat, with all periods displaying similar variance.

The first population level implication of environmental stochasticity is that the long-term growth rate of the population will be smaller than the expected average growth rate, and the higher the environmental variance the smaller the rate will be (Caswell, 2000; Lande et al., 2003). This also applies to populations structured by age and stages (Caswell, 2000; Tuljapurkar, 2013), which is the case with many soil animals (Walter and Proctor, 1999). However, while this is generally true for so called “white noise”, which is fluctuations without temporal autocorrelation (Figure 3), many time series of environmental variables display autocorrelation in their random component (Figures 3, 4). This is the so called “color” of noise, which has been shown to have different types and sometimes contrasting effects on the long-term rate of population growth (Ripa and Lundberg, 1996; Ruokolainen et al., 2009)—we argue this is particularly important in soil properties such as soil moisture. For example, in the time series we show in Figure 4, soil moisture and temperature, but especially soil moisture, clearly show fluctuations due to the day-night cycle and pulses of rain. After removing these fluctuations, the remaining random noise is positively autocorrelated, with highest variance displayed over short period of times. This is relatively easy to interpret: if soil becomes dry, it will stay so for a certain time until it rewets. As well as the rewetted soil may be subjected to another bout of drought but with loss of moisture that will take a certain time, meaning that the closer in time any two soil moisture measurements are, the more similar the measurements are likely to be (i.e., red noise). The implication is that both wet and dry conditions may persist longer than under purely white noise, and certain population models imply that red noise may reduce the risk of extinction (Ripa and Lundberg, 1996).

FIGURE 4

Figure 4. Two real world time series of soil moisture and temperature. In the time domain (top panels) both variables, but especially soil moisture, displays clear positive autocorrelation (over short time scales values tend to be similar). After detrending linear and periodic trends and isolating the stochastic component, the calculation of the autocorrelation coefficient phi confirms the presence of positive autocorrelation (i.e., red noise). With these parameters, one can now draw a desired number of stochastic trajectories, which can inform the modeling of population, community, and ecosystem level (e.g., energy fluxes and soil respiration) processes.

Regardless of the particular population model, however, an important tool to describe the structure of environmental stochasticity in time series analysis is the Fourier transformation (Bloomfield, 2004; Bush et al., 2017). The Fourier transform of a time series moves data from the time domain (x-axis) to the frequency domain (new x-axis); in the original time series, the y-axis is simply the measured variable (Figure 3), and in the Fourier transform the y-axis becomes variance (Figure 3). In practice, the Fourier transform allows detection of the frequencies or, conversely, time period mostly expressed in the data (i.e., that shows more variance). If all the frequencies are equally expressed, the pattern (or spectrum) of the Fourier transform is flat (white noise), and there is no relationship between frequency and variance. If low frequencies show more variance than high frequency there will be a negative correlation between frequency and variance (i.e., there will be positive autocorrelation in the time domain referred to as red noise). Fourier analysis of soil time series can thus quantify the autocorrelation structure of environmental stochasticity in soil, thereby elucidating the effect of this stochasticity on soil biota population and communities. This analysis is straightforward for abiotic variables especially for data collected from data loggers that allow the collection of time series of desired length and resolution (see example in Figure 4).

Unfortunately, the same is not possible for soil biota (Caruso and Bardgett, 2021). Yet, future modeling studies (see example in Figure 2) still benefit from incorporating forms of environmental stochasticity, that reflect the random structure observed in soil abiotic environments. For example, energetic food webs incorporate temperature as a key parameter via the fundamental metabolic scaling equation $I = i 0 M^{b} e x p (- \frac{E}{k T})$ . There are now soil food web models parameterized on field data (Potapov et al., 2021; Pettit et al., 2023) using the energetic approach (Moore and de Ruiter, 2012). Thus, in silico soil food web models (i.e., computer model estimations) can explore the effect of environmental stochasticity on energy fluxes using information on the stochastic structure that can be revealed by a Fourier analysis of soil temperature time series (Figure 4). These models could also simulate the impact of perturbation regimes, that changes the temporal structure of fluctuations in environmental variables.

3. Community level estimate of stochasticity

One major limitation in the study of stochastic fluctuations in soil biodiversity is the lack of time series (Bardgett and Caruso, 2020; Caruso and Bardgett, 2021). However, many ecological processes that happen over time potentially leave their signature in the spatial distribution of organisms, and there are some models that may allow estimate of environmental and stochastic variance from combined spatial and temporal time series (Lande et al., 2003; Botsford et al., 2019). We offer an example here, based on a publicly available time series of soil oribatid mites (https://doi.org/10.5061/dryad.tmpg4f4vt), which was the subject of previous work of ours (Caruso et al., 2020). The dataset consists of several 9-year time series replicated over multiple sampling locations. We fitted the model by Engen et al. (2002) to this dataset, and we provide an R script (see Data Availability Statement) that displays all the key step of the model fitting procedure. The key output of the model is a partitioning of the total community variance in terms of environmental variance, and a term that is the sum of demographic variation and overdispersion (which is relative to the underlying Poisson distribution assumed by the model). The environmental variance is further decomposed into two sources: general environmental variance, which is a forcing component that applies to all species, and species-specific environmental variance. The model takes the relative abundance of species, location and sampling time as input. It is fitted in stages: first, a bivariate Poisson distribution is fitted to the combined species relative abundance of all pairs of possible samples. Second, the variance and correlation parameters of each fit are combined with the spatial and temporal distance between each pair of samples. Third, an estimate is obtained for the parameters that control the dependency of the correlation parameter on the spatial and temporal distance between any two samples.

These parameters contain the main output, that is the estimate of environmental and demographic variance. Note however, as outlined in the previous section, because measurement errors are generally not well-incorporated into models and analyses, estimates of environmental and demographic stochasticity are likely overestimated. This model is based on the following main stochastic model, where the logarithm X of the abundance of any species i in the community has the following forms: $\frac{d X_{i}}{d t} = r_{i} - m X_{i} - A (X) + \frac{σ_{d}}{\sqrt{N_{i}}} D_{i} + σ_{e} S_{i} + σ_{E} E$ . Where the subscript i indicates a certain species, D, S, and E are Brownian motions and, respectively, represents demographic (D), species specific (S) environmental, and general environmental (E) effects. The sigma coefficient of each of this term (D, S, and E) is the variance we aim to estimate in order to quantify the relative importance of environmental (e+E) vs. demographic (d) stochasticity. The model is not neutral with respect to intrinsic growth rate r and the species level stochastic effect D and E, and it includes interactions between species A(X) as well as density dependence controlled by parameter m. With these simple assumptions, a crude but field-based estimate of environmental and demographic stochasticity is possible from the type of community data sets available to soil ecologists. Our calculations (see R script, Sup. Info.) show that for an old conifer forest (>100 yrs), the total variance observed in the temporal and spatial fluctuation of oribatid species could be partitioned as 85% environmental variance (σ_e+σ_E) and 15% demographic variance + overdispersion (that is σ_d+o). However, the method does not allow separation of demographic variance from overdispersion (the deviation in the relationship between mean and variance of the underlying Poisson distribution of the model), which is likely due to spatial aggregation observed at small scales (Engen et al., 2002). Regardless, the results empirically confirm the theoretical expectation for the main source of stochastic variation as environmental, but also non-neutral. Following, the 85% of total environmental variance can be partitioned into interspecific differences in the response to the general environmental variance (82% of the 85%, or 70% of the total) and local, species-specific response (18% of the 85%, or 15% of the total). The variance values were ${σ_{e}}^{2} = 0.36; {σ_{E}}^{2} = 1.37; {σ_{d + o}}^{2} = 0.29$ . Assuming no overdispersion we would have ${σ_{e}}^{2} = 0.36$ and ${σ_{d}}^{2} = 0.29$ , which we use to estimate an average critical population density of 12 individuals (using $N_{c} = 10 \frac{σ_{d}^{2}}{σ_{e}^{2}}$ ).

Most oribatid mite species in forest soils have densities above 50 individual m⁻². Using the formula for the critical equilibrium value N_c with a long-term stochastic rate λ = 1.6, ${σ_{e}}^{2} = 0.36$ and ${σ_{d}}^{2} = 0.29$ (see Caruso et al., 2020), N_c is <1, suggesting the risk of extinction is negligible. Although crude, the example shows that relatively short time series (9 years) with some spatial replication and community level data allows an estimate of environmental and demographic stochasticity in soil communities. The time scale of the time series might be tuned to the biological groups under consideration. For example, for bacteria a time series could occur over a single season. In the northern hemisphere at temperate latitudes, bacterial communities could be sequenced molecularly at 10 intervals from May to August (the growing season) in multiple locations with relatively affordability.

4. Emerging concepts and way forward

Stochasticity can be fully embraced in the study of the factors that control the structure and functions of soil biodiversity and aboveground-belowground linkages. Through modeling stochastic processes, we can better predict and identify populations, communities or ecosystems that are vulnerable to stochasticity in terms of fluctuations that can make soil biological activities uncertain and unstable (Bardgett and Caruso, 2020). This is important because instability in soil biological activities cause uncertainty and instability in fluxes that link aboveground and belowground components, and functions such as nutrient cycling. To obtain this goal, we need to improve data in terms of a quantification of measurement errors based on the theory of error propagation, and produce more time series of different soil biota across multiple locations. Experimental time series are the main tool to elucidate potential divergent (i.e., alternative stable state) soil-climate feedbacks through examining the resistance and resilience of soil biodiversity and function (Bardgett and Caruso, 2020). However, experimental tests and empirical data of this are lacking. In terms of new data and models, we need a quantification of the relative roles of different sources of stochasticity and explicit incorporation of stochastic autocorrelation structure into soil models (at population, community, and ecosystem levels). Improved experimental designs, theory and practical methods are already available to quantify and propagate measurement errors (for example see Caruso and Rillig, 2022). It is just a matter of applying them more extensively and critically, acknowledging errors at the various levels at which they occur (e.g., in energetic soil food webs). Time series remain challenging to obtain but there is an increasing number of contributions toward a temporal description of soil biota (e.g., Barreto et al., 2021), and especially for microbial groups the collection of time series within a single year or season is both feasible and informative of short terms dynamics.

We demonstrate that even a modest time series replicated over space allows a crude but reliable quantification of demographic and environmental stochasticity in soil communities and the quantification of the joint spatial and temporal autocorrelation observed both in biota and the environment, while a Fourier analysis of the stochastic components of time series can reveal the type of autocorrelation that characterizes stochastic fluctuations both in biota and the environment. We believe this approach will be very useful in future manipulative experiments, such as global change experiments that manipulate the intensity, frequency, and type of multiple perturbations (Rillig et al., 2019), which generate complex fluctuations.

Data availability statement

Publicly available datasets were analyzed in this study. This data can be found here: https://doi.org/10.6084/m9.figshare.22654375.v1.

Author contributions

TC, TB, and ZL all developed the ideas in this manuscript. TC coded and ran models and wrote the first version of the manuscript with ZL. All authors contributed to drafts and the final version of the manuscript with approval for publication.

Funding

Funding was provided by European Commission to TC (Structure and Ecological Niche in the Soil Environment; EC FP7-631399-SENSE) and by the Natural Sciences and Engineering Research Council of Canada (NSERC) Discovery Grants program to ZL (#418241-2012 and #05901-2019).

Conflict of interest

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

Publisher's note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

References

Adler, P. B. (2007). A niche for neutrality. Ecol. Lett. 10, 95–104. doi: 10.1111/j.1461-0248.2006.00996.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Bardgett, R. D., Bowman, W. D., Kaufmann, R., and Schmidt, S. K. (2005). A temporal approach to linking aboveground and belowground ecology. Trends Ecol. Evol. 20, 634–641. doi: 10.1016/j.tree.2005.08.005

PubMed Abstract | CrossRef Full Text | Google Scholar

Bardgett, R. D., and Caruso, T. (2020). Soil microbial community responses to climate extremes: resistance, resilience and transitions to alternative states. Philos. Trans. R. Soc. B Biol. Sci. 375, 20190112. doi: 10.1098/rstb.2019.0112

PubMed Abstract | CrossRef Full Text | Google Scholar

Bardgett, R. D., Freeman, C., and Ostle, N. J. (2008). Microbial contributions to climate change through carbon cycle feedbacks. ISME J. 2, 805–814. doi: 10.1038/ismej.2008.58

PubMed Abstract | CrossRef Full Text | Google Scholar

Bardgett, R. D., and Van der Putten, W. H. (2014). Belowground biodiversity and ecosystem functioning. Nature 515, 505–511. doi: 10.1038/nature13855

PubMed Abstract | CrossRef Full Text | Google Scholar

Bardgett, R. D., and Wardle, D. A. (2010). Aboveground-Belowground Linkages: Biotic Interactions, Ecosystem Processes, and Global Change. San Diego, CA: Oxford University Press.

Google Scholar

Barnes, A. D., Jochum, M., Lefcheck, J. S., Eisenhauer, N., Scherber, C., O'Connor, M. I., et al. (2018). Energy flux: the link between multitrophic biodiversity and ecosystem functioning. Trends Ecol. Evol. 33, 186–197. doi: 10.1016/j.tree.2017.12.007

PubMed Abstract | CrossRef Full Text | Google Scholar

Barnes, A. D., Jochum, M., Mumme, S., Haneda, N. F., Farajallah, A., Widarto, T. H., et al. (2014). Consequences of tropical land use for multitrophic biodiversity and ecosystem functioning. Nat. Commun. 5, 1–7. doi: 10.1038/ncomms6351

PubMed Abstract | CrossRef Full Text | Google Scholar

Barreto, C., Branfireun, B. A., McLaughlin, J., and Lindo, Z. (2021). Responses of oribatid mites to warming in boreal peatlands depend on fen type. Pedobiol. J. Soil Ecol. 89, 150772. doi: 10.1016/j.pedobi.2021.150772