Coevolutionary transitions from antagonism to mutualism explained by the Co-Opted Antagonist Hypothesis
General framework for eco-coevolutionary transitions from antagonism to mutualismWe develop a general framework in which we model interactions between host species i (density Hi) and its partner species k (density Fk), which are initially purely antagonistic. The model is general, but could be applied broadly to bacterial hosts and parasitic phages or plant hosts and animal or fungal partners, for example. The ecological dynamics of this community (without evolution) are given by:$$frac{d{H}_{i}}{{dt}}={g}_{i}{H}_{i}left(1-mathop{sum}limits_{j}{q}_{{ij}}{H}_{j}right)+{sum }_{k}{f}!_{ik}left[beta left({H}_{i},{F}_{k}right),alpha left({{H}_{i},F}_{k}right)right]$$
(4a)
$$frac{d{F}_{k}}{{dt}}=mathop{sum}limits_{i}{f}_{ki}left[beta left({H}_{i},{F}_{k}right),alpha left({{H}_{i},F}_{k}right)right]-{delta }_{k}{F}_{k}$$
(4b)
The first term of Eq. (4a) describes host population growth in the absence of partner species, where gi is its intrinsic per capita growth rate and qij is the competitive effect of host j on host i for other limiting factors. The general function fik describes the effects of interactions with partner k on host i: β(Hi, Fk) gives the potential mutualism and α(Hi, Fk) describes the antagonism. In Eq. (4b), the general function fki gives the effects of interactions with host i on partner k and δk is the partner’s per capita mortality rate.To derive an explicit eco-coevolutionary model, we apply Equation (4) to model interactions between a single host species and its exclusive partner species (for the sake of simplicity) in terms of host traits xi and partner traits yi (involved in interactions with host i); the ecological dynamics of which are given by:$$frac{1}{{H}_{i}}frac{d{H}_{i}}{{dt}}={g}_{i}left(1-{q}_{i}{H}_{i}right)+frac{bleft[{x}_{i}^{B}right]vleft[{x}_{i}^{V},{y}_{i}^{V}right]{F}_{k}}{{S}_{i}+vleft[{x}_{i}^{V},{y}_{i}^{V}right]{F}_{k}}-hleft[{x}_{i}^{H},{y}_{i}^{H}right]vleft[{x}_{i}^{V},{y}_{i}^{V}right]{F}_{k}$$
(5a)
$$frac{1}{{F}_{k}}frac{d{F}_{k}}{{dt}}=eleft[{y}_{i}^{V},{y}_{i}^{H}right]vleft[{x}_{i}^{V},{y}_{i}^{V}right]hleft[{x}_{i}^{H},{y}_{i}^{H}right]{H}_{i}-{delta }_{k}$$
(5b)
where b is the mutualistic benefits to the host, v is the visitation rate, Si is a saturation constant, h is the costs of antagonism to the host and its benefits to the partner, and e is the partner’s conversion efficiency. The mutualistic and antagonistic interactions are assumed to contribute additively to host population growth and multiplicatively to partner population growth, assumptions that may be valid for many types of interactions, but will not apply universally. To prevent unbounded population growth in the model, the effects of mutualism on host population growth are assumed to saturates with increasing partner density.The function b[xiB] gives the mutualistic benefits of the partner as a function of host trait xiB:$$bleft[{x}_{i}^{B}right]={b}_{{max },i}left(frac{2}{1+{e}^{-{B}_{i}^{{prime} }{x}_{i}^{B}}}-1right)$$
(6a)
where bmax,i gives the maximum mutualistic benefits and ({B}_{i}^{{prime} }) is a saturation constant. The interaction is purely antagonistic when xiB = 0. As xiB increases, the mutualistic benefits b[xiB] increase towards bmax,i.The function v[xiV, yiV] gives visitation rate as a sigmoid function of host trait xiV and partner trait yiV:$$vleft[{x}_{i}^{V},{y}_{i}^{V}right]=frac{{v}_{{max },i}}{1+{e}^{-{V}_{i}^{{prime} }left({x}_{i}^{V}+{y}_{i}^{V}right)}}$$
(6b)
where vmax,i is the maximum visitation rate and ({V}_{i}^{{prime} }) determines how rapidly visitation rate changes as host and partner traits change. As xiV or yiV increase, the visitation rate increases and approaches vmax,i when xiV + yiV → ∞. As xiV or yiV decrease, the visitation rate decreases and approaches zero when xiV + yiV → −∞. Negative values of xiV indicate that the host species is reducing its attraction of the partner species.The function h[xiH, yiH] gives the costs of antagonism to the host and its benefits to the partner, which is described via a sigmoid function of the difference between host trait xiH and partner trait yiH:$$hleft[{x}_{i}^{H},{y}_{i}^{H}right]=frac{{h}_{{max },i}}{1+{e}^{{H}_{i}^{{prime} }left({x}_{i}^{H}-{y}_{i}^{H}right)}}$$
(6c)
where hmax,i gives the maximum antagonism and ({H}_{i}^{{prime} }) determines how antagonism changes as the difference between host and partner traits increases. When xiH > yiH, antagonism declines and approaches zero when xiH – yiH → ∞, while when xiH < yiH, antagonism increases and approaches hmax,i when xiH – yiH → -∞ (unlike xiV, xiH cannot be negative).Partner traits yiV and yiH trade off with conversion efficiency via the function e[yiV, yiH] as defined by:$$eleft[{y}_{i}^{V},{y}_{i}^{H}right]={e}_{{max },i}{e}^{-left({c}_{I,i}^{V}{left({y}_{i}^{V}right)}^{2}+{c}_{I,i}^{H}{left({y}_{i}^{H}right)}^{2}right)}$$
(6d)
where emax,i is the maximum conversion efficiency when interacting with host i (when yiV = yiH = 0), and cI,iV and cI,iH determine how rapidly conversion efficiency declines as yiV or yiH increase, thus quantifying the costliness of traits yiV and yiH, respectively. This trade-off shape was chosen because it is unimodal and constrains conversion efficiency to always be positive. Host trade-offs are defined below (Eq. 8c).Host-partner coevolutionary dynamicsWe model coevolution via the adaptive dynamics framework17,18. Coevolution of a mutant host trait ximut and partner trait yimut (for any general traits xi and yi) is given by:$$frac{d{{x}_{i}}^{{mut}}}{dtau }={mu }_{x}{left.frac{partial {W}_{H}left({{x}_{i}}^{{mut}},{x}_{i},{y}_{i}right)}{partial {{x}_{i}}^{{mut}}}right|}_{{{x}_{i}}^{{mut}}={x}_{i}}$$
(7a)
$$frac{d{{y}_{i}}^{{mut}}}{dtau }={mu }_{y}{left.frac{partial {W}_{F}left({{y}_{i}}^{{mut}},{y}_{i},{x}_{i}right)}{partial {{y}_{i}}^{{mut}}}right|}_{{{y}_{i}}^{{mut}}={y}_{i}}$$
(7b)
where τ is the evolutionary timescale, μx and μy give, respectively, the rates of host and partner evolution, and WH(ximut,xi,yi) and WF(yimut,yi,xi) are the invasion fitness (per capita growth rate when rare) of a mutant host and partner species with trait ximut and yimut in a resident community with trait xi and yi, respectively. The partial derivatives ({left.partial {W}_{H}/partial {{x}_{i}}^{{mut}}right|}_{{{x}_{i}}^{{mut}}={x}_{i}}) and ({left.partial {W}_{F}/partial {{y}_{i}}^{{mut}}right|}_{{{y}_{i}}^{{mut}}={y}_{i}}) are the selection gradients.We model coevolution of mutualistic benefits from the focal partner species (via b), attraction (via v), and defense (via h). The invasion fitness of the mutant host and a mutant partner are given by:$${W}_{H}={g}_{i}left(1-qleft[{{x}_{i}}^{{mut}},{x}_{i}right]{{H}_{i}}^{ast }right)+frac{bleft[{x}_{i}^{B,{mut}}right]vleft[{x}_{i}^{V,{mut}},{y}_{i}^{V}right]{{F}_{k}}^{ast }}{{S}_{i}+vleft[{x}_{i}^{V,{mut}},{y}_{i}^{V}right]{{F}_{k}}^{ast }}-hleft[{x}_{i}^{H,{mut}},{y}_{i}^{H}right]{vleft[{x}_{i}^{V,{mut}},{y}_{i}^{V}right]{F}_{k}}^{ast }$$
(8a)
$${W}_{F}=eleft[{y}_{i}^{V,{mut}},{y}_{i}^{H,{mut}}right]vleft[{x}_{i}^{V},{y}_{i}^{V,{mut}}right]hleft[{x}_{i}^{H},{y}_{i}^{H,{mut}}right]{{H}_{i}}^{ast }-{delta }_{k}$$
(8b)
where Hi* and Fk* are species’ densities at the ecological equilibrium (of Eq. 5). The functions b, v, h, and e are given by Eq. (6a–d), respectively, where xi and yi are replaced with ximut in Eq. (8a) and yimut in Eq. (8b). The function q[ximut,xi] describes trade-offs between mutant host traits and mutant host competitive ability as defined by:$$qleft[{{x}_{i}}^{{mut}},{x}_{i}right]=1+{c}_{H,i}^{B}left({left({{x}_{i}}^{B,{mut}}right)}^{{s}_{i}^{B}}-{left({{x}_{i}}^{B}right)}^{{s}_{i}^{B}}right)+{c}_{H,i}^{V}left({left({{x}_{i}}^{V,{mut}}right)}^{{s}_{i}^{V}}-{left({{x}_{i}}^{V}right)}^{{s}_{i}^{V}}right)+{c}_{H,i}^{H}left({left({{x}_{i}}^{H,{mut}}right)}^{{s}_{i}^{H}}-{left({{x}_{i}}^{H}right)}^{{s}_{i}^{H}}right)$$
(8c)
If ximut > xi for any trait, the competitive effect experienced by the mutant host is increased by an amount taken to be proportional (for simplicity) to the difference between the trait values, ximut – xi, whereas if ximut < xi, the competitive effect experienced by the mutant host is decreased by that amount. The coefficients cH,iB, cH,iV, and cH,iH measure the costs associated with the trade-off for each trait, while the shape parameters siB, siV, and siH define whether the trade-offs are linear (si = 1), concave (si < 1), or convex (si > 1).Mutualism can evolve via the COA for all trade-off shapes (Supplementary Fig. 3). Parameter space plots show that the interaction transitions from antagonism to net mutualism when the costs associated with host traits underlying attraction (cH,iV) and defense (cH,iH) are within a range beyond which there is evolutionary purging of the partner (Supplementary Fig. 3a–c). Only with convex trade-offs can the net antagonism persist. The coevolution of mutualism also requires that the costs associated with partner traits underlying visitation (cFV) and antagonism (cFH) exceed a threshold (Supplementary Fig. 3d–f) below which there is evolutionary purging of the partner (linear or convex trade-offs) or the net antagonism persists (linear or concave trade-offs). Coevolution of mutualism occurs across greater parameter ranges when the trade-offs are linear or slightly concave because costs increase less rapidly than with convex trade-offs.Ecological model of plant-insect interactionsWe tailor the general model (Eq. 4) to model populations of D. wrightii (density Pw) and D. discolor (density Pd) interacting with M. sexta. We scale the model so that Pi = 1 in the absence of M. sexta: thus, Pi >1 indicates that pollination benefits exceed herbivory costs, and Pi < 1 indicates that herbivory costs exceed pollination benefits. The Datura species do not rely obligately on M. sexta and, consistent with ecology of the natural community (Box 1), the model incorporates the alternative host plant, Proboscidea parviflora (density Pp), and the alternative nectar source, Agave palmeri. The ecological dynamics of this community (without evolution) are given by:$$frac{1}{{P}_{i}}frac{d{P}_{i}}{{dt}}=left(1-{P}_{i}right)+frac{{b}_{i}{v}_{i}A}{H+{v}_{i}A}-{h}_{i}{L}_{i}$$
(9a)
$$frac{d{L}_{i}}{{dt}}=varepsilon {e}_{i}{v}_{i}{P}_{i}A-{m}_{i}{h}_{i}{L}_{i}-{d}_{i}{L}_{i}$$
(9b)
$$frac{{dA}}{{dt}}=mathop{sum}limits_{i}{rho }_{i}{m}_{i}{h}_{i}{L}_{i}-{d}_{A}A$$
(9c)
Equation (9a) describes the population dynamics of plant species i (D. wrightii, D. discolor, or P. parviflora). Equation (9b,c) give the dynamics of M. sexta: Li gives the larvae density on plant species i, which recruit into the adult population, A. Pollination is described by the term biviA/(H+viA), where bi is the per capita growth of plant species i due to pollination by the antagonist, vi is the visitation rate to plant species i per antagonist adult, and H is the saturation constant for pollination. Oviposition is given by εeiviPiA, where ei is the oviposition efficiency (number of eggs laid per floral visit) and ε is the fractional increase in egg production due to nectar-feeding at A. palmeri. Floral visits lead to both pollination and oviposition because these behaviors have been shown to be tightly linked in M. sexta19. Pollination and oviposition are given by saturating and linear functions, respectively, based on our data (Supplementary Data 1). Herbivory damage is given by the term hiLi, where hi is the herbivory rate per larvae on plant species i. Larvae mature at rate mihiLi, where mi is the maturation efficiency (fraction of larvae maturing on plant species i). Larval mortality on plant species i is di, adult mortality is dA, and ρi is pupae survival (due to data constraints, we include pupae survival in our estimates of maturation mi, set ρi = 1, and drop ρi from equations hereafter). Equation (9a) gives the dynamics of the alternative larval host plant, P. parviflora (bp = 0 and cannot evolve), which can coevolve attraction and defense. The alternative nectar source, A. palmeri, is incorporated within the model via the parameter ε.Model scalingWithout the antagonist, plant population growth is given by gi (1 – qiPi), where gi is the per capita growth rate of plant species i due to autonomous self-pollination or pollination by other species and qi is plant self-limitation. As qi is very difficult to quantify in nature, we scale the model so that Pi = 1 without the antagonist. We scale plant density ((hat{{P}_{i}}={q}_{i}{P}_{i})), larvae density ((hat{{L}_{i}}={q}_{i}{L}_{i})), herbivory rate ((hat{{h}_{i}}={h}_{i}/{q}_{i})), maturation efficiency ((hat{{m}_{i}}={q}_{i}{m}_{i})), and survival of pupae ((hat{{rho }_{i}}={rho }_{i}/{q}_{i})); where the hats denote scaled quantities and are dropped elsewhere for clarity. Thus, the model is scaled for parameterization, but is not non-dimensionalized. We then scale gi to 1 such that pollination benefits, bi, are estimated by the ratio of the seed set of moth-pollinated flowers to autonomously self-pollinated flowers. Parameter estimates are for scaled quantities.Interaction breakdown boundary for ancestral interaction in a one-plant species communityFor the ancestral insect to persist, its per capita growth rate must be positive when it is rare (i.e., at Pi* = 1, Li* = 0, A* = 0). In stage-structured models, the per capita growth rate is given by the dominant eigenvalue (λD) of the matrix:$$left[begin{array}{cc}-{m}_{i}{h}_{i}-{d}_{i} & varepsilon {e}_{i}{v}_{i}{{P}_{i}}^*\ {m}_{i}{h}_{i} & -{d}_{A}end{array}right]$$which is given by:$${lambda }_{D}=frac{1}{2}left(-{d}_{A}-{d}_{i}-{m}_{i}{h}_{i}+sqrt{{({d}_{A}+{d}_{i}+{m}_{i}{h}_{i})}^{2}-4({d}_{A}left({d}_{i}+{m}_{i}{h}_{i}right)-varepsilon {e}_{i}{v}_{i}{m}_{i}{h}_{i})}right).$$For the insect to persist, λD must have a positive real part, which occurs only when the second term in the square root of λD is negative; i.e., ({d}_{A}left({d}_{i}+{m}_{i}{h}_{i}right)-varepsilon {e}_{i}{v}_{i}{m}_{i}{h}_{i} , , 1). Applying ({f}_{i}=frac{varepsilon {e}_{i}{v}_{i}}{{d}_{A}}) and ({s}_{i}=frac{{m}_{i}{h}_{i}}{{{m}_{i}{h}_{i}+d}_{i}}), where fi is insect lifetime fecundity and si is the larval success (probability of larvae maturing rather than dying), yields Eq. (1).Interaction transition boundary in a one-plant species communityFor the interaction to transition from antagonism to mutualism, equilibrium plant density, Pi* must exceed one (see “Model scaling”). Setting Eq. (9b) to zero and solving for Pi* yields:({{P}_{i}}^{ast }=frac{{{m}_{i}{h}_{i}+d}_{i}}{varepsilon {e}_{i}{v}_{i}}left(frac{{{L}_{i}}^{ast }}{{A}^{ast }}right)). Setting Eq. (9c) to zero and rearranging terms then yields: (frac{{{L}_{i}}^{ast }}{{A}^{ast }}=frac{{d}_{A}}{{m}_{i}{h}_{i}}). Thus, ({{P}_{i}}^{ast }=frac{{{m}_{i}{h}_{i}+d}_{i}}{varepsilon {e}_{i}{v}_{i}}left(frac{{d}_{A}}{{m}_{i}{h}_{i}}right)) and (rearranging slightly) the condition for mutualism to arise is: ({{P}_{i}}^{ast }=left(frac{{d}_{A}}{varepsilon {e}_{i}{v}_{i}}right)left(frac{{{m}_{i}{h}_{i}+d}_{i}}{{m}_{i}{h}_{i}}right) , > , 1). Rearranging and applying ({f}_{i}=frac{varepsilon {e}_{i}{v}_{i}}{{d}_{A}}) and ({s}_{i}=frac{{m}_{i}{h}_{i}}{{{m}_{i}{h}_{i}+d}_{i}}) yields Eq. (2).Interaction breakdown boundary in a one-plant species communityIn the ancestral interaction, insect persistence is evaluated by whether or not it can increase from low density, which yields Eq. (1). Within the net mutualistic region, however, the insect cannot increase from very low density because it cannot buoy plant density sufficiently to maintain a positive per capita growth rate (mathematically, Eq. 1 cannot hold when Eq. 2 is satisfied). The mutualistic region is thus characterized by bistability (see Supplementary Figure 1), and the interaction breakdown boundary is determined by the conditions for the coexistence equilibrium to exist. At the coexistence equilibrium, the larval and adult densities are: ({{L}_{i}}^{ast }=frac{-B+sqrt{{B}^{2}-4{A}_{L}{C}_{L}}}{2{A}_{L}}) and ({A}^{ast }=frac{-B+sqrt{{B}^{2}-4{A}_{A}{C}_{A}}}{{2A}_{A}}), where ({A}_{L}=varepsilon {e}_{i}{{v}_{i}}^{2}{h}_{i}{d}_{A}), ({A}_{A}=varepsilon {e}_{i}{{v}_{i}}^{2}{m}_{i}{{h}_{i}}^{2}), (B=varepsilon {e}_{i}{{v}_{i}}^{2}{m}_{i}{h}_{i}left(frac{1}{{f}_{i}{s}_{i}}+frac{H}{{v}_{i}{m}_{i}}-left(1+{b}_{i}right)right)), ({C}_{L}=varepsilon {e}_{i}{v}_{i}H{d}_{A}left(frac{1}{{f}_{i}{s}_{i}}-1right)), and ({C}_{A}=varepsilon {e}_{i}{v}_{i}{m}_{i}{h}_{i}Hleft(frac{1}{{f}_{i}{s}_{i}}-1right)). For the coexistence equilibrium to exist, either CL and CA must be negative or B must be negative and Li* and A* must be real. CL and CA are negative when fi si > 1, which is Eq. (1) and cannot hold within the mutualistic region because Eq. (2) must be satisfied. However, B is negative when ({f}_{i}{s}_{i}left(left(1+{b}_{i}right)-frac{H}{{v}_{i}{m}_{i}}right) , > , 1), which is approximated by Eq. (3) when the last term is assumed to be small. For Li* and A* to be real, B2 – 4ALCL > 0 and B2 – 4AACA > 0. Assuming that the pollination saturation constant is small (i.e., H ≈ 0) yields CL ≈ CA ≈ 0 such that ({{L}_{i}}^{ast }approx frac{-B}{{A}_{L}}approx frac{{m}_{i}}{{d}_{A}{f}_{i}{s}_{i}}left({f}_{i}{s}_{i}left(1+{b}_{i}right)-1right)) and ({A}^{ast }approx frac{-B}{{A}_{A}}approx frac{1}{{h}_{i}}left({f}_{i}{s}_{i}left(1+{b}_{i}right)-1right)), which are both positive when fi si (1 + bi) > 1 as approximated by Eq. (3).Interaction transition and breakdown boundaries in a two-plant species communityThese boundaries are analytically intractable and are estimated by simulation (see codes provided online).Coevolutionary dynamics of plants and insectThe effects of plant traits xi and insect traits yi on the ecological dynamics of the interactions are given by:$$frac{1}{{P}_{i}}frac{d{P}_{i}}{{dt}}=left(1-{P}_{i}right)+frac{bleft[{x}_{i}^{B}right]vleft[{x}_{i}^{V},{y}_{i}^{V}right]A}{H+vleft[{x}_{i}^{V},{y}_{i}^{V}right]A}-hleft[{x}_{i}^{H},{y}_{i}^{H}right]{L}_{i}$$
(10a)
$$frac{d{L}_{i}}{{dt}}=varepsilon eleft[{y}_{i}^{V},{y}_{i}^{H}right]vleft[{x}_{i}^{V},{y}_{i}^{V}right]{P}_{i}A-{m}_{i}hleft[{x}_{i}^{H},{y}_{i}^{H}right]{L}_{i}-{d}_{i}{L}_{i}$$
(10b)
$$frac{{dA}}{{dt}}=mathop{sum}limits_{i}{m}_{i}hleft[{x}_{i}^{H},{y}_{i}^{H}right]{L}_{i}-{d}_{A}A$$
(10c)
We model coevolution of plant-insect interactions using the adaptive dynamics framework17,18 to link population dynamics and trait coevolution. The coevolution of mutant plant trait xmut and insect trait ymut (for general traits x and y) is given by Equation (7). We model the coevolution of pollination benefits from the antagonist, bi (via mutant plant trait xiB,mut), attraction (via mutant plant trait xiV,mut and mutant insect trait yiV,mut), and defense (via mutant plant trait xiH,mut and mutant insect trait yiH,mut). The invasion fitness of a mutant plant is given by:$${W}_{P,i}left({x}_{i}^{{mut}},{x}_{i},{y}_{i}right)=left(1-qleft[{x}_{i}^{{mut}},{x}_{i}right]{{P}_{i}}^{ast }right)+frac{bleft[{x}_{i}^{B,{mut}}right]vleft[{x}_{i}^{V,{mut}},{y}_{i}^{V}right]{A}^{ast }}{H+vleft[{x}_{i}^{V,{mut}},{y}_{i}^{V}right]{A}^{ast }}-hleft[{x}_{i}^{H,{mut}},{y}_{i}^{H}right]{{L}_{i}}^{ast }$$
(11a)
where Pi*, Li*, and A* are the densities of the plant, insect larvae per plant, and insect adults, respectively, at the ecological equilibrium (of Eq. 10). The functions (b[x_{i}^{B,mut}], v[x_{i}^{V,mut} , y_{i}^{V}],) and (h[x_{i}^{H,mut}, y_{i}^{H}]), describe the effects of mutant plant traits (x_{i}^{B,mut}), (x_{i}^{V,mut}), (x_{i}^{H,mut}), and (x_{i}^{H,mut}) on pollination benefits, attraction, and defense, respectively, which are defined by Eq. (6a–c), where xi is replaced with ximut (where the plant is the host species and the insect is the partner species). The function q[x,mut, xi] defines the trade-offs between mutant plant traits and the competitive ability of mutant plants, which is given by Eq. (8c) (with si = 1). At a coESS, ximut = xi for all traits such that q[ximut, xi] = 1 and the original definition of Pi >1 indicating that pollination benefits exceed herbivory costs is retained when pollination benefits evolve.Invasion fitness of a mutant insect is given by the dominant eigenvalue of its system of equations evaluated at the resident equilibrium. In a one-plant species community, the insect invasion fitness is:$${W}_{I,i}=frac{1}{2}left(-{d}_{A}-{d}_{i}-{m}_{i}{h}_{i}^{{mut}}+sqrt{{left({d}_{A}+{d}_{i}+{m}_{i}{h}_{i}^{{mut}}right)}^{2}-4left({d}_{A}left({d}_{i}+{m}_{i}{h}_{i}^{{mut}}right)-frac{varepsilon {e}_{i}^{{mut}}{v}_{i}^{{mut}}{h}_{i}^{{mut}}{d}_{A}left({d}_{i}+{m}_{i}{h}_{i}right)}{{e}_{i}{v}_{i}{h}_{i}}right)}right)$$
(11b)
where vimut, himut, and eimut are functions describing the effects of mutant insect traits on attraction, defense, and mutant oviposition efficiency, respectively, which are given by Eq. (6b–d), where yi is replaced with yimut. Invasion fitness of a mutant insect in a two-plant species community is given by the dominant eigenvalue of its system of equations evaluated at the resident equilibrium, which is analytically tractable, but sufficiently complicated that we do not include it here (see codes provided online).The curves where the selection gradients (see Eqs. 7) become zero give the evolutionary isoclines for the coevolutionary system. The points where the isoclines intersect give the coevolutionary singularities, which are coevolutionary stable states (coESSs) when they are stable for both plants and the insect. For tractability, the local stability of the coevolutionary singularities was assessed by carefully inspecting the selection gradient of each trait in the neighborhood of its coESS with all other traits held at their coESS as well as by simulating coevolutionary dynamics. Importantly, all three plant traits (xiB, xiV, and xiH) and both insect traits (yiV and yiH) all coevolve simultaneously in the model.Coevolution of the ancestral antagonistic interactionIn the ancestral interaction, pollination by the antagonist is impossible (bi = 0) and thus visitation only contributes to oviposition. From the plant perspective, the selection gradients for attraction and defense in the ancestral interaction are given by:$${left.frac{partial {W}_{P,i}}{partial {x}_{i}^{V,{mut}}}right|}_{{x}_{i}^{{mut}}={x}_{i}}=-{c}_{P,i}^{V}{{P}_{i}}^{ast }$$
(12a)
$${left.frac{partial {W}_{P,i}}{partial {x}_{i}^{H,{mut}}}right|}_{{x}_{i}^{{mut}}={x}_{i}}=frac{{h}_{{max },i}{H}_{i}^{{prime} }{e}^{{H}_{i}^{{prime} }left({x}_{i}^{H}-{y}_{i}^{H}right)}}{{left(1+{e}^{{H}_{i}^{{prime} }left({x}_{i}^{H}-{y}_{i}^{H}right)}right)}^{2}}{{L}_{i}}^{ast }-{c}_{P,i}^{H}{{P}_{i}}^{ast }$$
(12b)
Equation (12a) predicts that selection favors plant traits that reduce attracting the antagonist (e.g., reduced production of volatiles) and lower costs associated with competitive ability. We constrain xiV to be non-negative in the ancestral interaction so that xiV = 0 at the coESS; otherwise, xiV → –∞ and the plant always purges the insect given this model parameterization. Selection balances reduced herbivory damage (first term of Eq. 12b) with costs of reduced competitive ability (second term of Eq. 12b). Selection gradients for insect traits are sufficiently complicated that we do not include them here (see codes provided online); however, selection balances traits that increase visitation and overcome plant defenses with the costs associated with reduced oviposition. The ancestral coESSs are given in Supplementary Table 3.Coevolution of pollination benefits, attraction, and defenseThe evolution of mutant plant traits that allow the antagonist to pollinate it (bimut > 0) initiates the evolution of pollination benefits from the antagonist. The selection gradient for pollination benefits from the antagonist is given by:$${left.frac{partial {W}_{P,i}}{partial {x}_{i}^{B,{mut}}}right|}_{{x}_{i}^{{mut}}={x}_{i}}=frac{2{b}_{{max },i}{B}_{i}^{{prime} }{e}^{-{B}_{i}^{{prime} }{x}_{i}^{B}}vleft[{x}_{i}^{V},{y}_{i}^{V}right]{A}^{ast }}{{left(1+{e}^{-{B}_{i}^{{prime} }{x}_{i}^{B}}right)}^{2}left(H+vleft[{x}_{i}^{V},{y}_{i}^{V}right]{A}^{ast }right)}-{c}_{P,i}^{B}{{P}_{i}}^{ast }$$
(13a)
Equation (13a) shows that plants evolve traits to benefit from floral visits by the antagonist when selection for increased pollination benefits (first term of Eq. 13a) exceeds the costs associated with reduced competitive ability (second term of Eq. 13a).In the model, pollination benefits from the antagonist evolve via Eq. (13a) simultaneously with plant and insect traits affecting attraction and defense. The plant selection gradient for attraction is now:$${left.frac{partial {W}_{P,i}}{partial {x}_{i}^{V,{mut}}}right|}_{{x}_{i}^{{mut}}={x}_{i}}=frac{bleft[{x}_{i}^{B}right]{v}_{{max },i}{V}_{i}^{{prime} }{e}^{-{V}_{i}^{{prime} }left({x}_{i}^{V}+{y}_{i}^{V}right)}H{A}^{ast }}{{left(Hleft(1+{e}^{-{V}_{i}^{{prime} }left({x}_{i}^{V}+{y}_{i}^{V}right)}right)+{v}_{{max },i}{A}^{ast }right)}^{2}}-{c}_{P,i}^{V}{{P}_{i}}^{ast }$$
(13b)
The co-option of the antagonist has fundamentally changed selection on attraction (Eq. 13b vs. Equation 12a), which now balances traits affecting attraction (first term of Eq. 13b) with the costs of reduced competitive ability (second term of Eq. 13b). Co-option of the antagonist also modifies selection on defense (which is still given by Eq. 12b) by changing both trait values and equilibrium densities.Model parameterizationAll ecological parameters are estimated from empirical data. Here we parameterize the saturation constant H, maturation efficiency mi, larval mortality di, and adult mortality dA as well as the parameters for the alternative larval host plant and the alternative nectar source (see “Model validation” for other parameters).We cannot fit the saturation constant H to data because seed set saturates with even a single floral visit. We therefore estimate H as follows: D. wrightii flowers have a 91% chance of setting fruit30; thus, ({v}_{w}A/(H+{v}_{w}A))= 0.91 for a single visit (({v}_{w})A = 1). Solving (1/(H+1))= 0.91 for H yields: H = 0.1. H is assumed to be the same for D. discolor as pollination benefits saturate with a single visit for D. discolor. For maturation efficiency mi, only 0.5% of M. sexta larvae on D. wrightii survive through the final larval instar in nature34; thus, mw = 0.005. As M. sexta suffers 40% lower larval survival on D. discolor (5/8 larvae surviving to pupation) than on D. wrightii (10/10 larvae surviving to pupation) in our experiment19, we estimate that maturation efficiency is ~40% lower on D. discolor than on D. wrightii; i.e., md = (1 – 0.4)mw = 0.003. To estimate larval mortality, we note that larval survival is given by: ({m}_{i}={e}^{-{d}_{i}{D}_{i}}), where Di is development time. M. sexta has a larval stage of ~20 days on D. wrightii35 and there is no difference in development on D. wrightii and D. discolor, at least to the 5th instar19. Solving for di yields: dw ≈ 0.25 and dd ≈ 0.3. Finally, adults live ~5 days in the wild36. Assuming adult mortality is roughly the inverse of the lifespan: dA ≈ 0.2.For the alternative larval host plant, females lay similar numbers of eggs on D. wrightii and P. parviflora34; thus, visitation rate and oviposition efficiency are assumed to be the same as with D. wrightii; i.e., vp = vw and ep = ew. Because P. parviflora plants are of similar size and architecture as D. wrightii34, we assume that herbivory rate on P. parviflora is the same as on D. wrightii; i.e., hp = hw. (see “Model validation” for estimates of vw, ew, and hw). Only 1% of M. sexta larvae on P. parviflora survive through the final larval stage34; thus, mp = 0.01. As larvae have roughly the same development time on P. parviflora as on D. wrightii (~20 days37), solving ({m}_{p}={e}^{-{d}_{p}{D}_{p}}) yields an estimate of larval mortality on P. parviflora of: dp ≈ 0.25.For the alternative nectar source, A. palmeri provides M. sexta with copious amounts of nectar that females likely utilize for egg production38. M. sexta females lay 100–300 eggs/night39. If females foraging exclusively on D. wrightii lay the minimum 100 eggs/night and females that also forage at A. palmeri lay the maximum 300 eggs/night, then A. palmerii is estimated to increase oviposition by a factor of: ε = 3.Model validationPollination benefits (bi), visitation rate (vi), herbivory rate (hi), and oviposition efficiency (ei) all evolve simultaneously in the model. We independently validate the coESSs predicted by the models whenever possible by estimating these parameters using data that were not used to parameterize the models. We estimate bi via the ratio of the seed set of moth-pollinated flowers to autonomously self-pollinated flowers (autonomously self-pollinated seeds germinate as readily as do outcrossed seeds;30). Pollinated D. wrightii and D. discolor flowers set bw = 4.6 ± 0.2 and bd = 3.6 ± 0.1 times more seeds, respectively, than do autonomously self-pollinated flowers (D. wrightii: n = 21 fruit; D. discolor: n = 85 fruit). Moths averaged vw = 4.3 ± 0.6 floral visits to D. wrightii (n = 89 plants) and vd = 2.4 ± 0.4 floral visits to D. discolor (n = 33 plants) in our experiment19. Estimating the herbivory rate is very difficult in nature; however, we can make cursory estimates based on our data. A single M. sexta larvae can consume 1400–1900 cm2 of leaves, which is more than many D. wrightii plants in nature30. Assuming that an average D. wrightii plant supplies larvae with 1400 cm2 of leaves, the variation in leaf consumption (500 cm2) represents ~0.4 plants (=500/1400). Thus, M. sexta larvae are estimated to consume: hw ≈ 1 ± 0.4 D. wrightii plants. M. sexta larvae consumed roughly two times more D. discolor leaf biomass than D. wrightii leaf biomass based on our cursory estimates from our experiments; thus, hd = 2hw ≈ 2 ± 0.8. We estimate oviposition efficiency by the slope of a linear regression of the number of eggs versus the number of floral visits that each plant received from each female moth in our experiments19, which yields: ew = 0.6 ± 0.1 (n = 34 plants) and ed = 0.6 ± 0.2 (n = 24 plants) (Supplementary Data 1).Estimating evolutionary model parametersDirectly estimating evolutionary parameters with data is not possible. We therefore use theory to predict how key parameters affect eco-coevolutionary outcomes and to select reasonable parameter estimates. Our approach is as follows. We set the rates of plant and insect evolution to one (μx = μy = 1); these rates affect the speed of evolution, but not the coESSs. For each trait, we need to estimate the maximum value (bmax,i, vmax,i, hmax,i, and emax,i), the coefficient (({R}_{i}^{{prime} }), ({V}_{i}^{{prime} }), and ({H}_{i}^{{prime} })), and the associated costs (cP,iB, cP,iV, and cP,iH for plant i and cI,iV and cI,iH for the insect). Maximum trait values were chosen to constrain coevolution to a realistic range. We set the coefficients ({R}_{i}^{{prime} }), ({V}_{i}^{{prime} }), and ({H}_{i}^{{prime} }) to one for simplicity because the exact value of any trait x and y are themselves somewhat arbitrary. The costs associated with the traits therefore largely determine the coevolutionary outcomes in the model.We estimate the costs of each trait by systematically varying the costs of plant traits in the one-plant species community given reasonable values for the insect costs and then systematically varying the costs of insect traits while holding plant costs constant at their chosen values (Fig. 5). Parameter space plots show that the interactions transition from antagonism to net mutualism provided that the costs associated with insect traits underlying visitation (cI,iV) exceed a threshold below which the plant and insect engage in an evolutionary arms-race that results in the evolutionary purging of the antagonist (Fig. 5a, b). Only very rarely does the net antagonism persist. We assigned all insect traits a cost of 0.5 (black points in Fig. 5a, b) and then systematically vary the costs of plant traits associated with attraction and defense.Parameter space plots show that interactions transition from antagonism to net mutualism when the costs associated with defense are high relative to the costs associated with attraction (cP,iH > cP,iV); otherwise, coevolution drives evolutionary purging of the antagonist (Fig. 5c, d). When the costs associated with attraction and defense are both fairly high, the net antagonism persists. We assigned values of cP,iH and cP,iV to D. wrightii and D. discolor such that the parameters for D. discolor are closer to the threshold at which evolutionary purging occurs than are those of D. wrightii (Fig. 5d vs. 5c), reflecting the smaller range of ecological parameters over which M. sexta can persist with D. discolor versus with D. wrightii (Fig. 2b vs. 2a). Finally, the costs associated with pollination benefits from the antagonist (cP,iB) must be very high for the net antagonism to persist and we never observed evolutionary purging of the insect within the range of values used (see codes provided online). We assigned values of cP,iB so that pollination benefits to D. wrightii and D. discolor are well below their maximum values. Our estimates of evolutionary parameters are reported in Supplementary Table 2. Evolutionary parameters for P. parviflora are set equal to D. discolor because, in the absence of more information, both species are annual plants that may face broadly similar evolutionary constraints, at least relative to the perennial D. wrightii.Reporting summaryFurther information on research design is available in the Nature Research Reporting Summary linked to this article. More
