Simple measures to capture the robustness and the plasticity of soil microbial communities
Takashi Shimada Kazumori Mise Kai Morino Shigeto Otsuka


APS/123-QED





Department of Systems Innovation, School of Engineering, The University of Tokyo
Mathematics and Informatics Center, The University of Tokyo

Graduate School of Agricultural and Life Sciences, The University of Tokyo
Graduate School of Science, The University of Tokyo
Bioproduction Research Institute, National Institute of Advanced Industrial Science and Technology

Interdisciplinary Graduate School of Engineering Sciences, Kyushu University

Graduate School of Agricultural and Life Sciences, The University of Tokyo
Collaborative Research Institute for Innovative Microbiology, The University of Tokyo









 

§ ABSTRACT 
Soil microbial communities are known to be robust against perturbations such as nutrition inputs, which appears as an obstacle for the soil improvement.
On the other hand, its adaptable aspect has been also reported.
Here we propose simple measures for these seemingly contradicting features of soil microbial communities, robustness and plasticity, based on the distribution of the populations.
The first measure is the similarity in the population balance, i.e. the shape of the distribution function, which is found to show resilience against the nutrition inputs.
The other is the similarity in the composition of the species measured by the rank order of the population, which shows an adaptable response during the population balance is recovering.
These results clearly show that the soil microbial system is robust (or, homeostatic) in its population balance, while the composition of the species is rather plastic and adaptable.





Simple measures to capture the robustness and the plasticity of soil microbial communities
    Shigeto Otsuka
    September 9, 2024
==========================================================================================




 
Understanding of the stability of ecosystems has been a central question in ecology, often with emphasis on local stability, persistence, permanence, resilience, etc., and its relation to its diversity <cit.>.
It has also been inspiring broader field to consider the robustness of complex systems which consist of many components interacting each other more in general <cit.>.

However, empirical tests of various theoretical predictions on animal-plant ecosystems have been limited within natural experiment, i.e. comparing the systems under different environment and with different history. This is mainly due to the large scales in space and time scale, for conducting controlled experiments.

Microbial ecosystems are one of the best systems to resolve such situation, because they are easy to manipulate and their development can be observed in a shorter time scale.
Bacteria in soil are also known to be extremely abundant and diverse, with estimates generally ranging from 10^6 to 10^10 cells <cit.> and 6,400-38,000 taxa <cit.> per gram of soil. These bacteria play vital roles in biogeochemical cycling, plant growth, and the maintenance of terrestrial ecosystems. Reflecting the importance of soil functions for terrestrial ecosystems, there has been considerable effort invested in understanding the response of soil ecosystems to disturbance <cit.>. 
The possibility of such approach was brilliantly illustrated in a pioneering work on the microbial communities of natural soil and natural water <cit.>.
In that study, it was shown that the macroscopic characteristics of the micro-ecosystems such as those dry weight and pH can be modulated by an artificial natural selection process.
And the difference was inherited to the following systems even after quitting the natural selection.
This shows that the microbial ecosystem have both plasticity and robustness or stability in macroscopic sense, while the underlying microscopic mechanism remained largely unknown.


Advance of DNA sequencing technology has enabled microbiologists to obtain detailed information on dozens of microbial communities at one time <cit.>.
This has paved a way to elaborate time-series investigation of microbial ecosystem successions <cit.>.
While many have elucidated the succession of microbial communities in a descriptive manner, for example by cataloguing a plethora of microbial clades (taxonomic groups) that increase or decrease during the ecosystem development <cit.>, others summarize the dynamics of the systems using statistical features such as

co-occurrence network indicators or ARIMA modeling <cit.>.
Other characteristics such as the shape of species abundance distribution <cit.> and the network structures of interactions <cit.> have been also the target of investigation.
However, the interpretation of these features is by no means straightforward.

In this study, we propose simple and comprehensible measures to characterize the state of microbial communities, based on the balance and the composition of populations of the operational taxonomic units (OTUs).
These two measures are independent each other in the sense that there always exists a direction to modulate the community so that one measure changes while the other is kept constant.
It is shown that these measures well capture the conserved aspect (relating to the stability, robustness, etc.) and the adaptive aspect (plasticity) of the microbial community.



§ RESULTS



 §.§ System and Preparation

