- Email: [email protected]
Experimental evolution in silico: a custom-designed mathematical model for virulence evolution of Bacillus thuringiensis
Accepted Manuscript Title: Experimental evolution in silico: a custom-designed mathematical model for virulence evolution of Bacillus thuringiensis Author: Jakob Friedrich Strauß Philip Crain Hinrich Schulenburg Arndt Telschow PII: DOI: Reference:
S0944-2006(16)30017-4 http://dx.doi.org/doi:10.1016/j.zool.2016.03.005 ZOOL 25495
To appear in: Received date: Revised date: Accepted date:
30-11-2015 20-2-2016 17-3-2016
Please cite this article as: Strauss, Jakob Friedrich, Crain, Philip, Schulenburg, Hinrich, Telschow, Arndt, Experimental evolution in silico: a custom-designed mathematical model for virulence evolution of Bacillus thuringiensis.Zoology http://dx.doi.org/10.1016/j.zool.2016.03.005 This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.
Experimental evolution in silico: a custom-designed mathematical model for virulence evolution of Bacillus thuringiensis Jakob Friedrich Straußa, Philip Craina,b, Hinrich Schulenburgc, Arndt Telschowa,* a
Institute of Evolution and Biodiversity, Westfälische Wilhelms-Universität,
Hüfferstraße 1, D-48149 Münster, Germany b
DuPont Pioneer, 200 Powder Mill Rd, Wilmington, DE 19803, USA
c
Department of Evolutionary Ecology and Genetics, Christian-Albrechts-Universität
zu Kiel, Am Botanischen Garten 1-9, D-24118 Kiel, Germany
Corresponding author. E-mail address: [email protected] (A. Telschow).
1
Highlights
There is a theory–experiment mismatch regarding the evolution of virulence. Custom-designed mathematical models help analyse the results from experimental evolution. We propose a model customised to the virulence evolution of Bacillus thuringiensis. Virulence evolution is explained by context-dependent cost of toxin production. The model explains empirical results and makes testable predictions for future experiments.
Abstract Most mathematical models on the evolution of virulence are based on epidemiological models that assume parasite transmission follows the mass action principle. In experimental evolution, however, mass action is often violated due to controlled infection protocols. This “theory–experiment mismatch” raises the question whether there is a need for new mathematical models to accommodate the particular characteristics of experimental evolution. Here, we explore the experimental evolution model system of Bacillus thuringiensis as a parasite and Caenorhabditis elegans as a host. Recent experimental studies with strict control of parasite transmission revealed that one-sided adaptation of B. thuringiensis with non-evolving hosts selects for intermediate or no virulence, sometimes coupled with parasite extinction. In contrast, host–parasite coevolution selects for high virulence and for hosts with strong resistance against B. thuringiensis. In order to explain the empirical results, we propose a new mathematical model that mimics the basic experimental set-up. The key assumptions are: (i) controlled parasite transmission (no mass action), (ii) discrete host generations, and (iii) context-dependent cost of toxin production. Our model analysis revealed the same basic trends as found in the experiments. Especially, we could show that resistant hosts select for highly virulent bacterial strains. Moreover, we found (i) that the evolved level of virulence is independent of the initial level of virulence, and (ii) that the average amount of bacteria ingested significantly affects the evolution of virulence with fewer bacteria ingested selecting for highly virulent strains. These predictions can be tested in future experiments. This study highlights the usefulness of custom-designed mathematical models in the analysis and interpretation of empirical results from experimental evolution. Keywords: Bacillus thuringiensis; Evolution of virulence; Mathematical modelling; Experimental evolution; Mass action principle 2
1. Introduction Mathematical models are important tools that help define and elucidate hypotheses and theories regarding evolutionary processes (Otto and Day, 2007). Experimental evolution tests such hypotheses through controlled experiments, in which one or few variables are allowed to vary while all others are held constant (Kawecki et al., 2012). There is a growing number of studies in experimental evolution, ranging from work on the evolution of virulence (Bedhomme et al., 2012; Kubinak et al., 2012) to stress tolerance (Bettencourt et al., 1999; Crecy et al., 2009), mating behaviour (Leu and Murray, 2006), and host–parasite coevolution (Buckling et al., 2006; Paterson et al., 2010; Schulte et al., 2010; Masri et al., 2015). Yet, there is potential for misinterpretation of empirical results due to a “theory–experiment mismatch”. Mathematical models were often developed to describe natural processes, not experimental work. In order to overcome this mismatch, there is a need for (new) models that capture the essential structure of the experimental design (Day, 2001). We illustrate the “theory–experiment mismatch” by contrasting the results from models for virulence evolution with empirical results from experimental evolution. A key assumption in most host–parasite coevolution models is that the probability that an uninfected individual becomes infected is proportional to the number of infected individuals in a population, i.e., the mass action principle. This is an important assumption because transmission is critical for parasite spread, persistence and evolution. A classic example is the evolution of virulence in Myxoma viruses after their release into Australian rabbit populations for pest control (Anderson and May, 1982). An epidemiological model (Kermack and McKendrick, 1933; Baily, 1975) was modified to include a trade-off between parasite transmission and parasite virulence, and the reduction in virulence was explained as an optimal response of the parasite (for details see Anderson and May, 1982). Most of the present models on virulence evolution have a structure similar to that of epidemiological models (e.g., Nowak and May, 1994; Hethcote, 2000; Berngruber et al., 2013), and thus assume that parasite transmission follows the mass action principle. Parasite transmission and the different transmission mechanisms are known to be crucial for virulence evolution (e.g., Rafaluk et al., 2015b). As a consequence, transmission of parasites is strictly controlled in many experimental evolution studies (e.g., Schulte et al., 2010; Masri et al., 2015; but see Rafaluk et al., 2015a or Poullain et al., 2007 for examples without transmission control). Masri et al. (2015) performed 3
coevolution experiments with Bacillus thuringiensis (BT) as a parasite and Caenorhabditis elegans (CE) as a host. Parasite transmission was precisely controlled, in order to compare evolutionary scenarios where either both antagonists were allowed to reciprocally adapt to each other (i.e., coevolution) or where only one of the two was allowed to adapt, while the other was always re-introduced into the experiment from a stock culture (i.e., one-sided adaptation). One-sided adaptation of the parasite with non-evolving hosts favoured bacteria with intermediate or no virulence, sometimes coupled with parasite extinction. In contrast, host–parasite coevolution selected for bacteria with high virulence and for hosts with strong resistance against the bacteria. In these experiments, host generations were discrete and parasite transmission was controlled such that every new host generation was exposed to the same number of bacterial spores. Spores were obtained from the dead host bodies of the previous selection round (see Masri et al., 2015 for details). Obviously, parasite transmission in this experiment does not follow the mass action principle. Thus, most standard models for virulence evolution are unlikely to capture the dynamics in such experimental studies. In the experiment, BT is transmitted via spores between host generations. BT spore formation is likely to start after nutrition available from the host is depleted, which follows CE death and can thus be related to the time point of host mortality. In addition to spores, some BT produce Cry toxin crystals. Cry toxins damage host tissue, which ultimately provides access for the pathogens to previously indigestible nutrients and also ultimately kills the host (Agaisse and Lereclus, 1995). The presence and abundance of these toxins are the key determinants of BT virulence. Cry toxins are encoded on plasmids (Gonzales and Carlton, 1980), and toxinproducing BT carry these plasmids while non-toxin-producing BT do not. The crystals stay attached to the spore and are likely to provide a benefit in local competition when the spores germinate in the host gut (Cry toxins become active as soon as they are solubilised in the host gut). However, maintaining virulence, i.e. replicating plasmids as well as producing the toxin during spore formation, is likely costly and may reduce bacterial growth rates. While toxin-producing BT benefit while the host is still alive, non-toxin-producing BT can be assumed to benefit when the host is dead (for a review on BT biology see Nielsen-LeRoux et al., 2012). Here, we propose a new mathematical model to describe the evolution of virulence in the Masri experiment. Key assumptions are (i) controlled parasite transmission, (ii) 4
discrete host generations, and (iii) costly toxin production. We analyse two scenarios. First, we assume that toxin production reduces growth of the producing bacteria, but that the benefit of toxin is equally shared among toxin and non-toxin producers (scenario 1: context-independent cost of toxin production). Second, we assume that toxin producers have a benefit over non-toxin producers as long as the host is alive, but that non-toxin producers grow faster in the dead host (scenario 2: contextdependent cost of toxin production). The theoretical analysis revealed that only scenario 2 is able explain the experiments of Masri et al. (2015).
2. Materials and methods The mathematical model captures the basic structure of the one-sided parasite adaptation treatment as described in Masri et al. (2015). It consists of two basic parts: (1) the intra-host infection dynamics of two bacterial strains, and (2) the evolutionary dynamics including controlled bacterial transmission between hosts.
2.1. Intra-host dynamics A discrete time step model describes the intra-host infection dynamics. There are four variables. These are the number of toxin-producing bacteria non-toxin-producing bacteria amount of digestible nutrient indicated by ‘ ’ and ‘
, the number of
, the amount of indigestible nutrient
, and the
. The variables of subsequent time steps are
1’. Within each time step, the following events occur in
sequential order: (i) nutrient–toxin interaction, (ii) nutrient uptake and bacterial growth, and (iii) bacterial mortality. Intermediate values of the variables are indicated by + (after step i) and ++ (after step ii). The initial conditions are as follows. At the beginning of infection, the host is infected with a certain number of bacteria, the bacterial load b. At this point, there are toxin producers, and The amount of toxins
0
non-toxin producers, and it holds that
0 .
ingested by the host is purely determined by the number of
toxin-producing bacteria at the beginning of infection, i.e.,
0
. Note that
the amount of toxins does not change during the whole process of infection. Finally, the host consists at time step zero of 0
0
units of indigestible nutrients and
0 units of digestible nutrients.
The model consists of two phases that may differ with respect to the bacterial growth rates. Phase 1 starts with the bacterial infection of the host and lasts either until host 5
death or until a maximal number of time steps the fraction
is reached. Host death occurs when
of the indigestible nutrient nu has been converted into digestible
nutrient. The parameter
reflects the ability of the host to “resist” toxin-facilitated
nutrient digestion and is called level of host resistance. Phase 2 starts with host death and ends after
time steps. There is no second phase if the host does not die
in phase 1.
2.1.1. Phase 1: host alive (i) Nutrient–toxin interaction. Each unit of toxin converts one unit of indigestible nutrient into one unit of digestible nutrient. If the amount of indigestible nutrition exceeds the amount of toxin or is equal to it (1)
,
(2)
.
If
, then it holds that
, then it holds that 0,
(3) (4)
. and
The number of bacteria does not change during this step, i.e. .
(ii) Nutrient uptake and bacterial growth. The growth rates of toxin producers and non-toxin producers are denoted by
and
, respectively. Two cases need to be
distinguished. First, there is no food limitation inside the body if . Then, bacterial numbers and nutrient change to (5)
1
,
(6)
1
,
(7)
.
Here, c is the conversion rate from digestible nutrient into bacteria. Second, food limitation occurs if
. Then, digestible nutrient and
bacterial numbers after consumption and growth compute to (8)
1
,
(9)
1
,
(10)
0.
Indigestible nutrient does not change during this step, i.e.
. 6
(iii) Bacterial mortality. A fraction
of each type of bacteria dies. This step completes
the time steps. It holds that (11)
1
1
,
(12)
1
1
,
(13)
1
,
(14)
1
.
2.1.2. Phase 2: dead host The first phase ends with host death, i.e. when dnu units of indigestible nutrient have been converted into digestible nutrient. Note that high levels of host resistance (i.e., high values of
) result in dynamics with a prolonged first phase. The infection
dynamics of phase 2 are described by equations (1)–(14), but with the modification that bacterial growth rates change from
to
and from
to
. An important
aspect of phase 2 is that all nutrients are digestible. Therefore, the nutrient variables 0 and
are set to
1
, where
denotes the amount of
digestible nutrients at the end of phase 1.
2.2. Evolutionary dynamics The evolutionary model describes how bacterial frequencies change between subsequent selection rounds. For simplicity, we assume a uniform host population, i.e. all hosts in subsequent generations are characterised by the same parameter set (including
and
) and the infection dynamics in all hosts follows model (1)–(14).
Let us assume that hosts in the previous generation are infected by producers and
0
non-toxin producers (
bacterial numbers have changed to
and
). After
0
toxin
time steps, the
. Inter-host transmission
happens in the following way. Each host in the new generation is infected by bacteria. These consist of
toxin producers and
non-toxin
producers. Note that the number with which a new host is infected is kept constant at in the model, but the bacterial frequencies can vary. For the analysis, we iterated this process times, and plotted the frequency of toxin producers as a function of and for different parameter constellations. The parameter represents the number of selection rounds in Masri et al. (2015).
7
3. Results The model was used to investigate two biologically feasible scenarios of the evolutionary experiments. Scenario 1 assumes context-independent costs of toxin production (see Section 3.1.), while in scenario 2 toxin costs are context-dependent (see Section 3.2). Though both scenarios can explain the results of the one-sided parasite adaptation (OSPA) treatment in Masri et al. (2015) (i.e., loss of virulence, parasite extinction), only the second scenario captures the dynamics observed during host–parasite coevolution (i.e., increased virulence and resistance).
3.1. Scenario 1: context-independent cost of toxin production We assume that toxin-producing bacteria have a lower growth rate than non-toxin producers in both alive and dead hosts, i.e.
and
. The reduced
growth rates may reflect metabolic costs or delayed cell division, possibly caused by replication of plasmids involved in toxicity (Gonzales and Carlton, 1980) or toxin synthesis (Agaisse and Lereclus, 1995; Nielsen-LeRoux et al., 2012). We found that non-toxin producers outcompete toxin producers in all modelled scenarios. This result is expected because toxin production involves a cost and the benefit of toxins is shared equally among bacteria. Fig. 1 shows a typical example of the evolutionary dynamics. In every selection round, toxin producers decrease in frequency, which can be interpreted as evolution towards reduced virulence. However, this process comes to an end if toxin-producer frequency is below a certain critical threshold. If there are too few toxin-producing bacteria present in the initial infection to kill new hosts, then all hosts survive infection and the selection experiment comes to an end. Note that this qualitative result does not depend on the parameter choice. In particular, increase in virulence was not found for varying host resistance levels (results not shown). In conclusion, scenario 1 offers an explanation for the OSPA treatment, but not for host–parasite coevolution; it predicts the decrease or loss of virulence, but fails to explain selection for high virulence and how this relates to the observed increased levels of host resistance.
3.2. Scenario 2: context-dependent cost of toxin production Next, we analysed the model for context-dependent cost of toxin production. We assume that toxin-producing bacteria have a relative advantage in growth as long as the host is alive (phase 1), but have a disadvantage once the host is dead (phase 2). 8
This is modelled with growth rates for which
holds true.
Empirical studies support the view of context-dependent growth rates. In general, local access to nutrients is thought to be facilitated in alive hosts by the local presence of Cry toxins, which may not diffuse equally across space and are likely to remain in the vicinity of the toxin-producing cells (for a review on the biological mechanism see Raymond et al., 2010 and Nielsen-LeRoux et al., 2012). This suggests faster growth of toxin-producing bacteria in phase 1. However, once the host is dead, we assume that access to nutritional resources is not restricted. Because toxin producers may have a lower maximal growth rate due to plasmid replication or other factors (e.g., toxin synthesis), they are likely to be disfavoured relative to non-producing cells in dead hosts. Scenario 2 differs in comparison to scenario 1 because there is negative frequencydependent selection of virulence. The time until host death in the model is mainly determined by the amount of toxin with which the host is initially infected. Therefore, high frequencies of toxin producers result in rapid host death, intermediate frequencies in comparably late death, and if toxin producer frequencies are below a certain threshold, the host survives bacterial infection and phase 2 never occurs. These qualitative results are true for both scenario 1 and scenario 2; however, in scenario 2, the phase length determines the relative advantage of toxin producers over non-toxin producers. Negative frequency-dependent selection occurs because a high frequency of toxin producers results in quick death of the host (i.e., shorter phase 1) and toxin producers have a relative fitness disadvantage during the prolonged second phase. Conversely, low frequencies of toxin producers result in the slow death of the host and increased selective advantage. Fig. 2 illustrates the evolutionary dynamics for scenario 2 (Fig. 2D), and the intra-host dynamics in selection rounds 1, 17 and 60 (Fig. 2A, B and C, respectively). Given our benchmark parameter values, the frequency of toxin-producing bacteria increases with the number of selection rounds (Fig. 2D). In selection round 1 (Fig. 2A), phase 1 lasts for eight bacterial generations. After 17 selection rounds (Fig. 2B), phase 1 is shortened to seven generations. Finally, at selection round 60 (Fig. 2C), phase 1 is only six generations. The decrease in the length of phase 1 represents increased virulence of the pathogen, since the host lives for fewer bacterial generations relative to hosts in earlier selection rounds.
9
An important implication of the negative frequency-dependent selection of virulence is that the system can evolve to an equilibrium level of virulence. Fig. 3A illustrates this point: for the same parameter set, virulence increases from starting conditions of 11% and 25% of toxin-producing bacteria and decreases from starting conditions of 50%, 89% and 99%, with the system always evolving to the same equilibrium near 49%. Such an equilibrium is not reached for initial conditions of 10% or lower because the amount of toxin is not sufficient to kill the host, and the evolution experiment ends after the first selection round (indicated by an asterisk in Fig. 3A). Masri et al. (2015) reported reduced virulence and pathogen extinction under OSPA experimental conditions, but high virulence when host and pathogen were allowed to coevolve. The model for scenario 2 demonstrates both empirical observations theoretically. Fig. 3B shows that the frequency of toxin-producing bacteria is affected by host resistance. Increasing the level of host resistance selects for more virulent pathogens because toxin producers benefit from the prolonged phase 1. This finding offers a simple explanation for the coevolution experiment under the assumption that hosts evolve towards higher levels of resistance. Simulations with low host resistance can explain results from OSPA. For d = 0.1, virulence decreases in every selection round until it is too low to kill new hosts, thus ending the experiment (indicated by an asterisk). We further investigated the effect of bacterial load (i.e., the amount of bacteria that each host ingests in every selection round) on the evolution of virulence. Fig. 3C demonstrates that high bacterial load leads to low levels of virulence and low bacterial load to high levels of virulence. This effect can be understood against the background of frequency-dependent selection. In Fig. 3C, we vary the bacterial load, but start in all cases with 50% toxin producers. The bacterial load determines the amount of toxin, which is derived from the absolute number of initial toxin producers, . For example, bacterial loads of 100 and 200 result in 50 and 100 units of toxin, respectively. As a consequence, high bacterial load leads to a relatively shorter phase 1 and low bacterial load to a relatively longer phase 1, because more toxin results in faster host death. If phase 1 is comparatively long, this favours the toxin producers, subsequently leading to an overall higher virulence of the parasite population. Note that the length of phase 1 can be calculated analytically and computes to the smallest natural number larger than
/
.
10
Next, we investigated the effect of host body size on the evolution of virulence. Fig. 3D demonstrates that large body size leads to high levels of virulence and small body size to low levels of virulence. Analogous to Fig. 3C, the effect in Fig. 3D can be understood against the background of frequency-dependent selection. In Fig. 3D, we vary the host body size while the bacterial load stays constant. The body size determines the total amount of nutrition available, split by the host’s level of resistance
into
for phase 1 and 1
for phase 2. Increasing
increases the length of phase 1. It is important to note that in the model the toxin acts as a catalyst, providing a constant amount of nutrition to the bacteria, and that the amount of toxin can never exceed the bacterial load. As a consequence, large body size leads to a long phase 1 and small body size to a short phase 1, especially because the bacterial load is not coupled to body size. The effects of host body size and host level of resistance are counter-affected by the bacterial load.
4. Discussion The model presented here provides a mechanistic hypothesis for BT virulence evolution and makes specific predictions for future one-sided parasite adaption experiments. First, the analysis suggests that the parasite evolves towards an equilibrium level of virulence due to frequency-dependent selection. This hypothesis can be tested experimentally by varying the initial level of virulence of the parasite population and carefully monitoring how virulence and time until host death change in the course of evolution. From the model we expect (i) that different starting pathogen populations evolve to the same level of virulence and (ii) that the level of virulence is positively correlated to the average time from infection until host death. Second, the mathematical analysis suggests that virulence evolution is strongly affected by the bacterial load (i.e., the number of bacteria with which hosts are infected in every selection round). For a controlled experiment where the bacterial load is varied, the model predicts that high (low) bacterial loads select for low (high) levels of virulence. Third, we hypothesise that large (small) body size selects for high (low) levels of virulence. Note, however, that body size should not be interpreted as the actual weight of the host, but as the amount of nutrient accessible for the pathogen if the host is fully exploited. This hypothesis may therefore be difficult to test because “body size” in the sense of the model may be difficult to control experimentally.
11
Our model assumes that pathogen replication inside the hosts solely depends on pathogen growth rates, available nutrients and the host level of resistance. In real systems, pathogen growth is a complex process that depends on various factors, including bacteriocin production (Favret and Yousten, 1989), protective symbionts (Zug and Hammerstein, 2015), and the host immune memory (Sadd and SchmidHempel, 2006; Roth et al., 2009, 2010). Though we did not explicitly model these factors, we can hypothesise how they affect virulence evolution based on our theoretical analysis. Bacteriocins are chemical compounds produced by certain bacterial strains to suppress growth of other bacterial strains. Some BT strains are well known to produce bacteriocins (Favret and Yousten, 1989). The effect of bacteriocins on virulence evolution can be modelled in our theoretical framework by changing the growth rates of the toxin and the non-toxin producers relative to each other. A possible outcome may be evolution towards higher virulence. We expect this to happen if bacteriocins impede the growth of toxin producers because this will increase the relative advantage of toxin production due to a prolonged phase 1. Similar effects may occur in experiments that involve protective symbionts (i.e., symbionts that increase host resistance) or immune memory (i.e., previous challenges with pathogens increase host resistance). In particular, we expect protective symbionts to select for higher virulence of BT if they impede toxin and nontoxin producers equally. Though the main focus of the present study was to explain specific results from experimental evolution, the analysis also offers a new perspective on BT evolution in natural ecosystems. Let us assume the virulence evolution of BT follows the twophase model analysed above, i.e. toxin producers have an advantage in living hosts and non-toxin producers in dead hosts. In a natural community, however, BT coexists with other microbes and may experience severe competition for resources. For example, saprophagous microbes may quickly colonise invertebrates killed by BT, reducing the resources available. Based on our theoretical analysis, we expect that this resource reduction will select for higher BT virulence because the non-toxinproducing BT suffer relatively more from increased resource competition than the toxin producers. We therefore hypothesise that communities with high abundance of saprophages select for highly virulent BT. In particular, we expect that host species that are prone to quick colonisation by saprophages, e.g., due to certain foraging
12
behaviours or the microbiome composition, may suffer from especially virulent BT strains. These hypotheses may be tested in future empirical studies. How do our results relate to previous theoretical work? A standard approach to analysing virulence evolution is to determine evolutionarily stable levels of virulence within an epidemiological framework, often under the assumption of a trade-off between pathogen transmission and virulence (Anderson and May, 1982; Ewald, 1983). Several authors stressed the importance of the transmission mode for virulence evolution (Ewald, 1983; Day, 2001). Day (2001) further pointed out that the transmission mode in experimental evolution may deviate from standard models and that mathematical models which are applied to selection experiments need to take into account key characteristics of the experimental protocol. In our study, we implemented the experimental protocol of the Masri experiments (Masri et al., 2015) by differentiating between inter-host dynamics and intra-host dynamics. Further, we chose a difference equation model to mimic the discrete host generations in the serial passage experiments. Thus, instead of trying to generate a new theory of virulence evolution we propose a rather basic model that is useful for experimentalists to test future experiments. Nevertheless, this model may serve as the basis for more sophisticated theoretical analyses. One promising direction may be to include stochasticity, which is surely relevant in the experiments due to the small population sizes (Papkou et al., 2016, current issue). Another direction may be to analyse an epidemiological version of the model, in which parasite transmission follows the mass action principle. However, the analysis of these extended models goes beyond the scope of the current study and is left for future research.
5. Conclusion Most mathematical models on the evolution of virulence are based on epidemiological models that assume parasite transmission follows the mass action principle. In experimental evolution, however, mass action is often violated due to controlled infection protocols. This “theory–experiment’ mismatch” was investigated by focusing on experimental results from virulence evolution (Masri et al., 2015) that are not explained by standard mathematical models. We developed a new model that is not based on the mass action principle for parasite transmission and explains the empirical results. The present study suggests that there is a need for custom-
13
designed mathematical models to capture the dynamics of virulence changes during evolution experiments.
Acknowledgements We thank Sultan Beshir for technical assistance, Joachim Kurtz, Heiko Liesegang, Rebecca Schulte, Gerrit Joop and Esther Sundermann for fruitful discussions on host–parasite coevolution, and Florian Grziwotz and two anonymous reviewers for helpful comments on the manuscript. This study was funded by the German Science Foundation (SPP 1399, TE 976/2-1).
14
References Anderson, R.M., May, R.M., 1982. Coevolution of hosts and parasites. Parasitology 85, 411–426. Agaisse, H., Lereclus, D., 1995. How does Bacillus thuringiensis produce so much insecticidal crystal protein? J. Bacteriol. 177, 6027–6032. Baily, N.J.T., 1975. The Mathematical Theory of Infectious Diseases and its Applications. 2nd ed., Macmillan, New York. Bedhomme, S., Lafforgue, G., Elena, S.F., 2012. Multihost experimental evolution of a plant RNA virus reveals local adaptation and host-specific mutations. Mol. Biol. Evol. 29, 1481–1492. Berngruber, T.W., Froissart, R., Choisy, M., Gandon, S., 2013. Evolution of virulence in emerging endemics. PLOS Pathog. 9, e1003209. Bettencourt, B.R., Feder, M.E., Cavicchi, S., 1999. Experimental evolution of HSP70 expression and thermotolerance in Drosophila melanogaster. Evolution 53, 484–492. Buckling, A., Wei, Y., Massey, R.C., Brockhurst, M.A., Hochberg, M.E., 2006. Antagonistic coevolution with parasites increases the cost of host deleterious mutations. Proc. R. Soc. Lond. B 273, 45–49. Crecy, E., Jaronski, S., Lyons, B., Lyons, T.J., Keyhani, N.O., 2009. Directed evolution of a filamentous fungus for thermotolerance. BMC Biotechnol. 9, 74. Day, T., 2001. Parasite transmission modes and the evolution of virulence. Evolution 55, 2389–2400. Ewald, P.W., 1983. Host–parasite relations, vectors, and the evolution of disease severity. Annu. Rev. Ecol. Syst. 14, 465–485. Favret, M.E., Yousten, A.A., 1989. Thuricin: the bacteriocin produced by Bacillus thuringiensis. J. Invertebr. Pathol. 53, 206–216. Gonzalez Jr., J.M., Carlton, B.C., 1980. Patterns of plasmid DNA in crystalliferous and acrystalliferous strains of Bacillus thuringiensis. Plasmid 3, 92–98. Hethcote, H. W., 2000. The mathematics of infectious diseases. SIAM Rev. 42, 599– 653. Kawecki, T.J., Lenski, R.E., Ebert, D., Hollis, B., Olivieri, I., Whitlock, M.C., 2012. Experimental evolution. Trends Ecol. Evol. 27, 547–560.
15
Kermack, W.O., McKendrick, A.G., 1933. Contributions to the mathematical theory of epidemics. III. Further studies of the problem of endemicity. Proc. R. Soc. Lond. A 843, 94–122. Kubinak, J.L., Ruff, J.S., Hyzer, C.W., Slev, P.R., Potts, W.K., 2012. Experimental viral evolution to specific host MHC genotypes reveals fitness and virulence trade offs in alternative MHC types. Proc. Natl. Acad. Sci. U.S.A. 109, 3422– 3427. Leu, J.Y., Murray, A.W., 2006. Experimental evolution of mating discrimination in budding yeast. Curr. Biol. 16, 280–286. Masri, L., Branca, A., Sheppard, A.E., Papkou, A., Laehnemann, D., Guenther, P.S., Prahl, S., Saebelfeld, M., Hollensteiner, J., Liesegang, H., Brzuszkiewicz, E., Daniel, R., Michiels, N.K., Schulte, R.D., Kurtz, J., Rosenstiel, P., Telschow, A., Bornberg-Bauer, E., Schulenburg, H., 2015. Host–pathogen coevolution: the selective advantage of Bacillus thuringiensis virulence and its cry toxin genes. PLOS Biol. 13, e1002169. Nielsen-LeRoux, C., Gaudriault, S., Ramarao, N., Lereclus, D., Givaudan, A., 2012. How the insect pathogen bacteria Bacillus thuringiensis and Xenorhabdus/Photorhabdus occupy their hosts. Curr. Opin. Microbiol. 15, 220–231. Nowak, M.A., May, R.M., 1994. Superinfection and the evolution of parasite virulence. Proc. R. Soc. B 255, 81–89. Otto, S.P., Day, T., 2007. A Biologist’s Guide to Mathematical Modeling in Ecology and Evolution. Princeton University Press, Princeton, NJ. Papkou, A., Gokhale, C.S., Traulsen, A., Schulenburg, H., 2016. Host–parasite coevolution: why changing population size matters. Zoology 119, ##–##. Paterson, S., Vogwill, T., Buckling, A., Benmayor, R., Spiers, A.J., Thomson, N.R., Quail, M., Smith, F., Walker, D., Libberton, B., Fenton, A., Hall, N., Brockhurst, M.A., 2010. Antagonistic coevolution accelerates molecular evolution. Nature 464, 275–278. Poullain, V., Gandon, S., Brockhurst, M.A., Buckling, A., Hochberg, M.E., 2007. The evolution of specificity in evolving and coevolving antagonistic interactions between a bacteria and its phage. Evolution 62, 1–11.
16
Rafaluk, C., Gildenhard, M., Mitschke, A., Telschow, A., Schulenburg, H., Joop, G., 2015a. Rapid evolution of virulence leading to host extinction under host– parasite coevolution. BMC Evol. Biol. 15, 112. Rafaluk, C., Jansen, G., Schulenburg, H., Joop, G., 2015b. When experimental selection for virulence leads to loss of virulence. Trends Parasitol. 31, 426– 434. Raymond, B., Johnston, P.R., Nielsen-LeRoux, C., Lereclus, D., Crickmore, N., 2010. Bacillus thuringiensis: an impotent pathogen? Trends Microbiol. 18, 189–194. Roth, O., Sadd, B.M., Schmid-Hempel, P., Kurtz, J., 2009. Strain-specific priming of resistance in the red flour beetle, Tribolium castaneum. Proc. Biol. Sci. 276, 145–151. Roth, O., Joop, G., Eggert, H., Hilbert, J., Daniel, J., Schmid-Hempel, P., Kurtz, J., 2010. Paternally derived immune priming for offspring in the red flour beetle, Tribolium castaneum. J. Anim. Ecol. 79, 403–413. Sadd, B.M., Schmid-Hempel, P., 2006. Insect immunity shows specificity in protection upon secondary pathogen exposure. Curr. Biol. 16, 1206–1210. Schulte, R.D., Makus, C., Hasert, B., Michiels, N.K., Schulenburg, H, 2010. Multiple reciprocal adaptations and rapid genetic change upon experimental coevolution of an animal host and its microbial parasite. Proc. Natl. Acad. Sci. U.S.A. 107, 7359–7364. Zug, R., Hammerstein, P., 2015. Bad guys turned nice? A critical assessment of Wolbachia mutualisms in arthropod hosts. Biol. Rev. Camb. Philos. Soc. 90, 89–111.
17
Figure captions Fig. 1. Evolutionary dynamics for scenario 1 (context-independent cost of toxin production). Shown is the frequency of toxin- (red) and non-toxin- (blue) producing bacteria at the end of each selection round. During every selection round, non-toxin producers increase in frequency. If toxin producer frequency decreases to a certain threshold, too few toxins are available among spores to kill new hosts. This results in the end of the selection experiment at selection round eight. Parameters: 1.8,
2,
1000,
0.001,
30,
10,
0.001,
30,
10.
Fig. 2. Intra-host dynamics for scenario 2 (context-dependent cost of toxin production). (A–C) Intra-host dynamics of bacteria for different selection rounds. There are two phases that differ with respect to the cost and benefit of toxin production. Phase 1 (green) lasts until host death and is characterised by faster growth of toxin producers, whereas in phase 2 (red) toxin producers are outcompeted by non-toxin producers. The context-dependent cost results in frequency-dependent selection on bacteria because phase length is determined by toxin-producer frequency. Changes in phase length are shown in (A–C). The full evolutionary 1.5,
dynamics of the selection experiment is shown in (D). Parameters: 1,
2,
1000,
0.34,
0.001,
30,
10.
Fig. 3. Evolutionary dynamics of scenario 2 (context-dependent cost of toxin production) showing how parameter variation affects the evolution of virulence (measured as toxin-producer frequency at the end of each selection round): variation of (A) initial bacterial virulence, (B) host resistance, (C) number of bacterial spores used to initiate infection in each selection round, (D) host body size. For some parameter constellations in (A and B), virulence evolves to such low levels that hosts survive the infection. This results in the end of the selection experiment and is indicated by an asterisk. The legends depict the varied parameters in the same sequence as the plot (top-down). Other parameters and initial conditions: (A) 2, 10; 30, 0.001, 0.3,
1.1,
3.07,
1000,
1.5,
(B) 10; (C) 30,
1,
1.2,
1.5,
10; (D)
0.001,
30,
1.5,
50, 2,
0.3,
0.001,
1000,
1, 1,
2, 2,
50, 1000, 1000,
30, 0.001,
0.525, 50,
10. 18
Figure 1
19
Figure 2
20
Figure 3
21
Table 1. Parameters and symbols of the mathematical model. number of toxin- (non-toxin-) producing bacteria units of indigestible (digestible) nutrients initial number of toxin- (non-toxin-) producing bacteria bacterial load: total number of initial bacteria (
)
initial units of indigestible (digestible) nutrients growth rate of toxin- (non-toxin-) producing bacteria in first phase growth rate of toxin- (non-toxin-) producing bacteria in second phase conversion coefficient from digestible nutrient into bacteria level of host resistance (0…1) bacterial mortality number of time steps (bacterial generations within a single host) number of selection rounds
22
S0944-2006(16)30017-4 http://dx.doi.org/doi:10.1016/j.zool.2016.03.005 ZOOL 25495
To appear in: Received date: Revised date: Accepted date:
30-11-2015 20-2-2016 17-3-2016
Please cite this article as: Strauss, Jakob Friedrich, Crain, Philip, Schulenburg, Hinrich, Telschow, Arndt, Experimental evolution in silico: a custom-designed mathematical model for virulence evolution of Bacillus thuringiensis.Zoology http://dx.doi.org/10.1016/j.zool.2016.03.005 This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.
Experimental evolution in silico: a custom-designed mathematical model for virulence evolution of Bacillus thuringiensis Jakob Friedrich Straußa, Philip Craina,b, Hinrich Schulenburgc, Arndt Telschowa,* a
Institute of Evolution and Biodiversity, Westfälische Wilhelms-Universität,
Hüfferstraße 1, D-48149 Münster, Germany b
DuPont Pioneer, 200 Powder Mill Rd, Wilmington, DE 19803, USA
c
Department of Evolutionary Ecology and Genetics, Christian-Albrechts-Universität
zu Kiel, Am Botanischen Garten 1-9, D-24118 Kiel, Germany
Corresponding author. E-mail address: [email protected] (A. Telschow).
1
Highlights
There is a theory–experiment mismatch regarding the evolution of virulence. Custom-designed mathematical models help analyse the results from experimental evolution. We propose a model customised to the virulence evolution of Bacillus thuringiensis. Virulence evolution is explained by context-dependent cost of toxin production. The model explains empirical results and makes testable predictions for future experiments.
Abstract Most mathematical models on the evolution of virulence are based on epidemiological models that assume parasite transmission follows the mass action principle. In experimental evolution, however, mass action is often violated due to controlled infection protocols. This “theory–experiment mismatch” raises the question whether there is a need for new mathematical models to accommodate the particular characteristics of experimental evolution. Here, we explore the experimental evolution model system of Bacillus thuringiensis as a parasite and Caenorhabditis elegans as a host. Recent experimental studies with strict control of parasite transmission revealed that one-sided adaptation of B. thuringiensis with non-evolving hosts selects for intermediate or no virulence, sometimes coupled with parasite extinction. In contrast, host–parasite coevolution selects for high virulence and for hosts with strong resistance against B. thuringiensis. In order to explain the empirical results, we propose a new mathematical model that mimics the basic experimental set-up. The key assumptions are: (i) controlled parasite transmission (no mass action), (ii) discrete host generations, and (iii) context-dependent cost of toxin production. Our model analysis revealed the same basic trends as found in the experiments. Especially, we could show that resistant hosts select for highly virulent bacterial strains. Moreover, we found (i) that the evolved level of virulence is independent of the initial level of virulence, and (ii) that the average amount of bacteria ingested significantly affects the evolution of virulence with fewer bacteria ingested selecting for highly virulent strains. These predictions can be tested in future experiments. This study highlights the usefulness of custom-designed mathematical models in the analysis and interpretation of empirical results from experimental evolution. Keywords: Bacillus thuringiensis; Evolution of virulence; Mathematical modelling; Experimental evolution; Mass action principle 2
1. Introduction Mathematical models are important tools that help define and elucidate hypotheses and theories regarding evolutionary processes (Otto and Day, 2007). Experimental evolution tests such hypotheses through controlled experiments, in which one or few variables are allowed to vary while all others are held constant (Kawecki et al., 2012). There is a growing number of studies in experimental evolution, ranging from work on the evolution of virulence (Bedhomme et al., 2012; Kubinak et al., 2012) to stress tolerance (Bettencourt et al., 1999; Crecy et al., 2009), mating behaviour (Leu and Murray, 2006), and host–parasite coevolution (Buckling et al., 2006; Paterson et al., 2010; Schulte et al., 2010; Masri et al., 2015). Yet, there is potential for misinterpretation of empirical results due to a “theory–experiment mismatch”. Mathematical models were often developed to describe natural processes, not experimental work. In order to overcome this mismatch, there is a need for (new) models that capture the essential structure of the experimental design (Day, 2001). We illustrate the “theory–experiment mismatch” by contrasting the results from models for virulence evolution with empirical results from experimental evolution. A key assumption in most host–parasite coevolution models is that the probability that an uninfected individual becomes infected is proportional to the number of infected individuals in a population, i.e., the mass action principle. This is an important assumption because transmission is critical for parasite spread, persistence and evolution. A classic example is the evolution of virulence in Myxoma viruses after their release into Australian rabbit populations for pest control (Anderson and May, 1982). An epidemiological model (Kermack and McKendrick, 1933; Baily, 1975) was modified to include a trade-off between parasite transmission and parasite virulence, and the reduction in virulence was explained as an optimal response of the parasite (for details see Anderson and May, 1982). Most of the present models on virulence evolution have a structure similar to that of epidemiological models (e.g., Nowak and May, 1994; Hethcote, 2000; Berngruber et al., 2013), and thus assume that parasite transmission follows the mass action principle. Parasite transmission and the different transmission mechanisms are known to be crucial for virulence evolution (e.g., Rafaluk et al., 2015b). As a consequence, transmission of parasites is strictly controlled in many experimental evolution studies (e.g., Schulte et al., 2010; Masri et al., 2015; but see Rafaluk et al., 2015a or Poullain et al., 2007 for examples without transmission control). Masri et al. (2015) performed 3
coevolution experiments with Bacillus thuringiensis (BT) as a parasite and Caenorhabditis elegans (CE) as a host. Parasite transmission was precisely controlled, in order to compare evolutionary scenarios where either both antagonists were allowed to reciprocally adapt to each other (i.e., coevolution) or where only one of the two was allowed to adapt, while the other was always re-introduced into the experiment from a stock culture (i.e., one-sided adaptation). One-sided adaptation of the parasite with non-evolving hosts favoured bacteria with intermediate or no virulence, sometimes coupled with parasite extinction. In contrast, host–parasite coevolution selected for bacteria with high virulence and for hosts with strong resistance against the bacteria. In these experiments, host generations were discrete and parasite transmission was controlled such that every new host generation was exposed to the same number of bacterial spores. Spores were obtained from the dead host bodies of the previous selection round (see Masri et al., 2015 for details). Obviously, parasite transmission in this experiment does not follow the mass action principle. Thus, most standard models for virulence evolution are unlikely to capture the dynamics in such experimental studies. In the experiment, BT is transmitted via spores between host generations. BT spore formation is likely to start after nutrition available from the host is depleted, which follows CE death and can thus be related to the time point of host mortality. In addition to spores, some BT produce Cry toxin crystals. Cry toxins damage host tissue, which ultimately provides access for the pathogens to previously indigestible nutrients and also ultimately kills the host (Agaisse and Lereclus, 1995). The presence and abundance of these toxins are the key determinants of BT virulence. Cry toxins are encoded on plasmids (Gonzales and Carlton, 1980), and toxinproducing BT carry these plasmids while non-toxin-producing BT do not. The crystals stay attached to the spore and are likely to provide a benefit in local competition when the spores germinate in the host gut (Cry toxins become active as soon as they are solubilised in the host gut). However, maintaining virulence, i.e. replicating plasmids as well as producing the toxin during spore formation, is likely costly and may reduce bacterial growth rates. While toxin-producing BT benefit while the host is still alive, non-toxin-producing BT can be assumed to benefit when the host is dead (for a review on BT biology see Nielsen-LeRoux et al., 2012). Here, we propose a new mathematical model to describe the evolution of virulence in the Masri experiment. Key assumptions are (i) controlled parasite transmission, (ii) 4
discrete host generations, and (iii) costly toxin production. We analyse two scenarios. First, we assume that toxin production reduces growth of the producing bacteria, but that the benefit of toxin is equally shared among toxin and non-toxin producers (scenario 1: context-independent cost of toxin production). Second, we assume that toxin producers have a benefit over non-toxin producers as long as the host is alive, but that non-toxin producers grow faster in the dead host (scenario 2: contextdependent cost of toxin production). The theoretical analysis revealed that only scenario 2 is able explain the experiments of Masri et al. (2015).
2. Materials and methods The mathematical model captures the basic structure of the one-sided parasite adaptation treatment as described in Masri et al. (2015). It consists of two basic parts: (1) the intra-host infection dynamics of two bacterial strains, and (2) the evolutionary dynamics including controlled bacterial transmission between hosts.
2.1. Intra-host dynamics A discrete time step model describes the intra-host infection dynamics. There are four variables. These are the number of toxin-producing bacteria non-toxin-producing bacteria amount of digestible nutrient indicated by ‘ ’ and ‘
, the number of
, the amount of indigestible nutrient
, and the
. The variables of subsequent time steps are
1’. Within each time step, the following events occur in
sequential order: (i) nutrient–toxin interaction, (ii) nutrient uptake and bacterial growth, and (iii) bacterial mortality. Intermediate values of the variables are indicated by + (after step i) and ++ (after step ii). The initial conditions are as follows. At the beginning of infection, the host is infected with a certain number of bacteria, the bacterial load b. At this point, there are toxin producers, and The amount of toxins
0
non-toxin producers, and it holds that
0 .
ingested by the host is purely determined by the number of
toxin-producing bacteria at the beginning of infection, i.e.,
0
. Note that
the amount of toxins does not change during the whole process of infection. Finally, the host consists at time step zero of 0
0
units of indigestible nutrients and
0 units of digestible nutrients.
The model consists of two phases that may differ with respect to the bacterial growth rates. Phase 1 starts with the bacterial infection of the host and lasts either until host 5
death or until a maximal number of time steps the fraction
is reached. Host death occurs when
of the indigestible nutrient nu has been converted into digestible
nutrient. The parameter
reflects the ability of the host to “resist” toxin-facilitated
nutrient digestion and is called level of host resistance. Phase 2 starts with host death and ends after
time steps. There is no second phase if the host does not die
in phase 1.
2.1.1. Phase 1: host alive (i) Nutrient–toxin interaction. Each unit of toxin converts one unit of indigestible nutrient into one unit of digestible nutrient. If the amount of indigestible nutrition exceeds the amount of toxin or is equal to it (1)
,
(2)
.
If
, then it holds that
, then it holds that 0,
(3) (4)
. and
The number of bacteria does not change during this step, i.e. .
(ii) Nutrient uptake and bacterial growth. The growth rates of toxin producers and non-toxin producers are denoted by
and
, respectively. Two cases need to be
distinguished. First, there is no food limitation inside the body if . Then, bacterial numbers and nutrient change to (5)
1
,
(6)
1
,
(7)
.
Here, c is the conversion rate from digestible nutrient into bacteria. Second, food limitation occurs if
. Then, digestible nutrient and
bacterial numbers after consumption and growth compute to (8)
1
,
(9)
1
,
(10)
0.
Indigestible nutrient does not change during this step, i.e.
. 6
(iii) Bacterial mortality. A fraction
of each type of bacteria dies. This step completes
the time steps. It holds that (11)
1
1
,
(12)
1
1
,
(13)
1
,
(14)
1
.
2.1.2. Phase 2: dead host The first phase ends with host death, i.e. when dnu units of indigestible nutrient have been converted into digestible nutrient. Note that high levels of host resistance (i.e., high values of
) result in dynamics with a prolonged first phase. The infection
dynamics of phase 2 are described by equations (1)–(14), but with the modification that bacterial growth rates change from
to
and from
to
. An important
aspect of phase 2 is that all nutrients are digestible. Therefore, the nutrient variables 0 and
are set to
1
, where
denotes the amount of
digestible nutrients at the end of phase 1.
2.2. Evolutionary dynamics The evolutionary model describes how bacterial frequencies change between subsequent selection rounds. For simplicity, we assume a uniform host population, i.e. all hosts in subsequent generations are characterised by the same parameter set (including
and
) and the infection dynamics in all hosts follows model (1)–(14).
Let us assume that hosts in the previous generation are infected by producers and
0
non-toxin producers (
bacterial numbers have changed to
and
). After
0
toxin
time steps, the
. Inter-host transmission
happens in the following way. Each host in the new generation is infected by bacteria. These consist of
toxin producers and
non-toxin
producers. Note that the number with which a new host is infected is kept constant at in the model, but the bacterial frequencies can vary. For the analysis, we iterated this process times, and plotted the frequency of toxin producers as a function of and for different parameter constellations. The parameter represents the number of selection rounds in Masri et al. (2015).
7
3. Results The model was used to investigate two biologically feasible scenarios of the evolutionary experiments. Scenario 1 assumes context-independent costs of toxin production (see Section 3.1.), while in scenario 2 toxin costs are context-dependent (see Section 3.2). Though both scenarios can explain the results of the one-sided parasite adaptation (OSPA) treatment in Masri et al. (2015) (i.e., loss of virulence, parasite extinction), only the second scenario captures the dynamics observed during host–parasite coevolution (i.e., increased virulence and resistance).
3.1. Scenario 1: context-independent cost of toxin production We assume that toxin-producing bacteria have a lower growth rate than non-toxin producers in both alive and dead hosts, i.e.
and
. The reduced
growth rates may reflect metabolic costs or delayed cell division, possibly caused by replication of plasmids involved in toxicity (Gonzales and Carlton, 1980) or toxin synthesis (Agaisse and Lereclus, 1995; Nielsen-LeRoux et al., 2012). We found that non-toxin producers outcompete toxin producers in all modelled scenarios. This result is expected because toxin production involves a cost and the benefit of toxins is shared equally among bacteria. Fig. 1 shows a typical example of the evolutionary dynamics. In every selection round, toxin producers decrease in frequency, which can be interpreted as evolution towards reduced virulence. However, this process comes to an end if toxin-producer frequency is below a certain critical threshold. If there are too few toxin-producing bacteria present in the initial infection to kill new hosts, then all hosts survive infection and the selection experiment comes to an end. Note that this qualitative result does not depend on the parameter choice. In particular, increase in virulence was not found for varying host resistance levels (results not shown). In conclusion, scenario 1 offers an explanation for the OSPA treatment, but not for host–parasite coevolution; it predicts the decrease or loss of virulence, but fails to explain selection for high virulence and how this relates to the observed increased levels of host resistance.
3.2. Scenario 2: context-dependent cost of toxin production Next, we analysed the model for context-dependent cost of toxin production. We assume that toxin-producing bacteria have a relative advantage in growth as long as the host is alive (phase 1), but have a disadvantage once the host is dead (phase 2). 8
This is modelled with growth rates for which
holds true.
Empirical studies support the view of context-dependent growth rates. In general, local access to nutrients is thought to be facilitated in alive hosts by the local presence of Cry toxins, which may not diffuse equally across space and are likely to remain in the vicinity of the toxin-producing cells (for a review on the biological mechanism see Raymond et al., 2010 and Nielsen-LeRoux et al., 2012). This suggests faster growth of toxin-producing bacteria in phase 1. However, once the host is dead, we assume that access to nutritional resources is not restricted. Because toxin producers may have a lower maximal growth rate due to plasmid replication or other factors (e.g., toxin synthesis), they are likely to be disfavoured relative to non-producing cells in dead hosts. Scenario 2 differs in comparison to scenario 1 because there is negative frequencydependent selection of virulence. The time until host death in the model is mainly determined by the amount of toxin with which the host is initially infected. Therefore, high frequencies of toxin producers result in rapid host death, intermediate frequencies in comparably late death, and if toxin producer frequencies are below a certain threshold, the host survives bacterial infection and phase 2 never occurs. These qualitative results are true for both scenario 1 and scenario 2; however, in scenario 2, the phase length determines the relative advantage of toxin producers over non-toxin producers. Negative frequency-dependent selection occurs because a high frequency of toxin producers results in quick death of the host (i.e., shorter phase 1) and toxin producers have a relative fitness disadvantage during the prolonged second phase. Conversely, low frequencies of toxin producers result in the slow death of the host and increased selective advantage. Fig. 2 illustrates the evolutionary dynamics for scenario 2 (Fig. 2D), and the intra-host dynamics in selection rounds 1, 17 and 60 (Fig. 2A, B and C, respectively). Given our benchmark parameter values, the frequency of toxin-producing bacteria increases with the number of selection rounds (Fig. 2D). In selection round 1 (Fig. 2A), phase 1 lasts for eight bacterial generations. After 17 selection rounds (Fig. 2B), phase 1 is shortened to seven generations. Finally, at selection round 60 (Fig. 2C), phase 1 is only six generations. The decrease in the length of phase 1 represents increased virulence of the pathogen, since the host lives for fewer bacterial generations relative to hosts in earlier selection rounds.
9
An important implication of the negative frequency-dependent selection of virulence is that the system can evolve to an equilibrium level of virulence. Fig. 3A illustrates this point: for the same parameter set, virulence increases from starting conditions of 11% and 25% of toxin-producing bacteria and decreases from starting conditions of 50%, 89% and 99%, with the system always evolving to the same equilibrium near 49%. Such an equilibrium is not reached for initial conditions of 10% or lower because the amount of toxin is not sufficient to kill the host, and the evolution experiment ends after the first selection round (indicated by an asterisk in Fig. 3A). Masri et al. (2015) reported reduced virulence and pathogen extinction under OSPA experimental conditions, but high virulence when host and pathogen were allowed to coevolve. The model for scenario 2 demonstrates both empirical observations theoretically. Fig. 3B shows that the frequency of toxin-producing bacteria is affected by host resistance. Increasing the level of host resistance selects for more virulent pathogens because toxin producers benefit from the prolonged phase 1. This finding offers a simple explanation for the coevolution experiment under the assumption that hosts evolve towards higher levels of resistance. Simulations with low host resistance can explain results from OSPA. For d = 0.1, virulence decreases in every selection round until it is too low to kill new hosts, thus ending the experiment (indicated by an asterisk). We further investigated the effect of bacterial load (i.e., the amount of bacteria that each host ingests in every selection round) on the evolution of virulence. Fig. 3C demonstrates that high bacterial load leads to low levels of virulence and low bacterial load to high levels of virulence. This effect can be understood against the background of frequency-dependent selection. In Fig. 3C, we vary the bacterial load, but start in all cases with 50% toxin producers. The bacterial load determines the amount of toxin, which is derived from the absolute number of initial toxin producers, . For example, bacterial loads of 100 and 200 result in 50 and 100 units of toxin, respectively. As a consequence, high bacterial load leads to a relatively shorter phase 1 and low bacterial load to a relatively longer phase 1, because more toxin results in faster host death. If phase 1 is comparatively long, this favours the toxin producers, subsequently leading to an overall higher virulence of the parasite population. Note that the length of phase 1 can be calculated analytically and computes to the smallest natural number larger than
/
.
10
Next, we investigated the effect of host body size on the evolution of virulence. Fig. 3D demonstrates that large body size leads to high levels of virulence and small body size to low levels of virulence. Analogous to Fig. 3C, the effect in Fig. 3D can be understood against the background of frequency-dependent selection. In Fig. 3D, we vary the host body size while the bacterial load stays constant. The body size determines the total amount of nutrition available, split by the host’s level of resistance
into
for phase 1 and 1
for phase 2. Increasing
increases the length of phase 1. It is important to note that in the model the toxin acts as a catalyst, providing a constant amount of nutrition to the bacteria, and that the amount of toxin can never exceed the bacterial load. As a consequence, large body size leads to a long phase 1 and small body size to a short phase 1, especially because the bacterial load is not coupled to body size. The effects of host body size and host level of resistance are counter-affected by the bacterial load.
4. Discussion The model presented here provides a mechanistic hypothesis for BT virulence evolution and makes specific predictions for future one-sided parasite adaption experiments. First, the analysis suggests that the parasite evolves towards an equilibrium level of virulence due to frequency-dependent selection. This hypothesis can be tested experimentally by varying the initial level of virulence of the parasite population and carefully monitoring how virulence and time until host death change in the course of evolution. From the model we expect (i) that different starting pathogen populations evolve to the same level of virulence and (ii) that the level of virulence is positively correlated to the average time from infection until host death. Second, the mathematical analysis suggests that virulence evolution is strongly affected by the bacterial load (i.e., the number of bacteria with which hosts are infected in every selection round). For a controlled experiment where the bacterial load is varied, the model predicts that high (low) bacterial loads select for low (high) levels of virulence. Third, we hypothesise that large (small) body size selects for high (low) levels of virulence. Note, however, that body size should not be interpreted as the actual weight of the host, but as the amount of nutrient accessible for the pathogen if the host is fully exploited. This hypothesis may therefore be difficult to test because “body size” in the sense of the model may be difficult to control experimentally.
11
Our model assumes that pathogen replication inside the hosts solely depends on pathogen growth rates, available nutrients and the host level of resistance. In real systems, pathogen growth is a complex process that depends on various factors, including bacteriocin production (Favret and Yousten, 1989), protective symbionts (Zug and Hammerstein, 2015), and the host immune memory (Sadd and SchmidHempel, 2006; Roth et al., 2009, 2010). Though we did not explicitly model these factors, we can hypothesise how they affect virulence evolution based on our theoretical analysis. Bacteriocins are chemical compounds produced by certain bacterial strains to suppress growth of other bacterial strains. Some BT strains are well known to produce bacteriocins (Favret and Yousten, 1989). The effect of bacteriocins on virulence evolution can be modelled in our theoretical framework by changing the growth rates of the toxin and the non-toxin producers relative to each other. A possible outcome may be evolution towards higher virulence. We expect this to happen if bacteriocins impede the growth of toxin producers because this will increase the relative advantage of toxin production due to a prolonged phase 1. Similar effects may occur in experiments that involve protective symbionts (i.e., symbionts that increase host resistance) or immune memory (i.e., previous challenges with pathogens increase host resistance). In particular, we expect protective symbionts to select for higher virulence of BT if they impede toxin and nontoxin producers equally. Though the main focus of the present study was to explain specific results from experimental evolution, the analysis also offers a new perspective on BT evolution in natural ecosystems. Let us assume the virulence evolution of BT follows the twophase model analysed above, i.e. toxin producers have an advantage in living hosts and non-toxin producers in dead hosts. In a natural community, however, BT coexists with other microbes and may experience severe competition for resources. For example, saprophagous microbes may quickly colonise invertebrates killed by BT, reducing the resources available. Based on our theoretical analysis, we expect that this resource reduction will select for higher BT virulence because the non-toxinproducing BT suffer relatively more from increased resource competition than the toxin producers. We therefore hypothesise that communities with high abundance of saprophages select for highly virulent BT. In particular, we expect that host species that are prone to quick colonisation by saprophages, e.g., due to certain foraging
12
behaviours or the microbiome composition, may suffer from especially virulent BT strains. These hypotheses may be tested in future empirical studies. How do our results relate to previous theoretical work? A standard approach to analysing virulence evolution is to determine evolutionarily stable levels of virulence within an epidemiological framework, often under the assumption of a trade-off between pathogen transmission and virulence (Anderson and May, 1982; Ewald, 1983). Several authors stressed the importance of the transmission mode for virulence evolution (Ewald, 1983; Day, 2001). Day (2001) further pointed out that the transmission mode in experimental evolution may deviate from standard models and that mathematical models which are applied to selection experiments need to take into account key characteristics of the experimental protocol. In our study, we implemented the experimental protocol of the Masri experiments (Masri et al., 2015) by differentiating between inter-host dynamics and intra-host dynamics. Further, we chose a difference equation model to mimic the discrete host generations in the serial passage experiments. Thus, instead of trying to generate a new theory of virulence evolution we propose a rather basic model that is useful for experimentalists to test future experiments. Nevertheless, this model may serve as the basis for more sophisticated theoretical analyses. One promising direction may be to include stochasticity, which is surely relevant in the experiments due to the small population sizes (Papkou et al., 2016, current issue). Another direction may be to analyse an epidemiological version of the model, in which parasite transmission follows the mass action principle. However, the analysis of these extended models goes beyond the scope of the current study and is left for future research.
5. Conclusion Most mathematical models on the evolution of virulence are based on epidemiological models that assume parasite transmission follows the mass action principle. In experimental evolution, however, mass action is often violated due to controlled infection protocols. This “theory–experiment’ mismatch” was investigated by focusing on experimental results from virulence evolution (Masri et al., 2015) that are not explained by standard mathematical models. We developed a new model that is not based on the mass action principle for parasite transmission and explains the empirical results. The present study suggests that there is a need for custom-
13
designed mathematical models to capture the dynamics of virulence changes during evolution experiments.
Acknowledgements We thank Sultan Beshir for technical assistance, Joachim Kurtz, Heiko Liesegang, Rebecca Schulte, Gerrit Joop and Esther Sundermann for fruitful discussions on host–parasite coevolution, and Florian Grziwotz and two anonymous reviewers for helpful comments on the manuscript. This study was funded by the German Science Foundation (SPP 1399, TE 976/2-1).
14
References Anderson, R.M., May, R.M., 1982. Coevolution of hosts and parasites. Parasitology 85, 411–426. Agaisse, H., Lereclus, D., 1995. How does Bacillus thuringiensis produce so much insecticidal crystal protein? J. Bacteriol. 177, 6027–6032. Baily, N.J.T., 1975. The Mathematical Theory of Infectious Diseases and its Applications. 2nd ed., Macmillan, New York. Bedhomme, S., Lafforgue, G., Elena, S.F., 2012. Multihost experimental evolution of a plant RNA virus reveals local adaptation and host-specific mutations. Mol. Biol. Evol. 29, 1481–1492. Berngruber, T.W., Froissart, R., Choisy, M., Gandon, S., 2013. Evolution of virulence in emerging endemics. PLOS Pathog. 9, e1003209. Bettencourt, B.R., Feder, M.E., Cavicchi, S., 1999. Experimental evolution of HSP70 expression and thermotolerance in Drosophila melanogaster. Evolution 53, 484–492. Buckling, A., Wei, Y., Massey, R.C., Brockhurst, M.A., Hochberg, M.E., 2006. Antagonistic coevolution with parasites increases the cost of host deleterious mutations. Proc. R. Soc. Lond. B 273, 45–49. Crecy, E., Jaronski, S., Lyons, B., Lyons, T.J., Keyhani, N.O., 2009. Directed evolution of a filamentous fungus for thermotolerance. BMC Biotechnol. 9, 74. Day, T., 2001. Parasite transmission modes and the evolution of virulence. Evolution 55, 2389–2400. Ewald, P.W., 1983. Host–parasite relations, vectors, and the evolution of disease severity. Annu. Rev. Ecol. Syst. 14, 465–485. Favret, M.E., Yousten, A.A., 1989. Thuricin: the bacteriocin produced by Bacillus thuringiensis. J. Invertebr. Pathol. 53, 206–216. Gonzalez Jr., J.M., Carlton, B.C., 1980. Patterns of plasmid DNA in crystalliferous and acrystalliferous strains of Bacillus thuringiensis. Plasmid 3, 92–98. Hethcote, H. W., 2000. The mathematics of infectious diseases. SIAM Rev. 42, 599– 653. Kawecki, T.J., Lenski, R.E., Ebert, D., Hollis, B., Olivieri, I., Whitlock, M.C., 2012. Experimental evolution. Trends Ecol. Evol. 27, 547–560.
15
Kermack, W.O., McKendrick, A.G., 1933. Contributions to the mathematical theory of epidemics. III. Further studies of the problem of endemicity. Proc. R. Soc. Lond. A 843, 94–122. Kubinak, J.L., Ruff, J.S., Hyzer, C.W., Slev, P.R., Potts, W.K., 2012. Experimental viral evolution to specific host MHC genotypes reveals fitness and virulence trade offs in alternative MHC types. Proc. Natl. Acad. Sci. U.S.A. 109, 3422– 3427. Leu, J.Y., Murray, A.W., 2006. Experimental evolution of mating discrimination in budding yeast. Curr. Biol. 16, 280–286. Masri, L., Branca, A., Sheppard, A.E., Papkou, A., Laehnemann, D., Guenther, P.S., Prahl, S., Saebelfeld, M., Hollensteiner, J., Liesegang, H., Brzuszkiewicz, E., Daniel, R., Michiels, N.K., Schulte, R.D., Kurtz, J., Rosenstiel, P., Telschow, A., Bornberg-Bauer, E., Schulenburg, H., 2015. Host–pathogen coevolution: the selective advantage of Bacillus thuringiensis virulence and its cry toxin genes. PLOS Biol. 13, e1002169. Nielsen-LeRoux, C., Gaudriault, S., Ramarao, N., Lereclus, D., Givaudan, A., 2012. How the insect pathogen bacteria Bacillus thuringiensis and Xenorhabdus/Photorhabdus occupy their hosts. Curr. Opin. Microbiol. 15, 220–231. Nowak, M.A., May, R.M., 1994. Superinfection and the evolution of parasite virulence. Proc. R. Soc. B 255, 81–89. Otto, S.P., Day, T., 2007. A Biologist’s Guide to Mathematical Modeling in Ecology and Evolution. Princeton University Press, Princeton, NJ. Papkou, A., Gokhale, C.S., Traulsen, A., Schulenburg, H., 2016. Host–parasite coevolution: why changing population size matters. Zoology 119, ##–##. Paterson, S., Vogwill, T., Buckling, A., Benmayor, R., Spiers, A.J., Thomson, N.R., Quail, M., Smith, F., Walker, D., Libberton, B., Fenton, A., Hall, N., Brockhurst, M.A., 2010. Antagonistic coevolution accelerates molecular evolution. Nature 464, 275–278. Poullain, V., Gandon, S., Brockhurst, M.A., Buckling, A., Hochberg, M.E., 2007. The evolution of specificity in evolving and coevolving antagonistic interactions between a bacteria and its phage. Evolution 62, 1–11.
16
Rafaluk, C., Gildenhard, M., Mitschke, A., Telschow, A., Schulenburg, H., Joop, G., 2015a. Rapid evolution of virulence leading to host extinction under host– parasite coevolution. BMC Evol. Biol. 15, 112. Rafaluk, C., Jansen, G., Schulenburg, H., Joop, G., 2015b. When experimental selection for virulence leads to loss of virulence. Trends Parasitol. 31, 426– 434. Raymond, B., Johnston, P.R., Nielsen-LeRoux, C., Lereclus, D., Crickmore, N., 2010. Bacillus thuringiensis: an impotent pathogen? Trends Microbiol. 18, 189–194. Roth, O., Sadd, B.M., Schmid-Hempel, P., Kurtz, J., 2009. Strain-specific priming of resistance in the red flour beetle, Tribolium castaneum. Proc. Biol. Sci. 276, 145–151. Roth, O., Joop, G., Eggert, H., Hilbert, J., Daniel, J., Schmid-Hempel, P., Kurtz, J., 2010. Paternally derived immune priming for offspring in the red flour beetle, Tribolium castaneum. J. Anim. Ecol. 79, 403–413. Sadd, B.M., Schmid-Hempel, P., 2006. Insect immunity shows specificity in protection upon secondary pathogen exposure. Curr. Biol. 16, 1206–1210. Schulte, R.D., Makus, C., Hasert, B., Michiels, N.K., Schulenburg, H, 2010. Multiple reciprocal adaptations and rapid genetic change upon experimental coevolution of an animal host and its microbial parasite. Proc. Natl. Acad. Sci. U.S.A. 107, 7359–7364. Zug, R., Hammerstein, P., 2015. Bad guys turned nice? A critical assessment of Wolbachia mutualisms in arthropod hosts. Biol. Rev. Camb. Philos. Soc. 90, 89–111.
17
Figure captions Fig. 1. Evolutionary dynamics for scenario 1 (context-independent cost of toxin production). Shown is the frequency of toxin- (red) and non-toxin- (blue) producing bacteria at the end of each selection round. During every selection round, non-toxin producers increase in frequency. If toxin producer frequency decreases to a certain threshold, too few toxins are available among spores to kill new hosts. This results in the end of the selection experiment at selection round eight. Parameters: 1.8,
2,
1000,
0.001,
30,
10,
0.001,
30,
10.
Fig. 2. Intra-host dynamics for scenario 2 (context-dependent cost of toxin production). (A–C) Intra-host dynamics of bacteria for different selection rounds. There are two phases that differ with respect to the cost and benefit of toxin production. Phase 1 (green) lasts until host death and is characterised by faster growth of toxin producers, whereas in phase 2 (red) toxin producers are outcompeted by non-toxin producers. The context-dependent cost results in frequency-dependent selection on bacteria because phase length is determined by toxin-producer frequency. Changes in phase length are shown in (A–C). The full evolutionary 1.5,
dynamics of the selection experiment is shown in (D). Parameters: 1,
2,
1000,
0.34,
0.001,
30,
10.
Fig. 3. Evolutionary dynamics of scenario 2 (context-dependent cost of toxin production) showing how parameter variation affects the evolution of virulence (measured as toxin-producer frequency at the end of each selection round): variation of (A) initial bacterial virulence, (B) host resistance, (C) number of bacterial spores used to initiate infection in each selection round, (D) host body size. For some parameter constellations in (A and B), virulence evolves to such low levels that hosts survive the infection. This results in the end of the selection experiment and is indicated by an asterisk. The legends depict the varied parameters in the same sequence as the plot (top-down). Other parameters and initial conditions: (A) 2, 10; 30, 0.001, 0.3,
1.1,
3.07,
1000,
1.5,
(B) 10; (C) 30,
1,
1.2,
1.5,
10; (D)
0.001,
30,
1.5,
50, 2,
0.3,
0.001,
1000,
1, 1,
2, 2,
50, 1000, 1000,
30, 0.001,
0.525, 50,
10. 18
Figure 1
19
Figure 2
20
Figure 3
21
Table 1. Parameters and symbols of the mathematical model. number of toxin- (non-toxin-) producing bacteria units of indigestible (digestible) nutrients initial number of toxin- (non-toxin-) producing bacteria bacterial load: total number of initial bacteria (
)
initial units of indigestible (digestible) nutrients growth rate of toxin- (non-toxin-) producing bacteria in first phase growth rate of toxin- (non-toxin-) producing bacteria in second phase conversion coefficient from digestible nutrient into bacteria level of host resistance (0…1) bacterial mortality number of time steps (bacterial generations within a single host) number of selection rounds
22