The mode of host–parasite interaction shapes coevolutionary dynamics and the fate of host cooperation

Antagonistic coevolution between hosts and parasites can have a major impact on host population structures, and hence on the evolution of social traits. Using stochastic modelling techniques in the context of bacteria–virus interactions, we investigate the impact of coevolution across a continuum of host–parasite genetic specificity (specifically, where genotypes have the same infectivity/resistance ranges (matching alleles, MA) to highly variable ranges (gene-for-gene, GFG)) on population genetic structure, and on the social behaviour of the host. We find that host cooperation is more likely to be maintained towards the MA end of the continuum, as the more frequent bottlenecks associated with an MA-like interaction can prevent defector invasion, and can even allow migrant cooperators to invade populations of defectors.


INTRODUCTION
The maintenance of cooperation is an evolutionary conundrum: why invest resources if the benefits of the investment are returned to other individuals? A major line of explanation for the persistence of cooperation in the face of non-investing 'cheats' is that the benefits generated by cooperators return preferentially to cooperative individuals, as a result of non-random population structure [1][2][3][4][5][6]. Here, we extend the study of cooperation by turning our focus to a ubiquitous driver of population structuring: parasites.
Parasites shape host populations owing to the deleterious effects they have on their hosts [7,8]. Given that parasites rely on their hosts for resources, changes in the host population composition will also impact upon parasite demography and genetic structure [9]. A tight genetic interaction between hosts and pathogens can lead to ongoing host-parasite coevolution, defined as the reciprocal evolution of interacting hosts and parasites [10]. The interaction is usually antagonistic, and selects for hosts to evolve resistance and parasites to evolve infectivity. A key consequence of coevolution is the impact on genetic diversity of host and parasite populations, although this will critically depend on the type of coevolution dynamics. At one extreme, coevolution follows an arms race dynamics (ARD), with directional selection for increasing resistance and infectivity, purging diversity. At the other is fluctuating selection dynamics (FSD), where parasite genotypes specialize on host genotypes, potentially resulting in the maintenance of large amounts of host and parasite genetic diversity in time and space [11][12][13].
An important determinant of coevolutionary dynamics is the underlying specificity. The two most common representations of specificity are the matching alleles (MA) model and the gene-for-gene (GFG) model, although other variants exist, for example [14]. MA models are based upon a system of self/non-self recognition molecules where hosts can successfully defend against any parasite genotype that does not match their own [15,16]. Typical of many invertebrate immune systems, MA models assume that one parasite genotype will have a different subset of susceptible hosts than another parasite genotype. Infection is therefore determined by both the host and parasite genotypes, with such tight specificity sometimes leading to FSD [17]. The GFG model, favoured by plant pathologists, predicts whether infection is successful based on the interaction between resistance loci and virulence loci [13,18,19]. GFG models are often characterized by directional ARD [18,20,21], which can lead to the evolution of generalist parasites [22], although if there are costs associated with increased resistance/infectivity ranges, FSD can also arise [13,21].
As with coevolutionary dynamics themselves, MA and GFG models can be understood as the two extremes of a continuum of specificity, and the interactions between most hosts and parasites are likely to lie somewhere between the two extremes with some degree of specialization and some generalization. Agrawal & Lively [13] investigated the dynamics of hybrid MA -GFG models via the use of a single parameter which formed a continuum: MA at one extreme and GFG at the other. Accurate characterizations of specificity patterns have important consequences for predicting the evolution of virulence [23,24], patterns in local adaptation [25 -27] and the evolution of recombination [18,28]. Now we ask: how does genetic specificity shape host-parasite coevolution, host population genetic structuring and in turn, the maintenance of host cooperation?
A recent study on cooperation in a bacterial model system highlighted that strong selection on the host population (adaptation to the passaging environment and to an antagonistic phage virus) allowed a primarily cooperative host population to purge low-frequency cheats, presumably owing to cooperator alleles hitchhiking on beneficial resistance mutations [29]. Using a fully stochastic model of host-parasite dynamics, we dissect this bacteria-phage interaction, and extend to a coevolutionary time-scale, spanning GFG and MA mechanisms of interaction. We predict that the MA limit is the most favourable for the maintenance of cooperation because of increased genetic turnover allowing repeated purges of rare defector alleles.

