Making pathogens sociable: the [corrected] emergence of high relatedness through limited host invasibility.

Cooperation depends upon high relatedness, the high genetic similarity of interacting partners relative to the wider population. For pathogenic bacteria, which show diverse cooperative traits, the population processes that determine relatedness are poorly understood. Here, we explore whether within-host dynamics can produce high relatedness in the insect pathogen Bacillus thuringiensis. We study the effects of host/pathogen interactions on relatedness via a model of host invasion and fit parameters to competition experiments with marked strains. We show that invasibility is a key parameter for determining relatedness and experimentally demonstrate the emergence of high relatedness from well-mixed inocula. We find that a single infection cycle results in a bottleneck with a similar level of relatedness to those previously reported in the field. The bottlenecks that are a product of widespread barriers to infection can therefore produce the population structure required for the evolution of cooperative virulence.

Introduction

Pathogenic bacteria appear to be particularly dependent upon cooperation (West ). They rely upon the secretion of a wide range of signals and virulence factors such as pore-forming toxins, siderophores and quorum-sensing molecules (West and Buckling, 2003; Brockhurst ; Diggle ; Raymond ). While many retain the private benefits of virulence factors by injecting them locally into host cells using type IV and type III secretion systems (Bomberger ), a large number of pathogens release virulence factors that act at a distance (West and Buckling, 2003; Schmid-Hempel and Frank, 2007; Bomberger ; Raymond ). By releasing metabolites into extracellular space, the benefits of these virulence factors are exploitable by neighbours. These shared products have been commonly described as ‘public goods' (Sachs ). A key issue for the persistence of this form of cooperation is its vulnerability to exploitation by ‘cheater' mutants that reap the benefits of investment by their neighbours without paying the costs themselves (Frank, 1998). High relatedness, which increases the inclusive fitness of cooperators, can prevent cooperative strains from being out-competed by cheater mutants (Hamilton, 1964; Griffin ; Diggle ; Gilbert ). However, the population structure of pathogens has rarely been resolved at a fine enough scale to allow ecologically relevant calculations of relatedness (but see Raymond ).

Explaining how high relatedness arises is an area of much interest in microbial social biology (Nadell and 2013; Mitri ; Damore and Gore, 2012). Mechanisms such as kin discrimination promote high relatedness and social interactions among kin in higher organisms (West ; Gilbert ). In bacteria, kin-discrimination mechanisms are largely unknown. However, bacteria may be able to alter the relatedness of their local environment via antagonistic secretions such as bacteriocins, which can selectively impose mortality on non-kin (Garbutt ). Horizontal gene transfer can also actively increase relatedness in bacterial populations, as bacteria can ‘infect' neighbouring cells with genes required for cooperative activity (Rankin ). However, the broad impact of antagonism and horizontal gene transfer on bacterial population structure is still largely unknown. In general, stochastic processes may be very important in microorganisms. Stochastic sorting of different genotypes can occur during dispersal and during the division and growth of bacteria in biofilms. This can, theoretically at least, lead to population bottlenecks and the spatial structure required to maintain the production of public goods (Nadell and 2013; Korolev ; Mitri ; Datta ; Momeni ). Although spatial structure in biofilms clearly is important for maintaining high relatedness on a local scale, it is unclear how this translates to the maintenance of relatedness during transmission of bacteria from host to host. Here, we will explore the role of stochasticity during the invasion of hosts by bacterial pathogens and its impact on relatedness.

Mixed infections by genetically diverse bacteria of the same species are believed to happen less often than single strain infections in humans, but still happen in, on average, 21.7% of infections (Balmer and Tanner, 2011), and some studies show an even higher rate (Levert ). A number of studies have shown that the progress of infection through hosts is associated with population bottlenecks and the loss of pathogen diversity (Moxon and Murphy, 1978; Brown ; Grant ; Plaut ; Prajsnar ), so that a plausible general explanation for the low diversity (and therefore high relatedness) of bacterial parasites in hosts is the existence of strong bottlenecks during colonisation. However, many of these previous examples involve infections generated artificially in model hosts (Grant ; Prajsnar ), or diseases such as inhalational anthrax and meningitis (Moxon and Murphy, 1978; Plaut ), which do not represent a normal or successful mode of transmission for these pathogens. The importance of social interactions, and hence changes in diversity, has also not been explored for any of these pathogens.

