|Home | About | Journals | Submit | Contact Us | Français|
Coevolution is relentlessly creating and maintaining biodiversity and therefore has been a central topic in evolutionary biology. Previous theoretical studies have mostly considered coevolution between genetically symmetric traits (i.e. coevolution between two continuous quantitative traits or two discrete Mendelian traits). However, recent empirical evidence indicates that coevolution can occur between genetically asymmetric traits (e.g. between quantitative and Mendelian traits). We examine consequences of antagonistic coevolution mediated by a quantitative predator trait and a Mendelian prey trait, such that predation is more intense with decreased phenotypic distance between their traits (phenotype matching). This antagonistic coevolution produces a complex pattern of bifurcations with bistability (initial state dependence) in a two-dimensional model for trait coevolution. Furthermore, with eco-evolutionary dynamics (so that the trait evolution affects predator–prey population dynamics), we find that coevolution can cause rich dynamics including anti-phase cycles, in-phase cycles, chaotic dynamics and deterministic predator extinction. Predator extinction is more likely to occur when the prey trait exhibits complete dominance rather than semidominance and when the predator trait evolves very rapidly. Our study illustrates how recognizing the genetic architectures of interacting ecological traits can be essential for understanding the population and evolutionary dynamics of coevolving species.
Coevolution—reciprocal evolution in interacting species—has been considered a driving force for creating and maintaining biodiversity . Evolutionary ‘arms races’ in predator–prey, parasite–host and exploiter–victim systems [2,3] are thought to create an endless evolutionary chase called Red Queen dynamics [4,5]. In addition to accumulating empirical evidence of coevolution [1,5], recent studies have shown that evolution can be rapid enough to affect contemporary population dynamics [6–8]. Understanding ‘rapid coevolution’ is therefore important for understanding and predicting future ecological dynamics [9–11].
Theoretical studies on coevolution have mostly considered genetically symmetric traits in antagonistic interactions (i.e. coevolution either between continuous quantitative traits or between discrete Mendelian traits) . Continuous quantitative traits are generally modelled using phenotypically based approaches for trait evolutionary dynamics, partly because the genetic basis of ecologically important traits was largely unknown . The approaches include approximate quantitative genetics models for multilocus trait dynamics (e.g. [10–17]), mutation-limited asexual clonal models (Adaptive Dynamics) (e.g. [18–24]) and models for evolutionarily stable strategies (ESS) that cannot be invaded by other strategies (e.g. ). These studies have assumed that the coevolving traits are quantitative with a continuum of possible values . Coevolution between Mendelian predator and prey traits has also been studied (e.g. [27–31]). Indeed, there is a large theoretical literature going back to classical population genetics about exploiter–victim coevolution with Mendelian traits (reviewed by Seger ).
How do different genetic architectures affect coevolutionary dynamics? One important difference is that Mendelian models have an inherent tendency to cycle, first noted by Haldane . Multiple alleles, multiple loci or multiple species versions of these models readily generate complex cycles or chaotic dynamics [32,33]. A comparison of models with two versus more alleles suggests that cycles would be very likely when prey and predator have continuously varying quantitative traits, but in fact the reverse is true : coevolution of quantitative traits rarely produces cycles. The key difference is that in quantitative trait models, all phenotypes are near the population mean . Cycles in Mendelian models involve switches between extremes: high frequency of one predator type (feeding on its complementary prey genotype) causes increased frequency of the least complementary prey genotype. By contrast, in quantitative trait models, trait distributions change through gradual shifts in the mean of a Gaussian distribution. Coevolution-driven cycles can occur in quantitative trait models [14,19], but these require very strong external stabilizing selection for intermediate trait values . With external stabilizing selection weak or absent, the only possible type of ongoing dynamics is unbounded runaway selection for extreme traits [14,15].
Coevolution models with discrete Mendelian traits have frequently been used for parasite–host coevolution, whereas models for predator–prey coevolution mainly assume continuous quantitative traits . However, empirical studies have found many ecologically important traits controlled by Mendelian inheritance in prey (e.g. [34–36]) or predator (e.g. [37–39]) species . This suggests the potential importance of coevolution between traits with different genetic bases. To the best of our knowledge, no study has examined the outcome of coevolution between quantitative and Mendelian traits, despite empirical evidence for coevolution of this type . In Southeast Asia, some snakes can eat right-handed (dextral) land snails effectively but are not good at eating left-handed (sinistral) snails because of left–right asymmetry in the snakes [36,41]. Within the right-handed snake range, left-handed snail species frequently evolved by speciation from a right-handed ancestor, probably as a counter-adaptation against the biased predation [36,42]. Handedness of snails is determined by two alleles at a single nuclear locus [36,43], whereas the snake traits (teeth number on each side and foraging behaviour) are continuous quantitative traits  and therefore they seem to be controlled by many loci with small effects. These observations indicate antagonistic coevolution between quantitative and Mendelian traits. Will the dynamics resemble Mendelian–Mendelian coevolution, quantitative–quantitative coevolution, something in between, or are some new phenomena likely to occur?
In this study, we examine theoretically coevolution between a quantitative predator trait and a Mendelian prey trait. Our model is analogous to matching-allele interaction, in that predators are most able to capture prey whose phenotype matches their own (i.e. a bidirectional axis of vulnerability sensu ). We consider both purely coevolutionary dynamics (trait evolution that does not interact with population dynamics) and eco-evolutionary dynamics (traits and population sizes change on the same time scale). Analysis of our models shows that: (i) continuous-discrete coevolution can produce a rich array of dynamics including bistability (initial state dependence), anti-phase cycles, in-phase cycles and chaotic dynamics, (ii) semidominance of the prey trait is more favourable for predator persistence than complete dominance and (iii) increasing heritability of the predator trait can result in deterministic predator extinction. We also present some results for the reverse situation, Mendelian predator and quantitative prey traits.
We consider a predator–prey model with phenotype matching (a bidirectional axis) instead of a unidirectional axis of vulnerability (e.g. speed, toxin or armour: [10,12,13,16,17,44]). Phenotype matching means that predators are most successful at capturing prey when their trait value x is close in value to the prey trait value θ, and unsuccessful if |x−θ| is large. Examples of phenotype matching include habitat choice, pattern matching and lock-and-key mechanisms [9,13–15,18–20,45,46]. We assume that per capita attack rate on prey genotype i declines in a Gaussian manner when x differs from θi, Here αi is the maximum per capita attack rate, and τ determines how steeply the rate decreases as predator–prey mismatch increases .
The predator trait distribution is Gaussian, with mean and constant phenotypic variance σ2: We assume (i) that phenotypic variance has a genetic (heritable) component, and an environmental component, with and for simplicity, (ii) that the variance components are constant, which is commonly assumed in dynamic models for quantitative trait evolution [10–12,14,26].
Under our assumptions, the average attack rate by a predator on prey genotype i is
where i = 1, 2 and 3 correspond to prey genotypes A1A1, A1A2 and A2A2, respectively. The model represents complete dominance with two discrete prey trait values when θ1 = θ2 ≠ θ3. In semidominance, the heterozygotes have an intermediate trait θ2 = (θ1 + θ3)/2 so there are three distinct prey phenotypes. We assume that θ1 = –θ3 without loss of generality, so that θ2 = 0 under semidominance. The above models represent ends of a spectrum. There is a continuum between these extremes and probably some sort of smooth transition between their dynamics.
Our first model considers purely evolutionary dynamics driven by relative fitness, with total population sizes held constant as in . Let p be the A1 allele frequency. We assume prey genotype fitnesses (in discrete time, with discrete generations) to be where is prey generation time, r is the reproduction rate (equal for all genotypes), P is the predator density and defined by equation (2.1). Thus, the only advantage of one prey phenotype over another is predation risk. We assume that prey start each generation in Hardy–Weinberg equilibrium owing to random mating. To eliminate dynamic behaviours resulting only from the assumption of discrete generations, we pass to a continuous-time approximation in the usual way, by assuming small changes per generation ( 1) owing to weak selection (e.g. , §3.3). In the electronic supplementary material, appendix S1, we show that for small , the prey evolutionary dynamics are approximated by a continuous-time equation determined by the additive genetic variance and fitness gradient,
where is the per-predator attack rate or average predation rate,
We assume that predator fitness is proportional to per-predator total predation rate, which is determined by the prey abundance N, prey genotype frequencies, and the attack rate on each genotype, We then adopt a conventional fitness gradient model 
The speed of adaptive evolution in the predator quantitative trait is determined by the strength of selection and additive genetic variance . The assumption that the trait remains normally distributed was demonstrated to be a good approximation to trait dynamics in various multilocus models numerically . Equations (2.2) and (2.3) are our two-dimensional model for purely evolutionary dynamics. The attack rate is simplified by rescaling predator and prey traits relative to and defining With these, and dropping the prime, the attack rate is The predator's additive genetic variance on this scale is
Our second model allows evolutionary and ecological dynamics (changes in predator and prey abundance) to occur on the same time scale. Specifically, we put the predation rates from our evolutionary model into a Rosenzweig–MacArthur-type predator–prey model where the prey has logistic growth, and the predator has a linear (Holling type I) functional response. The dynamics of prey density, N, and predator density, P, are then
where r is the prey intrinsic growth rate, K is the carrying capacity, d is the per capita mortality and is the average predation rate, The combination of logistic growth and linear functional response means that with constant parameters, an equilibrium at which predator and prey coexist is always stable . The predator equation often includes conversion efficiency, but we can scale prey abundance so that this equals one.
The fitnesses implied by equation (2.4) are consistent with our evolutionary model (equations (2.2) and (2.3)). Prey reproduction (r in equation (2.4)) has no effect on equation (2.2) as long as it is the same for all genotypes, so the logistic density-dependence in equation (2.4) does not change the prey trait dynamics. Predator fitness (dP/dt)/P in equation (2.4) equals the per-predator attack rate, minus the constant d, so the fitness gradient and trait dynamics are still given by equation (2.3).
Assuming that all prey genotypes are equivalent except for the trait value θi (α = α1 = α2 = α3), we rescale the model as N′ = N/K, P′ = P/K and t′ = rt. With the new parameters d′ = d/r, and the (double) prime dropped, the full system becomes
with given by Thus, the model has four parameters (d, α, V, θ1) after rescaling.
With predator and prey abundance constant, we can eliminate more parameters by rescaling in the last two lines of equation (2.5), giving (with the asterisk dropped)
where and This leaves only two parameters, the scaled predator genetic variation (H) and prey trait difference (θ1).
We found that the complete dominance model is stable when the prey trait θ1 is smaller than the threshold (), and is unstable otherwise, by local stability analyses (figure 1a; the mathematical analysis is in the electronic supplementary material, appendix S2). The speed of predator evolution, H, does not affect local stability and the prey trait distance is the sole determinant of stability. At the stable equilibrium, the predator trait is between the prey homozygotes (), because the prey phenotypes are similar enough for predators to exploit both of them simultaneously.
By contrast, the parameter region in which the semidominance model is stable decreases gradually as θ1 increases (figure 1b). Thus, stability results from simultaneously increasing predator evolution speed and decreasing the prey trait difference (figure 1b). When the prey trait difference is smaller than about 1.12, slowing predator evolution can create limit cycles through a Hopf bifurcation (a black curve in figure 1b). When the predator evolution speed is smaller than about 0.628, increasing the prey trait difference also results in a Hopf bifurcation. This is illustrated in figure 2, using parameter values shown in figure 1b (black points with H = 0.2). We used an R  implementation of the pplane program (math.rice.edu/~dfield/dfpp.html) to compute and plot fixed points, nullclines, phase arrows, trajectories and local stable and unstable manifolds for saddle points (see the electronic supplementary material for our pplane R scripts). Increasing the prey trait distance first changes the equilibrium from a stable node to a stable focus (figure 2a). Then, a Hopf bifurcation occurs and the system shows a limit cycle (figure 2b). The cycle amplitude at first increases with the prey trait distance, but as prey allele frequency is bounded (0 ≤ p ≤ 1), trajectories spend longer periods of time near p = 0 and 1 (figure 2c,d). The dynamics in figure 2c,d look like heteroclinic (saddle-to-saddle) trajectories, but in fact they are more extreme versions of figure 2b: the unstable manifolds of the saddles on the boundary (red curves) quickly converge onto a limit cycle that passes very close to the saddles, approaching one saddle almost along its stable manifolds (the boundaries p = 0 and 1), then leaving almost along the unstable manifold and approaching the other saddle. This process repeats ad infinitum, so that the prey allele frequency has recurrent jumps between p = 0 and 1 (figure 2c,d).
As the prey trait distance increases further, the unstable focus changes to an unstable node (figure 2c), and then (when an eigenvalue passes through zero) the node becomes a saddle and at the same time splits into three equilibria (figure 2d; electronic supplementary material, appendix S2). The two new equilibria (figure 2d) are initially unstable when H is small, but they are locally stable when H is large and predators evolve rapidly. In that case, increasing the prey trait difference results in a bifurcation to bistable equilibria (figure 3a) and then to bistable limit cycles (figure 3b; electronic supplementary material, figures S1 and S2: a dashed black line and grey curve in figure 1b). Further increase in the prey trait difference causes the two bistable cycles to grow and then merge.
Bistability means that the outcome of coevolution depends on the initial condition (solid and dotted lines in figure 3). In the semidominance model (figure 1b), it is possible for rapid predator evolution to constrain predator trait dynamics between the heterozygote () and either one of homozygotes ( or –θ1) depending on the initial state, when predator evolution is rapid (H is large) and the prey trait difference is intermediate (θ1 is around 1.2), as follows. If the predator mean trait is near one of the homozygote phenotypes, heterozygotes increase in frequency in response to predation. Once heterozygotes are common enough, the predator trait evolves rapidly (owing to large H) towards the heterozygote phenotype (); as a result, the predator phenotype is closer to zero than the prey mean phenotype, and this reverses the direction of prey evolution back towards the original homozygote (figure 3). When the prey trait difference is small, predators can exploit the three prey phenotypes simultaneously and the system is stabilized. When the prey trait difference is large, the homozygotes are so different from the heterozygote that it becomes difficult for the predator trait to get ahead of the prey phenotype. Thus, large H and intermediate θ1 are necessary for bistability to occur.
By local stability analyses (electronic supplementary material, appendix S3) and numerical simulations, we found that coevolutionary dynamics in the four-dimensional eco-evolutionary system shows stable equilibria, limit cycles, bistable equilibria, bistable cycles or deterministic predator extinction depending on the prey trait (θ1) and predator mortality (d) (figure 4). Limit cycles are a consequence of coevolution, because the ecological subsystem is always stable with prey logistic growth and predator linear functional response (electronic supplementary material, appendix S3). The bifurcations of the complete dominance model are not affected by the speed of predator evolution, V, but increasing V stabilizes the semidominance model (electronic supplementary material, figure S3).
Deterministic predator extinction can occur when the predator equilibrium density becomes negative at the internal equilibrium (‘extinction’ regions in figure 4). Here the prey phenotypes are so different that predators cannot exploit them to invade with the intermediate trait () even when the prey density is at its equilibrium. In addition, coevolution can drive deterministic extinction with distinct prey phenotypes, as gradual predator evolution cannot keep up with jumping prey evolution (figure 5b: ‘extinction by coevolution’ regions in figure 4). This occurs because of the genetic asymmetry between prey and predator: while predator evolution is gradual, the prey trait distribution can shift from one extreme phenotype to the other by changing allele frequencies. Therefore, the complete dominance model (figure 4a) causes predator extinction more easily than the semidominance model (figure 4b) because the heterozygous phenotype works as a bridge between the two extreme phenotypes of homozygotes in the semidominance model.
Predator extinction depends on initial conditions as well as the predator trait heritability, Increasing the heritability (and the speed of predator adaptation, H and V) tends to stabilize purely coevolutionary cycles (figure 1; electronic supplementary material, figure S3), but it can also result in deterministic predator extinction (figure 5). When the speed of evolution is slow, the predator stays in between the prey homozygotes and behaves as a generalist that exploits both homozygotes simultaneously (figure 5a). However, fast evolution can cause overreaction of the predator trait to specialize on one of the homozygotes. After the prey approaches one extreme (homozygote) phenotype, the predator phenotype follows and gets close enough to the homozygote, and then the other prey allele starts to become common. The predator continues to evolve in the direction of the mean prey phenotype until the mean prey phenotype (evolving now in the opposite direction) crosses the predator phenotype. As the prey phenotype escapes, predator density drops until the predator no longer affects prey evolution and prey allele frequency stops changing. At this point, the homozygote that was previously common is rare, but it still exists. The predator is then trapped in a local basin around the now-rare homozygote (an evolutionary trap, sensu ) and cannot return to the intermediate trait value (figure 5b). Because predator abundance is a continuous variable it never becomes exactly zero in our model, so the predator trait continues to evolve even when predator density is extremely low (figure 5b). Therefore, increasing the speed of predator evolution expands the ‘extinction by coevolution’ regions in the phase diagrams (figure 4; electronic supplementary material, figure S3).
Adding population densities to the purely evolutionary model allows a greater range of dynamics to occur including anti-phase cycles (figure 6a), in-phase cycles (figure 6b) and chaotic dynamics (figure 6c). Anti-phase cycles are characterized by out-of-phase oscillations in predator and prey densities (i.e. predator maxima coinciding with prey minima and vice versa), and arise when the predator mortality and prey trait difference are large. In-phase cycles have a longer cycle period and occur when predator mortality and prey trait difference are both small (black points in figure 4b). Chaotic dynamics arises when the attack rate and heritability are large enough (electronic supplementary material, figure S4).
Our results show that genetic asymmetry is potentially important for understanding and predicting ecological and evolutionary dynamics. We explored the consequences of antagonistic coevolution between quantitative predator and Mendelian prey traits. As we did not include stabilizing selection towards an intermediate optima [14,45] or the Holling type II functional response, a continuous–continuous analogue of our model will result in either stable equilibrium or runaway escalation towards infinite trait values . By contrast, we found a wide variety of possible dynamics in discrete–continuous coevolution, including coevolutionary cycles (figure 2), bistable equilibria and bistable cycles (figure 3) in the purely evolutionary system (figure 1). When we couple population and evolutionary dynamics, deterministic predator extinction (figure 5), anti-phase cycles, in-phase cycles and chaotic dynamics (figure 6) are also possible (figure 4). Therefore, continuous–discrete coevolution is more similar to discrete–discrete coevolution in terms of stability, as it can produce coevolutionary cycles without stabilizing selection. Furthermore, the rich dynamics from simple models demonstrate that traits' genetic architectures can determine the dynamics of coevolution [15,45,54].
We found that semidominance can promote trait cycling when the prey trait difference is small, stabilize dynamics when the prey trait difference is large (figure 1), and prevent predator extinction (figure 4) compared to the complete dominance model, as heterozygotes work as a bridge between the two extreme homozygotes. Interestingly, predators mostly exploit heterozygotes in stable coexistence, and heterozygotes are maintained by subsidy from homozygotes via sexual reproduction (phenotypic subsidy) . As stable coexistence of three asexual genotypes sharing the same predator and resource is impossible , sexual reproduction is an important factor for coevolution with semidominance.
Increasing the predator trait heritability (h2) intuitively appears beneficial for the predator because it accelerates adaptive responses to prey phenotypic changes. Indeed, increasing heritability promotes population persistence under directional environmental change . However, larger heritability can cause the predator trait to overreact to prey frequency change, and eventually results in predator extinction (figure 5). This phenomenon was previously discussed in the context of periodic abiotic environmental change ; evolution of prey traits has a similar effect on the predator in discrete–continuous coevolution.
Using a different parameter rescaling (electronic supplementary material, appendix S4), we can examine the effects of heritability and trait variance separately. In the complete dominance model, the system is always stable when the prey trait difference is small (electronic supplementary material, figure S5a), but large trait variance stabilizes the system otherwise (electronic supplementary material, figure S5b). In the semidominance model, on the other hand, both heritability and variance are necessary for the system to be stable (electronic supplementary material, figures S5c,d). This underlines that the components of trait variation matter, not just the total amount of trait variation [47,55].
Anti-phase cycles were observed in predator–prey microcosm experiments with rapid prey evolution [7,58]. Our findings suggest that in-phase and chaotic cycles might also occur as a result of eco-evolutionary feedbacks in experimental studies. The possibility of in-phase cycles  and chaotic cycles [24,45] of (co)evolving predator–prey systems has been theoretically suggested, but our study illustrates another possible mechanism.
Our model's in-phase cycles (figure 6b) are similar to the canard cycles in models for prey evolution with multiple predators . A canard is a trajectory that spends a long time near an unstable object. In the in-phase cycles, prey allele frequency spends a long time near each extreme (p = 0 or 1) after the direction of selection has shifted to favour the opposite extreme, because of low trait variance, and then suddenly moves to the other extreme. As a result, the lag between traits is larger than that in anti-phase cycles (figure 6).
Our model is a simple representation of coevolution between quantitative and Mendelian traits, not a specific model for snake–snail coevolution mediated by handedness . For example, snail coiling direction is determined by their maternal genotypes , and this delayed inheritance can stabilize coevolutionary cycles (electronic supplementary material, appendix S5, and figure S6). In addition, the reproductive incompatibility between snails with opposite coils  will affect coevolutionary dynamics, because rare genotypes have low fitness ; even with complete dominance, reproductive incompatibility causes bistability . Therefore, the combination of delayed inheritance and reproductive incompatibility in land snails makes cyclic dynamics unlikely compared to our model.
Scale-eating cichlid fish and their prey are a classical example of negative frequency-dependent selection . There are right- and left-handed scale eaters, with ‘handedness’ determined by one locus with two alleles . Because prey fish behaviourally adapt to more frequent phenotypes  and learning can be described by a quantitative trait model , this may be regarded as coevolution between discrete predator and continuous prey traits. This reversed system also results in coevolutionary bistability with the reversed effects of parameters: smaller predator trait difference and larger prey evolution speed destabilize dynamics (electronic supplementary material, appendix S6 and figure S7). As behavioural dynamics are probably fast relative to evolution, the cichlid system may become unstable (electronic supplementary material, figure S7).
Quantitative genetics models for trait change along a fitness gradient are commonly used to model evolution, phenotypic plasticity and even learning by changing the time scale of adaptation . On the other hand, recent studies showed that different mechanisms of rapid adaptation can produce distinct ecological dynamics. For example, phenotypic plasticity, rapid evolution and evolution of plasticity can cause distinct eco-evolutionary dynamics . The number of loci that control ecologically important traits also affects ecological speciation , evolutionary rescue  and coevolution [14,28,45]. Our study suggests that different genetic architectures can result in distinct eco-evolutionary dynamics. Functional genomic approaches (e.g. [34,35]) and further studies on coevolution with different forms of genetic asymmetry such as ploidy , epigenetics  and recombination  may produce insights on contemporary eco-evolutionary dynamics with diverse genetic architectures.
We thank A. Sasaki, N. G. Hairston Jr and M. Hoso for discussion and valuable comments on the earlier manuscript. We also thank two anonymous reviewers and members of the Ellner lab for their helpful comments.
Supplementary R scripts: Dryad http://dx.doi.org/10.5061/dryad.7jq44.
M.Y. conceived the study. M.Y. and S.P.E. constructed and analysed the models. M.Y. wrote the initial manuscript and S.P.E. substantially contributed to manuscript revisions.
We declare we have no competing interests.
M.Y. was supported by the Japan Society for the Promotion of Science (JSPS) Postdoctoral Fellowship for Research Abroad (24-869) and is supported by Hakubi Center for Advanced Research and John Mung Program of Kyoto University. S.P.E. was supported by National Science Foundation (NSF) DEB-1256719 and DEB-1353039.