Eco‐evolutionary feedbacks between prey densities and linkage disequilibrium in the predator maintain diversity

Diversity occurs at multiple scales. Within a single population, there is diversity in genotypes and phenotypes. At a larger scale, within ecological communities, there is diversity in species. A number of studies have investigated how diversity at these two scales influence each other through what has been termed eco‐evolutionary feedbacks. Here we study a three‐species ecological module called apparent competition, in which the predator is evolving in a trait that determines its interaction with two prey species. Unlike previous studies on apparent competition, which employed evolutionary frameworks with very simple genetics, we study an eco‐evolutionary model in which the predator's trait is determined by two recombining diallelic loci, so that its mean and variance can evolve, as well as associations (linkage disequilibrium) between the loci. We ask how eco‐evolutionary feedbacks with these two loci affect the coexistence of the prey species and the maintenance of polymorphisms within the predator species. We uncover a novel eco‐evolutionary feedback between the prey densities and the linkage disequilibrium between the predator's loci. Through a stability analysis, we demonstrate how these feedbacks affect polymorphisms at both loci and, among others, may generate stable cycling.

In recent years, there has been a growing body of literature devoted to understanding the feedbacks between ecological and evolutionary processes. Such work has demonstrated that these feedbacks can have important reciprocal effects. For example, the evolution of traits affects classical ecological questions, such as when and how species coexist or when interspecific diversity is promoted, and the ecology of species interactions affects classical evolutionary questions, such as how intraspecific diversity is maintained (Yoshida et al. 2003;Duffy and Sivars-Becker 2007;Lankau and Strauss 2007;Sanchez and Gore 2013;Frickel et al. 2016). Up to now, eco-evolutionary feedbacks have predominantly been studied in models developed with a very simple genetic framework, for instance when a trait of a focal species is determined by a single genetic locus (e.g., Pimental 1968;Levin and Udovic 1977;Wilson and Turelli 1986;Gokhale et al. 2013;Schreiber et al. 2018), or by a quantitative trait, for which only the mean responds to selection (e.g., Schreiber et al. 2011;Patel and Schreiber 2015;Cortez and Patel 2017), or by adaptive-dynamics approaches that assume asexual inheritance (e.g., Abrams and Kawecki 1999;Rueffler et al. 2005;Abrams 2006aAbrams , 2006bVasseur and Fox 2011).
Although these and other studies demonstrate the reciprocal effect between ecological and evolutionary dynamics, and provide insight into the resulting consequences, they constrain the possible evolutionary responses by ignoring the complex genetic basis of most quantitative traits (Falconer and Mackay 1996;Charlesworth and Charlesworth 2010). Population genetics studies demonstrate that these complexities can have important impacts on the evolution and fitness of the species by affecting the distribution of the focal trait in multiple ways (for reviews, see Bürger 2000;Johnson and Barton 2005). Thus, it is natural to suspect that these evolutionary processes within one focal species affect its interactions with other species in the community. Hence, to uncover a full picture of the implications of eco-evolutionary feedbacks, we need to dig deeper into details of the population-genetic processes within a species and analyze their impacts on inter-and intraspecific diversity. In our work, we begin this effort by analyzing an eco-evolutionary apparent competition model, in which the predator is evolving in a trait that is governed by two diallelic, recombining loci and determines the interaction with the two prey.
Apparent competition (Holt 1977) is a well-documented community ecological module consisting of two prey species that negatively affect each other through a shared predator ( Fig. 1; for a recent review, see Holt and Bonsall 2017). In particular, by nurturing the shared predator, each species suppresses the other. This phenomenon has been widely observed in nature among a variety of ecosystems, including marine (Schmitt 1987;Menge 1995), tropical (Morris et al. 2004), agricultural (Settle and Wilson 1990;Muller and Godfray 1997), and host-parasite systems (Bonsall and Hassell 1997;Greenman and Hudson 2000).
Furthermore, a number of empirical examples show that intraspecific variation in heritable traits of the predator can influence consumption rates of different prey. For example, in the well-studied three-spine stickleback species (Gasterosteus aculeatus), Robinson (2000) showed that heritable traits such as number of gill raker, mouth width, and body size, influence whether an individual preferentially consumes brine shrimp (Artemia nauplii) or amphipods (Hyallela). In another example, Ehlinger (1990) demonstrates that the length of the pectoral fin in bluegill sunfish (Lepomis macrochirus) determines whether an individual consumes prey living in open or prey living in vegetative habitats. When phenotypes are heritable and influence important ecological interactions, they give rise to the potential for eco-evolutionary feedbacks.
Like many important ecological modules, such as exploitative competition, predator-prey, or intraguild predation, apparent competition has been the subject of theoretical studies to examine the mechanisms and potential impacts of eco-evolutionary feedbacks, when one or more species are evolving in a trait affecting species interactions. Previous studies have developed and analyzed a model in which an evolving predator has a trait governed by a single locus with two alleles that determines its consumption of the prey (Wilson and Turelli 1986) and, additionally, how efficiently it uses what it consumes (Schreiber et al. 2018). Both studies show that, generally, genetic polymorphisms are more likely to be maintained with eco-evolutionary feedbacks, and thus intraspecific diversity is promoted.
On the other extreme, and considering a slightly different scenario, Schreiber et al. (2011) analyzed a model in which the shared predator evolves in a quantitative trait according to Lande (1976) classical Gaussian model in which only the mean phenotype responds to selection, but the genetic and phenotypic variance, as well as the higher moments of the genotype and the phenotype distribution, remain constant. They found that eco-evolutionary feedbacks enable coexistence of the three species, that is, interspecific diversity. In each of these studies, researchers attribute the maintenance of diversity to the following mechanism: as one prey increases in abundance, traits that strengthen the predator's ability to consume that prey are selected for, and this releases the less abundant prey from predation, which in turn generates an opposing selection pressure.
These and other studies were useful for building an understanding of basic mechanisms of eco-evolutionary feedbacks and their impacts. However, their simplified genetics ignores potentially important features of ecologically relevant traits that may often be polygenic. What makes the full dynamics of polygenic traits so difficult to study analytically is the interaction between selection and recombination. Selection generates statistical associations between alleles at different loci, so called linkage disequilibria, and recombination breaks up these associations. The latter may have the advantage of increasing the genetic variance, and thus accelerating the response to selection, but it also imposes a load by breaking up haplotypes carrying favorable allele combinations, thus reducing fitness. Additionally, epistasis, that is, nonadditive interactions among loci on traits or fitness, and differences in locus effects are common (Phillips 2008) and may sometimes elevate or generate intraspecific variation (e.g., Gimelfarb 1989;Nagylaki 1989;Hastings and Hom 1990;Gavrilets and Hastings 1993;Bürger 2000;Hermisson et al. 2003;Turelli and Barton 2006;Pontz et al. 2018).
By assuming a constant variance, the quantitative genetics approach of Schreiber et al. (2011) does not allow for the variation to change as a consequence of evolution. In reality, allele frequencies and linkage disequilibria change in response to selection and therefore also the genetic variance. In fact, numerous population-genetic studies show that these changes can be complex and are difficult to predict analytically (reviewed in Bürger 2000). Although single-locus models admit changes in the variance, they do not account for the complexities of polygenic traits; in particular, the trait can assume only three values, which imposes serious constraints on the evolutionary dynamics.
The goal of this work is to uncover some of the mechanisms and impacts of eco-evolutionary feedbacks when ecologically important traits are governed by multiple loci, and to explore how the underlying genetic architecture influences these mechanisms. In particular, we extend the previous eco-evolutionary work in apparent competition to account for more complex genetic processes by examining the roles of epistasis, differences in locus effects, recombination, and linkage disequilibrium. Due to the analytical complexity of multilocus models, we restrict our study to two loci with two alleles. This model is sufficiently complex to capture these components of genetic architecture, while simple enough to enable analytical results and a mechanistic understanding. We begin by introducing our new eco-evolutionary model that couples Lotka-Volterra population equations with classical two-locus population genetics equations. With this model, we first ask under what conditions is there coexistence among the species (intraspecific diversity)? Second, when does the predator population maintain polymorphisms (intraspecific diversity)? In these analyses, we explore the eco-evolutionary dynamics and show that it exhibits surprising features.

