Multiple host use and the dynamics of host switching in host–parasite systems

The link between multi‐host use and host switching in host–parasite interactions is a continuing area of debate. Lycaenid butterflies in the genus Maculinea, for example, exploit societies of different Myrmica ant species across their ranges, but there is only rare evidence that they simultaneously utilise multiple hosts at a local site, even where alternative hosts are present. We present a simple population‐genetic model accounting for the proportion of two alternative hosts and the fitness of parasite genotypes on each host. In agreement with standard models, we conclude that simultaneous host use is possible whenever fitness of heterozygotes on alternative hosts is not too low. We specifically focus on host‐shifting dynamics when the frequency of hosts changes. We find that (i) host shifting may proceed so rapidly that multiple host use is unlikely to be observed, (ii) back and forth transition in host use can exhibit a hysteresis loop, (iii) the parasites' host use may not be proportional to local host frequencies and be restricted to the rarer host under some conditions, and (iv) that a substantial decline in parasite abundance may typically precede a shift in host use. We conclude that focusing not just on possible equilibrium conditions but also considering the dynamics of host shifting in non‐equilibrium situations may provide added insights into host–parasite systems.


Introduction
Host-parasite interactions are ubiquitous across ecosystems and are widely recognised as important drivers of dynamic evolutionary change, a perspective that inspired the Red Queen hypothesis (van Valen, 1973), the idea that parasitism promotes recombination and sex (e.g. Hamilton, 1980;Bell & Smith, 1987;Salathé et al., 2008) and the concept of the geographic mosaic of coevolution (Thompson, 2009). The range of strategies and counterstrategies that have evolved in antagonists is impressive, and, in many cases, parasites can overcome the defence of a potential host only by being highly specialised. Yet in other cases, e.g. leaf mining or galling insects, communities are dominated by generalist species utilising a range of host species (Askew & Shaw, 1974;Bailey et al., 2009) with frequent shifting from one host to another (Hardy & Cook, 2010). Also, where invasive host species reach enemy free space, they recruit natural enemies from the local parasite species pool, either expanding the latter's host range or stimulating switching to the new resource all together (Cornell & Hawkins, 1993;Schönrogge et al., 2012).
The genus of lycaenid Maculinea (syn. Phengaris; Fric et al., 2007) butterflies provides a system where understanding the constraints for host switching would contribute greatly to our understanding of its current patterns of host use. The five European species have been studied extensively as social parasites infiltrating colonies of Myrmica ants. The impact of Maculinea larvae on their host colonies can indeed be dramatic: depending on whether they follow a cuckoo or predatory strategy (Thomas et al., 2005b); Maculinea infections can result in the complete demise of an infected ant colony. Adoption of young caterpillars into ant colonies (the butterfly larvae begin their life cycles by feeding on specific host plants; Thomas et al., 2005a) involves a number of complex and intricate interactions with their (potential) hosts including chemical deceit, acoustic signalling, and behavioural mimicry (Thomas et al., 2005b;Barbero et al., 2009;Fürst et al., 2012;Casacci et al., 2013).
The complexity of these adaptive traits seemed consistent with high levels of host specificity, and early field sampling suggested that each Maculinea species survived only with a single Myrmica host species despite larvae being carried to the nest by up to five different Myrmica species (Thomas et al., 1989). Notwithstanding artefactual errors in the literature (Thomas et al., 2005a), this view had to be revised. It is now clear that most Maculinea species successfully exploit colonies of different Myrmica species in different parts of their range (Thomas et al., 2005a,b,;Nash et al., 2008;Tartally et al., 2008). The predatory species Ma. arion and Ma. teleius show generally lower intra-population host specificity than the cuckoo species (Thomas et al., 2005a) but nonetheless cannot establish viable populations on non-primary hosts 2009). This even holds for Ma. teleius (primary host is My. scabrinodis), where rather high survival rates on a secondary host (My. rubra) have been reported (Thomas et al., 2005b). Yet, a few populations are apparently supported by My. rubra in otherwise My. scabrinodis-exploiting regions in central Europe (Witek et al., 2008;. The cuckoo species, Ma. alcon and Ma. rebeli, are typically highly specific to a single host antlocally and uniformly across wide regions, with just 0%-2% survival in the nests of non-host Myrmica (e.g. Thomas et al., 2005b;Witek et al., 2008). Host switches only occur at greater distances across Europe, involving mimicry of ants with very different chemical recognition profiles (cuticular hydrocarbon profiles, CHC) (Elmes et al., 2002), i.e. My. scabrinodis and My. rubra or My. ruginodis in the case of Ma. alcon and My. sabuleti and My. schencki in the case of Ma. rebeli Thomas et al., 2013). Remarkably, multiple host use has been reported (My. rubra and ruginodis) for Ma. alcon for some but not all the study sites in Denmark where both hosts were present. My. scabrinodis, an important host in other regions was not utilised, however, despite being quite abundant on some sites (Als et al., 2002). In some parts of Europe, notably in Austria (Schlick-Steiner et al., 2004), the Italian Apennines  and the Carpathian Basin (Tartally et al., 2008) near the range boundaries of different ant-exploiting types, Ma. rebeli is recorded as exploiting more than one host species on the same site. It is unclear whether these unusual sites support hybrids, co-occurring cryptic (sub-) species, polymorphisms, true generalists, or whether the Myrmica populations are unusually receptive to intruders in these regions (Thomas et al., 2005a). A recent and extensive review demonstrated that, wherever multiple hosts occurs on a site, a non-proportional use of host ant species is observed in the majority of cases (Tartally et al., 2019). And where multiple host use did occur, the CHC profiles of alternative host populations seem to be very similar . In the case of the Austrian sites, CHC profiles of pre-adoption larvae seem to match those of two different host species, initially suggesting a generalist strategy of Ma. rebeli in that region (Schlick-Steiner et al., 2004). Moreover, other explanations are not precluded since groups of five Ma. rebeli larvae at a time were used to extract semio-chemicals; furthermore, later studies of Ma. rebeli from other regions showed that the critical host-specific mimetic semio-chemicals are secreted by the parasites only after several days spent in host nests, while individuals adopted by other species of Myrmica suppress these secretions and rely (usually unsuccessfully) on camouflage from semio-chemicals acquired through contact with their hosts Thomas et al., 2013). In summary, with the few exceptions for Ma. alcon, Ma. teleius and Ma. rebeli mentioned, Maculinea populations are typically not successful in locally exploiting different hosts even if alternative putative host species co-occur, in some cases in greater abundance than the primary host (Thomas et al., 1989;2005a;Steiner et al., 2003;Nash et al., 2008;Tartally et al., 2008;Andersen et al., 2014). Further, some sites have been monitored over many years (up to 45 generations; e.g. Thomas et al., 2009;Ugelvig et al., 2011;Andersen et al., 2014), but to our knowledge, none of these studies provides direct evidence for a temporal host-shift, despite the inference from the study by Nash et al. (2008).
The failure to observe multiple local host use is puzzling as hostparasite interactions, as a matter of principle, introduce an element of frequency-dependent selection into the evolutionary arms race (Dybdahl & Lively, 1998;Thompson, 2009) and the potential for rapid turn over in community composition. Rare host species are more likely to evade specialised parasitism providing possible benefits in competition with other host species (Carius et al., 2001;Rolff & Siva-Jothy, 2003;Yoder & Nuismer, 2010). In multispecies communities, this should favour coexistence of multiple host-parasite pairs (Chaianunporn & Hovestadt, 2011) or introduce frequent turnover of antagonistic species (pairs) in local communities (Carvalho & Crisp, 1987;Savill & Hogeweg, 1999;Chaianunporn & Hovestadt, 2012;Rabajante et al., 2015), or promote coevolutionary alternation. We have to recognise, however, that shifting host use within a species of interbreeding parasites, a process called coevolutionary alternation (see Nuismer & Thompson, 2006), may not progress in the same way as, e.g. the invasion of a new parasite species into a community that utilises a different (and more abundant) host.
One reason why host shifting may in fact be difficult has long been recognised: heterozygote disadvantage (see Clarke & O'Donald, 1964;Ayala & Campbell, 1974). Whenever host utilisation requires a specificand at least in partgenetically based strategy or profile to match the defensive means of a specific host, genetic heterozygotes emerging from the mating of homozygotes that are each well adapted to alternative hosts may perform poorly on either host. Because the probability for homozygotes is proportional to p 2 , rare alleles (with proportion p in the population) will occur much more frequently in heterozygotes. Under such conditions, it may be extremely difficult for a new mutant or immigrant to establish in a parasite population despite the fact that its host may be more abundant than that of the resident population.
The argument of heterozygous disadvantage (which by definition does not apply if we consider non-interbreeding species in a community) is not a new one. Yet previous analysis of the phenomenon has primarily focused on generating statements about equilibrium conditions. Although, in natural landscapes or metapopulations that are dynamic and stochastic, where immigration of host and parasites from other sites occurs continuously, and where the abundances of the alternative host species may continuously change, it may be more informative to investigate how the presence of a heterozygote disadvantage affects the timing and dynamics of transitions from using one host to another. To explore these questions in more detail, we provide a simple model of a local host-parasite population that is linked to an external world by immigration.

General description
We develop a model for the population and genetic dynamics of a single population of a diploid parasite that may undergo coevolutionary alternation. Later on we introduce continuous immigration of parasite individuals from an external source, e.g. other surrounding populations; assuming ongoing mutations would have a comparable effect. We assume a simple genetic system where alleles h or c code for the specialisation on either of two possible host species (H or C); the role of the relative performance of heterozygotes w hc on either host's population dynamics is one of the critical parameters to explore. Our approach here differs from the study by Nuismer and Thompson (2006) in three important attributes: (i) association with different hosts occurs at random, i.e. parasites do not show a preference for one or the other host; consequently, after mating eggs or offspring are randomly associated with hosts in proportion to the hosts' abundance. Yet, model behaviour should be similar if parasites would need to search longer for a rare suitable host and thus loose net-fertility compared to a parasite phenotype that can infect a more abundant host. (ii) Traits affecting fitness on different hosts are discrete and cannot be expressed in full at the same time, and (iii) we assume diploid organisms with a single locus genetics and random mating between adult parasites. These assumptions appear justified in the case of Maculinea as these butterflies deposit eggs on specific host plants. Only later is each young larva carried into the nest of the first foraging worker of any Myrmica species that encounters it, regardless of whether it is a host ant or one of up to four non-host Myrmica species that commonly forage beneath the initial foodplants (Thomas et al., 2005b). This assumption is not undermined by recent research that indicates that Ma. arion uses host-plant defence chemicals to root disturbance as cues to bias oviposition to plants growing in the vicinity of any Myrmica nest (Patricelli et al., 2015). Therefore, phenotypes of Maculinea teleius and Ma. nausithous may bias oviposition to their shared food plants that grow near their respective host species by selecting growth forms that coincide with the host's optimum niches (Thomas et al. 2005a). The effect is likely not strong and there is no indication that Maculinea genotypes specialised on different hosts would differ intra-specifically in their egg-laying behaviour; such an effect could easily be incorporated, however, in the model that follows.
With the assumptions provided above and given that the relative proportion f of host H remains constant in time, the temporal dynamics of gene h in the parasite's population is defined by the following equation derived in full in Supplement A of the Supporting Information: Here, term p t is the proportion of allele h, q t = 1 − p t is the proportion of allele c in the parasite population at time t, w X,xy provides the density-independent probability of survival of an egg or young larvae of parasite of genotype xy on host X (either H or C), and s t, X is the density-dependent survival of larvae competing on host X that survived previous density-independent selection. We do not define an explicit population model for the hosts; later on, we consider, however, what would happen if, under the impact of parasitism or for any other reason, the relative abundance of different hosts would change in time. All further mathematical details of model derivation are outlined in Supplement A of the Supporting Information, but the model's properties can be summarised in the schematic Fig. 1: We can recognise three possible situations (I, II, and III) under the assumption that the proportion of hosts f remains fixed (more on this below): (I and I*) The overall fitness expectation w of one homozygote (yy = hh or cc) parasite is larger while that of the other homozygote is below fitness expectation of heterozygotes w yy (f ) ≤ w xy (f ) < w xx (f ). In this case, only one evolutionary outcome is possible independent of starting conditions: The homozygote with higher fitness expectation will ultimately dominate the gene pool. The gene can always invade a population dominated by the other allele (use of the alternative host) because it has a higher fitness even when rare (as a heterozygote). (II) If host proportions are so that the fitness expectations for both homozygotes are higher than that of the heterozygotes w xy (f ) < w yy (f ) ≤ w xx (f ) the evolutionary outcome will depend on the starting conditions: Neither allele could invade a population dominated by the other allele because the rare allele would mostly occur in heterozygotes while the dominant would mostly occur in fitter homozygotes. (III) The fitness of heterozygotes is above that of both homozygotes w xx (f ) ≤ w yy (f ) < w xy (f ). In this case, a stable coexistence in host use is predicted as both alleles could invade a population dominated by the other. A necessary condition for this to happen is thatat least in a certain range of values for fthe fitness of heterozygotes is above the intersection point of the two lines for the homozygotes' fitness (see Supplement B of the Supporting Information for formal proof).
Our model differs from simple textbook models only with regard to the additional effect of the relative proportion f of host H. The relative abundance of hosts operates as a weighting factor, especially relevant for the expected overall fitness of homozygotes. For this reason, it depends on the specific combination of values of f on the one hand and the expected fitness of different genotypes on different hosts (the set w X,xy of fitness values) on the other, whether (i) stable coexistence is possible, i.e. whether a classical heterozygote disadvantage or a benefit emerges ( Fig. 1b), and (ii) how a transition from using one host to using another might proceed. More trivial, situations are also conceivable, e.g. that heterozygotes always have a fitness benefit whatever the relative host abundance f. Note, however, that across a gradient in f we can only combine case (I) and (II) or case (I) and (III), but all three cases never simultaneously emerge for a given set of fitness values w X, xy . Figure 1c shows an extreme case where the heterozygotes never have a larger fitness than either of the homozygotes, i.e. a situation where only case (II) applies whatever the value of f.

Numerical simulations of dynamic system
As outlined, one of the two alleles would go extinct in a closed system if the relevant parameters were not to create a type III zone; this is only possible if fitness of heterozygotes is rather large on at least one of the possible hosts and relative host abundance f falls into an intermediate range. If regional genetic diversity in host use exists, however, immigration could maintain rare alleles at low levels in any local population (as could mutations). Further, the proportional abundance of hosts (f t ) or the susceptibility of hosts may change in time, not least due to the impacts of the established parasite type on its host (e.g. Ashby et al., 2019); corresponding evidence exists for the Myrmica ants infected by Maculinea (Hochberg et al., 1994;Elmes et al., 1996;Thomas et al., 1997). In general, however, there is little reason to believe that Maculinea parasitism has a relevant effect on Myrmica hosts' population dynamics or evolution . This is one reason why we abstain from implementing a full (e.g. Cortez & Weitz 2014) eco-evolutionary host-parasite model here that would allow us to determine the degree to which the host-parasite interaction itself might be responsible for inducing a cycle in host abundancewe will thus not provide answers to this question. Note that such a model would also need to be more complex and parameter rich with at least 2 × 3 additional parameters for describing the interaction of the three genotypes with the two possible hosts and 2-4 parameters for regulation of the two hosts' population dynamics. Our more limited goal is to describe the temporal dynamics of host switching in the parasite provided that host frequencies do change over time as a possible explanation for why multiple host use might be observed so rarely even though the opportunity (availability of alternative hosts) seems to be there. The reason for such changes in host abundance might indeed be the impact of the parasites (eco-evolutionary feedback), but it could also be due to changes in the abiotic environment, interactions between hosts, or the impact of other predators that attack the same host species. For the validity of the arguments that follow it is not important to know which factor may be responsible for generating the dynamics in host abundance or availability.
In the following, we will provideagain at a very basic levela scenario for the dynamics of host shifting in systems where host abundance f t fluctuates over time in a sinusoidal pattern with phase φ and where an equal and small number I hh = I cc of nonmated parasites of both homozygotes immigrates at any timestep. As said, we do this without defining a proper dynamics for the host population (but see discussion); we just investigate how a shift in host use would proceed provided that f t gradually changes in time. The details of implementation of this dynamic system are provided in the Supplement C of the Supporting Information.

Data evaluation
We iterate the dynamics of the system according to the equation system defined in the Supporting Information appendices A (a) (b) (c) Fig. 1. Schematic graph of expected fitness for homozygotes hh (red line), homozygotes cc (blue line) and heterozygotes hc (dark grey line) in dependence of the proportion f of host of type H and under the assumption that parasite genotypes are randomly paired with hosts. (a) Fitness of heterozygotes is never larger than that of both homozygotes at the same time, (b) fitness of heterozygotes may be larger than that of both homozygotes, and (c) fitness of heterozygotes is always lower than that of either homozygote. Three principal zones can be identified: Fitness of one homozygote is above and that of the other homozygote below that of heterozygotes (I and I*), fitness of both homozygotes is above fitness of heterozygotes (II), and fitness of both homozygotes falls below that of heterozygotes (III). The exemplary arrows in paled colours show how, in each region, in which direction fitness of allele h (red arrow) or c (blue arrow) would change with increasing proportion -arrows always start at the value for heterozygotes as rare alleles will only occur in heterozygotes. A classical heterozygote benefit only emerges in region III where the rare gene can always invade and stable coexistence is possible; outcome in region II depends on starting conditions (more details in text and Supplement B of the Supporting Information). and C over several complete cycles of host frequency f t with start values f 0 = f max and p 0 = 0.99 and by applying the modification for allele frequency and adult population size as specified in Supplement C of the Supporting Information, eq. (S9). We keep the value for f t fixed at f 0 for the first 50 generations to allow equilibrating of p t and only then implement the cyclic change in f t as defined by Supplement C (Supporting Information) eq. (S8); results for these first 50 generations are removed from data set before data evaluation. For presentation, we follow the temporal dynamics of the proportion p t of allele h in the population (measured at the end of the cycle, i.e. in the surviving and emerging larvae but before immigration) and parasite population size N t over the next 800 generations. For each simulation run, we then define the proportion of the total simulation time (generations) where the frequency of both alleles (h and c) falls into a coexistence range defined by the arbitrarily selected detection threshold d = 0.025; we assume that typical field investigations would fail, e.g. due to limited sample size, to provide evidence for genes that make up less than 2.5% of the gene pool. Choosing other values for d has no principal effect on results. We thus count any situation as coexistence where d ≤ p t ≤ 1 − d during the last two full cycles in f t . We further define the minimum values for p t and N t reached during the final simulation cycle.
We generate corresponding results for combinations of parameter values coming from the range w H,hc = w C, hc 2 [0… 1], net fertility (surviving egg number) R 2 [1.5…6.5], I hh / K = I cc /K 2 {0.001, 0.005, 0.025}, and the length φ of a full cycle through f max ! f min ! f max taken from the values φ 2 {200, 100, 50} generations; note that a transition from the maximum to the minimum value of f t takes half as long as these values, i.e. just 25 generations in the fastest cycle. The range of value for R covers a fair spectrum of net-fertility values reported for insects (Hassell et al., 1976), w X,xy values cover the whole possible spectrum and values for immigration rate, and I xx /K cover a range as it may occur in metapopulations without being so large as to create panmictic populations.

Results
All numerical simulations rapidly settle into stable cycles driven by the temporal dynamics in host abundance f t : after a maximum of four cycles, values for p t and N t estimated for f t values in the same phase (f t = f t + ϕ ) become identical, i.e. p t = p t + ϕ and N t = N t + ϕ . In Fig. 2, we provide some exemplary graphs for the scenarios where the fitness of homozygotes in the wrong host is always zero (w X,yy = 0) showing the dynamics in allele frequency and parasite population size as the proportion f t of host H fluctuates through time. We recognise that with low heterozygote fitness (w X,xy = 0.1) the parasite population is typically dominated by just one allele; transitions from using one host to using the other occur in very brief periods, i.e. within a few generations. This implies that field biologists would often enough fail to report coexistence of parasites specialised on different host even though both hosts are available at all times. Note also that the dynamics of these transitions provides a case of hysteresis: The proportion of host H has to drop to very small values before the transition to host C comes about, whereas for the backtransition to occur f t has first to increase to rather large values again. This further implies that at most times, the abundance of allele h does not reflect the relative abundance of host H (values in the f t -p t phase plot are typically far from the main diagonal; Fig. 2d.1). These effects are subdued as we enlarge values for w X,xy , that is, the transition from using one host to using the other becomes more gradual. With values of w X,yy = 0.7, the frequency of allele h more or less corresponds to the actual frequency of host H. We further note that periods of transition in host use may be associated with a substantial reduction in the parasite's population size, especially if w X,xy is small (Fig. 2c and Supporting Information Fig. A1).
The scenario presented in Fig. 2 corresponds to that shown in Fig. 1a for the case w X,yy < 0.5, i.e. the case where a stable equilibrium does not exist in the static case, and with that in Fig. 1b w X,xy > 0.5. In Supporting Information Fig. A2, we provide results for similar simulations but with the assumption that w X, xy = w X,yy . It is noteworthy that in these latter simulations the strong cycles in allele frequency p t and the hysteresis described above are maintained also for the high values w X,xy = w X,yy = 0.7. Further, a host shift occurs the later in the cycle, the larger w X,xy , whereas in Fig. 2b, we recognise the opposite trend.
In Fig. 3, we summarise the effects of heterozygote fitness w X,xy , parasite net growth rate R, proportional number of parasite immigrants (I/K of type hh and cc), and length φ of a full cycle through host abundance f for the fraction of generations where coexistence would be detected at a threshold of d ≥ 0.025. Coexistence is affected by all variables investigated with an interesting uni-modal response to changes in net fertility R: For values of R < 2.5, the fraction of generations with detectable coexistence declines with increasing R but beyond that value the fraction of generations with reported coexistence would gradually increase again. Massive immigration of homozygotes, i.e. c. 5% of maximum adult population size (lower row of panels in Fig. 3) also greatly increase the proportion of generations where both alleles would be found in the population (note that the proportion is measured at the end of the larval life cycle, i.e. after selection but before immigration). Finally, with faster host cycles, a larger proportion of generations will fall into the coexistence range. Overall, the preconditions to find coexistence (with the liberal definition applied here) are high heterozygote fitness in combination with either very low or rapid parasite growth rates and substantial parasite immigration (see discussion).
In the alternative scenario with w X,xy = w X,yy , we find a more complex pattern; note that this scenario corresponds to that shown in Fig. 1c with a zone II across the whole range of host abundance (unstable equilibrium). Especially with low immigration, we find large and separate domains in parameter space where coexistence of species would rarely be detected (Supporting Information Fig. A3). One domain (low heterozygote fitness) is in principle similar to that obtained in the previous simulations but a second with intermediate values for fitness w X,xy = w X,yy and rather large values of R emerges because a host transition does not occur at Fig. 3. Effect of heterozygote fitness w X,xy , parasite net fertility R, immigration, and host cycle length φ on the fraction of generations with coexistence identified at a detection level of d = 0.025. Immigration is identical for homozygotes of type hh and cc and expressed in proportion to the maximum habitat carrying capacity (K = 1000). w X,yy = 0 in all simulations. Coexistence was identified over the last two full host cycles, i.e. over the last 400, 200, or 100 generations out of the 850 generations iterated. all noticeable by the fact that p t values (initialised with p 0 = 0.99) never fall to values below 0.9 (results not shown). Finally, as w X,xy = w X,yy approaches the value of 1, all fitness values become similar and, with symmetric immigration as assumed here, stable coexistence emerges with p t values close to 0.5.

Discussion
The genetic model presented here is very simple, neither reflecting the specificities of Maculinea genetics nor accounting for many aspects of its ecology; in particular, we ignore possible eco-evolutionary feedback effects as was done in the model for haploid organisms published by Ashby et al. (2019). We believe, however, that we can take some messages from the model that should not critically depend on the simplifying assumptions we made here; in Supplement D of the Supporting Information, we discuss the effect of relaxing some of these assumptions. In addition, the argumentation provided will also apply if alternative hosts were not different host species but just distinct host genotypes, e.g. genotypes with different CHC profiles; such profiles play a great role in the nest-mate recognition of social insects (Howard & Blomquist, 2005). In this case, local defensive coadaptation in host species would be possible only if gene flow between genotypes is sufficiently low, as seems to be the case in Danish My. rubra populations that are host to Ma. alcon .
Clearly, the most critical component in our model is the fitness matrix defining performance of different genotypes on different hosts and thus the trade-off between specialising on one host and the fitness when utilising alternative hosts. For the standard scenario assumed our results provide a clear message. To observe stable coexistence in alternative host use, the fitness of heterozygotes must not be too low. Otherwise, the parasite population will at most times be dominated by one or the other strategy as stable coexistence is impossible. If the frequency of hosts (f) itself remains stable in time, a coexistence of strategies is only possible, if the parameter combination defines a zone III situation, as shown in Fig. 1b, i.e. if 1 − w < f < w with w H,hc = w C, hc = w and w > 0.5, in other words, if host abundances and fitness parameters are so that heterozygotes have an overall fitness benefit compared to either homozygous parasite. This may, for example, be the situation in case of Danish Ma. alcon populations that infect the very similar (in terms of CHC profiles) My. rubra and My. ruginoides but not the more dissimilar but also common My. scabrinodis that is the primary host in other regions (Als et al., 2002;Tartally et al., 2019). Similarly, in some Italian sites, ant species simultaneously exploited are quite similar in their CHC profiles, whereas Ma. rebeli did not survive in colonies of another principle host (in other European regions), My. scabrinodis . Another interesting case of multiple host use on a local site is the simultaneous use of My. scabrinodis and My. vandeli where the latter ant species is itself believed to parasitise My. scabrinodis colonies and thus to have evolved CHC profiles similar to that of My. scabrinodis (Tartally et al., 2019). In these cases, the heterozygote disadvantage may be low enough to allow for a stable mixed host use strategy.
More interesting, however, are the results from our dynamic simulations where we assume a gradual and cyclic shift in host relative abundance f t ; note that the change in f t might also reflect evolutionary changes in host accessibility or defence itself driven by the exploitation by parasites (Nuismer et al., 2005;Ashby et al., 2019). Periods of multiple host use may be observable in these scenarios, but typically the transition from using one host to using the other should be fast so that short-term (field) studies are unlikely to provide evidence for the coexistence of alternative host use strategies. Mechanisms that would promote the likelihood of (detectable) coexistence are fast parasite growth, fast host cycles (trivially so as transitions occur more frequently) and substantial gene flow into the population (see below). Although, for transitions in host use to occur, host abundance must show a considerable change in relative abundance. In the special case of the Myrmica-Maculinea interaction, it is unlikely that the population-wise effect of Maculinea on Myrmica populations is sufficient to create such dramatic changes : even though individual infected colonies may suffer severely (e.g. Thomas & Wardlaw, 1992;Wardlaw et al., 2000), only a small proportion of colonies is usually affected. Observed changes in Myrmica populations were rather a consequence of, e.g. land-use changes (Elmes & Thomas, 1992). Further, whenever coexistence of strategies is observed, it is likely that hosts are not utilised in proportion of their relative abundance (Tartally et al., 2019), an observation that agrees with the model behaviour shown in Fig. 2 and Supporting Information Fig. A2. Disproportional host use was, for example, observed for the Danish populations of Ma. alcon with multiple host use (Als et al., 2002) and in some Italian sites .
Of special interest in this context is the non-linear effect of parasite net fertility R on resulting coexistence (Fig. 3): at very low values of R, coexistence is possible because transition from one strategy to the other progresses just too slowly in relation to the host cycle in f, i.e. host cycles are too fast for the parasite system ever to settle into a stable equilibrium. Otherwise, with increasing R coexistence is more likely because high values of R benefit the rare gene more than the resident gene, as resource competition on the rarely affected host is lower.
In a similar vein, we note that the symmetric immigration assumed in our model disproportionately favours the rare gene immigration works like a numerical boost to fertility R with an increasing per capita effect the smaller the resident population size. This allows host switching at larger population size and makes observation of coexistence more likely. This may be the case in some areas with a regional mosaic in CHC-profiles of one host, My. rubra, infected by Ma. alcon in Denmark (Als et al., 2002;Nash et al., 2008). On the contrary, in Maculinea usually one host-use strategy dominates over very large regions (Sielezniew & Stankiewicz, 2002;Settele et al., 2005;Thomas et al., 2005a;Tartally et al., 2008;Witek et al., 2008;Patricelli et al., 2010). Immigration is then likely biased in favour of the dominant strategy and consequently rather diminishes the ability of a rare invading strategy to spread in a local population even if the corresponding host species was more abundant than the host used by the majority. To express this differently, in a landscape dominated by a single host-use strategy, establishment of an alternative strategy is less likely due to the fact that the majority of immigrants will follow the resident strategy, too. Emergence of a landscape level mosaic in host use could be promoted, however, if the parasites themselves impose a strong effect on their host population as compared to, e.g. the role of habitat Chaianunporn & Hovestadt, 2011); corresponding empirical evidence has been suggested for parasitic plants (Pennings & Callaway, 1996).
We must stress the difference between the coexistence of different host-parasite combinations in the community context as contrasted to multiple hosts use by a single parasite species (genetic diversity). Multiple host use within a parasite species leads to the formation of heterozygotes, whose fitness we extensively examined in our scenarios; we have demonstrated that the formation of heterozygotes may undermine the coexistence of multiple host-use strategies. Where multiple species pairs coexist, this question obviously never arises, because heterozygotes are never formed (see Rabajante et al., 2015). We can consequently conclude that mechanisms like assortative mating, whereby mating primarily occurs between individuals using the same host, or genetic dominance (Durinx & Van Dooren, 2009) would make local coexistence in host use more likely (Sorenson et al., 2003).
The population genetic dynamics principles we have presented here are as such simple and presumably apply also under different circumstances. Yet, we do not provide cues on how relevant the mechanism we provide might be in any particular system, for example whether host cycles and thus host shifting may be driven by the eco-evolutionary feedbacks resulting from the host-parasite interaction itself (Ashby et al., 2019). Note, however, that many of these models have been developed with haploid organisms respectively a quantitative genetics in mind that cannot be applied to the dynamics of discrete genotypes as we have assumed here (for overview, see Govaert et al., 2019;Lion, 2018), yet Lively (2012) provides an example of an simple ecoevolutionary model assuming a diploid two-allele genetics. Therefore, the population genetic effects described here do as such not depend on the mechanism underlying a cycle in hosts' abundance and/or susceptibility. Indeed, if a host woulde.g. due to changing habitat conditionsshift permanently from high to low frequency, the potential, subsequent host-shifting dynamics in a parasite population would approximately be described by just a half-cycle of the dynamics shown in Fig. 2 or Supporting Information Fig. A2.
Our modelling exercise was inspired by empirical observations of the carefully studied Myrmica-Maculinea association. We thus want to evaluate empirical observations in the light of our model. In particular, we think that the model provides a framework for interpreting the many and sometimes diverging observations of the Maculinea-Myrmica system (e.g. Thomas et al., 1998;2005a;Sielezniew & Stankiewicz, 2002;Settele et al., 2005;Pech et al., 2007;Tartally et al., 2008).
The most general statements of our model are (i) observations of local coexistence of alternative host-use strategies should be rare provided that the fitness of heterozygous Maculinea genotypes is rather low (and equally that of homozygotes on the alternative host); for this expectation to hold, it is not important whether the underlying genetics is as simple as assumed here as long as offspring of matings between alternative host-use strategies develop intermediate phenotypes. (ii) Where several hosts are successfully infiltrated, host use should typically not be in proportion to their relative abundance. In this context, we are not aware of any studies were Maculinea individuals have been genotyped, but both of the above predictions are supported by the comprehensive review recently published by Tartally et al. (2019): on about 70% of 419 Maculinea sites, only one host species was utilised despite the nearly universal presence of alternative hosts. And where alternative hosts were used, utilisation was rarely in proportion to host abundance, with a number of instances where the rare host was over-exploited (cf. fig. 5 in the study by Tartally et al., 2019). Further, most of the records of more than one host ant species producing adults Maculinea on the same site are of the predatory Ma. arion and Ma. teleius, and at least someprobably most of theseoccur because the main host was parasitised initially, but all the brood were eaten in the final spurt of growth causing the ant colony to desert its nest leaving the brood parasite behind. A neighbouring colony of another Myrmica spp may then bud off and occupy the vacant nest site, which contains Maculinea pupae (Thomas & Wardlaw, 1992;Thomas et al., 2005a).
In principle, parasites might evolve the ability to exploit several hosts without a relevant trade-off (Joshi & Thompson, 1995) as may be the case for many insect herbivores. In the specific case of Maculinea, and especially in the case of the very well integrated and interacting cuckoo species 2005a;Als et al., 2004), the route to successful adoption involves many steps of deceit including presentation of the right CHC profile (e.g. Thomas et al., 2005a;Nash et al., 2008), and behavioural and acoustic simulation (Barbero et al., 2009;Schönrogge et al., 2016). It seems unlikely that the full suite of these adaptations could be expressed phenotypically in two or more sets so that larvae could invade alternative hosts equally well, in particular if the alternative hosts are rather different in their attributes. Thus, to the best of our knowledge, this has not been investigated for Maculinea by testing the host-specificity of offspring bred from males and females that originate from populations using different host species either between or within sites. Moreover, hosts themselves would be under pressure to evolve contrasting recognition systems that cannot be mimicked at the same time Yoder & Nuismer, 2010;Chaianunporn & Hovestadt, 2011).
Interestingly, the predatory Ma. teleius is the Maculinea species with the strongest evidence for the local exploitation of multiple Myrmica host species, at least on a minority of local sites. On current knowledge Ma. teleius is also the least specialised of the Maculinea species due to its specific life style and developed defence mechanisms (Thomas et al., 2005a;. In the light of our model, this observation could be explained if the heterozygote disadvantage is lower in this species, resulting in a situation allowing for stable coexistence (type III zone in our Fig. 1). This is not necessarily an intrinsic attribute of the social parasite but could also emanate from the fact that its principal secondary host, My. rubra, is the most highly polygynous of the known ant host species, a situation that tends to result in a less precise chemical recognition system, making it more readily invaded by close but imperfect chemical mimics  (Gardner et al., 2007;Fürst et al., 2012). Further, adoption of Ma. alcon from sites where they use different hosts took longer than that of larvae from populations using only one hostand adoption of larvae taken from the same site as the ants used in the tests was fastest (Als et al., 2001); these observations indicate local adaptation of butterfly populations and a cost to a generalist strategy or a disadvantage for heterozygotes. The expression of a heterozygote disadvantage (and the disadvantage of homozygotes on wrong hosts) may also depend on the general well-being of individual host ant colonies and thus on environmental conditions. For example under laboratory conditions, where ants are both unstressed and supplied with food ad libitum, Myrmica ants are much more tolerant of infiltration by social parasites; it was only when ants were stressed that the Maculinea larvae were destroyed in colonies of host species to which they were not adapted Schönrogge et al., 2004;Thomas et al., 2013). Thus, unusually favourable conditions for the Myrmica host ants, in particular regions or in certain years, might also promote coexistence of alternative host use strategies.
As stated earlier, any mechanism that promotes genetic separation between host races, like assortative mating, would promote local coexistence in host use. In Maculinea, assortative mating might occur if, for example, individuals using different hosts tend to emerge on different dates or if different hosts were spatially separated, so that individuals using the same host are more likely to mate with each other. It is noteworthy in this context that Maculinea butterflies are known for their low mobility (Nowicki et al., 2005;Ugelvig et al., 2012) with evidence for the formation of home ranges even within rather small sites . Nonetheless, the fact that evidence for local coexistence of host races in Maculinea is rare tends to undermine the proposition of the existence of true but cryptic species' (see Als et al., 2004) that utilise alternative hosts.
Given the fact that Maculinea larvae typically inflict strong damage on host ant colonies (Thomas et al., 1997;2009), it remains a separate question as to why, even though theoretically expected, evidence for massive population cycles and coexistence of different strategies at the local or regional scale have rarely been reported. Several factors may be responsible. First, Myrmica colonies on a given site are not susceptible to Maculinea infiltration if their foraging ranges (>1-3 m) do not include the host plants that Maculinea need in their early life stages 2005a); as explained already before, a high proportion of host colonies are thus often located in safe sites and the overall effect of Maculinea on its host's population is usually small. Second, typical habitats may show spatial structure with a mosaic of habitat conditions favouring one or the other of the possible host species . Third, in cases where egg deposition on host plants is largely random, there would be an implicit decline in Maculinea net-fertility as soon as the spatial cover of host declines, as larvae are less likely to be found by host ants. This density-dependent effect may in turn reduce the pressure on hosts as hosts become rarer.
In summary, with our model, we want to exemplify that an analysis of the dynamic attributes of interactive systems may provide additional insights beyond those provided by equilibrium analyses alone. Here, we show that, given that certain assumptions are met, host shifting may rarely occur and only at a moment when an alternative host becomes (much) more abundant than the host currently used. Furthermore, transitions may then occur at such a fast rate that the likelihood for simultaneously observing use of alternative hosts may indeed be small. Rural Affairs (U.K.); Agence nationale de la recherche (France); Formas (Sweden); and Swedish Environmental Protection Agency (Sweden) through the FP6 BiodivERsA Eranet. We thank two anonymous reviewers for their valuable comments. We express our deepest gratitude to Graham Elmes, who contributed so much to our general understanding of ant ecology. This manuscript owes much to his persistent enthusiasm and he will be sadly missed.