MATERIAL AND METHODS
We follow a stochastic framework widely used in population dynamics [30]. Numerical simulations were performed using an exact method, the Gillespie algorithm [31]. The analysis takes place within the framework of an ecological model of the infection of bacteria by lytic phages, using stochastic population dynamics. The parameter values used are based on biologically realistic values [32].
(a) Simple host-parasite model To start with, we consider a single type of host (a bacterium) and its parasite (a virus-in particular, a lytic phage). The changes in the discrete numbers of bacterial hosts n H and free-living viruses n V are due to four processes: host birth, competition among hosts, virus death and lysis of a bacterium by a phage (resulting in the instantaneous death of the bacterium and the release of y copies of the virus; we assume that the latent period is negligible). They can be captured using the following reactions: From these microscopic reaction rates, we can derive equations for the average concentrations of hosts and parasites (see electronic supplementary material, S1), H and V, which follow the deterministic trajectories (in a mean-field approximation): The system is similar to a predator-prey system [33], differing solely in the interpretation and magnitude of the quantity y which stands for the burst size and is of order 10-100-while in a predator -prey system the quantity (y -1) would be termed ecological efficiency and takes a value less than 1.
(b) Host-parasite coevolution model The simplest scenario that allows for host -parasite coevolution is that between two types of host and two types of parasite. Using the notation of a single locus and two alleles, we will call the hosts H a and H A and allow for mutations among them, and similarly for the two types of parasite V a and V A .
Depending on the model of specificity, the infectivity and susceptibility ranges of host and parasite (i.e. the parameters of the four different lysis reactions that can in principle take place) will vary. As in Agrawal & Lively [13], we will use a single parameter q to translate across the MA-GFG continuum (MA: q ¼ 0, GFG: q ¼ 1). The parameter p represents the binding efficiency of a specialist virus (which can only bind to one host type) and p g ( p) the binding efficiency of a generalist virus (which can bind to two different host types). We have chosen the form of these interactions to be a linear interpolation from MA to GFG, as summarized in table 1.
When q ¼ 0 (at the MA end of the continuum), there are two specialist host-parasite pairs: (V a , H a ) and (V A , H A ). When q ¼ 1 (at the GFG end), one virus is a specialist (V a can only infect H a ) but the other one is a generalist (V A can infect both H a and H A ), albeit at the cost of a smaller binding efficiency. From the point of view of the hosts, H A carries a resistance allele (to V a ) and incurs in a cost of resistance in the form of a death rate z that scales with q. Therefore, the reactions that we must add to the ones considered in the simple host-parasite model are Starting from the stochastic formulation and proceeding as in the previous simple host-parasite model, we can derive the following deterministic equations for the evolution of the average concentrations of hosts and parasites (see electronic supplementary material, S1): Note that there is no direct competition between parasites, only an indirect one due to sharing common prey (hosts) when q . 0. The effect of q can be interpreted as a modification of the structure of the host-parasite network Table 1. Interaction parameters (lysis rates) for the four host-parasite pairs of the coevolution model ( §2b).
Coevolution and cooperation B. J. Z. Quigley et al. 3743 of interactions (see electronic supplementary material, S6). This coevolution model can in fact be taken to be a minimal motif within a generic host -parasite food web. The change from q ¼ 0 (MA) to q ¼ 1 (GFG) corresponds to the change from a modular to a nested network structure (diagonal and lower-triangular interaction matrices, respectively).
For a recent study about the structure of experimental bacteria -phage networks, see Flores et al. [34]. In our model, there is no clear time-scale separation between the ecological and the evolutionary dynamics (between the dynamics on and of the host -parasite network). Evolutionary dynamics here refers to the mutations allowed within the definition of the model (not to the evolution of the strength of the host -parasite links i.e. of q itself-that is left for future work). The stochastic population dynamics will take the ecosystem through a series of network configurations, as nodes are populated and links activated-due to mutations-or removed-due to extinctions. For an example of a model of community assembly and for further references, see Capitán & Cuesta [35].