It is therefore unclear how natural population processes, over the entire natural cycle of infection, affect relatedness at a spatial scale that is meaningful for microbial social interactions and on-going transmission. For Bacillus thuringiensis, a widespread pathogen of invertebrates (Raymond ), the importance of social interactions and high relatedness is very clear. These bacteria produce proteinaceous crystal (Cry) toxins at sporulation that determine virulence and host-range and that are required for pathogenicity (Schnepf ). Cry toxin production is cooperative at the scale of a single host, and natural populations are characterised by high levels of relatedness (Raymond ), although the population processes leading to high relatedness have not yet been clearly established. Here, we hypothesise that the successive bottlenecks involved in transmission can generate levels of relatedness (at the host scale) similar to that seen in field populations of the entomopathogen, B. thuringiensis (Raymond ). Thus, general within-host population processes should preserve the high levels of population structure required for effective cooperation and maintenance of virulence.

Materials and methods

We explored experimentally and theoretically whether within-host population dynamics were sufficient to explain observed high levels of relatedness at the level of the host. We specifically considered the dynamics of mixed infections and explored four parameters that can affect the structuring of co-infecting genotypes. The theoretical model was parameterised and tested using in vivo experiments with competing genotypes of B. thuringiensis, marked with different spontaneous antibiotic resistance (to rifampicin and nalidixic acid). Insects were inoculated orally by feeding on droplets containing a well-defined dose (300 cfu (colony-forming unit) in total, approximately LD90) with varying proportions of each genotype.

Model derivation

The probability density function of the final strain abundance was derived based on an invasion-growth model. Time to successful invasion for each strain was exponentially distributed:

with Λ the successful invasion rate of strain i and t0 the time of invasion. After invasion, each strain grew deterministically and at a constant rate. The density of strain i(n) at time t can then be defined as follows:

with γ, the growth rate within the host. When considering two strains (a and b), it is possible to rescale the parameters without loss of generality as follows: the number of infectious pathogens successfully entering a host per unit time (λ=(Λ+Λ)/γ), the proportion of the strain a in the inoculum of co-invading genotypes (between 0 and 1;p), the relative growth rate of strain b (γ=γ/γ) and the time since infection (τ=γt). On the basis of these assumptions, we were able to derive the full probability distribution of the proportion (x) of the two strains after a certain time (f(x|p,λ,γ,τ); see Supplementary Appendix S1 in Supporting Information for more details). The final proportion of each strain was highly dependent on which strain invaded first, because that strain has more time to grow and become dominant. We also explored how sensitive the model was to the choice of growth model and show that the results will hold up for a Lotka-Volterra competition model as long as both strains are sufficiently similar (see Supplementary Appendix S2 in the Supporting Information).

Relatedness was used as a measure of the population structure after infection. We calculate the relatedness coefficient as defined by Queller and Goodnight, 1989 from our model to further explore the effect of the infection barrier. The resulting relatedness coefficient is between 0 and 1. Relatedness can then be interpreted as the probability to pick (with replacement) two individuals of the same strain from one subpopulation (insect) minus the probability to pick two individuals of the same type from the population as a whole, divided by 1 minus the probability to pick two individuals of the same type from the populations as a whole. If the subpopulation was dominated by one strain, but not the whole population, then the relatedness is close to 1. Note though, that if the population as a whole is also dominated by the same strain, then relatedness coefficient is low, because the genetic similarity in the sub-population reflects that across the entire population. We calculated the relatedness using the most likely parameter values. See Supplementary Appendix S1 in Supporting Information for more details on the model derivation and fitting.

Model validation

A likelihood-based method was used to find the most likely parameter values given the data. The likelihood of different parameter values was calculated based on the abundance counts collected during the experiment. We then used an adaptive Differential Evolution-based algorithm to find the parameters which resulted in a maximum likelihood value (Brest ). Next, we bootstrapped the data to find the upper and lower limits of the parameter values using a 95% confidence interval. The relative importance of the underlying processes could be inferred from these results.

We compared the findings of the model presented above with an alternative model, where we assumed that the main limiting process in the infection was the survival of the inoculum to initial infection. In this case, a low survival rate would lead to large variation in the final proportion (x), because the final proportion is dependent on the number of bacteria of each strain to survive till successful founding of the infection. Comparing the results of the fit of our model with this alternative model was performed to quantify whether our model provides a better explanation of the observed experimental results. See Supplementary Appendix S3 in the Supporting Information for more details on this alternative model.

Experimental design

