General considerations of the model
We consider the simple case of a host-microbiome system in which every host is associated with a single species of bacteria that is transmitted to the host's offspring with perfect fidelity. We take the generation time of a host to be much larger than for bacteria and we probabilistically determine the survival of each host and bacterium, according to their state at the end of the respective generation time (as detailed below). Each surviving host and bacterium gives rise to one offspring that inherits the traits of its parent, subject to a small random modification depending on a constant mutation rate, μ (no epigenetics is considered). The host and its bacterial community are jointly exposed to a toxin of concentration T, thus creating a stress that impacts the survival probability of the host and each of its bacteria. This stress depends both on their intrinsic traits and on how they interact with one another. To investigate whether and how the coupling between the survival of host and bacteria could support non-traditional modes of adaptation, we consider broadly applicable types of interactions between host and bacteria. The mathematical representations of these interactions was chosen to simplify the identification and analysis of general effects which apply to many host-microbiome systems (as opposed to a model designed to fit a specific system).
We start by defining the toxic stress experienced by individual host and bacterium. Since this stress depends on the level of toxin, T, and on the individual’s sensitivity to the toxin, x, we define the instantaneous toxic stress for host and bacteria as SH = xH T and SB = xB T, respectively. Accordingly, this stress can be alleviated by cell-intrinsic reduction in sensitivity and/or by secretion of a detoxifying agent, “detox”, which reduces the toxic challenge (with or without associated cost).
Unlike in a germ-free system, a host in a composite host-microbiome system is influenced (and/or dependent) on bacterial-derived nutrients and various other factors [43,44,45]. Exposure to toxin may therefore lead to physiological stress to the host due to a significant loss of bacteria. An indirect stress to the host can also be induced by factors that promote a significant excess of bacteria. We model these effects by considering a physiological stress, Sph, which depends on deviations from a “preferred” size of the bacterial population. From the bacterial perspective, on the other hand, the host provides a niche of a particular size (carrying capacity for bacteria). In the simpler case of free-living bacteria in a fixed environment, the carrying capacity is typically modelled by a constant parameter, representing the amount of extractable resources. The fixed niche assumption, however, does not necessarily hold when the bacteria are accommodated inside a host which can modulate the size of the niche under stress [23]. Since we do not know in advance whether and how a host’s stress influences the number of bacteria that can be accommodated, we constructed a population model in which this influence is determined by natural selection.
Altogether, the model considers host-microbiome interactions that are mediated by: (i) mutual alleviation of toxic challenge via secretion of a detoxification agent (“detox”), (ii) dependence of the hosts’ well-being on the size of the bacterial population and (iii) modulation of the bacterial niche size based on the toxic stress experienced by the host.
Model formulation
For each host and bacterium, we assign a probability of survival to reproduction, PH and PB, defined respectively as:
$$ {\boldsymbol{P}}_{\boldsymbol{H}}=\left(\mathbf{1}-{\boldsymbol{N}}_{\boldsymbol{H}}/\mathbf{2}{\boldsymbol{K}}_{\boldsymbol{H}}\right)\mathbf{\exp}\left[-\left({\widehat{\boldsymbol{S}}}_{\boldsymbol{H}}+{\widehat{\boldsymbol{S}}}_{\boldsymbol{P}\boldsymbol{h}}\right)\right] $$
(1)
$$ {\boldsymbol{P}}_{\boldsymbol{B}}=\left(\mathbf{1}-{\boldsymbol{N}}_{\boldsymbol{B}}/\mathbf{2}{\boldsymbol{K}}_{\boldsymbol{B}}\right)\mathbf{\exp}\left(-{\boldsymbol{S}}_{\boldsymbol{B}}\right) $$
(2)
Here, NH and NB are the population sizes of hosts and bacteria per host, respectively, KH is the maximal number of hosts that can be supported by the external environment (carrying capacity for hosts) and KB is the number of bacteria that can be accommodated in the host (carrying capacity for bacteria). The toxic and physiological stress to the host, Ŝ H = <SH >t and Ŝ Ph = ln(<NB>t/KB0) + (1 - <NB>t/KB0), are defined respectively in terms of time averages of SH and NB over a host generation (interval between host reproduction events; recall that the probability of survival is calculated only at the end of each generation). The physiological stress vanishes when the time-averaged bacterial population, <NB>t, reaches a size determined by the fixed parameter, KB0. The latter also sets an inverse scale (1/KB0) for the negative impact of losing too many bacteria or having to support excess numbers of bacteria [44].
To test if selection might favor hosts that react to toxic stress by modulating the niche available for bacteria [46,47,48,49], we consider a population of hosts, each with a distinct dependence of the carrying capacity on the toxic stress of the host. For that, we define KB as:
$$ {\boldsymbol{K}}_{\boldsymbol{B}}={\boldsymbol{K}}_{\boldsymbol{B}}^{\mathbf{0}}\left(\mathbf{1}+\boldsymbol{\delta} \cdot {\boldsymbol{S}}_{\boldsymbol{H}}\right) $$
(3)
where δ is an evolvable trait, determining how the bacterial niche in the host is affected by the toxic stress it experiences. Since bacteria can affect this stress by secreting detox on a timescale shorter than a host generation, KB is jointly influenced by the host and the bacteria. To enable unbiased analysis of how KB changes in response to selection under exposure to toxin, we considered a starting population of hosts with a broad distribution of δ’s, symmetric around zero.
We assume that all the hosts and their bacteria are exposed, at time t, to the same influx of active toxin, θ(t), applied instantaneously (i.e., in one bacterial generation, Δt). This toxin can be neutralized by release of detox from the host and each of its bacteria [44, 49,50,51]:
$$ \boldsymbol{T}\left(\boldsymbol{t}+\boldsymbol{\Delta} \boldsymbol{t}\right)=\boldsymbol{T}\left(\boldsymbol{t}\right)\ \mathbf{\exp}\left(-{\boldsymbol{\uplambda}}_{\mathbf{B}}\mathbf{\sum}{\mathbf{y}}_{\mathbf{B}}-{\boldsymbol{\uplambda}}_{\mathbf{H}}{\mathbf{y}}_{\mathbf{H}}\right)+\boldsymbol{\theta} \left(\boldsymbol{t}\right) $$
(4)
where yH and yB are the instantaneous amounts of detox secreted inside the host (by resident bacteria and the host itself) and λH and λB are the respective detoxification capacities of host and bacteria. We assume that all the bacteria of a given host benefit equally from the total amount of detox, regardless of their individual contributions to this total detox. The effect of having a cost associated with the secretion of detox by the bacteria is investigated in an extended version of the model (Additional file 1).
The evolvable traits of the model (x, y and δ) are initially drawn from trait-specific distributions and are modified by the joint actions of mutation and selection. Surviving bacteria divide at every time step of the simulation (Δt), while the surviving hosts reproduce every τ generations of bacteria (so that the host generation time is τ Δt). We consider the simplest reproduction model in which each of the surviving hosts and bacteria gives rise to one offspring that inherits the traits of its parent, subject to a small random modification depending on a constant mutation rate, μ:
$$ {\boldsymbol{Z}}_{\mathbf{offspring}}={\boldsymbol{Z}}_{\mathbf{parent}}+\boldsymbol{\eta} \mathbf{\surd}\boldsymbol{\mu } -{\boldsymbol{\beta}}_{\boldsymbol{z}}\boldsymbol{\mu} \left({\boldsymbol{Z}}_{\mathbf{parent}}-{\boldsymbol{Z}}_{\mathbf{0}}\right) $$
(5)
Here Z corresponds to any of the evolving traits x, y and δ, η is a standard Gaussian deviate with zero mean, and the parameters, Z0 and βz, are trait-specific coefficients controlling the peak and width of the steady state distributions (specified in Methods). Note that 1/βz sets a characteristic time for the distribution of a trait Z to return to steady state, following an initial perturbation. The values of βy and βδ were chosen to support broad distributions of y and δ, respectively. To prevent a trivial solution in which all the individuals are completely insensitive to toxin, the sensitivity distribution (for Z = xH and xB) is truncated at x = 0. We also avoid negative values of detox secretion by setting negative y values in Eq. 5 to zero. The remaining dynamic variables are updated in every generation of bacteria (NB, T, SH, SB and KB) and host (NH). The current study was based on an initial population of 32000 hosts (NH = KH = 32000) with 100 bacteria per host (NB = KB0 = 100). The host generation time was set to τ = 100 bacterial generations and all the mutation rates were μ = 10− 3 per generation (for both host and bacteria).
Stress-dependent adjustment of bacterial niche size
We examined the effects of exposure to a single pulse of toxin, T0, applied at t0 (i.e. θ(t0) = T0). On timescales smaller than one host generation (100Δt), the bacterial community undergoes selection for less sensitive bacteria, accompanied by a drop in the bacterial population size (Figs. 1a,b). In a system with only one level of selection (e.g. free-living bacteria), this would be the only adaptive change. However, when the bacterial population is symbiotically coupled to a host, the survival of each host and bacterium depends also on the amount of detox secreted by the bacteria (Fig. 1c). The secretion is higher for hosts which react to the toxic stress by increasing their carrying capacity for bacteria (i.e. hosts with δ > 0; Additional file 1: Figure S1A). This leads to stress-dependent selection of hosts which provide a larger bacterial niche KB (Fig. 1d), thus increasing the number of resistant bacteria beyond KB0 (Fig. 1a). The benefit from this increase is two-fold: Alleviation of the negative impact of losing bacteria (by assisting recovery of the bacterial population; Fig. 1e, Additional file 1: Figure S1B) and elevation of the total amount of secreted detox (Fig. 1f, Additional file 1: Figure S1C). However, when <NB>t is larger than KB0, the benefit from higher detox secretion is counteracted by the negative impact of bacterial overload. The overall effect of these positive and negative factors adjusts the bacterial population size in a stress-dependent manner which tends to maximize the probability of survival of the host.
Stress-dependent adaptation within a host generation
Microbiomes that are modified by the stress in one host generation can be transmitted to the host’s offspring, potentially increasing its stress tolerance. In order to evaluate the possibility and magnitude of this outcome, we introduce a new measure, termed the “Lamarckian”. It quantifies the change in the survival probability of the host's offspring due to (stress-dependent) microbial variations that were induced during the lifetime of the parental host. To take into account only those changes that were induced by the environmental stress, we compared the survival of offspring hosts to the survival of their parents as determined by the initial state of these parents. To implement this analysis in the simulation, we identify the hosts which survived a generation of exposure, revert them to the initial state of their microbiome and apply a new simulation to the reverted hosts (denoted “cloned parents”) and their offspring (Fig. 2a). We then compare the survival rates of the offspring (SROffs) to that of their cloned parents (SRCP) and define the Lamarckian, L, as:
$$ \boldsymbol{L}={\boldsymbol{SR}}_{\boldsymbol{Offs}}/{\boldsymbol{SR}}_{\boldsymbol{CP}}-\mathbf{1}, $$
(6)
so that it is positive if the average survival increases due to transfer of changes acquired during the previous host generation. The use of the initial state of the parental host and its microbiome allows us to distinguish the increase of tolerance due to selection of initially better fit parents, from the gain of tolerance due to transmission of changes acquired during a host generation (not present in the initial parental clones).
For a given λB, we find that L is an increasing function of the injected amount of toxin, vanishing only at low T0 (Fig. 2b). For a given T0, on the other hand, the Lamarckian has a non-monotonic dependence on λB. This is manifested by a nearly constant L > 0 over a range of small λB, followed by an increase to a maximum at intermediate values of λB and lastly, a decline at sufficiently large λB (Fig. 2c).
The positive Lamarckian is the result of transgenerational transfer of a bacterial population that was selected for lower toxin sensitivity during the parental host generation (Fig. 2d). To determine how these bacteria increase the probability of survival of the hosts’ offspring, we analyzed the toxic and physiologic stress in the offspring vs. their cloned parents. For small enough λB, the benefit from bacterial secretion of detox is negligible and the positive Lamarckian is primarily due to alleviation of the physiological stress in the offspring (Additional file 1: Figure S2). This is due to inheritance of bacteria that are less sensitive to toxin (Fig. 2d), so that the population size of bacteria in the exposed offspring remains closer to the preferred value (KB0) compared to the bacterial population size in their cloned parents. At intermediate values of λB, the offspring have an additional benefit due to the detox secreted by their toxin-resistant bacteria, thus making a second contribution to the Lamarckian (Fig. 2e,f). However, when λB is large enough to support substantial neutralization of toxin during a single host generation (Additional file 1: Figure S3), the selection pressure on both hosts and their microbiomes is weakened and the Lamarckian decreases because of the diminished difference between parents and offspring (Additional file 1: Figure S2B).
Selection of hosts based on collective traits of their bacterial community (‘Microbiome selection’)
When the toxic pressure persists over timescales larger than one host generation (Fig. 3a), the selection favors hosts with bacterial communities that secrete higher amounts of detox per bacterium, <yB >P (Fig. 3b). Since this selection is determined primarily by the microbiome as a whole and not by individual bacteria, we will refer to it as Microbiome selection. When the secretion of detox comes at a cost to the individual, the microbiome selection for detox is weakened, but it is still apparent over a broad range of cost levels (Additional file 1: Figure S4A,B). The negative effect of the cost on the survival probability of bacteria (Supplementary Information, Eq. 2′) aggravates the initial loss of bacteria and increases the physiological stress to the host (Additional file 1: Figure S4C). This promotes selection of hosts that can partially alleviate this stress by accommodating larger numbers of bacteria (Additional file 1: Figure S4D). The cost on bacterial detox therefore strengthens the selection of hosts which accommodate more bacteria at the expense of weakening the selection for increased detox per bacterium. The Lamarckian effect, on the other hand, is not compromised by the cost of detox (Additional file 1: Figure S5A,B), because the increase of physiological stress in parental hosts is larger than the corresponding increase in their offspring (Additional file 1: Figure S5C).
In the current model, Microbiome Selection occurs only at the time of host reproduction. If the toxin persists over a period longer than μ− 1 bacterial generations and the elimination of mutations is sufficiently slow (i.e. small enough βy), the selection is accompanied by significant accumulation of bacterial mutations. Such accumulation enhances the selection for higher <yB>, thus increasing the detoxification rate (Fig. 3a, inset) and expediting host adaptation (Fig. 3c). This is accompanied by extended persistence of high detox levels (Fig. 3b) and by elevated detox variability across host-microbiome systems (Fig. 3b, inset). Additional increase of variability under stress is observed in the carrying capacity for bacteria and in the size of the bacterial population (Additional file 1: Figures S6A,B).
Following the neutralization of toxin, the selected bacterial mutations persist over a characteristic timescale of 1/μ = 10 host generations, thus providing a ‘memory’ of the previous exposure. To evaluate the influence of this ‘memory’ on the tolerance to new exposures, we analyzed the response to repeated pulses of injected toxin, separated by time intervals shorter than 10 host generations. These re-exposures led to repetitive microbiome selections occurring at a rate that is sufficient to oppose the relaxation of <yB>P to its (lower) equilibrium value (Fig. 4a vs. Fig. 3b). The resulting enhancement of detoxification (Fig. 4b), reduced the selection pressure on the host (Fig. 4c) and enabled the survival of intrinsically less resistant hosts and bacteria (Fig. 4d). Progressive reduction in the intrinsic resistance of the host due to repetitive selections of higher bacterial detox, is reminiscent of Genetic Assimilation by successive selections of host-intrinsic alleles [52, 53]. In the case of Microbiome selection, however, the gradual change in the population of hosts is caused by recurrent selection of variations in the bacterial population (analogously viewed as “Bacterial Microbiome Assimilation”). Bacterial variations emerge on faster timescales compared with germline mutations in the host genome, but they are considerably less stable than host-intrinsic mutations. Nonetheless, when the repertoire of host-intrinsic alleles available for selection is limited, the hosts’ population may become more strongly dependent on variations that emerge within the host’s lifetime (e.g. bacterial and epigenetic variations).
Potential strategies for Lamarckian estimation in experimental settings
Quantification of the Lamarckian in the model was done by reverting a subset of host-microbiome systems to their initial state and re-subjecting them to toxin. Since we cannot apply this procedure to experimental data, the Lamarckian of a real system would need to be approximated by other means, which may be context-dependent. In organisms such as flies and worms, where the bacteria can be removed without a significant impact on survival (e.g. by egg dechorionation followed by rearing on a sufficiently rich diet [54,55,56]), the Lamarckian can be approximated in steps that are conceptually similar to the simulation procedure: first, the hosts are exposed to a challenge and their offspring are cleared of bacteria and separated into two subpopulations. One of these subpopulations is re-colonized with (‘naïve’) microbiota from untreated hosts (as in refs. [35, 57]), while the other is colonized by (‘experienced’) microbiota from a group of hosts which survived exposure to a challenge. The Lamarckian would then be estimated from the ratio between the survival rates of hosts with experienced vs. naïve microbiota (i.e. L ≈ SR Exp. microb. / SR Naïve microb. - 1). This evaluation, however, neglects other types of changes that may have been acquired and transmitted to offspring (e.g. transfer of small RNAs [10], altered deposition of maternal RNA [58], persistent chromatin modifications [8], horizontal transfer of biochemical signals [59] and/or other modes of local niche construction [14]). Additional consideration that may affect the evaluation is horizontal transmission of bacteria to bystander hosts and/or to offspring of other hosts. The above effects can be taken into account by removing the bacteria from two untreated subpopulations, re-colonizing hosts with ‘naïve’ and ‘experienced’ microbiota, respectively, and estimating the Lamarckian from the survival of these colonized populations under challenge. More generally, it should also be possible to obtain a relative measure of the Lamarckian by manipulating the microbiome (or any other factor) in a subpopulation of hosts and evaluating the relative difference in offspring adaptation compared to offspring of non-manipulated parents (taken from the same distribution of hosts).