An eco‐evolutionary feedback loop between population dynamics and fighter expression affects the evolution of alternative reproductive tactics

Abstract Surprisingly, little is known about how eco‐evolutionary feedback loops affect trait dynamics within a single population. Polymorphisms of discrete alternative phenotypes present ideal test beds to investigate this, as the alternative phenotypes typically exhibit contrasting demographic rates mediated through frequency or density dependence, and are thus differentially affected by selection. Alternative reproductive tactics (ARTs), like male fighters and sneakers, are an extreme form of discrete phenotype expression and occur across many taxa. Fighters possess weapons for male–male competition over access to mates, whereas sneakers are defenceless but adopt tactics like female‐mimicking. Because fighters in some species mortally injure conspecifics, this raises the question whether fighter expression can feed back to affect population size and structure, thereby altering the selection gradient and evolutionary dynamics of ART expression in an eco‐evolutionary feedback loop. Here, we investigated how the eco‐evolutionary feedback loop between fighter expression and population size and structure affects the evolution and maintenance of ARTs. We introduced intraspecific killing by fighters in a two‐sex, two‐ART population model parameterized for the male dimorphic bulb mite (Rhizoglyphus robini) that includes life‐history differences between the ARTs and a mating‐probability matrix analogous to the classic hawk–dove game. Using adaptive dynamics, we found that the intraspecific killing by fighters can extend the range of life‐history parameter values under which ARTs evolve, because fighters that kill other fighters decrease fighter fitness. This effect can be nullified when benefits from killing are incorporated, like increased reproduction through increased energy uptake. The eco‐evolutionary feedback effects found here for a dimorphic trait likely also occur in other fitness‐related traits, such as behavioural syndromes, parental care and niche construction traits. Current theoretical advances to model eco‐evolutionary processes, and empirical steps towards unravelling the underlying drivers, pave the way for understanding how selection affects trait evolution in an eco‐evolutionary feedback loop.