(c) Cooperation and coevolution model
We introduce a third model by adding a cooperationdefection dilemma to the coevolution model described in §2b. Now every host has two traits: the original one (that determines resistance to parasites) labelled by a and A, plus a social trait labelled C for cooperators and D for defectors. So we consider four types of host in total: C a , D a , C A and D A , as well as the two types of viruses V a and V A . We capture the social dilemma posed by the interaction between cooperators and defectors by cooperators paying a fixed cost (w), while providing benefits to every individual in the local population, be they a cooperator or a defector. Specifically, we include an extra birth rate for all hosts, proportional to the fraction of cooperators in the population (see electronic supplementary material, S1). The reactions we must add are the cost of cooperation and mutations in the cooperation-defection trait:

RESULTS AND DISCUSSION
To analyse the influence of host-parasite coevolution on the maintenance of host cooperation, we take a step-wise approach. First, we consider the impact of host evolution on host genetic structure, in the face of a static environmental challenge (following Morgan et al. [29]). We then introduce coevolutionary dynamics, without explicit cooperation and defection, to explore the basic impacts on host structure of reciprocal evolutionary antagonism. Finally, we analyse the full model, to track the fate of cooperation in the context of antagonistic coevolution.
(a) Host evolution only We begin by simplifying the full model ( §2c) by removing parasite evolution (starting with V a only and setting m V ¼ 0), in order to mimic the system studied in Morgan et al. [29]. In agreement with Morgan et al. [29], we find that strong selection on a non-social trait (here parasite resistance) generates positive frequencydependent selection for cooperation (see figure 1 and the electronic supplementary material, S2). The simple heuristic model in Morgan et al. [29] offered an interpretation for this effect: because beneficial non-social mutations are more likely to arise in the numerically dominant cooperator population, the cooperation allele is likely to hitchhike with these beneficial mutations that sweep through the population. However, the model paid for its simplicity with a number of stringent assumptions: no cost of cooperation, no demographic dynamics, deterministic gene frequency change and no coevolution. Introducing a cost of cooperation in our explicitly dynamical and stochastic framework illustrated that, as the relative cost of cooperation increases, cooperators are required to be initially present at higher frequencies in order for cooperation to be maintained (figure 1). When there is no cost of cooperation, the positive frequency-dependent effect of strong selection favours the maintenance of the more common genotype (more than 50%). As the cost of cooperation is increased, this 'break-even' threshold is shifted towards higher frequencies of cooperators, as the intrinsic cheater advantage to defectors becomes increasingly important. This result aids the interpretation of Morgan et al.'s experimental results, where the break-even point ranged from about 90 per cent cooperators (environment and phage adaptation) to about 98 per cent cooperators (environment adaptation only). These biases from the purely symmetrical 'survival of the commonest' result (figure 1, w ¼ 0) reflect an increasing cost of cooperation, relative to the presumed benefits of adaptation to the environmental challenge.
The Morgan et al. [29] model predicted that the frequency-dependent effect highlighted in figure 1 would break down at low densities (as neither cooperator nor defector lineages are likely to gain resistance) and also at high densities (as both are likely to gain resistance, again removing the relative advantage to either lineage).
In figure 2, we revisit these predictions, and find in contrast that the frequency-dependent effect is robustly maintained at higher densities. Our explicitly dynamical and stochastic model allows resistance mutations to be separated in time, so even if at high densities both lineages are likely to gain a critical resistance mutation, there is still a decisive advantage to gaining the resistance mutation first-and this is more likely to be the case in the dominant (higher frequency) lineage.