As a model system of community dynamics in response to perturbation, we here use time-series data of soil microbial communities receiving nutrient input <cit.>. This dataset consists of soil microbial community structures (populations of OTUs) at five time points, i.e. 0, 3, 10, 17 and 24 days after initial preparation. As shown in TABLE <ref>, the nutrient input condition consists of the combination of nitrogen source: ammonium chloride (NH4Cl, denoted by A in the labeling of the conditions) or urea (U), carbon source: glucose (G) or cellobiose (C). We also have two types of input schedule: either once in the beginning of the observation (day 0) or addition of one quarter of the amount four times every week (days 0, 7, 14, and 24), denoted by 1 and 4 in the labeling of the condition, respectively. Negative controls without nutrient input were also prepared (“Control”).

This nutrient amendment schemes were designed to cover various intensities and types of perturbations.
Urea and cellobiose represent more slowly acting nutrition sources than ammonium and glucose, since the former two require enzymatic degradation before being catabolized by microbes.
Regarding the input schedules, the one-time and the four-time schemes represent large abrupt perturbation and small but continuous perturbation, respectively.

Since we have three samples for each condition, a soil sample α is specified by the combination of all these parameters as:
α = (t, d, s), where n ∈{ Control, AG1,AG4, AC1, AC4, UG1, UG4, UC1, UC4},  d ∈{ 0, 3, 10, 17, 24 }, and s ∈{ 1, 2, 3 } represent the treatment condition, observation day, and the sample number for each condition, respectively. Because the data d=0 is only for Control, the total number of the soil samples is (1 × 5 × 3) + (8 × 4 × 3) = 111.







 §.§ Species Abundance Distributions

We first confirm the basic characters of the control communities by its relative abundance distribution (RSA) of operational taxonomy units (OTU), which we shall call OTU-RSA.

For this purpose, we compare the fittings by classical theoretical distributions for species abundance of ecological systems, namely, log-normal distribution, power-law (Zipf/Pareto) distribution, and the negative binomial distribution.
An example of fitting is shown in the left panel of Fig. <ref>.
OTU-RSA of control soils are found to be better fit by log-normal distribution and power law distributions than by the negative binomial distribution, with a distinctive difference in the log-likelihood and equivalently in AIC, since all these distributions have two parameters.
The obtained fitting parameter values consistently indicate that the OTU-RSA of the soil microbe ecosystem before any disturbance is broad (see Materials and Methods and SI for detail).
This feature is consistent with the systems under natural environment<cit.>, implying that the control soils are keeping such an intact state.

The abundance distributions of the soils under or after nutrition input are also fat-tailed, with typically narrower shape comparing to that of control soils (Fig. <ref>, right panel).
A typical response of RSA to the nutrition input is characterized by an intial decrease of diversity
which is measured by the number of detectable OTUs or by Shannon entropy of OTU abundances

    S = - ∑_r p(r) ln p(r),

and an increase of the dominance of top tens of OTUs.
This decrease is followed by the recovery to a smoother distribution which is nearer to that of original control soils (TABLE <ref> and Supplemental Information).
Such changes in RSA clearly illustrates a homeostatic aspect of the soil microbial ecosystem: Fertilization typically first leads to more oligopolistic abundance distribution but the system later shows a resilience toward the original population balance.



 §.§ Similarity in RSA shape

The above observation about the typical change in the shape of OTU-RSAs naturally motivates us to define a similarity measure solely based on that.
A simple and natural choice is to take the linear correlation coefficient between the relative OTU abundances of soil α and soil β in the same rank r,

    B_αβ^R_B≡∑_r=1^R_B( p^α_r - p^α) ( p^β_r - p^β)/√(∑_r=1^R_B( p^α_r - p^α)^2)√(∑_r=1^R_B( p^β_r - p^β)^2 ),

where p^χ_r represents the relative abundance of the OTU whose abundance rank is r in soil χ, and p^χ = . ∑_r=1^R_B p^χ_r / R_B denotes the average abundance up to the maximum rank R_B to be taken into account. Since OTU-RSA has broad shape, the maximum rank R_B in the following will be set at 1,000 to include the information from as diverse OTUs as possible.

