Host phenology regulates parasite–host demographic cycles and eco‐evolutionary feedbacks

Abstract Parasite–host interactions can drive periodic population dynamics when parasites overexploit host populations. The timing of host seasonal activity, or host phenology, determines the frequency and demographic impact of parasite–host interactions, which may govern whether parasites sufficiently overexploit hosts to drive population cycles. We describe a mathematical model of a monocyclic, obligate‐killer parasite system with seasonal host activity to investigate the consequences of host phenology on host–parasite dynamics. The results suggest that parasites can reach the densities necessary to destabilize host dynamics and drive cycling as they adapt, but only in some phenological scenarios such as environments with short seasons and synchronous host emergence. Furthermore, only parasite lineages that are sufficiently adapted to phenological scenarios with short seasons and synchronous host emergence can achieve the densities necessary to overexploit hosts and produce population cycles. Host‐parasite cycles also generate an eco‐evolutionary feedback that slows parasite adaptation to the phenological environment as rare advantageous phenotypes can be driven extinct due to a population bottleneck depending on when they are introduced in the cycle. The results demonstrate that seasonal environments can drive population cycling in a restricted set of phenological patterns and provide further evidence that the rate of adaptive evolution depends on underlying ecological dynamics.

predator-prey, herbivore-plant, and parasite-host systems have demonstrated the importance of seasonal activity patterns on interspecies interactions (Abbott & Dwyer, 2007;Greenman et al., 2004;Kamo & Sasaki, 2002;Taylor et al., 2013). Seasonal activity patterns determine the temporal abundance of a population, which modifies the strength of interspecies interactions (Barber et al., 2016;Bewick et al., 2016;Burkett-Cadena et al., 2011;Miller-Rushing et al., 2010;Paull & Johnson, 2014). Here, we demonstrate the consequence of seasonal activity patterns on parasite-host population dynamics and how seasonal patterns can result in parasite-host population cycles.
Additionally, we explore how parasite-host population cycles can alter the rate of parasite virulence evolution.
Population cycling generally starts with an overexploitation of resources followed by a population crash that allows resources to rebound (Myers & Cory, 2013). In the classic lemming demographic cycles, lemmings overconsume plant resources resulting in dramatic declines in lemming population sizes in subsequent years due to plant scarcity (Krebs, 2013). The plant populations are released from lemming herbivory and increase in abundance, providing sufficient resources for lemming population growth and a restart of the demographic cycle. Intrinsic, delayed density-dependent drivers such as these can account for the periodic or quasiperiodic oscillatory population dynamics observed in many ecologically coupled systems (Myers, 2018). Seasonal activity patterns, or phenology, influence the impact of interspecies interactions on demographic dynamics (van Asch & Visser, 2007;Yang & Rudolf, 2010). That is, seasonal activity patterns determine the number and type of interspecies interactions by altering the proportion of a population that is active throughout the year. For example, measles transmission is tightly linked to school terms such that transmission peaks when children are in school and crashes during vacation periods (Fine & Clarkson, 1982;Finkenstädt & Grenfell, 2000). Similarly, variation in demographic dynamics can impact species evolution, for example, resource-driven changes in host abundance are predicted to impact parasite virulence evolution (Hite & Cressler, 2018). Prior theoretical research demonstrated that the total number of parasite infections, which determines the parasite population size, varied dramatically among different host phenological patterns (MacDonald et al., 2021). Furthermore, the virulence strategies that maximize parasite fitness also differed among phenological patterns due to the differences in the temporal distribution of new infections. However, this work restricted host demographic feedbacks such that the potential for population cycles subsequent effect on evolutionary dynamics could not be investigated.
Changes in the population sizes of interacting species that result from ecological interactions can also influence the rate or direction of evolutionary change (Govaert et al., 2019). These eco-evolutionary Here, we explore eco-evolutionary feedbacks driven by parasite infection in a seasonal environment. We extend a previously published modeling framework (MacDonald et al., 2021)