(b) Host-parasite coevolution
We now turn our focus to the effects of coevolution on the host population dynamics, first in the absence of any social conflict (as described in §1b). Note that the reciprocally antagonistic nature of the bacteria-virus interaction can readily lead to instability, and loss of the pairing in a spatially unstructured population [36]. In order for the system to be stable (i.e. for hosts and parasites to coexist, without global extinctions), a substantial cost of generalism must be imposed (see the electronic supplementary material, S4), both on the resistant host and on the generalist parasite. A cost of wider viral infectivity range is expressed in our model in terms of a reduced binding efficiency, consistent with experimental data [37] and a trade-off mechanism based on antagonistic pleiotropy [38], i.e. mutations that allow the virus to infect a broader host range trade off against viral productivity. A cost of generalized resistance is captured in our model via an intrinsic per capita cost to the bacterial host, consistent with experimental data [39 -43].
Following the establishment of a robust coexistence regime (see the electronic supplementary material, figure S4), we turn to an analysis of host population bottlenecks as a function of the specificity of the genetic interaction between host and parasite. Consistent with Agrawal & Lively [13], we observe a reduced number of host population bottlenecks as q is increased-i.e. moving from MA to GFG (figure 3). Underlying each bottleneck is in fact an attempt of invasion by the A genotype (originated by a mutation) into the resident a population. The structure of the host-parasite network of interactions, fixed by the value of q (see electronic supplementary material, S6), determines the outcome. Invasion attempts fail repeatedly at the MA end, as the invading H A genotype causes an explosion in its specialist opponent V A , which in turn abruptly depletes the H A population in a negative feedback loop (so we observe many bottlenecks, as the a/A battle continues).
In contrast, a positive value of q stabilizes the system, preserving a balanced polymorphism in host resistance alleles. Indeed, at the GFG end H A manages to successfully invade once and then coexist with the resident H a due to the extra link in the network (electronic supplementary material, S6, its virus V A is devoting part of its resources to attack H a too) and because the costs of generalism come into play: the binding efficiency of V A to H A is lower; the intrinsic growth rate of H A is also lower due to the cost of generalism, paradoxically raising its chances of a successful invasion. It seems that these costs of generalism increase the period of the H A -V A host-parasite oscillation, thus enabling the successful establishment of H A (see electronic supplementary material, S7).

(c) Coevolution and cooperation
We now study the full model of §2c, where there is a cooperation-defection dilemma for the hosts, as well as host-parasite coevolution. We can now ask: what is the Coevolution and cooperation B. J. Z. Quigley et al. 3745 fate of cooperation as we move along the host-parasite specificity continuum from MA (q ¼ 0) to GFG (q ¼ 1) limits?
The MA model is associated with a diversification of both host and parasite genotypes [44] and so at first sight is the most inimical to the preservation of cooperation, which is typically favoured by local genetic homogeneity [45]. However, we find that cooperation is most robust to cheater invasion in the MAs limit ( figure 4). The probability of defector takeover grows with q as we move from MA to GFG (figure 4). In other words, cooperation is maintained for longer periods if host-parasite specificity is MA-like. This result is robust across a large region of parameter space (although system stability may be limited-see §3b).
The non-intuitive association between diversifying MAtype interactions and the preservation of cooperation can be understood when it is recognized that the diversifying selection is not acting on the cooperation locus, but only on the resistance locus. The result of the continual strong alternating selection on resistance is a continual cycle of genetic bottlenecks purging diversity at linked loci, which here means the cooperation locus. In the MA limit, recurrent a/A bottlenecks prevent rare defectors from invading-so cooperation is maintained (any D A defectors that may have arisen are wiped out, just like the majority of C A , due to the burst in V A ). On the contrary, in the GFG limit, a single switch from a to A hosts paves the way for subsequent defector takeover on a static A allele background. Examples of the corresponding typical time series are shown in the electronic supplementary material, figure S3.
We show that the genetic structure of host-parasite interaction influences the outcome of their coevolutionary dynamics. These dynamics in turn determine the structuring of the host population via genetic bottlenecks, which we find to be more numerous towards the MA end of the MA -GFG continuum. Population bottlenecks are known to influence the social behaviour of the host [46], favouring cooperation.