The average B similarities of the soil communities under each treatment condition t observed on day d to the control soil,

    B̅(t, d) = ∑_t_α = t, d_α = d;  t_β =  ControlB_αβ/∑_t_α = t, d_α = d;  t_β =  Control 1,

are shown in Fig. <ref> (bars).
We can see that the temporal evolutions of the soils in this measure against various inputs (different combinations of nutrition and the different schedules) universally show a homeostatic response, i.e. the similarity to the original control community first drops and that recovers later on.
The response of the soils against the continuous inputs is particularly interesting, because the recovery takes place during the successive nutrition input. This is why we regard it homeostatic, rather than a mere relaxation back to the original state under absence of the perturbation.






 §.§ Similarity in OTU abundance rank orders

We next seek a measure to capture a totally different aspect of the soil microbial community.
Because the systems have been confirmed to show robustness in its population balance (the shape of RSA), the remaining degree of freedom in OTU abundance distribution is in its composition.
Suppose that two communities have the same OTU-RSA, then yet two communities can be different if OTUs occupying each rank position are different each other.
In other words, the order of OTUs in relative abundance (whether a particular OTU is major to another particular OTU) is independent of the shape of OTU-RSA distribution.

Our similarity measure for the composition of OTUs, with respect to the OTU-RSA information, is defined by the Kendall's rank correlation between the abundance ranks of top N_C major OTUs of the two communities in comparison, α and β:

    𝒞_αβ^N_C≡#_c - #_d/N_C (N_C -1)/2,

where #_c and #_d represent the numbers of concordant OTU rank pairs, i.e.
( r^α_i - r^α_j ) ( r^β_i - r^β_j ) > 0,
and discordant OTU rank pairs, i.e.
( r^α_i - r^α_j ) ( r^β_i - r^β_j ) < 0,
respectively.
Where r^χ_k represents the abundance rank of OTU k in the soil χ. If a certain OTU is not detected in the sample, we replace r^γ_k by a uniform number r_max which is larger than the total number of OTUs. Therefore, for example, if both OTUs i and j are missing in a soil α, any combination of OTU abundance ranks in the other soil β (r_i^β and r_j^β) will not be counted as concordant or discordant.
The parameter N_C is chosen to be large enough to take as much information as possible and not too large to avoid an error from the detection limit.
To determine proper N_C, we first check how many OTUs are commonly observed in the soils under different conditions.
For this purpose, we evaluate the overall popularity (majority) of each OTU by aggregating the all OTU abundances across the available conditions.  
We find that almost all of top 300 major OTUs in the aggregated data are found to be shared in all soils (shown in SI, Fig. <ref>).
Some of the OTUs in the following overall majority rank range becomes to be absent (undetectable) in some soils.
Therefore, in the following we choose N_C = 100, OTUs up to this are expected to provide the majority-minority information well above the detection limit of OTU abundances in the current experiment.
One can confirm that the following results are not sensitive to the precise choice of N_C (for example, for N_C = 30 and N_C = 300).

The present measure 𝒞 takes 1 if the rank orders between OTUs (major-minority relations) is kept for all pairs, and it takes a value near 0 if there is no correlation.
It can take negative value if the reversal in the majority order is dominant, ultimately down to -1 when all the majority order pair are flipped.
Because this measure is independent of the shape of RSA distribution, and also because the swap in the ordering of OTU abundances does not change the RSA shape, the previously defined “population-balance measure” B and the currently introduced “composition measure” C can be said to be independent, or orthogonal, to each other.

The OTU composition similarities 𝒞 measured using the top N_C = 100 major OTUs are shown in Fig.<ref>, where bars again represent the average similarity of the samples with treatment condition t and the observation day d to the control soils:

    C̅(t, d) = ∑_t_α = t, d_α = d;  t_β =  ControlC_αβ/∑_t_α = t, d_α = d;  t_β =  Control 1.

The similarity among the control soil is confirmed to be high in this metric. Compared to that, the average similarities of the soils 3 days after nutrition inputs to the control soil is low, meaning that the composition of OTUs are changed greatly from the control soils.
What is remarkable is that the similarity of the perturbed soils to the control soils shows further decrease in the later time.
This tendency is universally observed for different conditions, and for the the choice of N_C as shown in the Table I and SI.
Recalling that the ℬ similarity shows an increase in the same moment of later time stage, the observed further change in OTU composition takes place during the recovery process of the population balance.
This implies that the response in the composition measure 𝒞 is not just due to a slower response, but is capturing the plasticity of the soil community to the nutrition input.