Model
In our eco-evolutionary model, we consider the joint dynamics of the densities N 1 = N 1 (t) and N 2 = N 2 (t) of two prey species (or two genotypes of a single clonally reproducing species), the density P = P(t) of their common predator, and the frequencies of the diploid predators' genotypes producing a quantitative (phenotypic) trait ( Fig. 1). In the absence of predation, prey k experiences logistic growth with intrinsic growth rate ρ k and carrying capacity K k . The predator population growth rate depends on its consumption of prey and each individual's quantitative trait determines the attack rates on the two prey species. Usually, we suppress the dependence of the variables on time t.
The trait value of an individual predator is determined by two recombining diallelic loci. The effects of the alleles for locus one, labeled A 1 and A 2 , are − 1 2 γ 1 and 1 2 γ 1 , respectively, and for locus two, labeled B 1 and B 2 , are − 1 2 γ 2 and 1 2 γ 2 , respectively. We label the four possible gametes A 1 B 1 , A 1 B 2 , A 2 B 1 , A 2 B 2 by i = 1, 2, 3, 4, respectively, and denote the frequency at time t of gamete i by p i (t). The genotype i j consists of the gametes i and j, and the loci contribute additively to the genotypic value g i j , that is, without dominance and epistasis. Thus, for instance, an individual with genotype A 1 B 1 /A 1 B 2 has genotypic value g 12 = (− 1 2 γ 1 − 1 2 γ 2 ) + (− 1 2 γ 1 + 1 2 γ 2 ) = −γ 1 . Here, γ 1 and γ 2 are both positive and, without loss of generality, we assume γ 1 + γ 2 = 1, and so each individual predator has a genotypic value between −1 and 1 (Fig. 1). Because we ignore environmental contributions to the trait, the phenotypic values equal the genotypic values g i j and are given by The attack rate of an individual with phenotypic value g on prey k is a function of g. We assume specifically that it is a quadratic function given by where we impose the constraints θ 1 ≤ −1, θ 2 ≥ 1, and α k > s((−1) k + θ k ) 2 for k = 1 and k = 2, so attack rates are positive and the optimum feasible phenotype for attacking prey 1 or 2 is −1 or 1, respectively (Fig. 1). The parameters θ k , α k , and s determine the shape of the attack rates and the induced trade-offs in efficiency of attacking one or the other prey. In particular, α k is the maximum attack rate on prey k, which occurs only at the optimum θ k (which may be outside the range of feasible phenotypes) and s is the selection intensity. For fixed θ k and α k , smaller s indicates a smaller range of possible attack rates and decreases the trade-off (Fig. 2, solid curves). If θ 2 = −θ 1 is increased while s is decreased such that the range of attack rates stays the same, then smaller s decreases the concavity of the trade-offs (Fig. 2, dashed curves). We use θ k and s to tune the shape of the trade-off and study the effects on the eco-evolutionary feedbacks. Following Schreiber et al. (2011), we assume the fitness of an individual predator depends linearly on the attack rates and Trade-offs between attack rates on each prey. Each curve is a parametric plot of attack rate on prey 1 (a 1 ) versus attack rate on prey 2 (a 2 ) over the genotypic range from -1 to 1. The maximum attack rate is α 1 = α 2 = 1.5, and selection intensity is s = 0.375, 0.3125, 0.25, 0.1785, for red, orange, green, and blue, respectively. For solid curves, the optimum phenotype for attacking prey 1 is θ 1 = −1 and that for prey 2 is θ 2 = 1. Note that as s decreases, the range of attack rates decreases. For dashed curves, −θ 1 = θ 2 = 1, 1.2, 1.5, and 2.1, for red, orange, green, and blue, respectively, so that the range of attack rates stays constant and only the shape of the curve changes. For the dashed curves, trade-offs become less concave as s decreases.
densities of the two prey. Hence, the fitness of a predator with phenotype g is w(g) = e 1 a 1 (g)N 1 + e 2 a 2 (g)N 2 − d , where d is the intrinsic death rate, and e k is the conversion efficiency of prey k, which measures the efficiency of transforming prey consumption into reproductive output. Because the attack rates are concave functions of g, so is their linear combination w(g). This induces negative epistasis in fitness between the loci as well as dominance in fitness within loci (see Supporting Information S1.1 for details). Because the fitness of each genotype depends on the density of the two prey so does the strength of selection. We model the evolutionary dynamics of the predator species by following the standard population genetics approach (e.g., Bürger 2000). We assume that the predator species is diploid, mates randomly with respect to the trait, and the trait exhibits no sex differences. We also ignore random drift and mutation. Therefore, zygotes are in Hardy-Weinberg proportions and it is sufficient to follow gamete frequencies. We write w i j for the fitness of genotype i j, where w i j = w(g i j ) is computed from (3) with g i j as described above.
Approximating the standard population genetics discretetime dynamics by a continuous-time dynamics, as is appropriate if both selection and recombination are weak and time is rescaled (e.g., Nagylaki and Crow 1974), the evolution of gamete frequencies is given byṗ Here, w i = j w i j p j is the marginal fitness of gamete i,w = i, j w i j p i p j is the mean fitness of the predator population, ξ 1 = ξ 4 = −ξ 2 = −ξ 3 = 1, r ≥ 0 is the recombination rate (and not the recombination probability as in a discrete-time model), and D = p 1 p 4 − p 2 p 3 denotes linkage disequilibrium.
For our purpose, it will be more convenient to describe the evolutionary dynamics by the allele frequencies p = p 1 + p 2 of A 1 , q = p 1 + p 3 of B 1 , and linkage disequilibrium D instead of the gamete frequencies. The reverse transformation is (This follows from the definition of D by utilizing that p i = p i j p j .) As is well known, the constraints p i ≥ 0 and Altogether, this yields the coupled eco-evolutionary dynamicsṄ is the predator population mean fitness and is the mean attack rate on prey k. Full expressions with partial derivatives evaluated are given in Appendix A. The first term in the expression forṗ describes direct selection on locus one, which is frequency-and densitydependent because ∂w/∂ p depends on population sizes as well as on p and q (see Appendix A). The second term arises from indirect selection on locus one caused by the association (measured by linkage disequilibrium D) with locus two. The third term is due to epistasis in fitness because with nonadditive interactions among loci, linkage disequilibrium affects mean fitness, whereas without nonadditive interactions, ∂w/∂ D = 0. Forq, an analogous description applies. ForḊ, the first two terms describe the effects of the additive component of selection, the third term is due to epistasis, and the last term quantifies the erosion of linkage disequilibrium by recombination. Importantly, mean fitness w depends on the prey population densities (8) and so selection on the predator's trait feeds back with the ecology.
In this framework, the mean phenotypic value of the trait in the population isḡ and the genetic variance is (e.g., Bürger 2000, p. 49) This shows that for given allele frequencies p and q, positive linkage disequilibrium inflates the variance, and negative linkage disequilibrium deflates it. In our model (1), selection acts on the allele frequencies and linkage disequilibrium to drive the mean phenotype toward an ecology-dependent optimum and then, to reduce the genetic variance. Explicit expressions for the mean attack rates and the mean fitness in terms ofḡ and V g are given in Appendix A.
To understand how eco-evolutionary feedbacks affect interspecific diversity, we first derive conditions to ensure the coexistence of the two prey along with the shared predator (section "Results: Species Coexistence"). Then, we ask how eco-evolutionary feedbacks affect the maintenance of polymorphisms. To do so, we partition our analysis into two parts: in part one (section "Loci in or near linkage equilibrium"), we assume that the recombination rate is large relative to the strength of selection (r s) and in part two (section "Linked loci"), we assume that the recombination rate is small relative to the strength of selection (r s).