| INTRODUC TI ON
Understanding the mechanisms of how a phenotypic character distribution of individuals within a population changes over time is a first step towards understanding how the joint dynamics of ecological and evolutionary processes affect populations (Smallegange & Coulson, 2013), of which we still know remarkably little (Hendry, 2016). One type of characters, or phenotypes, that are particularly likely to be influenced by both ecological and evolutionary processes are polymorphisms of discrete alternative phenotypes; examples of which include mating phenotypes in males (major vs. minor), protective phenotypes (armed vs. unarmed) or life cycle phenotypes (single vs. multiple reproductive bouts) (Oliveira, Taborsky, & Brockmann, 2008). The alternative phenotypes typically exhibit contrasting demographic rates: For example, mating and life cycle phenotypes differ in reproductive strategy and output, growth and even survival rates. When the contrasting demographic rates of alternative phenotypes are mediated through frequency and/or density dependence (Oliveira et al., 2008), which can be differentially affected by selection (Sinervo, Svensson, & Comendant, 2000), eco-evolutionary feedbacks are likely to occur. Alternative reproductive tactics (ARTs) are extreme forms of such phenotypic variation within single populations, and arise over evolutionary time when there is high competition for mates (Oliveira et al., 2008). Typically, ARTs occur in one of the sexes, often the male, and are discrete phenotypes with distinct mating tactics (Oliveira et al., 2008). There are usually two phenotypes: fighters and sneakers. In many male dimorphic species, male morph expression is a conditional strategy (Oliveira et al., 2008), so that the specific (permanent) morph a male develops into depends on its condition during a critical point in ontogeny. Whereas the environment to a large extent determines the condition of a male, the threshold response of what morph a male of a certain condition develops into, is in many species under polygenic control and heritable (Oliveira et al., 2008). Fighters typically are large, may mature slowly and possess weapons that are used to obtain and guard mates, whereas sneakers are small and have no weapons, and resort to alternative methods of gaining access to females such as circumventing fighters or maturing early, and consequently mate earlier in life than fighters from the same age cohort.
Currently, the environmental threshold (ET) model is the leading theory to explain the evolution and maintenance of conditionally expressed ARTs (Hazel, Smock, & Johnson, 1990;Hazel, Smock, & Lively, 2004). The ET model is based on the premise that ART expression depends on whether or not an individual reaches a critical threshold, or switch point, during ontogeny. This threshold is assumed to be based on a continuously distributed, polygenic trait, called the "liability," which can be a hormone profile. ART expression occurs in response to a cue (such as body size) that reliably informs on an individual's status within the mating environment. Individuals with a cue value above the threshold express one phenotype, while those below the threshold express the alternative. The ET model assumes that, in response to environment-specific individual-level selection, ARTs have evolved different fitness functions, through which selection can affect the distribution of individual liabilities.
Because ART frequency depends on the distribution of individual liabilities and the cue distribution, both are taken into account in determining how ART fitness influences the evolution of liabilities and hence the evolution of the threshold. Predictions from the ET model on the evolution of the threshold, and thereby ART expression, have been tested successfully in experimental evolution studies (Tomkins, Hazel, Penrose, Radwan, & LeBas, 2011). However, when one such experimental study was repeated, but this time allowing for population feedback by keeping generations overlapping when applying selective mortality regimes (as opposed to taking a non-random sample of individuals of the current generation to start the next), some evolutionary responses to selective mortality were diametrically opposite to the predictions from the ET model (Smallegange & Deere, 2014). If we understand why this mismatch occurred, we would be one step closer to understanding how the joint dynamics of ecological and evolutionary processes, and their interactions, affect population size, growth and persistence. A logical way forward is therefore to investigate how the ecological dynamics of shifts in population size and structure affect the evolution of ART expression by altering the selection gradient, that is the ecology-to-evolution pathway, and how evolutionary change in ART expression in turn affects population size and structure, that is the evolution-to-ecology pathway, within an eco-evolutionary feedback loop (Figure 1a).
This raises the question whether fighter expression can feed back to affect population size and structure, thereby altering the selection gradient and evolutionary dynamics of ART expression in an ecoevolutionary feedback loop. In this study, therefore, we explored to what extent the evolution and maintenance of ARTs is affected by such eco-evolutionary interaction between evolution of fighter expression and (ecological) change in population size and structure. To do so, we employed the adaptive dynamics framework (Geritz et al., 1998;Metz, Nisbet, & Geritz, 1992), which has been instrumental in understanding how eco-evolutionary feedback influences the evolutionary trajectories of trait dynamics (e.g. Lion, 2017). Because adaptive dynamics is built on the premise that trait evolution is the result of invasion by rare mutants (and not of short-term shifts in genotype frequencies), the eco-evolutionary feedback loop (Figure 1a) is played out at different time-scales: There is a temporal sequence in K E Y W O R D S adaptive dynamics, alternative reproductive phenotype, cannibalism, coexistence, environmental threshold model, population dynamics which ART expression evolves, followed by equilibration of population size and structure, followed by ART evolution, and so on, until an ESS of ART expression is reached. We modelled the population dynamics using an existing two-sex, two-ART demographic model, which comprises the demographic rates of adult females, fighters, sneakers and their offspring, zygotes ( Figure 1b). The proportion of male zygotes that develops into fighters is denoted by β ( Figure 1b); thus, 0 < β < 1 represents male dimorphism. The parameter β can be interpreted as being determined by the threshold of ART expression in a conditional strategy analogous to the ET model (Smallegange & Johansson, 2014), which can thus evolve in response to selection (Hazel et al., 1990(Hazel et al., , 2004. Specifically, we aimed to understand (a) how the probability u i that an individual of a targeted life stage i is killed when encountered by a fighter (which determines the proportion Q of the targeted life stage that is killed by fighters: Figure 1b) affects the evolutionarily stable strategy (ESS) of β, (b) how the probabilities that individuals of two different, targeted life stages are killed when encountered by fighters interactively influences this ESS, and (c) how intraspecific killing affects the effect of survival of the male morphs and competition parameters on this ESS. Intraspecific killing by fighters has mainly been linked to the killing of competitor males in conflicts over access to mates; the intraspecific killing of juveniles and females can also occur for other reasons, such as the early elimination of competition, or as an additional food source (Fox, 1975). Regardless of why fighters kill conspecifics, it results in increased mortality of the life stage targeted by the killing fighters, and it is probable that the demographic consequences of this intraspecific killing can affect the evolution of ARTs through the reduced survival of the targeted life stage (ecology-to-evolution pathway) and that the evolutionary trajectory of ART expression can in turn result in changes in population structure (the evolutionto-ecology pathway) ( Figure 1a). We also varied the competition parameters probability of engaging in competition, sneaker advantage, fighter costs and reward, as these have been shown to influence the evolution and maintenance of ARTs (Smallegange & Johansson, 2014). By including these parameters in our analyses, we could assess whether interactive effects between intraspecific killing and any of these parameters on the evolution of male morph coexistence occur. In addition, we explored how variation in survival affects the evolution and maintenance of ARTs in response to variation in intraspecific killing, in order to assess whether the effect of intraspecific killing depends on background mortality rates. For our study, we extended the two-sex, two-ART population model of Smallegange and Johansson (2014) by including intraspecific killing by fighters, and parameterized the model for bulb mites (Rhizoglyphus robini). Sneaker male bulb mites are known as scramblers, and, in our model analyses, we therefore refer to males as either being a fighter or a scrambler.
F I G U R E 1 (a) The eco-evolutionary feedback loop of Alternative reproductive tactics (ART) expression where shifts in population size and structure affect the selection gradient of ART expression (ecology-to-evolution pathway), and where evolutionary shifts in ART expression affect population size and structure (evolution-to-ecology pathway). (b) Stylized life cycle of a male dimorphic species in which fighters can kill conspecifics. Females (x), adult fighters (f) and adult scramblers (s) produce zygotes (z) at rates F i (i = x, f, s). Zygotes develop into females, fighters or scramblers at rates that are set by the sex ratio (ρ) and male morph ratio (β). Terms P i and G i (i = z, x, f, s) denote survival and growth rates, respectively. In this life cycle, fighters can kill individuals of each life stage, which reduces the survival rate Pi  (i = z, x, f, s). Zygotes are undifferentiated with respect to sex or male morph. The proportion of zygotes that develops into males equals ρ, and the proportion of males that develops into fighters equals β. A biological mechanism underlying the fraction β can be found in the ET model (Hazel et al., 1990(Hazel et al., , 2004. The development into an ART is determined The rate (individuals per time step) of surviving and staying in a stage equals P i (i = z, x, f, s), which, for adults, is calculated as P j = jj , with jj as the survival rate of morph j (i = x, f, s), and for zygotes as the following: The fertility functions of females and males of morph j (j = f, s) are given by F x = (B x,s + B x,f )∕(2n x ) and F j = 1∕(2n x )B x,j , respectively, which both depend on the number of births, B x,j, resulting from matings between females and males of morph j. The number of births from matings between females and morph j is defined as B x,j = ke x,j p x,j , and depends on (a) the clutch size (k), which decreases with female density, k = k 0 ∕(1 + n x ), (b) the encounter rate between a male of morph j and a female e x,j = e 0 n x n j ∕(n x + n s + n f ), with e 0 as the number of individuals of any life stage encountered by the focal individual per time step, and (c) the probability that an encounter results in a successful mating (p x,j ). The probability that an encounter results in a successful mating depends upon the strength of male-male competition (c m ), and the probability that a male of morph j gains access to a female when competing with a male of morph j is given by the following hawk-dove gamelike payoff matrix: In this matrix, V is the probability of accessing a female without costs, C are the fighter costs in terms of the probability of gaining access to a female, and ε is the probability of sneaking successfully.
We assume ε < V. This results in the probability that an encounter leads to a successful mating: The resulting population projection matrix includes all of the demographic rates of the four stage classes: The two-sex, two-ART population model is defined as where n t is the population vector at time t and A t is the projection matrix A at time t that is determined by the population vector at time t.

| Adjusting the baseline structure: including intraspecific killing
In the baseline structure of the model, encounter rates are calculated using only the densities of the adult stages (females, fighters and scramblers). Because we use this model to also investigate the effect of intraspecific killing on zygotes (see below), the encounter rate with zygotes has to be included. To achieve this, we calculated the encounter rate between individuals of stage j and individuals of stage i as: The focus of this study was to investigate the effect of intraspecific killing on the evolution of β, which is the fraction of male zygotes that mature into fighters. Let u i be the probability that an individual of a targeted life stage i (i = z, x, f, s) is killed when encountered by a fighter. We refer to the probability u i as the killing success, which ranges from zero to one. When u i equals zero, (1) Intraspecific killing is incorporated into the model by multiplying the survival rate without intraspecific killing (P i ) (default survival rate) with the proportion of individuals of the targeted morph that is not killed by fighter males (1 − Q i ). This means that we assume that individuals that mature, that is move from the zygote to an adult stage, are not vulnerable to intraspecific killing. This is probable, because during maturation, R. robini enters an immobile, quiescent stage during which individuals hide, and we assume that it is unlikely that these individuals are attacked by fighters. Table 1 gives an overview of all of the parameters and parameter values used in the model.

| Population model
Based on all of the demographic rates and intraspecific killings, we derive the following population projection matrix N: The population model is defined as n t+1 = N t n t , where n t is the population vector at time t and N t is the projection matrix at time t (Equation 8) that is determined by the population vector at time t. When a parameter was varied in one of the analyses, the range in which the parameter was varied is also given. All values, except for those under "Intraspecific killing," are taken from Smallegange and Johansson (2014).

| Evolutionary dynamics
The adaptive dynamics approach (Geritz et al., 1998)  . However, in this study we took a more system-specific approach (without losing sight of the more general questions) and incorporated differences in survival rates between the morphs and sexes (Equation 8). For this reason, we were unable to solve the model analytically, and instead, we ran simulations. To this end, we first used the population model n t+1 = N t n t with the resident trait value β to create a time series of 1,000 time steps t, in order to arrive at the equilibrium densities at t = 1,000. Whether equilibrium was reached was verified by checking whether the dominant eigenvalue was equal to one.
We then replaced the resident trait value of β with the mutant value β' and calculated the dominant eigenvalue of the resulting matrix N' to arrive at the invasion fitness. Therefore, it is assumed that the invading mutant experiences the population structure of the resident population.
The direction of evolutionary change is determined by the selection gradient, which is defined as the slope of the invasion fitness with regard to the variant trait at β' = β. When the selection gradient is positive (negative), a mutant with a slightly higher (lower) trait value can invade the population and replace the resident. In adaptive dynamics, a candidate evolutionary endpoint can be found when the selection gradient equals zero. When this point is also resistant to invasion and is an attractor for gradual evolution, it can be considered an ESS (McGill & Brown, 2007). We assessed these ESS criteria using pairwise invasibility plots (PIPs). The invasion fitness is always equal to zero when β' = β and represents the primary isocline in the PIP. Situations in which the invasion fitness is zero while β' ≠ β represents the secondary isocline. The intersects between the primary isocline and the secondary isocline represent candidate evolutionary endpoints (β*). We assessed whether each candidate evolutionary endpoint was convergence and evolutionarily stable using three visual criteria (Geritz et al., 1998): (a) if the vertical line through β* is entirely in a negative invasion fitness area, this indicates that β* cannot be invaded, and is therefore evolutionarily stable; if this is not the case, β* is an evolutionary repeller, and the value of β will evolve away from this value; (b) any resident strategy can be invaded by a mutant closer to the singular strategy, β*, if the area to the left of the secondary isocline and immediately above the primary isocline and the area to the right of the secondary isocline and below the primary isocline are positive invasion fitness areas; in this case, β* is convergence-stable; (c) if the horizontal line through β* is entirely within a positive invasion fitness area, this indicates that β* can always spread through the population when initially rare.
When all of these criteria are met, β' can be considered an ESS (β ESS ).
In some cases, the singular strategy is an evolutionary repeller (see also Results); in all other cases, the selection gradient vanishes and singular strategies are convergence-stable. However, once adopted by the resident population, these singular strategies are evolutionarily neutral, indicating that mutant strategies have the same fitness as the resident population. The implications of fitness equality are common in classical game matrices, including the hawk-dove game.
In the adaptive dynamics framework, such evolutionary neutrality represents a limit case between an ESS and an evolutionary branching point (Geritz et al., 1998). It has been shown that these points can turn into an invasion-resistant strategy through only slight adjustments of the model structure (Dieckmann & Metz, 2006), and we will refer to them henceforth as ESSs.
Using this procedure, we explored how β  Table 1. size (Smallegange, 2011). Male final instars above a size threshold are more likely to mature into fighters; below the size threshold, they are more likely to mature into scramblers (Smallegange, 2011). All parameter values are taken from Smallegange and Johansson (2014), except values for u i (i = z, x, f, s) that we set to zero as default, but in our model, analyses varied each value of u i between zero and unity (Table 1). In R. robini, the sex of individuals is determined by sex chromosomes (XX females, X0 males) (Oliver, 1977), thus, we assume no sex ratio adjustment by females (unlike, e.g., in Alonzo & Sinervo, 2001 and set ρ = 0.5. There are deviations in adults from a 1:1 sex ratio due to differences in life-history trajectories between females and males (Smallegange, 2011), and due to adult fighters killing males (Smallegange & Deere, 2014). These differences and their effects on β ESS , however, are taken into account in our model as we take a life cycle approach.

| RE SULTS
Firstly, we explored how the success of intraspecific killing by fighter males influences the ESS of β, β ESS . If the probability that a fighter kills a zygote (u z ) or another fighter (u f ) increases, the value of β ESS decreases but is always higher than zero (Figure 2), indicating that fighters and scramblers always coexist if fighters kill zygotes or other fighters. The probability that a fighter kills a female (u x ) has no influence on β ESS , and hence no effect on male morph coexistence ( Figure 2). If the probability that a fighter kills a scrambler (u s ) is very low (u s < 0.005), evolutionary bistability occurs: β ESS1 is be- and ε (sneaker benefit)) on β ESS . Varying u z and each of the competition parameters interactively affects the value of β ESS (first column of Figure 6; interactive effects are inferred from non-parallel isoclines in each panel), but the region of parameter space across which male morph coexistence occurs is only slightly affected, and only at low values of u z and ε (Figure 6 m). Varying u x and each of the competition parameters shows that u x has little or no effect on β ESS, and neither is the region of parameter space across which male morph coexistence occurs affected (second column of Figure 6).
Varying u f and each of the competition parameters interactively affects the value of β ESS , except in the case of the costs of fighting for fighter C (third column of Figure 6; interactive effects are inferred from non-parallel isoclines). Only when varying u f and ε, the sneaker advantage, is the region of parameter space across which  Figure 6).

| D ISCUSS I ON
One of the key questions in evolutionary biology is: What maintains genetic and phenotypic variation? Prime examples of such diversity are male dimorphic species (Oliveira et al., 2008). Recent empirical results on the evolution and maintenance of ARTs emphasize how population feedback within an eco-evolutionary feedback loop can influence the evolution of ARTs; that is, allowing for population feedback in experimental evolution studies on ART expression in the bulb mite shows that the observed direction of evolution of male morph expression can be different, and even diametrically opposed to the direction predicted by evolutionary theory (Smallegange & Deere, 2014). Negative population feedback is also the main driver of the expression of the genetically determined, alternative reproductive strategies of female side-blotched lizards (Uta stansburiana) (Sinervo et al., 2000). Orange-throated lizard females have few and large offspring, whereas yellow-throated females have many and small offspring. In years of low population density, yellow-throated females are favoured, but the many offspring cause an increase in population growth that overshoots carrying capacity. These new conditions favour orange-throated females, and their offspring cause a decrease in population growth that undershoots carrying capacity. Again, only by taking the population feedback into account, can we understand the evolution and maintenance of these genetically determined, alternative morphs (Sinervo et al., 2000).
Here, we studied how the eco-evolutionary interaction between evolution of ART expression and shifts in equilibrium population size and structure affects theoretical predictions of the ESS for male morph expression in the bulb mite, in relation to the extent to which fighters killed conspecifics of different life stages. We found that the killing of females had no effect on the evolution of male morph coexistence. This is probably because, although the killing of females reduces the number of reproducing females, each reproducing female produces larger clutches, as clutch size is only dependent on the number of females. In population-dynamic equilibrium, larger clutch sizes therefore compensate for the smaller number of reproducing females. In reality, there will be a limit to the number of offspring that a female can produce in a single clutch, and when this is taken into account, we expect that the killing of females will affect fighter expression. Particularly, if there would be no (genetic) constraints on sex ratio expression (not in case of R. robini which has chromosomal sex determination), in which case ART coexistence depends on a complex interplay between ART frequency, female choice and sex allocation (Alonzo & Sinervo, 2001 (Smallegange, 2011). According to Claessen et al. (2004), the combination of these processes can lead to the stabilization of population dynamics when mortality weakens competition. We indeed found stable population dynamics when fighters killed zygotes or fighters. In addition, the combination of these processes can lead to evolutionary bistability when killing has an indirect positive effect on the cannibal, such as reduced competition. In our model, fighters directly reduce mating competition by killing scramblers, and we indeed found bistability when fighters killed scramblers. We did not account for any nutritional benefits that could be obtained from an intraspecific killing event. Therefore, to explore this scenario further, we adjusted the model such that intraspecific killing is now a cannibalistic event, and the energy gained from such an event is directly allocated to reproduction (Supporting information Appendix S1). This analysis shows that when fighters cannibalize other fighters, increased reproduction through cannibalism can compensate for the aforementioned negative fitness effects of the killing of other fighter males (Supporting information Figure   S1). It is important to note that cannibalism only has a positive fitness effect on fighter expression if the gain in reproductive effort is greater than the reduction in population growth rate through the killing of fighter males. As, to our knowledge, it is unknown how much energy mites actually gain from feeding on conspecifics, we cannot compare this prediction against empirical data. In addition, it has been observed that bulb mites of other life stages feed on in- To what extent this would affect the evolution of fighter expression, and hence male morph coexistence, remains to be explored.
Recently, a third male morph termed the mega-scrambler has been identified in R. robini (Stewart, van den Beuken, Rhebergen, Deere, & Smallegange, 2018). Systems comprising three ARTs are not uncommon (Rowland & Emlen, 2009;Shuster & Sassaman, 1997;Sinervo & Lively, 1996), which raises the question of how the model could be extended to address three-ART systems. Friedman  In conclusion, this study shows that intraspecific killing can influence the evolution and maintenance of male ARTs. Interestingly, we showed that when fighters kill other fighters, this reduces their fitness and thus increases the relative fitness benefits of scrambler males. Perhaps in this way, selection can still favour scrambler expression, even if scramblers would suffer increased mortality, like Smallegange and Deere (2014) observed in their experimental evolution study, where selection against scramblers resulted in increased (and not decreased) scrambler expression. This "hidden" scrambler benefit could mean that ARTs are more likely to evolve in species with a similar life history. However, if fighters gain reproductive benefits from the killing of conspecifics through cannibalism, the detrimental effects on fighter fitness can be nullified.
Alternatively, we have some evidence to suggest that maternal selection, that is the mother's phenotype affecting morph expression of her male offspring, could operate in bulb mites (Smallegange, 2011). In general, under certain conditions, such maternal effects on selection can cause a population to evolve maladaptively away from a fitness peak (Kirkpatrick & Lande, 1989), which could also explain the experimental evolution finding that selection against scramblers increased scrambler frequency (Smallegange & Deere, 2014). When investigating the evolution and maintenance of ARTs, one should therefore not only consider the effects of male-male competition and life-history differences between phenotypes (cf. Smallegange & Johansson, 2014) or maternal selection, but also how the effects of population-dynamical processes such as intraspecific killing and cannibalism feedback to influence the evolution of ARTs. The influence of such eco-evolutionary feedbacks on trait evolution is not necessarily limited to polymorphisms like ARTs, but likely also plays a role in the evolution of other traits directly linked to demography, such as behavioural syndromes.
The eco-evolutionary feedback loop that we studied is one where ecology and evolution are played out at different time-scales.
Eco-evolutionary interactions can also emerge when ecology and evolution occur on similar time-scales (Sinervo et al., 2000), complicating the eco-evolutionary process. But is the short-term or the long-term eco-evolutionary selective process most decisive in the evolution of a trait? Short-term, non-adaptive plasticity can, for example, potentiate evolution by increasing the strength of directional selection (Chevin, Lande, & Mace, 2010), whereas shortterm adaptive plasticity can constrain long-term evolution (Price, Qvarnström, & Irwin, 2003). Current empirical steps towards unravelling drivers of the eco-evolutionary process (e.g. Sinervo et al., 2000), in concert with theoretical advances to model shortterm eco-evolutionary processes (Coulson et al., 2017) that can be incorporated within the adaptive dynamics framework, pave the way for understanding how selection on different time-scales affects trait evolution. This is particularly pertinent now with ongoing changes in the environment and climate, as eco-evolutionary dynamics are most common during transient periods, for example, following an environmental change (Hiltunen & Becks, 2014).