| Within-season dynamics
The model describes the transmission dynamics of a freeliving, obligate-killer parasite that infects a seasonally available host (Figure 1). The size of the emerging host cohort in season n, ŝ(n), is determined by the number of hosts that reproduced in season n − 1. ŝ(n) enters the system at the beginning of the season over a period given by the function g(t,t l ). Hosts have non-overlapping generations and are alive for one season. The parasite (v) infects hosts and must kill the host to release new infectious progeny. We assume that the parasite is monocyclic and completes one generation per season. The monocyclic constraint is enforced by assuming that only the first generation of parasites in a season, v 1 , has enough time to release the second generation of parasites, v 2 . This transmission scenario occurs in many natural parasites (e.g., univoltine insects parasitized by ichneumonid wasps (Campbell, 1975;Delucchi, 1982;Kenis & Hilszczanski, 2007)). Parasites may effectively complete only one round of infection per season if the second parasite generation does not have enough time in the season to release new parasites in short-lived hosts or if the susceptible host stage is present for such a short period of time each season that there are no susceptible host stages available when the first generation of parasites kills infected hosts.
We refer to the generation of parasites that infect the susceptible host stage, s, as v 1 and the parasite progeny released from infected hosts as v 2 . τ is the delay between infection by v 1 and host death when v 2 is released. We ignore the progression of s hosts to later life stages as it does not impact transmission dynamics. The initial conditions at the start of each season are is the number of parasites at the beginning of season n as determined by the number of parasite progeny produced in n − 1. The transmission dynamics in season n are given by the following system of delay differential equations: where μ is the host death rate, δ is the decay rate of parasites in the environment, α is the transmission rate, β is the number of parasites produced upon host death, and τ is the delay between host infection and host death (Table 1). We make the common assumption for freeliving parasites that the removal of parasites through transmission (α) is negligible (Anderson & May, 1981;Caraco & Wang, 2008;Dwyer, 1994), i.e., (1b) ignores the term − s(t)v 1 (t). As virulence is the lifetime reduction in host fitness due to infection, we assume that parasites with shorter times between infection and host death (short incubation periods) are more virulent. Thus, τ is equivalent to the inverse of virulence where low virulence parasites have long τ and high virulence parasites have short τ. All parameters with their respective values are described in Table 1.
The function g(t, t l ) is a probability density function that captures the per-capita host emergence rate by specifying the timing and length of host emergence. We use a uniform distribution (U( ⋅ )) for analytical tractability, but other distributions can be used.
where t l denotes the length of the host emergence period, and T denotes the season length. The season begins (t 0 = 0) with the emergence of the susceptible host cohort, ŝ(n). The host cohort emerges from 0 ≤ t ≤ t l . v 2 parasites remaining in the system at t = T give rise to the initial parasite population in the following season (v 2 (T) =v(n + 1) = v 1 (0)). Parasites that have not killed their host by the end of the season do not release progeny. Background mortality arises from predation or some other natural cause. We assume that infected hosts that die from background mortality do not release parasites because the parasites are either consumed or the latency period corresponds to the time necessary to develop viable progeny (Wang, 2006;White, 2011). We solve Equations 1a-c analytically, Appendix S1.

| Between-season dynamics
We investigate the impact of the feedback between host demography and parasite fitness on parasite evolution by allowing the size of the emerging host cohort be a function of the number of uninfected hosts remaining at the end of the prior season using a difference , F I G U R E 1 Diagrammatic representation of the infectious cycle within each season. The host population (ŝ(n)) at the start of season n is the offspring of uninfected hosts that survived and reproduced at the end of the prior season. The parasite population at the start of season n(v 1 (0) =v(n)) is derived from infected hosts killed by the parasite prior to the end of season n − 1 and survived in the environment until the end of the season (v 2 (T)). All parasites emerge at the beginning of the season (t = 0) while all hosts emerge at a constant rate between 0 ≤ t ≤ t l . The rate of new infections is density-dependent resulting in the majority of infections occurring near the beginning of the season when susceptible host and free parasite densities are high. Parasite-induced host death at time postinfection releases parasite progeny (v 2 ) into the environment where they decay in the environment at rate . The monocyclic parasite progeny (v 2 ) do not infect uninfected hosts within the same season. Parasite progeny that survives in the environment to the end of the season comprises the parasite population that emerges in the following season (v 2 (T) =v(n + 1)) In Appendix S1, we find analytical solutions for both ŝ(n + 1) and v(n + 1). However, we primarily explore the between-season dynamical behavior of the model numerically as analytical solutions cannot be used in parameter ranges that lead to population cycles. We discuss the stability analysis in more detail in Appendix S1.