§ DISCUSSION

In this study, we have proposed simple similarity measures to characterize the response of the soil microbial communities to nutrition inputs.
The first one, ℬ, is based on the balance of the populations (shape of OTU-RSA), and the other one, 𝒞, is based on the composition of OTUs.
The both similarities among the control soil communities are kept during our observation time of 24 days, meaning that the stationary character of the communities under no external inputs is well captured by these macroscopic measures.

When the nutrition inputs are applied, the immediate response of the communities appears as the significant drops in both similarities.
The typical temporal evolutions of the similarities in the later time are, however, different each other.
The  ℬ similarities of the perturbed soils to the original communities universally show a recovery in time. The recovery is observed even under the continuous nutrition input.
Therefore this measure universally characterizes the robust (or homeostatic) nature of the soil communities.
In 𝒞, on the other hand, an opposite character of the soil's response is detected. The 𝒞 of the perturbed soils generally show further decrease in time, in the middle of the recovery period in ℬ similarity (10, 17, and 24 days after the onset of nutrition inputs).
This further decrease is universal among different nutrients and input schedule, while the continuous input tends to results in more monotonic decay.
Therefore the plastic (or adaptive) aspect of the soil communities is also well captured by the 𝒞 similarity.

Most of the attempts which have so far been made to disentangle the dynamics of complex microbial community structures rely on beta-diversity metrics between spatially/temporally distant communities. The concept of beta-diversity appears to be two-fold. While it was originally proposed as a between-community difference in the presence/absence of each species <cit.>, currently popular metrics such as Bray-Curtis <cit.> and weighted UniFrac distances <cit.> bear the information on the quantitative balances between community members.
The two indices we have proposed here are approximately congruent with these beta-diversity-based distances, in the sense that both ℬ and 𝒞 are reflected.
This can be seen in the present data that the drastic changes the soil bacterial communities underwent after the initial perturbations (i.e. between days 0 and 3) is visible by the beta-diversity measures (Fig. <ref> in SI). However, in these measures, the rapid recovery in species balances (i.e. change in ℬ) has overwhelmed and masked the long-lasting effect on the change in microbial community members (i.e. change in 𝒞). It is also notable that the contrast between two schemes of nutritional input schedule, one-time and four times, are obscure in ℬ (and in canonical beta-diversity) but clearly highlighted by 𝒞.
Therefore, by decomposing the canonical beta-diversities into two features by ℬ and 𝒞, otherwise mingled differences between incubation conditions are better illuminated.

For the understandings of the robustness and plasticity of interacting communities in more general sense, present discovery can be regarded as a good guiding principle.
Let us consider trying the construction of a dynamical community model, under a constraint that which can explain the coexistance of observed robustness and plasticity. As is evident from the large change in OTU-RSA shape, the homeostatic response to the nutrition input we have observed by ℬ is beyond linear (local) stability argument. Therefore, the observed trajectory should be explained as a transition from an attractor (a stable state) to another attractor triggered by the change of nutrition parameter or a large disturbance in population distribution. A challenge in such approach is that one should reproduce the order of the response: the recovery in OTU balance takes place first and the dynamics of the OTU composition change is slower.
Another possibility to explain the robustness in OTU-RSA, with keeping the degree of freedom of change in the OTU composition, would be a more statistical (entropic) mechanism under some hidden microbial constraints. The neutral theory is a good candidate for this, though it is not straightforward to test that for the observed response. This is because the environment in this study is dynamically changing and hence the key parameters in the framework, such as the birth and death rate, dispersion, and so on of each OTU are not kept. And even if the effect of nutrition input is nicely modeled by the abrupt change in such parameters, to reproduce the responses in the later time remains as a difficult task.

