Darwinian Selection Induces Lamarckian Adaptation in a Holobiont Model

Current models of animal evolution focus on selection of individuals, ignoring the much faster selection of symbiotic bacteria. Here we take host-symbiont interactions into account by introducing a Population Genetics-like model of holobionts exposed to toxic stress. The stress can be alleviated by selection of resistant individuals (host and bacteria) and by secretion of a detoxification agent ("detox"). By defining a new measure, termed the Lamarckian, we show that selection of resistant bacteria over one generation of hosts leads to stress-dependent increase in the tolerance of the host's offspring. This benefit is mediated by co-alleviation of toxic and physiologic stress. Prolonged exposure leads to further adaptation by 'group selection' of bacterial communities with higher detox per bacterium. These findings show that Lamarckian adaptation can arise via interactions between two levels of Darwinian selection within a holobiont system. The conclusions and modelling framework are applicable to diverse types of holobiont systems.


Introduction
Evolutionary adaptations are commonly thought to be driven by genetic mutations occurring on a timescale of many generations. Selection of individuals with rare beneficial mutations can then enable the adaptation of a large population without the emergence of new variations within a single lifetime. This view has recently been brought into question (1-7) by evidence for non-Mendelian epigenetic phenomena (8)(9)(10)(11), genome editing and mobility systems (12,13), niche construction (14) and contribution of symbiotic micro-organisms to heritable variation of the whole unit (holobiont) (3,5,6,15,16). The latter may be of particular interest because it violates the fundamental Neo-Darwinian assumption of a single level of selection per individual. Considering the individual as an interacting community with multiple levels of selection and generation times may have transformative implications on the dynamics and outcomes of evolution (3,5,6). The genetic content and species composition of symbiotic bacteria can change substantially within a single generation of the host and are inherited by vertical and/or horizontal transmission of bacteria (3,5,6,(17)(18)(19). Since the symbiotic microbiome is an integral part of host development and physiology (16), variations in the microbiome can influence the state and phenotypes of the host. The resulting changes are not necessarily limited to somatic tissue of the host. They can also extend to the host germline (20) and can contribute to non-Mendelian inheritance of environmentally-induced phenotypes (17). In line with these considerations, it has been suggested that rapid diversification and transmission of symbiotic bacteria may contribute to rapid evolution of their host (3,5,6,15,21). This possibility, however, has not yet been functionally confirmed by direct experimentation and the modes of adaptation that are available for a population of holobiont communities have not been investigated in a mathematical framework. Current models of evolutionary change are typically limited to population genetics with a single species of individuals (22) or to ecological (23)(24)(25) and evolutionary game theory models (26) (27) of interactions between free-living or symbiotic (28,29) species. However, these models do not consider the evolution of a holobiont population in which every individual is an interacting community of host and bacteria.
Here, we construct a general modeling framework of host-symbiont evolution under exposure to toxic stress. The holobiont in this 'Population Symbio-Genetics' (PSG) model is considered as a single unit undergoing two levels of Darwinian selection occurring on different timescales for host and bacteria. We show that this structure can exhibit Lamarckian offspring adaptation due to stress-induced during the lifetime of the parental holobiont. The Lamarckian adaptation is enabled by selection and transfer of resistant bacteria and is mediated by context-dependent alleviations of distinct types of stress (toxic and physiologic). Beyond this Lamarckian effect we show that Darwinian selection over timescales larger than a host generation, promote 'group selection' of bacterial communities with higher average detoxification per bacterium (in contrast to selection of better fit individual bacteria which takes place within a host generation). This form of group selection enhances the adaptation independently of the Lamarckian effect and is accompanied by stress-dependent variability across holobionts (despite the fixed mutation rate).

Formulation of the PSG Model
The model considers a population of holobionts under exposure to active toxin of concentration T.
Interactions between each host and its bacterial symbionts 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 by the stress state of the host.
For each host and bacterium, we define the probability of survival to reproduction (PH and PB, respectively), as follows:

(2) PB = (1 -NB /2KB) exp(-SB)
Here, NH is the number of hosts and KH is the maximal number of hosts that can be supported by the external environment (carrying capacity for hosts). KB is the number of bacteria that can be accommodated in the host (carrying capacity for bacteria) and NB is the number of bacteria per host.
Representing host-and bacterial-specific sensitivities to toxin by xH and xB, respectively, we define the instantaneous toxic stress to host and bacteria as SH = xH T and SB = xB T. Ŝ H = <SH>t represents the average of SH over a host generation time (interval between host reproduction events). Additionally, Ŝ Ph (NB) = ln(<NB > t /KB 0 ) + (1 -<NB > t /KB 0 ) corresponds to a 'physiological stress' over a host generation due to deviations from the bacterial population size that is preferred in toxin-free conditions. This size is determined by KB 0 , which also provides an inverse scale (1/KB 0 ) for the negative impacts of losing bacterial-derived metabolites (when NB < KB 0 ) or having to support excess numbers of bacteria (NB > KB 0 ) (30). Accordingly, the physiological stress vanishes when NB =KB 0 . The averages over a host generation (i.e. <SH> t and <NB> t) are introduced as practical approximations, enabling us to simplify the model by evaluating the survival of the host only at the time of reproduction.
Modulation of the bacterial niche by the state of the host (31)(32)(33)(34) is modelled by changing KB as a function of the instantaneous host stress: where  is a host-intrinsic property, determining the direction and extent of the change in the bacterial niche size due to the stress of the host. 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 an unbiased analysis of how KB might evolve under different exposures to toxin, we consider a starting population of hosts with an 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 in a holobiont-specific manner by release of detox from the host and each of its bacteria (30,(34)(35)(36):

) T( t+Δt) = T(t) exp(-λB ∑yB -λH yH) + θ(t)
where yH and yB are the amounts of detox secreted inside a given holobiont at time t, and λH and λB are constant detoxification capacities. While all the individuals within a given holobiont equally benefit from the total amount of secreted detox, they have individual-specific sensitivities to toxin.
The evolving traits of the host and its bacteria (x, y and ) are initially drawn from trait-specific distributions. Surviving bacteria divide at every time step of the simulation (Δt) and living hosts reproduce every  generations of bacteria (so that 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, µ (Supplementary Methods).
Here z corresponds to any of the evolving traits x, y and ,  is a standard Gaussian deviate with zero mean, and z0 and z are trait-specific coefficients, controlling the peak and width of the steady state distributions (specified in Supplementary information). Values of y and  were chosen to support wide distributions of y and , respectively. For the sensitivity traits (z = xH and xB), the distribution is truncated at x = 0, so as to prevent a trivial solution in which all the individuals are completely insensitive to toxin. We 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). In this study, we consider an initial population of 32000 hosts (NH = KH = 32000) with 100 bacteria per host (NB = KB 0 = 100).
We set the host generation time to  = 100 bacterial generations and all mutation rates to µ = 10 -3 per generation (for both host and bacteria).

Stress-dependent adjustment of bacterial niche size
We first 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 (Fig. 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 coupled to a live host, the survival of the holobiont depends also on the amount of detox secreted by the bacterial community (Fig. 1C). This secretion is higher for hosts which react to the toxic stress by increasing their carrying capacity for bacteria (i.e. hosts with  > 0; Supplementary Fig. S1A).
This leads to stress-dependent selection of hosts which provide a larger bacterial niche, KB (

Stress-dependent Lamarckian adaptation
Bacteria that are modified by the stress in one host generation can be transmitted to the next generation, potentially increasing the stress tolerance of the offspring. To quantitatively evaluate this possibility we introduce a new measure, termed the "Lamarckian", representing the gain in offspring tolerance due to acquisition of traits by its parent. To evaluate the contribution of acquired (as opposed to initial) traits, we identify the subpopulation of hosts that survived the first generation of exposure and apply a new simulation to these parents at their initial state ("cloned parents") and to their offspring ( Fig. 2A). We then compare the survival rates of the offspring (SROffs) to that of their "cloned parents" (SRCP). The use of a parental subpopulation with the original state of microbiota allows us to distinguish increase of tolerance due to selection of initially better fit parents from a gain of tolerance due to transmission of changes acquired during a host generation (not present in the initial parental clones). We then define the Lamarckian, L, as: (6) so that it is positive if the average tolerance increases due to transfer of changes acquired during a host generation.
For a given λB, we find that L is an increasing function of the injected amount of toxin and vanishes at low T0 (Fig. 2B). Analysis of L for different choices of λB (at a given T0) reveals a non-monotonic dependence on λB, manifested by an essentially 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 transferring bacterial population which acquired toxin resistance during the parental host generation (Fig. 2D). To determine how this bacterial transfer confers a gain of toxin tolerance in the hosts' offspring, we compared the offspring and their cloned parents with respect to the toxic and physiologic stress under exposure to the toxin. For small enough λB, the reduction of toxic stress by bacteria is negligible and the positive Lamarckian is primarily due to alleviation of the physiological stress in the offspring (Supplementary Fig. S2A). At intermediate values of λB, the alleviation of bacterial loss enhances the neutralization of toxin, providing an additional contribution to the Lamarckian (Fig. 2E,F).
However, when λB is large enough to support substantial neutralization of toxin during a single host generation ( Supplementary Fig. S3), the Lamarckian decreases due to the reduction of stress in both the parents and the offspring (Supplementary Fig. S2B).

Stress-dependent group selection and 'Bacterial Assimilation'
When the toxic pressure persists over timescales larger than one host generation (Fig. 3A), the selection favors holobionts whose bacterial communities have higher average levels of detox per bacterium, <yB> ( Fig. 3B). This 'group selection' is based on a collective property of the bacterial community (as opposed to selection of individual bacteria). In the current model, this 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 (small enough y), the selection is accompanied by significant accumulation of bacterial mutations. This enhances the group selection for higher <yB>, thus leading to toxin-dependent increase in detoxification rate (Fig. 3A, inset) and expedited holobiont adaptation (Fig.   3C). Group selection for higher detox is also accompanied by extended persistence of high detox levels (Fig.   3B) and by elevated detox variability across holobionts (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 ( Supplementary Fig. S4A,B). Bacterial mutations that were selected under stress persist over a characteristic timescale of 1/µ = 10 host generations after neutralization of the toxin, thus providing 'memory' of the previous exposure. To evaluate the influence of this 'memory' on tolerance to new exposures, we analyzed the response to repeated pulses of injected toxin, separated by time intervals shorter than 10 host generations. The reexposures lead to successive selections which oppose the relaxation of <yB> to its (lower) equilibrium value ( Fig. 4A vs. Fig. 3B), resulting in enhanced detoxification (Fig. 4B) and reduced impact on the host population (Fig. 4C). This enhancement reduces the environmental pressure and enables the survival of intrinsically less resistant hosts and bacteria (Fig. 4D). The progressive increase in host sensitivity due to successive selections for higher detox per bacteria, <yB>, demonstrates a hitherto unrealized mode of assimilation. This mode of 'Bacterial Assimilation' is analogous to 'Genetic Assimilation' due to successive selections of host-intrinsic alleles (37,38). However, assimilation by bacteria might be more effective because of the faster generation of a larger reservoir of variations (variations are generated by many bacteria per host on time scales that are much shorter compared to generation time of new genetic alleles in the host). This advantage of bacterial assimilation is particularly critical when the host population is small, or alternatively, when the repertoire of host-intrinsic alleles has been trimmed by previous selections. The absence of faster group selection reflects a lack of mechanism (in our model) for changing <yB> during a single host generation (with the possible exception of rare cases of rapid changes in <yB> due to amplification of very small numbers of resistant bacteria). This limitation can be removed by allowing the stress of the host to influence bacterial phenotypes (e.g. by subjecting y to stress-dependent dynamics similar to that of KB). This scenario is, in fact, expected in every holobiont due to the numerous options for 2-way interactions between the host and its symbionts. Such an influence could make a substantial The Lamarckian was evaluated by reverting a subset of holobionts to their initial state and re-subjecting them to toxin. Since this procedure cannot be realized experimentally, the Lamarckian of real holobionts should be approximated by other context-dependent means. For example, when symbiotic bacteria can be removed without compromising the survival of their toxin-free hosts (e.g. flies and worms that develop from dechorionated eggs on good diet conditions), the Lamarckian can be approximated in steps that are conceptually similar to the simulation procedure. First, the holobionts are exposed to a challenge and their offspring are then 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. (17,20)), while the other is colonized by ('experienced') microbiota from a group of hosts which survived exposure to a challenge. The populations of hosts with naïve and experienced microbiota are then compared with respect to their ability to survive the challenge and the Lamarckian is approximated from the respective survival rates, L  SR Exp. bact / SR Naïve bact) -1. This approximation, however, neglects other types of changes that may have been acquired and transmitted to offspring (e.g. small RNAs (10), maternal RNA (39), persistent chromatin modifications (8), horizontal transfer of biochemical signals (40) or other modes of local niche construction (14), etc). An additional source of deviation arises when the transmitted change is picked up by bystander offspring, which did not inherit the acquired change from their own parents. The extent to which these additional modes influence the evaluation of the Lamarckian may be assessed with a different ordering of the steps in the above procedure. The modified ordering starts with removal of bacteria from two untreated subpopulations and continues with re-colonizing these subpopulations with 'naïve' and 'experienced' microbiota, respectively, exposing the colonized subpopulations to the challenging condition and re-estimating the Lamarckian from their differential survival. 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).
The above modelling framework does not aim to describe a particular holobiont, but rather to provide a general paradigm for considering adaptation to stress in a system with interactions between two levels of selection, each occurring on a different timescale. Combining this structure with pure Darwinian selection and constant rate of mutations gives rise to Lamarckian adaptation, group selection and stress-induced variability. These outcomes are not pre-specified but rather depend on the strength of toxic and physiological stress. Since the model relies on properties that are shared by diverse types of holobionts, the above conclusions are expected to be broadly applicable to a variety of holobionts and may also apply to other forms of equivalent organizations. Additionally, the modelling framework can be readily adjusted to consider additional factors, such as having multiple species of symbionts and pathogens (with interspecies competition and/or cooperation), asynchronous reproduction modes, epigenetic effects, ecological influences and transfer of bacteria (and/or toxin) between hosts.

Supplementary Methods:
The simulation starts with a population of hosts, each carrying a population of 100 bacteria. Host and bacterial properties (phenotypes) are initially drawn from defined distributions (steady state of Eq. 5 without toxin) with the parameters x0 = 0.25, x = 10, y0 = 0, y = 0.1, 0 = 0 and  = 0.1.
In every time step of the simulation (one bacterial generation), each bacterium reproduces if its survival probability (Eq. 2) is larger than a random number (between 0 and 1) drawn from a uniform distribution. Each of the surviving bacterial (parent) persists at its current state and gives rise to a modified bacterium (offspring), while dead bacteria are discarded. At the end of one host generation (100 time steps), the reproduction of hosts is determined based on the survival probability in Eq. 1.
Non-surviving hosts are discarded and each of the surviving hosts gives rise to a parent and offspring host as follows: -The parent retains its current state (sensitivity, detox, delta) and the state of its bacterial population.
-Following 99 bacterial generations, an offspring host is created with properties defined by Eq. 5.
Negative values of the sensitivity and detox are prevented by taking the absolute value of the outcome in Eq.1. Each offspring receives a copy of the bacterial population of its parent. These populations are then iterated forward one bacterial generation, the surviving bacteria reproduce so as to define the initial state of the bacterial populations in the next host generation of the parent and its offspring.