Results: Species Coexistence
We begin by analyzing when the two prey species coexist with one another despite the shared predator. Throughout all of our analyses, we assume that the carrying capacities of each prey, K k , are large enough to support the persistence of the predator (d < 2 k=1 e kāk K k ) for all values ofā k possible given the genotypes. We make this assumption because, when there is no predator to evolve, the community dynamics are uninteresting because each prey approaches its own carrying capacity. We find that the two prey coexist if For derivations of these conditions, see Appendix B.
Heuristically, when prey 1 is rare and prey 2 is abundant, alleles A 2 and B 2 will be selected for in the predator, driving the entire population toward extreme genotypes close to 1, that is, allele frequencies p and q to zero. Hence, the first condition ensures that prey 1 has positive growth rate when rare. Analogously, the second condition ensures that prey 2 has positive growth rate when rare. Together these give that the two species coexist provided both prey are initially present.
When the carrying capacities are large, the conditions (12) can be approximated as Hence, the two prey will coexist if the ratio of the intrinsic growth rate of prey 2 to 1 is between the ratios of the predator's respective attack rates at the extreme genotypes. Condition (13) is always satisfied when growth rates are equal and attack rates are symmetric with respect to the predator genotypic values (as in Fig. 1). Furthermore, when the phenotypic optima for attacking prey 1 and 2 are θ 1 = −1 and θ 2 = 1, respectively, then (13) simplifies to which illustrates the dependence on the selection intensity s (Fig. 3). Increasing s strengthens the trade-off between attacking one prey versus the other (solid lines in Fig. 2) and allows for coexistence. Note also that coexistence does not depend on recombination. In the following sections, we assume that condition (12) is met and ask what the effects are on the maintenance of intraspecific diversity, that is, polymorphisms within the predator.
. Stronger selection, that is, stronger concavity of the attack rates and a stronger trade-off, facilitates coexistence. Parameters for maximum attack rates are α 1 = α 2 = 1.

