## Abstract

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.

## 1. Introduction

Coevolution—reciprocal evolution in interacting species—has been considered a driving force for creating and maintaining biodiversity [1]. 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) [5]. 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 [3]. 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. [25]). These studies have assumed that the coevolving traits are quantitative with a continuum of possible values [26]. 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 [32]).

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 [32]. 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 [32]: coevolution of quantitative traits rarely produces cycles. The key difference is that in quantitative trait models, all phenotypes are near the population mean [32]. 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 [15]. 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 [32]. 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 [40]. 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 [36]. 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 [41] 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* [3]). 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.

## 2. Model

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

*α*is the maximum

_{i}*per capita*attack rate, and

*τ*determines how steeply the rate decreases as predator–prey mismatch increases [47].

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
2.1where *i* = 1, 2 and 3 correspond to prey genotypes A_{1}A_{1}, A_{1}A_{2} and A_{2}A_{2}, 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.

### (a) Evolutionary dynamics model

Our first model considers purely evolutionary dynamics driven by relative fitness, with total population sizes held constant as in [14]. Let *p* be the A_{1} 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. [48], §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,
2.2where 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 [26]
2.3

The speed of adaptive evolution in the predator quantitative trait is determined by the strength of selection and additive genetic variance [49]. The assumption that the trait remains normally distributed was demonstrated to be a good approximation to trait dynamics in various multilocus models numerically [50]. 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

### (b) Eco-evolutionary dynamics model

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
2.4where *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 [51]. 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 (d*P*/d*t*)/*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 2.5with given by Thus, the model has four parameters (

*d*,

*α*,

*V*,

*θ*

_{1}) after rescaling.

## 3. Results

### (a) Purely evolutionary dynamics

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)
3.1where 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 1*a*; 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 1*b*). Thus, stability results from simultaneously increasing predator evolution speed and decreasing the prey trait difference (figure 1*b*). 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 1*b*). 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 1*b* (black points with *H* = 0.2). We used an R [52] 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 2*a*). Then, a Hopf bifurcation occurs and the system shows a limit cycle (figure 2*b*). 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 2*c*,*d*). The dynamics in figure 2*c*,*d* look like heteroclinic (saddle-to-saddle) trajectories, but in fact they are more extreme versions of figure 2*b*: 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 2*c*,*d*).

As the prey trait distance increases further, the unstable focus changes to an unstable node (figure 2*c*), and then (when an eigenvalue passes through zero) the node becomes a saddle and at the same time splits into three equilibria (figure 2*d*; electronic supplementary material, appendix S2). The two new equilibria (figure 2*d*) 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 3*a*) and then to bistable limit cycles (figure 3*b*; electronic supplementary material, figures S1 and S2: a dashed black line and grey curve in figure 1*b*). 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 1*b*), 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.

### (b) Eco-evolutionary dynamics

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 5*b*: ‘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 4*a*) causes predator extinction more easily than the semidominance model (figure 4*b*) 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 5*a*). 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* [53]) and cannot return to the intermediate trait value (figure 5*b*). 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 5*b*). 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 6*a*), in-phase cycles (figure 6*b*) and chaotic dynamics (figure 6*c*). 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 4*b*). Chaotic dynamics arises when the attack rate and heritability are large enough (electronic supplementary material, figure S4).

## 4. Discussion

### (a) Comparing discrete–continuous with continuous–continuous coevolution

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 [15]. 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].

### (b) Semidominance versus complete dominance

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) [55]. As stable coexistence of three asexual genotypes sharing the same predator and resource is impossible [56], sexual reproduction is an important factor for coevolution with semidominance.

### (c) Effects of increasing heritability

Increasing the predator trait heritability (*h*^{2}) 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 [8]. 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 [57]; 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 S5*a*), but large trait variance stabilizes the system otherwise (electronic supplementary material, figure S5*b*). In the semidominance model, on the other hand, both heritability and variance are necessary for the system to be stable (electronic supplementary material, figures S5*c*,*d*). This underlines that the components of trait variation matter, not just the total amount of trait variation [47,55].

### (d) Complicated cycles

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 [59] 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 6*b*) are similar to the canard cycles in models for prey evolution with multiple predators [60]. 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).

### (e) Comparison of empirical systems: snake–snail and cichlid systems

Our model is a simple representation of coevolution between quantitative and Mendelian traits, not a specific model for snake–snail coevolution mediated by handedness [36]. For example, snail coiling direction is determined by their maternal genotypes [42], 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 [43] will affect coevolutionary dynamics, because rare genotypes have low fitness [42]; even with complete dominance, reproductive incompatibility causes bistability [42]. 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 [37]. There are right- and left-handed scale eaters, with ‘handedness’ determined by one locus with two alleles [38]. Because prey fish behaviourally adapt to more frequent phenotypes [37] and learning can be described by a quantitative trait model [61], 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).

## 5. Conclusion

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 [61]. 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 [62]. The number of loci that control ecologically important traits also affects ecological speciation [42], evolutionary rescue [63] 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 [31], epigenetics [64] and recombination [65] may produce insights on contemporary eco-evolutionary dynamics with diverse genetic architectures.

## Data accessibility

Supplementary R scripts: Dryad http://dx.doi.org/10.5061/dryad.7jq44.

## Authors' contributions

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.

## Competing interests

We declare we have no competing interests.

## Funding

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.

## Acknowledgements

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.

- Received December 8, 2015.
- Accepted February 24, 2016.

- © 2016 The Author(s)

Published by the Royal Society. All rights reserved.