| Parasite evolution
Evolutionary invasion analysis (Geritz et al., 1998;Metz et al., 1992) was used to study parasite adaptation to different seasonal host activity patterns. We first extend system (1) to follow the invasion dynamics a rare mutant parasite: where m subscripts refer to the invading mutant parasite and its corresponding traits. The initial conditions at the beginning of each season where v * and ŝ * are end of season equilibrium densities for parasite and host, respectively. See Appendix S2 for details of the time-dependent solutions for Equations (2a-2e).
The invasion fitness of a rare mutant parasite depends on the density of v 2m produced by the end of the season (v 2m (T)) in an environment with a resident parasite at equilibrium density v * . When system dynamics are equilibrial, the mutant parasite invades in a given host phenological scenario if the density of v 2m produced by time T is greater than or equal to the initial v 1m (0) = 1 introduced at the start of the season (v 2m (T) ≥ 1). When < T − t l , mutant invasion fitness can be found using When > T − t l , mutant invasion fitness can be found using To study the evolution of virulence traits in equilibrial environments, we assume that resident and mutant strains are identical at all other traits (e.g., = m ). Note that because there is no trade-off between β and τ, the parasite growth rate in the host is  (3a) and (3b) incorporates the effect of the resident on the population state (number of susceptibles over one season). This means that v 2m (T) is not a measure of R 0 , which by definition assumes a nondisease environment. Thus, we can use v 2m (T) as defined in (3a) and (3b) as a maximand in evolutionary dynamics (Lion & Metz, 2018).
In the present study, cycling can occur when host carryover is included in the model for some parameter ranges. When parasite-host dynamics are cycling (3) no longer reliably predicts the outcome of parasite evolution as periods of low host density can drive adaptive mutants to densities less than 1. From a purely mathematical standpoint, the criterion v 2m (T) ≥ 1 correctly predicts which mutants can invade in cycling populations. However, the invasion criterion does not account for the possibility that a mutant parasite that invades in its first season can drop below 1 in a later season. We thus conduct simulation analysis to verify that the evolutionary stable level of virulence is qualitatively the same as previous results.
The simulation analysis was done by first numerically simulating system (1) with a monomorphic parasite population. A single mutant parasite is introduced at the beginning of the 100th season when the system dynamics have settled on their attractor. The mutant's virulence strategy is drawn from a normal distribution whose mean is the value of τ from the resident strain ( m = r +  (0, 0.1)). System (2) is then numerically simulated with the resident and mutant parasite.
New mutants arise randomly after 1000 seasons have passed since the last mutant was introduced, at which point system (2)