Spontaneous antibiotic resistance mutants of B. thuringiensis kurstaki HD-1 were isolated from the commercial biopesticide preparation, DiPel WP (Valent Biosciences, Libertyville, IL, USA), by plating high numbers of cells (108 + cfu) on 100 μg ml−1 rifampicin or 15 μg ml−1 nalidixic acid. Antibiotic-resistant mutants with reduced fitness costs were isolated after one to three rounds of host passage in Plutella xylostella (Garbutt ), and identified by rapid growth on selective plates. Lack of cross resistance between two mutants R4K1 RifR and P1G NalR was confirmed by thickly streaking sporulated cultures across LB agar plates containing appropriate antibiotics. Although these strains had a recent common ancestor, we also investigated the pleiotropic consequences of resistance mutations with several life history experiments. These included pathogenicity assays at three doses with 30–40 insects per dose, as per the droplet assay described below and measures of spore productivity in infections with a single strain (with at least N=40 cadavers per strain).

Sporulated cultures of both strains were produced by growing dense lawns of bacteria on HCO sporulation media (Lecadet ) at 30 °C for 1 week. Spores and Cry toxins were recovered from plates and washed twice in sterile saline (0.85% NaCl), before being diluted into 10 ml of saline and stored at −20 °C in 0.5 ml aliquots. Defrosted spores were enumerated by plating serial dilutions; counts were made within 48 h of infecting insects.

P. xylostella larvae were reared on artificial diet (without antibiotics) as described previously (Raymond ). Late second instars were infected with 300 cfu of Bt in droplet assays. The final droplet mix contained 10 mM sucrose, 7.5% v/v green food dye (Dr Oetker, Bielefield, Germany) and 0.4% w/v agar (Oxoid, Basingstoke, UK), and 40% v/v cabbage extract (filtered liquid from boiled cabbages). The cabbage juice, sucrose and food dye were filter-sterilised before being used to dilute the spores; this mixture was then combined (50:50) with molten 0.8% agar at 60 °C. The resultant inoculum (with a final concentration of 300 cfu μl−1) was briefly held at 50 °C in a heat block (for no more than 10 min) while 1 μl droplets were dispensed into each well of 48-well plates using pre-warmed pipette tips. A single larva was added to each well, and plates were tightly sealed with damp tissue paper. Larvae were allowed to feed for 18 h. When droplets were at least 75% consumed and green dye was visible within larvae (Figure 1), insects were moved onto artificial diet for 5 days. Mortality rates were typically greater than 90% in these assays. Cadavers were transferred to 2 ml homogenisation tubes containing 10 μl of sterile water and incubated at 30 °C for at least 7 days. This procedure ensures that all Bt within cadavers had fully sporulated before homogenisation and final enumeration. Before homogenisation, 500 μl of sterile saline (0.85% NaCl) was added to each homogenisation tube and cadavers were pasteurised at 65 °C for 20 min to kill any remaining vegetative cells. Cadavers were homogenised in a beadbeater (Qiagen Tissue Lyzer, Manchester, UK) using 4-mm diameter steel ball bearings. Homogenates were serially diluted in saline and plated out on rifampicin (100 μg ml−1) and nalidixic acid (15 μg ml−1) LB agar plates in 15 μl droplets before being incubated overnight at 30 °C. As a general rule, dilutions providing colony counts between 5 and 50 per droplet were included in the data set. Replicated droplets counts (within cadavers) confirmed the repeatability of this enumeration technique.

Figure 1

Second instar diamondback moth larvae were infected with B. thuringiensis spores in agarose droplets containing green food dye. This allowed us to readily identify insects that had consumed a defined dose.

Results

The model predicts the exact probability of finding a certain final proportion of competing bacteria in a host given different parameter values (Figure 2). To better understand the implications of the possible probability density functions, we classified them according to three different attributes. Firstly, the final population can be ‘well-mixed', where there is a high probability of the abundances of the two strains being about equal in each cadaver. In contrast, the population can also be ‘bimodal', where there is a high probability that one strain will dominate the other strain in a cadaver. Although, in this case, each sample is dominated by one strain, both strains have an equal probability of becoming the dominating strain. Finally, the population can be ‘skewed', where most cadavers are dominated by one particular strain, that is, the population as a whole is dominated by one strain. Next, we studied in detail how the different parameters influence these different attributes of the probability density function and the resulting relatedness of the population as a whole.

Figure 2

Three typical outcomes of the model, which can be classified as follows: (i) well-mixed, i.e., there is a high probability of finding both strains at roughly equal abundances, (ii) bimodal, where one strain dominates the other, but both strains have an equal probability of being the dominant strain and (iii) skewed, where one strain is more likely to be the abundant strain.