Results: Maintenance of Polymorphisms
To investigate how eco-evolutionary feedbacks affect the maintenance of the polymorphism, we consider two assumptions: in section "Loci in or near linkage equilibrium" the loci are close to or in linkage equilibrium (D ≈ 0), and in section "Linked loci" they are linked with some recombination rate, typically 0 < r < s.
In each case, we begin by assuming both prey are symmetric. By this, we mean that the two prey are equivalent in the absence of the predator and the attack rates of the predator on the two prey are symmetric with respect to the predator genotypes, that is, As a consequence, the attack rates satisfy a 1 (−g) = a 2 (g) (Fig. 1). We make this assumption to avoid prey differences so we can focus on the effects of the predator's genetic dynamics. This assumption also enables the derivation of clear analytical results. Then, we drop this symmetry assumption and study how the results change when the prey have different growth rates.
With the assumption of symmetric prey, we find that there is at least one equilibrium point at which all three species coexist and at which both loci are polymorphic. Such an equilibrium is obtained by setting allele frequencies p = q = 1 2 , prey population densities N 1 = N 2 , and solving the equilibrium conditions for N 1 , predator population density P, and linkage disequilibrium D under the assumptions N 1 > 0, P > 0, and − 1 4 ≤ D ≤ 1 4 . The (uniquely determined) expressions for the equilibrium valuesN 1 andP are complicated and given in Supporting Information 1, S1.2, and in the Mathematica notebook in Supporting Information 2, and that for the linkage disequilibriumD is given below in equation (17b). We call this the fully polymorphic, symmetric equilibrium.
To understand the maintenance of diversity, we first analyze the stability of this fully polymorphic, symmetric equilibrium and then investigate the eco-evolutionary dynamics numerically using a numerical solver in Mathematica. The stability analysis provides insight into the eco-evolutionary feedbacks that affect the maintenance of intraspecific diversity.
Local stability of equilibrium points is determined by the stability modulus of the associated Jacobian matrix, s(J ), which is the largest of the real parts of the eigenvalues. When s(J ) < 0, the equilibrium is asymptotically stable, and when s(J ) > 0, it is unstable. Typically, evaluating stability of six-dimensional systems is analytically cumbersome, and only really feasible through computational methods. However, we identified a way to greatly simplify the problem by a change of variables. In particular, we define Therefore, instead of looking at the two prey densities individually, we will examine the total prey density N + and the difference N − between the two prey densities. With this change of variables and the ordering (N + , P, D, p, q, N − ), the symmetric equilibrium becomes andN + is given in Supporting Information 1, equation (S.32). This linkage disequilibrium is generated by epistasis, whose strength is proportional to the total prey density: −2γ 1 γ 2 seN + (eq. S.26 and S.27c in Supporting Information 1). It is negative because the concavity of the fitness function tends to favor genotypes with intermediate phenotypes. Because recombination acts to reduce the magnitude of linkage disequilibrium,D becomes more negative as s/r increases. In particular,D converges to − 1 4 as s/r → ∞ andD converges to zero as r/s → ∞.
With the change of variables, the Jacobian at this equilibrium now takes the simpler form of a block diagonal matrix: Here, J N+ P D is a submatrix of the total prey density, predator density, and linkage disequilibrium, and J pq N− is a submatrix of the two allele frequencies and the difference in prey densities. The exact expressions of the submatrices are given in Supporting Information 1, S1.3, and the derivation in the Mathematica notebook in Supporting Information 2. This change of variables is enlightening for two reasons. First, the block diagonal structure of the Jacobian spotlights two sets of critical eco-evolutionary first-order feedbacks (i.e., how components affect other components) that affect stability: a feedback between (1) total prey density, predator density and, intriguingly, linkage disequilibrium, and (2) allele frequencies at both loci and the difference between the two prey densities (Fig. 4). This suggests that, locally around the equilibrium point, the ecoevolutionary dynamics can be partitioned into these two uncoupled groups of variables and examined separately.
Upon further examination of the submatrices, we observe higher order feedbacks (that is, how components affect the feedbacks between other components). Interestingly, the total prey density, predator density, and the linkage disequilibrium (N + , P, D) influence the feedbacks between the allele frequencies and the prey density difference ( p, q, N − ; Fig. S1). In contrast, the feedbacks between the predator density, total prey density, and the linkage disequilibrium are unaffected by the allele frequencies and the prey density difference. In a subsequent section, we discuss how these higher order feedbacks affect the dynamics.
Second, by reducing a six-dimensional Jacobian matrix into two, more tractable, three-dimensional submatrices (J N+ P D and J pq N− ), we are able to derive analytical results for the stability conditions. Namely, we get two sets of conditions, one associated with each set of feedbacks (s(J N+ P D ) < 0 and s(J pq N− ) < 0). Full expressions of the conditions are given in Supporting Information 1, S1.3. For local stability, both conditions must be met. Additionally, when these two conditions are met, numerical solutions suggest that the equilibrium is globally stable. Conversely, if ei-ther condition is not satisfied, then the equilibrium is unstable and we can determine which of the two eco-evolutionary feedbacks is driving instability. Furthermore, in this case, the bifurcation patterns of J N+ P D and J pq N− provide insight into the ensuing eco-evolutionary dynamics when the equilibrium is unstable.