| RE SULTS
Parasites with high fitness in some seasonal environments can drive dynamic parasite-host cycles resembling classical consumer-resource cycles ( Figure 2). In the present model, parasites that can achieve sufficiently high densities infect and sterilize a substantial levels below those that can destabilize host-parasite dynamics and cause demographic cycles.
The parasite densities that can be attained in each host phenological scenario determines whether the system reaches stable equilibrial inter-annual dynamics or quasiperiodic parasite-host population cycles (Figure 3). In the majority of scenarios in which cycling occurs, the discrete dynamics form a closed invariant curve in the phase plane in which the phase is incommensurate, and thus, the asymptotic trajectory fills the invariant curve by never repeating itself ( Figure 4a). That is, the population sizes of both the host and parasite do not repeat across seasons, resulting in quasiperiodic cycles that are likely generated by a Neimark-Sacker bifurcation (Strogatz, 2018; see Appendix S1).
Additional environmental factors that promote high parasite density, such as low environmental decay rates, can increase the parameter region where cycling occurs (  (a), to destabilize demographic dynamics resulting in a bifurcation that drives quasiperiodic parasite-host dynamics. The bifurcation diagram shows end of season parasite densities for parasites with different virulence phenotypes (τ) during seasons 800-900 in a system where the host season is short (T = 4) and hosts emerge synchronously (t l = 1). The most fit parasites (2.75 < < 3.26) achieve densities that can disrupt dynamics and cause cycling. Parasites with virulence phenotypes that are too high ( < 2.75) or too low ( > 3.26) do not cause parasite-host cycles in this host phenological environment. (b-c) The population dynamics of hosts (b) and parasites (c) in a system experiencing quasiperiodic population cycles ( = 2.8, T = 4, t l = 1, other parameters found in Table 1) after reaching the quasiperiodic attractor. High parasite densities (ex. season 3-4) infect and sterilize a large proportion of the host population resulting in a dramatic host population decline (ex. seasons 4-5). The limited number of susceptible hosts causes a subsequent decline in parasite populations (ex. seasons 5-7). Host density rebounds once relieved from infection pressure (ex. seasons 6-8) allowing the parasites to exploit the host population again, driving a continuation of quasiperiodic cycling. In both panels: T = 4, t l = 1, all other parameters found in Table 1 advantageous mutant parasite strain. This eco-evolutionary feedback results in the extinction of many advantageous mutants and a reduced rate of evolution toward the virulence strategy that optimizes parasite fitness. That is, many mutant parasites with adaptive phenotypes that arise in a cycling system will not increase in frequency and ultimately be lost from the population.

| DISCUSS ION
The observed parasite-host population cycles emerge from a delayed density-dependent mechanism characteristic of consumerresource feedbacks (Turchin, 2013). In this system, the discrete host activity period introduces a delayed carryover effect in which the number of infected hosts in one season governs the host population size in the next. Although consumer-resource interactions can drive cycling in continuous time models, cycles are less likely to occur without an externally imposed delay (Keeling & Rohani, 2011 The stable parasite-host dynamics observed in some host phenological patterns differs from seminal theoretical studies demonstrating chaotic dynamics at all population growth rates of lethal parasites (May, 1985). Our results suggest that host phenology can stabilize host-parasite dynamics and provide one potential explanation for why chaotic dynamics are often not observed in natural obligate-killer parasite systems. Other model parameters such as natural host mortality rate and parasite decay rate also modulate parasite population sizes and thus also alter which phenological scenarios can lead to periodic population cycles. Furthermore, several factors that have been shown to impact the probability of dynamic population cycles not explored in this model could also modulate the phenological scenarios in which cycling could be expected (Hilker et al., 2020;Koella & Doebeli, 1999). For example, higher infectedhost fecundity would likely stabilize the dynamics for a greater range of phenological patterns.
Population cycles resulting from a consumer-resource ecological feedback precipitates an eco-evolutionary feedback that affects the rate of adaptive evolution. In this model, parasites with advantageous mutations always invade non-cycling systems. That is, advantageous mutants displace residents both in systems where the parasite is not sufficiently adapted in a host phenological environment that could support high parasite densities as well as in systems where the host phenological pattern cannot support densities sufficient to cause population cycling even for optimally adapted parasites. By contrast, only a fraction of parasites with adaptive mutations introduced into cycling systems can invade, effectively reducing the rate of adaptive evolution. These results suggests that adaptation in cycling seasonal disease systems is likely to proceed more slowly. Parasites with adaptive mutations that do not invade fail to increase sufficiently to prevent their extinction due to the F I G U R E 3 Parasite-host cycles occur in some, but not all, host phenological patterns. Boundary plot shows host phenological patterns where dynamics are stable ("endemic equilibrium") or cycling for parasites possess the optimal virulence trait for their phenological environment. Parasites are more likely to achieve the densities necessary to drive cycles when host emergence periods are short (small values of T) and host emergence is synchronous (small values of t l ). All other parameters found in Table 1 3 cling predator-prey model (Yamamichi et al., 2014). An assumption of the current model is that mutants are introduced at the beginning of random seasons, regardless of parasite population size. However, mutants are less likely to arise when parasite population size is small (Crow & Kimura, 1970) suggesting that the true impact of population cycling on the evolutionary rate is likely greater than estimated here. That is, the proportion of advantageous mutants lost in cycling populations in nature is likely greater than found here as mutants are more likely to arise at points in the cycle when parasite populations are large and on the precipice of crashing.
Our results extend previous theory on the interaction between cycling and the rate of adaptive evolution. Many previously published investigations focus on the impact of temporal fluctuations on long-term evolutionary outcomes (e.g., Donnelly et al., 2013;Ferriere & Gatto, 1995;Ferris & Best, 2018;Grunert et al., 2021;Metz et al., 1992). The results presented here suggest that temporal cycling can slow the rate of adaptive evolution by constraining when adaptive mutants can successfully invade, even if the long-term evolutionary outcome remains constant. That is, while prior work revealed the most advantageous long-term evolutionary strategies (Ferriere & Gatto, 1995;Metz et al., 1992), our approach identified the demographic conditions leading to the extinction of advantageous mutants. In addition, our results support previously published conclusions showing that evolutionary adaptions can not only drive demographic cycles (Ferriere & Gatto, 1993;Metz et al., 1995) but also extends this work through our result that cycling slows the rate of adaptive evolution.
These results suggest that spatial variation in host phenology could drive differences in demographic dynamics observed across geographic space. For example, parasite-host systems in more extreme latitudes and at higher altitudes are more likely to cycle than conspecifics in less extreme environments (Baltensweiler & Fischlin, 1988;Klemola et al., 2002;Schott et al., 2010). The activity periods in the more extreme environments tend to be shorter, and hosts may emerge more synchronously  in line with the predictions from the current model. These predictions could be tested empirically by studying the population dynamics of disease systems with forest F I G U R E 4 Mutant parasites with more adaptive virulence phenotypes often fail to invade when resident parasite and host dynamics are cycling. The phase plane (a) shows the discrete time limit cycle for host (ŝ) and resident parasite (v) densities (T = 4; t l = 1; = 2.8 in this example). The blue section denotes the phase of the parasite-host cycle when rare adaptive parasites can invade; the same mutant fails to invade when introduced at all other time points despite having the same selective advantage. The line (a) depicts the same iteration (six seasons) of the quasiperiodic dynamics of this system as illustrated in (b1) and (b2). An advantageous mutant fails to invade (red line, b3) if introduced in seasons when host density will decrease (red point, b1) and resident parasite density is moderate or high (red point, b2). The same advantageous mutant can invade (blue line, b3) and eventually replace the resident parasite if it is introduced when host density will increase (blue point, b1) and resident parasite density is low (blue point, b2). m = 2.81, all other parameters found in  (Best, 2018;Ferris et al., 2020). Future theoretical and empirical investigations into the impact of parasite-host cycles on the evolution of host resistance alleles, as seen in Gypsy moth populations (Elderd et al., 2008), could determine whether parasite-host co-evolution would stabilize population dynamics for a greater range of host phenological patterns. Similarly, the strength and possibly direction of selection on hosts will fluctuate as the system cycles, potentially favoring alternative host phenological patterns that in turn select for parasite traits with lower impacts on host fitness. Another interesting extension is the role genetic drift could play for parasite adaptation in stable versus cycling dynamics (Kennedy & Dwyer, 2018). The impact of drift on parasite evolution in cycling populations is highly complex and difficult to predict a priori. We will extend the current model to incorporate neutral evolution in future studies.
Relaxing some of the assumptions in this model is unlikely to qualitatively alter the major conclusions. For example, relaxing the monocylic parasite life cycle assumption will likely not change the result that cycles occur more readily in environments with short seasons and synchronous host emergence. Polycyclic parasites may even drive cycles for a larger range of phenological patterns as multiple infection cycles within a season can exacerbate decreases in host densities. Similarly, relaxing the obligate-killer assumption will likely decrease but not eliminate the range of phenological patterns that drive cycles by decreasing the impact on host fitness. Although the model as presented applies to only a narrow range of parasites in nature, many more parasite-host systems conform to models that include these extensions such as soil-borne plant pathogens, demicyclic rusts, postharvest diseases, and many diseases infecting univoltine insects (Crowell, 1934;Gaulin et al., 2007;Holuša & Lukášová, 2017;Zehr, 1982).
Environmental conditions such as phenology impact the frequency of interspecies interactions and thus the ecological importance of the interaction on population demography. Here, we show that short host seasons and synchronous host emergence allow parasites to reach densities sufficient to destabilize population dynamics and cause demographic cycling. The rate of adaptive parasite evolution in a cycling population is substantially slower than in an equilibrial population as beneficial mutations are more likely to go extinct when host population sizes are small or parasite population sizes are large. These results demonstrate that externally imposed environmental conditions such as host phenology can be important determinants of population cycling. It is important to consider ecological dynamics when predicting evolution by natural selection.

ACK N OWLED G EM ENTS
We would like to thank Erol Akçay for useful discussions. We would also like to thank Bryce Morsky and the members of the Mideo Laboratory for looking at early drafts of the manuscript and two anonymous reviewers for constructive comments that greatly improved the manuscript.
F I G U R E 5 Adaptive evolution proceeds more slowly during parasite-host demographic cycles than in stable equilibrial systems. Parasite evolution toward an intermediate optimal virulence strategy occurs more slowly when population demography is cycling (a) than in an equilibrial dynamic system (b). (a) Increases in parasite density as parasites evolve drive demographic cycling for 2.75 < < 3.26 . Population cycling delays adaptive evolution as rare mutants fail to invade the system when introduced at many phases of the dynamic cycle despite their selective advantage. By contrast, rare advantageous mutants always invade systems with stable dynamics (b). Plots show twelve independent simulations for each set of parameters-six runs starting at a virulence level lower than optimum and six runs starting at a virulence level higher than the optimum-where an adaptive mutant is introduced into the population no more than once every 1000 seasons. Evolutionary time represents the cumulative number of adaptive mutants sequentially introduced into each population. The average evolutionary time needed to reach the optimal virulence strategy is higher in the cycling system ((a) 21 mutants, range: 6-42 mutants) than in the stable system (b. 14 mutants, range: 6-27 mutants). Population cycling could not occur in (b) as the host cohort size remained constant across seasons (ŝ = 10 8 ); host cohort size in (a) (ŝ(n)) was determined by the number of hosts that reproduced in season n − 1. T = 4, t l = 1, all other parameters found in Table 1 (a) (b)