Finally, present finding is also meaningful in the agricultural context. It is important in agriculture to put organic matter such as compost into the soil to maintain the quality of the soil. Fertilizers as well are basically essential for producing crops on agricultural field soil. And, giving carbon source and/or nitrogen source to soil in this way can be regarded as one of the disturbance of soil ecosystem. It was reported that fertilization had a greater effect on soil bacterial community structure than crop rotation and growth stage (Guo et al. 2020). It is also shown that soil organic carbon and total nitrogen are two of the major determinants of community composition (Li et al. 2017). The disturbance of the shape of the OTU-RSA curve for the soil amended with carbon and nitrogen sources (glucose and ammonium chloride, respectively) (Fig. 1) clearly indicate a rapid changes in bacterial community structure. In addition, when combined with traditional analyses such as non-metric multidimensional scaling (NMDS) (SI Fig. 4), the overall picture of the changes that occur in the community structure becomes easier to understand.



§ METHODS



 §.§ Sample Preparation and Primary Sequence Analysis

Nucleotide sequence files registered in DDBJ/ENA/GenBank under the accession numbers DRR157393-157518 (DRA007564) and DRR169428-169475 (DRA008081) were retrieved. These sequences were generated by amplicon sequencing targeting 16S rRNA genes in soil microcosms amended with labile carbon and nitrogen sources (partly described in <cit.>). These genes are the commonly-used "marker" for the identification of prokaryotes. The soil microcosm incubation experiment was designed in a full-factorial manner: besides control samples without nutritional amendment, eight treatment groups, differing in carbon source (glucose or cellobiose), nitrogen source (ammonium chloride or urea), and timing(s) of nutritional amendment (weekly amendment over the experiment period or abrupt amendment of quadruple amount at the beginning of incubation), were prepared (Table 1). The total amount of carbon and nitrogen input was the same in all groups. Twelve microcosms were constructed for each group, three of which were destructively sampled after 3, 10, 17, and 24 days of incubation. Additionally, the pre-incubation (i.e. day 0) soil was examined in triplicate. Sampled soil was subjected to total DNA extraction, purification, amplification of partial 16S rRNA gene sequence, and Illumina high-throughput sequencing as described previously <cit.>. 

Low-quality (i.e. expected errors of 0.5 bases or more) sequences were removed using USEARCH v9.2.64 <cit.>. Subsequently singleton sequences were removed, and the remaining sequences were clustered into operational taxonomic units (OTUs) at a similarity threshold of 97% or more using UPARSE <cit.>. The OTU representative sequences (i.e. most frequently-observed sequence within those clustered into the OTU) were subjected to UCHIME2 <cit.> to filter out chimera-like OTUs. Finally, all the quality-filtered sequences, including singletons, were mapped to OTU representative sequences using USEARCH v9.2.64. To detect organelle-like OTUs, the OTU representative sequences were taxonomically annotated by RDP classifier <cit.> trained with Greengenes 13_8 clustered at 97% <cit.>. OTUs annotated as "family mitochondria" or "class Chloroplast" were precluded from further analyses.



 §.§ Species Abundance Distributions

To characterize the SAD, we compare the fittings by log-normal distribution

    P^LN (μ, σ; x)
    =
    1/σ x √(2π)exp{- (ln x - μ)^2/ 2σ^2},

power-law (Zipf/Pareto) distribution

    P^PL (α, θ; x) = αθ^α/(x+θ)^α+1,

and the negative binomial distribution

    P^NB (s, p; x) = Γ(x+s)/Γ(x+1)  Γ(s) p^s (1-p)^x,

where x denote the OTU abundance.

The obtained fitting parameter values;
lnμ∼lnσ∼ 1.5 for the log-normal distribution,
shape parameter 1.0 < α < 1.5 for Pareto distribution,
and the size 0 < s < 0.75 and the mean p ∼ 0.01 for the negative binomial distribution,
consistently indicate that the OTU-SAD of the soil microbe ecosystem before any artificial nutrition input is very broad, i.e. critical or fat-tailed (see SI for detail).




§ ACKNOWLEDGEMENT

TS was partly supported by JSPS KAKENHI grant number JP18K03449 and JP23K03256.
TS and SO were partly supported by JSPS KAKENHI grant number JP17K19220.
KMo was partly supported by JSPS KAKENHI Grant Number JP20K19885.











§ SUPPLEMENTARY INFORMATION




 §.§ NMDS





 §.§ ℬ similarity to control soils





 §.§ Similarity of the control soils in OTU composition






 §.§ Characteristics of RSA distributions


















 §.§ Sample Accession ID