Conflict of interest
The authors have no competing or conflicting interests to report.

R
net fertility of parasites f t relative abundance of host H (abundance of host C is 1 − f) p, q = 1 − p proportion of alleles h and c in parasite population N t number of adult (mating) parasites in population Z t number parasite zygotes (eggs) of genotypes hh, hc, and cc generated (eq. S1) Z t, X number parasite zygotes associated with host X, X 2 (H, C) (eq. S2) L t, X number parasite larvae adopted by host X (eq. S3) L 0 t, X number parasite larvae surviving densitydependent competition on host X (eq. S5) K carrying capacity of parasite population (defined by hosts' abilities) w X,xy adoption probability (fitness) of parasite with genotype xy 2 (hh, hc, cc) on host X φ length of full host cycle from f max ! f min ! f max I xx number of homozygous immigrants of genotype xx

Supporting information
Additional supporting information may be found online in the Supporting Information section at the end of the article.

Figure A1
Effect of heterozygote fitness, net growth rate, immigration, and host cycle length on the minimum parasite population size reached during a full cycle of host abundance. Results from the simulations shown in Figure of main text. Maximum carrying capacity for parasites is K = 1000. Figure A2. Dynamics of host transition under fluctuating relative abundance f of host H. All parameters and settings as in Figure of main text but here with w X,xy = w X,yy throughout. Blue dots show starting point in phase diagrams, red dots end-points. Figure A3. Effect of heterozygote fitness w X,xy , net growth rate R, immigration I cc /K, and host cycle length on the fraction of generations with coexistence reported at the detection level 0.025. All settings identical to Figure except that here w X, xy = w X,yy .