LOCI IN OR NEAR LINKAGE EQUILIBRIUM
If the ratio r/s of recombination rate and selection intensity is large, as will often be the case for unlinked loci, a population starting in linkage disequilibrium (D = 0) approaches linkage equilibrium quickly. After this short initial phase, it remains close to linkage equilibrium in so-called quasi-linkage equilibrium (D of order s/r ), and its evolutionary dynamics is closely approximated by the much simpler dynamics obtained by imposing D = 0 (Nagylaki et al. 1999). In particular, the stability properties of equilibria for sufficiently large r/s are the same as for D = 0. Under the assumption D = 0, equation (7f) forḊ becomes void, so that the eco-evolutionary dynamics (1) reduces to a five-dimensional system, and the equations for allele frequency change, (7d) and (7e), simplify drastically. This is the case we study in this subsection. As a consequence of the results in Nagylaki et al. (1999), the results derived below for D = 0 apply qualitatively whenever selection is sufficiently weak relative to recombination.
If additionally the prey are symmetric, as we assume in the sequel and is given in (15), then there is a fully polymorphic, symmetric equilibrium point at which all species coexist. It is given by (N + ,P, 0, 1 2 , 1 2 , 0); see Supporting Information 1 and 2 for full expressions. Furthermore, numerically solving for equilibrium points suggests that this is the unique fully polymorphic equilibrium with all three species coexisting. We find that ifD = 0, the stability modulus s(J pq N− ) is always positive, and hence, this equilibrium is always unstable due to the feedbacks between allele frequencies and the difference in prey densities (Supporting Information 1, S1.3.1).
Trivially, there are exactly four equilibria in which both loci are monomorphic, that is, p = q = 0, p = 1 − q = 0, 1 − p = q = 0, and p = q = 1, and in which all three species might coexist (see Supporting Information 2 for full expressions). In addition to these equilibria, there may be other equilibria in which one locus is polymorphic, which we find numerically.
To infer the eco-evolutionary dynamics, we use the following visualization on the two-dimensional allele-frequency plane (( p, q); Fig. 5). First, we use previous results that in the absence of evolutionary trait changes, the solutions of the ecological equations in (1) always approach a unique globally stable equilibrium (Takeuchi and Adachi 1983). For each pair of allele frequencies (and associated average attack rates), we determine the unique globally stable population equilibrium, and note whether it supports all three populations (coexistence) or just one prey and one predator, which are the only possibilities under our assumptions. This divides the allele frequency plane into three distinct regions, which indicate the community that is reached if allele frequencies remain within these bounds (shading in Fig. 3 indicates coexistence).
We then solve for the equilibrium population densities and determine the evolutionary dynamics (arrows in Fig. 5) when the population densities are at equilibrium. From this, we predict the eco-evolutionary dynamics. Although this is a projection of five dimensions onto two, hence may not always give a full picture of the dynamics, numerical solutions of the full dynamics agree with these predicted dynamics. Indeed, as suggested by the twodimensional visualization, we observe numerically that there is a bistability with these two boundary equilibria. In other words, the predator evolves to one of the boundary equilibria. Importantly, this shows that when loci are in linkage equilibrium, polymorphism cannot be maintained at both loci. As outlined above, this extends to weak linkage disequilibrium. This result is not surprising. In the (near) absence of linkage disequilibrium, the evolutionary dynamics follows selection, which acts to increase mean fitness by driving the mean phenotypeḡ to its optimum θN − /N + and, then, reducing the genetic variance V g (see Appendix A eq. A3). At the symmetric equilibrium, the mean phenotype is at its optimum of zero. As shown by (11), this double polymorphism maintains a high genetic variance, especially with weak or no linkage disequilibrium. Therefore, selection will drive allele frequencies away from this equilibrium (and at least one locus will go to fixation), while still maintaining a mean phenotype close to zero. This phase occurs in the shaded regions of Fig. 5A, B).
Whether a polymorphism is maintained at one locus depends on the locus effects. If the loci have equal or similar effects (i.e., γ 1 ≈ γ 2 ), the two stable equilibria are such that both loci are monomorphic and alleles with opposing effects are fixed (Fig. 5A). At these equilibria, both prey have equal or similar density and the optimum phenotype, at or close to zero, is closely achieved by the corresponding double homozygous genotypes A 1 B 2 /A 1 B 2 and A 2 B 1 /A 2 B 1 .
If the locus contributions are sufficiently unequal, the two stable equilibria are such that the locus with the major effect is polymorphic and that with the minor effect is monomorphic (Fig. 5B). The reason is similar as above, except that with sufficiently unequal effects, a mean phenotype close to zero can be maintained only if the locus with major effect remains polymorphic. These results resemble the evolutionary study in Hastings and Hom (1990).
Next, we drop the symmetric prey assumption and through numerical analysis, we examine the maintenance of polymorphisms when prey are asymmetric in intrinsic growth rate (but still within the bounds for prey coexistence given in (12)). We find that there is still one fully polymorphic equilibrium and that it is unstable. The predator evolves to the equilibrium point on the boundary with mean phenotype better suited for consuming the prey with higher growth rate. When the two loci have equal effects, this promotes polymorphisms in one of the two loci (Fig. 5C). Here, which locus maintains polymorphism depends on the initial allele frequencies as well as population densities. If loci contribute unequally to the trait, then another scenario is possible: the predator either evolves to a monomorphic equilibrium or a single-locus polymorphic equilibrium (Fig. 5D).