(d) Spatial structure
The assumption that individuals encounter one another at random in well-mixed populations is made primarily for mathematical simplicity and tractability [47,48]. In natural populations, however, this assumption is misplaced. Biological populations are largely broken up into subpopulations linked loosely by migration [49,50], known as a metapopulation. Within a single patch, cheats have a selective advantage and cooperation inevitably breaks down [51][52][53]. However, in the context of a metapopulation, cooperating demes grow to higher densities than cheating demes, enabling cooperation to persist across a spatially structured population [51][52][53].
In order to interpret our results in a metapopulation context, we explore how the dynamics within a patch is altered after a single migration event in a hostevolution-only scenario (as in §3a). Each combination of arriving and resident genotype is tested, for different values of the migrant-to-resident ratio and of the time of arrival. See the electronic supplementary material, S5 for details and figures. In contrast to the classical invasion hierarchy between cooperators and defectors, we find that rare cooperators can successfully invade locally abundant defectors if they carry the favourable resistance allele. Of course, the converse also holds, defectors can readily invade cooperators if they carry an additional resistance advantage (see the electronic supplementary material, figure S5). If, however, both are susceptible, resident cooperators have a high chance of resisting an invasion (within our observation time window) unless the migrant defectors are very numerous. Regarding relative timings, the later the resistant migrants arrive, the higher the chance that the residents have already evolved a resistant mutation themselves-which neutralizes the competitive advantage of the migrants.
A complete study involving repeated migration and coevolution in a full metapopulation setting is outside the scope of this article. We have, however, examined the two key issues, invasion hierarchies and the longevity of cooperator patches. In figure 4, we demonstrate how the combination of interaction genetics and coevolution can define the longevity of cooperator patches, in the face of defector challenge. Specifically, we find that in the MA limit, defector invasion is much impaired and so cooperator patches are able to persist for longer, increasing the global prevalence of cooperation across a metapopulation. Expanding our theoretical framework to incorporate more complex mechanisms of interaction (e.g. inverse GFG, multi-locus interactions, [14,21]) and non-random host-parasite interactions will present important challenges for future work. Alongside theoretical development, further experimental tests are vital, and also offer important applied considerations.
(e) Applied context Our theoretical analysis has largely been grounded in the biology of bacteria-phage interactions [29]. This . Probability of defector takeover across the MA -GFG continuum (MA q ¼ 0, GFG q ¼ 1) in the full coevolution and cooperation model (see § §2c and 3c) for different costs of generalism z. We plot the fraction of runs where defectors grow to be more than 99 per cent of the initially cooperator-only host population, within a 2000 h time window of observation and subject to non-extinction of the system (more than 85% of runs are stable in this case). Parameter values are runs per set, other parameters as in figure 2. Initial conditions: V a ¼ 10 6 ml 21 , C a ¼ 10 2 ml 21 , V A ¼ C A ¼ D a ¼ D A ¼ 0. Cooperation is maintained for longer times under a MA model of specificity.
interaction has enormous public health interest [54], due to the clinical importance of bacterial pathogens and the increasing interest in the use of 'phage therapy' as a novel mechanism of pathogen control [55,56]. Electronic supplementary material, figure S4 highlights a concern for phage therapy-the administration of a lethal phage virus does not always result in the eradication of the target bacterial population, even in the absence of spatial structure. For much of the parameter range, the bacterial population evolves resistance to the phage and continues to persist, despite the ability of the phage to co-evolve and chase the bacteria through genotype space.
We go on to show that the nature of this coevolutionary persistence has great significance for the fate of bacterial cooperation, defining the ability of cooperators to resist local replacement by cheats. Bacterial cooperation typically involves the secretion of shared extracellular enzymes, toxins and polymers-important virulence factors when released within a human host [45,57]. Therefore, the coevolutionary interaction with phages may have an important impact on the virulence of bacterial pathogens. In the MA limit, diversifying selection on bacterial resistance alleles will pose a formidable barrier to the invasion of cheats from rare, increasing the resistance of resident wild-type bacteria to recently proposed mechanisms of pathogen control via the introduction of 'Trojan Horse' cheater lineages [58,59].
A recent study has demonstrated that coevolutionary arms races between bacteria and phage decelerate over time [43], giving way to fluctuating selection. Moreover, bacteria-phage coevolution within natural soil communities follows an FSD [60]. These results suggests that MA-type dynamics, characteristic of fluctuating selection, are likely to be important in determining the phenotypic properties of parasite and host populations. Given the positive effect on cooperation at the MA extreme, this highlights the potential of our model to help explain the maintenance of cooperation in natural populations in a broader context.