Relatedness is high when the genetic similarity within host communities is high relative to the population at large (Queller and Goodnight, 1989). Bimodality in part reflects the predicted relatedness of pathogens after one round of infection, as under highly bimodal distributions competing genotypes will be assorted into different cadavers. We measured the bimodality of the probability density function by calculating the coefficient of bimodality, which is a function of the third and fourth moment of the distribution and ranges between 0 and 1 (Ellison, 1987). Low invasion rates (λ) resulted in a high coefficient of bimodality and therefore in high relatedness (Figures 3a and c). In contrast, increases in the difference between the growth rates (γ) or the difference in initial abundance in the inoculum (p) also result in high bimodality, but have a less pronounced effect on relatedness, because they also increase skewness in the population (Figures 3b and c). Skewness ensures that a single strain dominates both the within- and between-host communities and therefore has a less pronounced effect on relatedness. These results explain why changes in the invasion rate have an especially profound effect on the population structure and relatedness

Figure 3

Effect of different parameter values on: (a) bimodality, (b) skewness and (c) relatedness. Under these parameter values, skewness is not affected by the relative invasion rate or time since infection, so both parameters are represented by the same line in b. Default parameter values used: proportion in the inoculum (P=0.5), relative invasion rate (λ=1) and time since infection (τ=10).

Mixed infections in insect hosts using near isogenic B. thuringiensis strains revealed strongly bimodal patterns of genotype frequencies in infectious spores produced at the end of infection (Figures 4a and c), a pattern consistent with low invasion rates. Maximum likelihood-based model fitting produced a close fit between the theoretical model and each of three independent experiments initiated with different genotype frequencies. Using the maximum likelihood, we calculated the Akaike Information Criterion (AIC; Akaike, 1973) and compared this model with an alternative model, which assumed that the main cause of this variation is due to small founding populations. We found essentially no empirical support for the alternative model when compared with our model (Burnham and Anderson, 2002; see Supplementary Appendix S3 in the Supplementary Information for more details). Results were consistent for all data fitted simultaneously and for the three independent data sets. Parameterisation of the model in the three experiments was also consistent; in each case, the relative growth rates were similar (Table 1), and the overall relative growth rate of the NalR strain was 1.09, indicating that the fitness of the two mutants was very similar. Given the estimated parameter values, relatedness is consistently high (0.7–0.8) regardless of inoculum proportion (Figure 4d). This highlights that the low infection rates measured in this system result in high relatedness even when starting from a well-mixed population in the inoculum.

Figure 4

Model fitting results. Blue corresponds to the model when fitted to all the data, while red corresponds to the model fitted to the data with that proportion in the inoculum. (a–c) The lines show the predicted probability density function. The dotted lines show the same probability distribution, but taken over intervals. This makes visual comparison between the model and the data easier. The grey histogram shows the data. Each sub-figure shows the data and results grouped by inoculum proportion. (d) Relatedness as dependent on the initial inoculum population using the parameters from the fitted models.

Table 1

Results of the parameter fitting

Inoculum (p)	Parameter	Max L	5th percentile	95th percentile
All data	λ	0.388255	0.369545	0.433066
All data	γ	1.08953	1.04781	1.12001
All data	τ	12.8511	12.3539	12.8601
0.278	λ	0.45007	0.403914	0.525843
0.278	γ	1.04707	1.03256	1.09277
0.278	τ	12.8601	12.3224	12.8601
0.429	λ	0.398937	0.353042	0.546645
0.429	γ	1.13228	1.08419	1.27297
0.429	τ	12.3659	10.9993	12.3659
0.806	λ	0.401246	0.37108	0.487739
0.806	γ	0.9846	0.963852	1.03481
0.806	τ	12.5151	11.9057	12.5233

We show the maximum likelihood parameter value and the uncertainty in the parameter value. The percentile interval has been measured using a resampling approach. Results between the different data sets are consistent, with the invasion rate between 0.35 and 0.55. Relative growth rates refer to the growth of the NalR strain/ RifR and are close to 1; the time of measurement is between 10 and 13.

More direct calculations of relative fitness validate the results of the model fitting. The change in proportion of a particular mutant (P1G NalR) can be used to calculate fitness (Ross-Gillespie ). After dividing this estimate by log2 (total spore counts), we can estimate the relative growth rate per generation at 1.08 (N=332 larvae, 95 % confidence intervals 1.06–1.10). Additional life history differences confirmed that there exist small, but significant differences in pathogenicity (χ2=7.43, df =1, P=0.0064, Supplementary Information S4) which would equate to approximately a 15% difference in mortality rates at our experimental dose. The faster growing and less pathogenic strain (NalR) also produced slightly more spores per insect host in single strain infections (F1,83=7.87, P=0.0063; Supplementary Information S4).

Discussion