LINKED LOCI
When recombination is weaker than selection, the fully polymorphic symmetric equilibrium exhibits considerable negative linkage disequilibrium, given by (17b). If the prey are symmetric, that is, assumption (15) holds, then all three species coexist at the fully polymorphic equilibrium (17). We analyze the stability of this equilibrium using both sets of feedbacks and determine the impacts of recombination and of locus effects. Additionally, we compare these stability conditions to those in the absence of eco-evolutionary feedbacks to determine the role of feedbacks on the maintenance of diversity.
We find that when the locus effects are equal (i.e., γ 1 = γ 2 ), the fully polymorphic symmetric equilibrium is always unstable (Fig. 6; Supporting Information 1, S1.3.2), by the same reasoning as in section "Loci in or near linkage equilibrium". Furthermore, the symmetric equilibrium is also unstable if locus effects are unequal but recombination is sufficiently strong relative to selection. Again, the reasoning is similar to that when loci are in linkage equilibrium (section "Loci in or near linkage equilibrium"; Supporting Information 1, S1.3.3) because at the symmetric equilibrium, the optimum phenotype of zero is realized, but with considerable genetic variance. Selection acts to deplete genetic variance, while maintaining a mean phenotype close to zero.
In contrast, if the locus effects are unequal and the recombination rate is below a threshold, the symmetric equilibrium is stable (Fig. 6; Supporting Information 1, S1.3.3). With lower recombination, the average fitness of A 1 B 2 /A 2 B 1 offspring is higher than that of A 1 B 1 /A 2 B 2 offspring, because of negative epistasis due to the concavity of the fitness function w(g) in (3). This leads to the coadapted gene combination A 1 B 2 /A 2 B 1 becoming more prevalent in the population than A 1 B 1 /A 2 B 2 . In this case, strong negative linkage disequilibrium reduces the genetic variance in (11) sufficiently much to allow for stability of this double polymorphism. These findings are analogous to those for the well-understood, purely evolutionary two-locus model of constant quadratic stabilizing selection on a quantitative trait (see Supporting Information 1, S1.1; or, in more detail, Bürger 2000, Chapter 6.2).
The above analysis gives little insight into the ecoevolutionary dynamics when this equilibrium is unstable. Instability occurs via destabilization of either one (or both) of the two sets of feedbacks (i.e., the two conditions s(J N+ P D ) < 0 and s(J pq N− ) < 0). By examining the bifurcation patterns associated with J N+ P D and J pq N− , we determine the eco-evolutionary dynamics close to the bifurcation points at which destabilization is driven by either the feedbacks in total prey density, predator density and linkage disequilibrium (N + , P, D) or allele frequencies and difference between prey densities (p, q, N − ), respectively.
If s(J N+ P D ) > 0 and s(J pq N− ) < 0, then feedbacks between total prey density, predator density, and linkage disequilibrium drive instability (red shaded region in Fig. 6). From the expression (S.34) of J N+ P D in Supporting Information 1, we find that the transition from stability to instability when s(J N+ P D ) = 0 always occurs via a Hopf bifurcation: the focal polymorphic equilibrium is unstable and stable cycles emerge (Supporting Information 2).
Surprisingly, but as suggested by the feedbacks driving the instability, it is the total prey density, predator density, and the linkage disequilibrium that cycle around the equilibrium (Fig. 7). In these cycles, as total prey density N + increases, predator density P follows, as in classical predator-prey dynamics, and additionally, linkage disequilibrium D within the predator population also follows. By examining the phase lags in our numerical solutions, we find that the predator density is a quarter-phase lagged from the total prey density and the linkage disequilibrium is a half-phase lagged from the predator density (Fig. 7). The amplitude of these cycles increases and then decreases with increasing recombination rates in an interesting pattern (Fig. S2).
If parameters are sufficiently near the bifurcation point, the cycles of total prey density, predator density, and linkage disequilibrium (N + , P, D) are sufficiently small that the higher order feedbacks do not impact the stability of the equilibrium for allele frequencies and prey density difference: p, q, N − approach their equilibrium values. Conversely, when the cycles in N + , P, and D are large enough, these higher order feedbacks become destabilizing and drive cycles in p, q, and N − as well (Figs. 8 and S2).
If s(J N+ P D ) < 0 and s(J pq N− ) > 0, then the feedbacks between allele frequencies and the prey densities difference drive instability (white region below red curve in Fig. 6). Numerically, we find that the transition from stability to instability of the polymorphic equilibrium point occurs via a pitchfork bifurcation when parameters are such that s(J pq N− ) = 0. In particular, by numerically solving from various initial conditions, we observe that the fully polymorphic symmetric equilibrium becomes unstable and two new stable equilibria emerge, which are also polymorphic. Furthermore, these two new equilibria remain polymorphic only over a narrow range of parameters. Beyond this narrow range of parameters, the eco-evolutionary dynamics approach one of two alternative equilibrium points in which the locus with smaller effects fixes for one of the two alleles and all three species coexist (Fig. S2). Similar to the case of loci in linkage equilibrium (Fig. 5), the initial state of the community determines the alleles that fix and the ultimate densities of the three species.
Finally, we explore how asymmetries in the prey's intrinsic growth rate affect doubly polymorphic equilibria. Without loss of generality, we assume that prey 1 has a higher growth rate than prey 2. We can already see from the coexistence conditions (1) that when prey 1 has a much higher growth rate, it will exclude prey 2. For small differences, we numerically solve for the nearly symmetric doubly polymorphic equilibrium (N 1 ,N 2 ,P,p,q,D) and determine the critical recombination rate below which this equilibrium is stable.  Figure 6A (γ 1 = 0.7 and r = 0.06). For values of recombination rate r near the Hopf bifurcation, feedbacks cause cycles between total prey density, predator density, and linkage disequilibrium, while allele frequencies approach an equilibrium value ( p → 1 2 , q → 1 2 ) as t → ∞. The difference between the two prey densities converges to zero, that is, N − → 0 (not shown). Cycles are due, in part, to the fluctuating epistasis, which is proportional to the total prey density N + . (B) Numerics show that the predator density and linkage disequilibrium are roughly one-quarter and three-quarters lagged, respectively, from the prey density.
We find that as prey 1 intrinsic growth rate ρ 1 increases, the allele frequency at the major locus,p, increases, and the frequency at the minor locus,q, decreases (Fig. S3a). Predator density P increases as well (Fig. S3d), and linkage disequilibrium D becomes weaker, that is, less negative (Fig. S3c). For small recombination rates, the equilibrium density of prey 1 is larger than that of prey 2 (N 1 >N 2 ), as one might expect when prey 1 has larger intrinsic growth rate. Interestingly, for slightly higher recombination rates, this pattern is reversed: prey 2 is more abundant than prey 1 at equilibrium (N 2 >N 1 ; Fig. S3a). Finally, the nearly symmetric equilibrium points remain stable for higher recombination rates than when both prey are symmetric (Fig. S3a).

Discussion
In our work, we develop a model to further our understanding of how eco-evolutionary feedbacks affect the maintenance of diversity at both the ecological and evolutionary level. Our model couples a two-locus population-genetic model to a Lotka-Volterra ecological model of three interacting species. We study how these feedbacks affect ecological diversity (the coexistence of species) and evolutionary diversity (the maintenance of polymorphism). We show that accounting for multilocus genetics of traits uncovers novel mechanisms for eco-evolutionary feedbacks between the genetics of the evolving species and the species that it interacts with.
Even with this fairly simple model, the nature of feedbacks between evolution and ecology are complex. Our approach to linearize around a double polymorphic equilibrium point enables us to identify the first-order feedback loops, which affect stability of this equilibrium. As our study demonstrates, this is a potentially powerful tool in identifying feedbacks. When prey are symmetric, our block diagonal Jacobian indicates that there are two decoupled groups of feedbacks that affect stability: (1) the total prey density (N + ), the predator density (P), and linkage disequilibrium (D); and (2) the allele frequencies at both loci ( p, q) and the difference between the prey densities (N − ). To understand the feedback in the first group, suppose the total prey density increases from equilibrium. Then, there is a two-fold response (Fig. S4). On the one hand, the predator density increases, which then puts increased pressure on both prey to cause their densities to decrease, as in classical predator-prey cycles. On the other hand, as the total prey density N + increases, the relative fitness difference between individuals with intermediate phenotypes (close to 0) and individuals with extreme phenotypes (close to 1 or -1) increases (right upper loop in Fig. S4). Hence, selection for intermediate gametes (A 1 B 2 and A 2 B 1 ), and therefore negative epistasis gets stronger (epistasis is proportional to total prey density: −γ 1 γ 2 seN + ; Supporting Information 1, S1.1). This leads to more negative linkage disequilibrium D, especially if recombination is low. As a consequence, and shown by (11), the genetic variance decreases and the overall fitness of the predator increases (as is shown by substituting (11) into (A3)); this amplifies the negative impact the whole predator population has on the prey populations. Decreasing N + has the opposite effects and is illustrated in the left lower loop of Figure S4.
In summary, fluctuations between periods of strong and weak selection, caused by changes in total prey density, alter the ge-netic composition of the predator population, which in turn drive changes in prey density. We refer to this as the fluctuatingstrength-of-selection feedback (Fig. S4). This novel feedback highlights the importance of considering genetic architecture in identifying the mechanisms of eco-evolutionary feedbacks, because it is mediated by the linkage disequilibrium induced by epistatic fitness effects of linked loci and frequency-dependent selection acting on them.
The mechanism behind the feedback between the allele frequencies and the difference in prey densities is as follows: suppose we perturb the system at equilibrium by increasing the density of prey 2, hence N − = N 2 − N 1 (Fig. S5). Consequently, the fitness curve shifts: predator individuals with higher phenotypic values are more fit and there is selection for increased phenotypic values. This drives the allele frequencies p and q at both loci to decrease. As they decrease, the predator becomes better at consuming prey 2 and worse at prey 1, which reduces the difference N − in the prey densities. As N − becomes negative, selection for low phenotypic values is generated. We refer to this feedback as the alternating-direction-of-selection feedback (Fig. S5).
Our model predicts that both types of eco-evolutionary feedbacks affect the maintenance of diversity by (1) promoting coexistence of species, (2) promoting the maintenance of genetic polymorphisms, and (3) driving novel eco-evolutionary cycles. In particular, here, we highlight that the alternating-direction-ofselection feedback promotes coexistence: the predator evolves to genotypes better suited for consuming the more abundant species, and this alleviates competitive pressures on the rarer prey species. Previous theoretical and experimental work has recognized this mechanism for promoting coexistence in other types of ecological interactions, including species competing via a shared predator, as in apparent competition (Schreiber et al. 2011(Schreiber et al. , 2018, via direct competition (Lankau and Strauss 2007;, and intraguild predation (Ellner and Becks 2011;Patel and Schreiber 2015).
Although eco-evolutionary feedbacks have been recognized to promote the coexistence of species, there has been significantly less attention to how they influence the maintenance of intraspecific genetic variation. Consistent with previous work on polymorphisms in two-locus models with constant stabilizing selection (Supporting Information 1, S1.1), we show that if selection is weak relative to recombination, the double polymorphism, which always has negative linkage disequilibrium, is unstable (cf. Hastings and Hom 1990). If selection is strong relative to recombination, the double polymorphism is stable, as in the model with constant selection (cf. Gavrilets and Hastings 1993;Bürger 2000). Although selection for intermediate phenotypes when both prey are present facilitates the maintenance of polymorphism, recombination hampers it by breaking apart coadapted gene combinations. In summary, if either selection is weak or strong relative to recombination, the eco-evolutionary feedbacks in our model have no qualitative effects on the stability of the fully polymorphic symmetric equilibrium compared to the purely evolutionary model with constant stabilizing selection.
For intermediate recombination rates, we have identified that the two feedbacks affect the stability of the double polymorphic equilibrium in our model in different ways, which we infer from their effects on the Jacobian. In addition to promoting the coexistence of the two prey (discussed in a previous paragraph), the alternating-direction-of-selection feedback supports stability of the symmetric equilibrium, hence of a full polymorphism, for higher recombination rates than for the purely evolutionary dynamics with fixed population densities (compare the blue curve with the grey dashed curve in Fig. 5A).
This alternating-direction-of-selection feedback resembles a previously recognized mechanism in one-locus and multilocus models that enables the maintenance of polymorphisms: certain forms of externally driven temporal fluctuations in the environment (for example, through seasons), in which the optimal genotype changes in time, have been shown to maintain genetic vari-ation, especially in concert with epistasis or dominance (e.g., Haldane and Jayakar 1963;Ellner and Sasaki 1996;Bürger and Gimelfarb 2002;Wittmann et al. 2017).
In contrast, the fluctuating-strength-of-selection-feedback generates a novel type of eco-evolutionary cycling in apparent competition between total prey density N + , predator density P, and linkage disequilibrium D. Separately, neither the ecological equations of our model (Takeuchi and Adachi 1983) nor the evolutionary equations with either constant fitnesses (Gavrilets and Hastings 1993;Bürger and Gimelfarb 1999) or with negative frequency-dependent selection (Bürger and Gimelfarb 2004) exhibit cycling. The eco-evolutionary cycling we observe emerges from a Hopf bifurcation and is a consequence of the fluctuating epistasis induced by the fluctuating-strength-of-selectionfeedback. Interestingly, cycles occur even if the difference in prey densities N − and the allele frequencies p and q are constant. Cycling of D with constant allele frequencies was observed by Gandon and Otto (2007) in a model of host-parasite evolution going back to Nee (1989). Also in their model, the cycles are caused by fluctuating epistasis. As shown in Fig. 8, in our model all six quantities can also cycle.
In principle, the two types of eco-evolutionary feedbacks our study highlights can be tested empirically. An appropriate empirical system would involve populations whose (1) densities were measurable over time, (2) interactions were influenced by some key phenotype, and (3) there was an identified and tractable genetic basis of this key phenotype, such as a few loci of major effect. With their novel approaches, a number of previous empirical studies have paved the way for testing the predictions from our theoretical study. For example, ecological communities involving species of rotifers and algae in chemostats have been successfully developed and used to identify eco-evolutionary feedbacks, by tracking population densities and phenotypic changes over time (Yoshida et al. 2003;Hiltunen et al. 2014;Kasada et al. 2014). Other studies have investigated how genotype frequencies change over time in response to ecological (rotifer-algae system; Meyer et al. 2006) or environmental fluctuations (Drosophila melanogaster; Bergland et al. 2014). Additionally, in evolutionary genetics, studies are continuing to identify the genetic basis of phenotypes in a variety of organisms (Albert et al. 2008;Lamichhaney et al. 2016;Lawson and Petren 2017). Developing an appropriate and practical empirical system is key to being able to test our predictions.
Two previous theoretical studies examined the ecoevolutionary dynamics in apparent competition, in which the predator was evolving in a single locus. Our results are consistent with results from Schreiber et al. (2018) when one locus has no effect on the genotype (γ 1 = 0 or γ 2 = 0). If, additionally, the intrinsic growth of both prey is sufficiently faster than the predator's consumption, then we recover the model considered in Wilson and Turelli (1986), who pointed out that feedbacks enable the maintenance of a single-locus polymorphism, via a similar mechanism as our alternating-directions-of-selection feedback. There are many real biological examples in which predator and prey have similar generation times, implying that their respective biological rates are on similar time scales. Furthermore, many ecologically relevant traits are complex and the outcome of multiple loci with different effects. Accounting for these additional biological complexities may uncover new mechanisms of ecoevolutionary feedbacks and cyclic dynamics not previously recognized.
Previous theoretical (Rueffler et al. 2005;Schreiber et al. 2011;Patel and Schreiber 2015) and empirical (Meyer et al. 2006;Bolker et al. 2009;Kasada et al. 2014;Frickel et al. 2016) work highlighted that the shape of evolutionary trade-offs has important implications for eco-evolutionary dynamics. Indeed, the trade-offs influence the two feedbacks we identify here. With less concave trade-offs (smaller selection intensity s and larger θ, that is, more extreme optimum phenotypes for attacking prey), changes in the allele frequencies result in a greater change in the difference between the two prey densities. This has a destabilizing effect on the equilibrium via the alternating-direction-of-selection feedback. In contrast, less concave trade-offs have a stabilizing effect via the fluctuating-strength-of selection feedback. By lowering the fitness of the predator, there is a lower potential for relative fitness differences, which leads to only weak selection for intermediate phenotypes and diminishes the feedback on linkage disequilibrium. Altogether, for less concave trade-offs, feedbacks between the allele frequencies and the prey difference become more destabilizing, whereas feedbacks between the predator density, total prey density, and linkage disequilibrium become more stabilizing. Hence, there is an optimal trade-off for stability (Fig. S6).
In this work, we study the simplest possible multilocus model, with two diallelic loci, to account for more complex genetic architecture than in previous studies. We find novel ecoevolutionary feedbacks between genetics of the evolving species and the community that previously were not recognized but are important for the stability of polymorphisms. As many ecologically relevant phenotypes are governed by more than two genetic loci (Albert et al. 2008;Lawson and Petren 2017), an important avenue for future work is to understand how eco-evolutionary dynamics with more than two loci influence the coexistence and the maintenance of polymorphisms. Do we recover similar types of feedbacks when we consider more than two loci? Are there additional novel types of eco-evolutionary feedbacks? Due to the complexity of multilocus dynamics, we suggest that gaining insight into these questions will require computational techniques.

AUTHOR CONTRIBUTIONS
SP and RB designed the study, analyzed the model, and wrote the paper. D = sÑ + [γ 1 (1 − 2 p) + γ 2 (1 − 2q)][γ 1 (1 − 2 p) + γ 2 (1 − 2q) The equations forṗ andq show that sÑ + is a measure of the strength of density-dependent selection. In addition, the last term in the second line forḊ is always negative due to (6b). This reflects the well-known fact that stabilizing selection induces negative linkage disequilibrium through its epistatic effects.
Associate Editor: C. Bank Handling Editor: D. W. Hall

Supporting Information
Additional supporting information may be found online in the Supporting Information section at the end of the article.
Supporting Information 1: Stability analysis of fully polymorphic symmetric equilibrium. Supporting Information 2: Mathematica notebook. Figure S1. Higher-order feedbacks impacting stability of the symmetric double polymorphic equilibrium. Figure S2. Bifurcations from the symmetric equilibrium as a function of the recombination rate. Figure S3. Bifurcation of equilibria for asymmetric prey. Figure S4. Schematic of fluctuating-strength-of-selection feedback between total prey density N + , predator density P, and linkage disequilibrium D. Figure S5. Schematic of alternating-direction-of-selection feedback between allele frequencies p and q and prey-density difference N . Figure S6. Effects of shape of trade-off on stability of the fully polymorphic symmetric equilibrium.