Here, we have shown that the simple natural process of infection can produce strong between-host heterogeneity in pathogen distribution from an initially well-mixed inoculum containing competing genotypes. Clonal mixing, or low relatedness, in reservoirs of naturally distributed bacteria therefore presents no real challenge to socially based explanations for maintenance of secreted virulence factors with moderate fitness costs, because a single round of infection can produce the high levels of relatedness required to maintain cooperation. Mixed infections of pathogenic bacteria are relatively infrequent in many species infecting human hosts (Balmer and Tanner, 2011). Rarity of multiple infections can be explained by bottlenecks in the host colonisation process as well as barriers to super-infection, the invasion of a host already colonised. For example, Staphylococcus aureus actively prevents super-infection in a rat model (Margolis ), and this species has relatively low rates of mixed infection (<10%) (Cespedes ; Balmer and Tanner, 2011). Where there are no barriers to super-infection, as in Haemophilus influenzae (Margolis ), levels of reported mixed infections are correspondingly higher (53%) (Balmer and Tanner, 2011).

It could be argued that strong bottlenecks are a feature of infection processes in a small invertebrate host. However, the spatial scale at which social interactions occur among microbes can be even smaller (Slater ) and several pathogens display patterns of declining diversity as infections proceed (Moxon and Murphy, 1978; Brown ; Grant ; Plaut ; Prajsnar ). Low diversity may even be the product of strategies to maintain very low densities within phagocytes to promote immune evasion (Prajsnar ). In this study, the spatial scale of study is appropriate for the major virulence factors (Cry toxins) of our chosen pathogen, because cooperation and cheating based on these solubilised toxins occurs in the midgut of single insects (Raymond ). Insect cadavers are a suitable spatial scale for assessing relatedness because of their importance for spore production and ongoing disease transmission, which can occur orally via cannibalism or consumption of patches of foliage contaminated by cadavers (Knell ; Raymond , 2012). A caveat here is that transmission can also occur indirectly, when bacteria from a soil reservoir colonise plants (Raymond ). Note also, that infections in this study used a relatively high dose of inoculum, capable of killing 80–90% of hosts. This density of spores (>102 cfu cm2) is rare in natural environments (Collier ), and low dose infections may further reduce invasibility and increase relatedness (Moxon and Murphy, 1978). Despite this expected difference in doses in the field and the laboratory, measurements of population structure in this study were high and comparable with field measurements of relatedness at the Cry toxin locus of B. thuringiensis in the field (Raymond ).

It has been argued that fluctuating population structure poses a challenge to the generality of inclusive fitness theory in microbes (Damore and Gore, 2012). The reasoning is that real-world complexities such as non-linear fitness benefits and changes in relatedness are likely to lead to dynamical fluctuations in b, the benefit deriving from cooperation. This makes the ‘disentangling of fitness effects from population structure often impossible' (Damore and Gore, 2012). Our data suggest that this problem is not as severe as they suggest. We have previously made some attempt to link population processes to the maintenance of cooperation (Raymond ) and would argue that population dynamic processes can act to stabilise relatedness in this pathogen. In the field, relatedness at the locus of Cry toxin production in B. thuringiensis can fluctuate within a single patch, but is very stable in the whole population, despite variation in bacterial density and the overall proportion of toxin producers (Raymond ). In this study, population bottlenecks ensure that relatedness at nearly neutral loci remained high, despite variation in the initial frequencies of competitors. In the field, similar bottlenecks in infection and transmission (Zhou ) or in the colonisation of plants from spore reservoirs in the soil (Raymond ) may also stabilise high relatedness.

Previous observations of extreme bottlenecks during infection have been interpreted as evidence for the ability of bacteria to act independently during infection, as infections appear to be established from single cells, in some circumstances (Moxon and Murphy, 1978; Rubin, 1987). The ‘independent action hypothesis' proposes that the probability of infection resulting from a single cell is independent of the presence of conspecifics (Druett, 1952). The independent action hypothesis has widespread importance in epidemiological modelling and risk assessment (Haas, 1983; Buchanan ). Bottlenecks have therefore been interpreted as being evidence which supports the independent action hypothesis (Moxon and Murphy, 1978; Rubin 1987). However, the fact that pathogen populations pass through clonal or near-clonal bottlenecks does not mean that cells are acting independently. Shared exploitation of public goods is critical in the virulence of B. thuringiensis (Raymond ), and this may be the case in many other bacterial pathogens (Schmid-Hempel and Frank, 2007). Public-goods-based virulence therefore challenges the key assumption of independent action, as cells must work together to successfully overcome barriers to infection, as recent experiments with B. thuringiensis have demonstrated (Cornforth ). The fact that diversity is lost as infections proceed does not undermine the idea that bacteria can cooperate; rather it illustrates how the restricted invasibility of pathogens can facilitate a population genetic structure that favours cooperative interactions.