The maintenance of genetic diversity under host–parasite coevolution in finite, structured populations

As a corollary to the Red Queen hypothesis, host–parasite coevolution has been hypothesized to maintain genetic variation in both species. Recent theoretical work, however, suggests that reciprocal natural selection alone is insufficient to maintain variation at individual loci. As highlighted by our brief review of the theoretical literature, models of host–parasite coevolution often vary along multiple axes (e.g. inclusion of ecological feedbacks or abiotic selection mosaics), complicating a comprehensive understanding of the effects of interacting evolutionary processes on diversity. Here we develop a series of comparable models to explore the effect of interactions between spatial structures and antagonistic coevolution on genetic diversity. Using a matching alleles model in finite populations connected by migration, we find that, in contrast to panmictic populations, coevolution in a spatially structured environment can maintain genetic variation relative to neutral expectations with migration alone. These results demonstrate that geographic structure is essential for understanding the effect of coevolution on biological diversity.

work within a meta-population framework (e.g. Gandon, 2002;Gandon et al., 1996;Gandon & Michalakis, 2002;Gavrilets & Michalakis, 2008). Specifically, coevolution has been repeatably demonstrated to result in host local (mal)adaptation, meaning hosts are more (or less) fit in their local environment than in non-local environments (Nuismer, 2006;. As local (mal)adaptation requires the presence of genetic variation, these results suggest that there remain open questions as to how genetic diversity is maintained across geographic scales. Here we explore if, and how, spatial structure maintains host genetic diversity by comparing the transient and long-term dynamics of expected heterozygosity (Gillespie, 2004, p. 15), a measure of diversity at a single locus, across multiple geographic scenarios.
Geographic structure has a long-appreciated role in shaping coevolutionary dynamics (Lion & Gandon, 2015). Interactions between co-evolutionary dynamics and spatial structure are central to the geographic mosaic theory of coevolution (Gomulkiewicz et al., 2007;Thompson, 2005). Although the majority of co-evolutionary models have assumed well-mixed populations (Buckingham & Ashby, 2022), numerous theoretical models have explored how coevolution is affected by spatial structure in the form of a meta-population 'island model' (Frank, 1991b, Gomulkiewicz et al., 2000, Lively, 1999, Nuismer et al., 1999, in population lattices (Gandon et al., 1996) or across continuous space (Ashby et al., 2014;Lion & Gandon, 2015).
Geographic structure tends to promote fluctuating selection (Gómez et al., 2015), greater host resistance (Ashby et al., 2014) and lower parasite infectivity (Best et al., 2011). Spatial structure -particularly in the form of a meta-population -also results in local adaptation (Gandon et al., 1996;Nuismer, 2006;Thrall et al., 2016) and thus the maintenance of genetic variation. Nevertheless, the dynamics of heterozygosity in spatially structured co-evolving systems has yet to be rigorously compared to neutral expectations for the dynamics of genetic diversity in the absence of reciprocal natural selection. Here we perform such a comparison to characterize the contribution of host-parasite interactions on genetic diversity. Specifically, we ask In contrast, allele frequency dynamics in the MA model are characterized by classic Red Queen allele frequency cycles (Woolhouse et al., 2002), resulting in a 'dynamic polymorphism' associated with (indirect) negative frequency-dependent selection (Clarke, 1979;Takahata & Nei, 1990;Tellier et al., 2014). These Red Queen cycles are mathematically characterized as 'neutrally stable', making them particularly sensitive to model assumptions. As such, it is difficult to directly compare existing theoretical results between models with and without spatial structure, namely MacPherson et al. (2021a, 2021b)) and Gandon et al. (1996), Gandon, 2002), as these models differ along multiple axes known to impact co-evolutionary dynamics. To understand the interplay between multiple evolutionary processes, we begin by summarizing six evolutionary axes, described in detail in the Theoretical background section, that have previously been shown in the theoretical literature to qualitatively affect the outcome of coevolution. Due to varying model assumptions, a robust analysis or tabulation of how each of these axes impacts diversity is challenging to obtain from existing literature. We therefore focus our own models on two specific axes that have qualitative impacts on antagonistic coevolution, namely spatial structure (i.e. single patch vs. two patch vs. island-mainland) and the synchrony of life history (i.e. continuous vs. discrete time). Varying and contrasting spatial structure allows us to examine co-evolutionary diversity in the presence (two-patch and island-mainland models) and absence (single-patch model) of migration and how the effect of migration depends on the genetic composition of the migrants. Namely, whether the co-evolutionary dynamics shaping migrant genetic diversity are coupled (two-patch model) or uncoupled (island-mainland model) from the genetic diversity in the focal patch. As discussed in more detail below, the In our summary of all six axes, we describe how our models fit into this literature and illustrate how the model comparisons developed below are designed to disentangle the role of geographic structure on the maintenance of genetic diversity.

| The genetic basis of coevolution
As noted earlier, the transient and long-term outcome of coevolution differs between models of the genetic interaction. The GFG model resembles the asymmetric interaction between resistance and virulence genes found in many plant-pathogen systems and was first formulated to reflect flax-rust interactions (Flor, 1956). In this model, hosts carry either 'susceptible' or 'resistant' alleles. Susceptible hosts can be infected by 'avirulent' and 'virulent' parasite genotypes, whereas resistant hosts can be infected by only the virulent parasite genotype. The maintenance of genetic variation in GFG models is well characterized as requiring a trade-off or cost associated with parasite infectivity and host resistance (Salathé et al., 2005). In the absence of such costs, the most infectious parasite and most resistant hosts will fix in the population. The maintenance of genetic variation is supported by various factors, including genetic factors (e.g. overdominance in diploids); asynchrony in host and pathogen allele frequency dynamics (e.g. spatial asynchrony from population structure or temporal asynchrony arising from seed dormancy) and eco-evolutionary or epidemiological evolutionary feedbacks .
The matching alleles (MA) model is akin to a lock-and-key molecular interaction thought to underlie many animal-pathogen immune interactions (Dybdahl et al., 2014). In the MA model, host and parasite genotypes are paired such that hosts are susceptible to matching parasite genotypes and experience reduced susceptibility to mismatching ones. The symmetry of the MA model makes the coevolutionary dynamics sensitive to model assumptions and interacting evolutionary and ecological processes, hence it is more difficult to characterize than the GFG model and thus it will be our focus here. Irrespective of particular model assumptions, the host and parasite allele frequencies in the single locus-two allele MA model generally exhibit five equilibria (Tables 2.4-2.6 in Nuismer, 2017). First there is one polymorphic equilibrium where the host (p H ) and parasite (p P ) allele frequencies are both at 0.5. The other four equilibria characterize the fixation or loss of the focal allele in the host and parasite [(p H = 0, p P = 0), (p H = 1, p P = 0), (p H = 0, p P = 1), (p H = 1, p P = 1)].
However, the stability of these equilibria, and hence the evolutionary dynamics, do depend on model assumptions. Below we summarize previous results of assumption-dependent MA model dynamics and their implications for the maintenance of genetic diversity.

| Discrete versus continuous time
Depending on the life history of the focal system, models may be specified in either discrete or continuous time (Nuismer, 2017) with discrete-time models reflecting synchronized reproduction and casualty events and asynchronous life histories most adequately captured by a continuous-time model. The qualitative behaviour of the MA model differs between discrete and continuous time. The discretetime MA model is characterized by an unstable polymorphic equilibrium from which the system cycles away, resulting in the deterministic loss of genetic variation in the absence of mutation (MacPherson et al., 2021a;M'Gonigle et al., 2009) or other stabilizing forces (MacPherson et al., 2021b). The continuous-time MA model results in neutrally stable cycles (Woolhouse et al., 2002) that neither decrease nor increase genetic variation (as measured by expected heterozygosity) in the absence of genetic drift, predicting that heterozygosity at these loci should be lost at the same rate as under neutrality. As we lay out below, however, finite population size can qualitatively alter co-evolutionary dynamics. Given the qualitative dependence of evolutionary dynamics on life history, our model comparison contrasts the maintenance of diversity in discrete versus continuous time.

| Finite population size
Modelling coevolution in a finite population, recent work by MacPherson et al. (2021a) found that, compared to the neutral expectation, MA coevolution depletes genetic variation. Indeed, previous work suggests that indirect, NFDS does not maintain genetic variation in the absence of a high degree of mutation or migration (Agrawal & Lively, 2002;Borghans et al., 2004;Ejsmond & Radwan, 2015;Frank, 1991a;Frank, 1993;Lively, 1999;M'Gonigle et al., 2009). In addition, genetic drift in finite populations eventually results in the fixation of an allele in either the host or parasite.
Although selection in the MA model is symmetrical when both species are polymorphic, following fixation in the parasite (host), the MA model results in directional selection favouring the mismatching host (or matching parasite) type, rapidly depleting variation.

Models of coevolution in spatially structured populations have
also emphasized the importance of finite populations. When the local population is small (and selection is relatively weak), selection cannot counter drift, resulting in the loss of variation; conversely, when the population size is large, drift does not overwhelm selection and diversity is maintained (Gandon & Nuismer, 2009;Nuismer et al., 2010). As in MacPherson et al. (2021aMacPherson et al. ( , 2021b, here we focus on the dynamics of heterozygosity in finite populations relative to the corresponding expectations under neutral drift.

| Population dynamics
In addition to host-parasite coevolution, MA models may include additional population dynamics, such as ecological or epidemiological feedbacks, which change the number or density of individuals in a given population across time. When population dynamics are included, selection is both frequency and density dependent, resulting in feedbacks between evolutionary and population dynamics (Ashby et al., 2019). Although the inclusion of ecological or epidemiological evolutionary feedbacks may have little effect on the population dynamics (MacPherson & Otto, 2018;Nuismer, 2017), there are often large qualitative and quantitative impacts on host-parasite coevolution (Frank, 1991a), stabilizing polymorphisms (Gandon et al., 1996;MacPherson & Otto, 2018) and hence increasing heterozygosity in infinite (Ashby et al., 2019) and finite (MacPherson et al., 2021b) populations. Despite their potential importance, here we focus on evolution in spatially structured environments in the absence of ecological feedbacks. This analysis provides the necessary foundation with which future models can examine eco-evolutionary feedbacks.
Since migration is expected to have a stabilizing effect, we would expect the qualitative change in stability due to migration in the non-ecological model, which is neutrally stable in the absence of migration, to be similar to that of an ecological model, which is already weakly stable in the absence of migration.

| Evolutionary potential
Variation in the extent of local adaptation in co-evolving populations is often attributed to the relative evolutionary lability of the interacting species. Determinants of lability include the relative strengths of natural selection (interaction specificity and virulence; Gandon, 2002), the generation of novel genetic variation (mutation and recombination rates; Gandon & Michalakis, 2002), and relative generation time (life span and recovery rate; Gandon & Michalakis, 2002). In the context of finite populations, generation times impact the rate of drift between the interacting partners. In the models presented below, generation time is determined by the natural death rates of the host and pathogen.

| Abiotic selection mosaics
An important consideration in the context of the geographic mosaic theory of coevolution is abiotic selection. Spatial heterogeneity in abiotic selection can increase global and local polymorphism (i.e. genetic variation; Frank, 1991b. When parasite migration is low relative to the host, such abiotic selection mosaics favour host local adaptation (Kirkpatrick & Barton, 1997). Previous work by Nuismer (2006) suggests that for local adaptation to be permanent, both hosts and parasites must experience a selection mosaic; however, it is unclear from this work if local adaptation is permanent in the absence of mutation and in the presence of genetic drift (Gandon & Nuismer, 2009;Kirkpatrick & Barton, 1997). As with population dynamics, while we do not consider abiotic selection mosaics here, the results below provide a foundation on which interactions between biotic and abiotic drivers of diversity can be situated.

| ME THODS AND RE SULTS
For our work, we choose to focus on two particular axes highlighted by the theory summarized above to have qualitative impacts on antagonistic coevolution: spatial structure and the synchrony of life history. In particular, it remains unclear how geography and time (discrete vs. continuous) contribute to the maintenance of genetic diversity relative to neutrality and how these two factors interact.
We test the role of these axes for the maintenance of diversity under host-parasite coevolution using six models ( Figure 1). We consider three forms of population structure: (1) a single patch, (2) two patches connected via migration and (3) an island-mainland framework. We examine each of these population structures in both continuous and discrete time.
Understanding the consequences of host-parasite coevolution for the maintenance of genetic diversity requires defining a quantifiable metric of genetic diversity. Here, we consider the expected heterozygosity at a single haploid, biallelic locus. If the frequencies of the two alleles are p and q = 1 − p, the expected heterozygosity, simply referred to as heterozygosity here, is given by H = 2pq.
Heterozygosity is maximized when allelic diversity is maximized and p = q = 0.5. By using heterozygosity as our measure of genetic diversity, we are by extension considering only 'static diversity', i.e. the coexistence of multiple genotypes at a given point in time, and the comparison to the analogous neutral model. This exists in contrast to temporal diversity, which captures the cycling through different predominant genotypes across time.
For the sake of simplicity, we focus our attention solely on host heterozygosity with the expectation that parasite diversity should behave similarly given the symmetry of the MA model. Below, we discuss each of the six models in turn, with a focus on answering two specific questions: under what conditions is the equilibrium host heterozygosity non-zero, and when does coevolution slow the transient loss of heterozygosity relative to neutral drift?
For each model, we begin by analysing the co-evolutionary dynamics in the infinite population size limit. The stability of the F I G U R E 1 Overview of the six matching alleles models considered and their main results for the maintenance of host genetic diversity. For a given formulation of time-discrete or continuous-we consider three forms of population structure: one patch, two patches connected via migration and an islandmainland connected by migration. Bolded text in the top half of the figure indicates results in which coevolution maintained diversity, either transiently (relative to neutrality) or in the long term. Note that the two-patch, continuous-time model maintains diversity only transiently under the conditions stated. The bottom half of the figure highlights model features, including host-parasite coevolution, spatial set-up and time scale. In the co-evolutionary model, hosts and parasites undergo both extrinsic birth and death as well as death and subsequent birth resulting from successful infection. Geography is defined as one of three options: one patch, two patches with migration and island-mainland with migration. Time scale is either continuous or discrete, with models correspondingly specified as either Moran or Wright-Fisher models. In the Moran model, a single generation is defined by the time needed for N death events, where N is host population size. In the Wright-Fisher model, a single generation occurs as each of N individuals simultaneously die and are replaced. In both the Moran and Wright-Fisher model schematics, grey circles are reflective of individual hosts, where crosses indicate a host death. Purple and green arrows represent extrinsic birth-death and infection induced birth-death events for the host. corresponding deterministic models provides useful analytical insights into stochastic behaviour and hence can illuminate the maintenance versus loss of diversity in finite populations. As illustrated by MacPherson et al. (2021aMacPherson et al. ( , 2021b), deterministic models characterized by stable polymorphic equilibria are expected to lead to maintenance of diversity in the corresponding finite population (i.e. stochastic) models. Similarly, neutrally stable polymorphic equilibria are expected to behave like genetic drift and unstable polymorphic equilibria are expected to deplete genetic diversity (i.e. expected heterozygosity).
Each model below is specified in terms of biological parameters (e.g. the strength of selection, migration rates), however, as we illustrate with the one-patch model, the effect of these parameters on the stability properties are highly confounded. Hence, we analyse the behaviour of each model below in two complementary ways. First, we take a mathematical approach, deriving a set of 'composite' parameters that uniquely characterize the stability properties of the deterministic model but have no direct biological interpretation. By varying these composite parameters, we are able to capture the full range of co-evolutionary dynamics in the system and formulate a comprehensive description of the system's behaviour. We then take a biological approach varying evolutionary rates and hence map the mathematical behaviour to an underlying biological mechanism.

| One-patch model
We begin by reviewing the dynamics of host genetic variation in a single patch originally presented in MacPherson et al. (2021a) as it provides important introduction to the modelling design and contrast to the spatially structured models presented below. The onepatch Moran model (Moran, 1958) tracks the frequency of hosts (p H ) of type 1and parasites (p P ) of type 1, both with a single biallelic locus in a panmictic population with constant, finite host and parasite population sizes. Note the frequencies of type 2 hosts and parasites are q H = 1 − p H and q P = 1 − p P respectively. The evolutionary dynamics are characterized by three processes (reflected by the individual terms in Equation 1a, 1b below): genotype-dependent matching alleles infection (occurring at rates X and Y for matching and mismatching host and parasite genotypes respectively), natural host mortality (occurring at rate ) and parasite mortality (occurring at rate ). The first of these processes results in reciprocal natural selection driving coevolution, whereas the second determine the relative generation time of the two species.
A natural time scale on which to analyse the dynamics of genetic variation is 'host generations' as it allows us to directly compare heterozygosity in the co-evolutionary model to the expectation under neutral genetic drift. In a continuous-time stochastic Moran model, a generation can be defined as the expected time until H birth/death events, where H is the host population size. Given the frequency-dependent nature of coevolution, however, this time will depend not only on the rate of natural host and parasite death ( and ), but also on the rate of virulent infections and hence on the frequency of (non-)matching hosts and parasites. For additional details on rescaling time into host generations, see Appendix S1.2 (Supplemental Materials).
Although we are interested in the stochastic co-evolutionary dynamics in a finite population, the characterization of the dynamics in the infinite population size limit can provide valuable analytical insights into the behaviour of the system. In our model below, α gives the probability of host death given infection, X and Y give the probabilities of successful infection of a host by matching and nonmatching parasite genotypes, respectively, and D 1 is the parameter rescaling time into units of host generations, which depends on both selective (infection) and neutral (extrinsic birth-death) events (see Appendix S1.2). The terms in parentheses reflect the natural death and infection events (see Figure 1). To keep the total population size constant, all birth events are coupled with a corresponding death and vice versa. Natural turnover occurring at rate ( ) in the host (parasite) leads to a coupled death/birth event in the host (parasite) with no net effect on the allele frequency. Matching (mismatching) infection occurs at rate X (Y) with four subsequent coupled events of host death/birth and parasite birth/death occurring at rate .
Such coupled events can lead to either a net gain or net loss of the focal allele depending on the host and parasite allele frequencies.
Formulating the analogous deterministic system, we track the frequency of type 1 hosts p H and type 1 parasites p P : After simplification, and noting that terms related to natural mortality (i.e. those containing and ) exactly cancel out, we arrive at The term (X − Y)  Figure S2). In addition, coevolution transiently depletes host diversity faster than drift alone due to strong directional selection following allele fixation in the parasite (  Figure S2).
Although the internal equilibrium exhibits neutrally stable cycles, the boundary equilibria are unstable (i.e. possess positive leading eigenvalues; see Appendix S1.7). We find the extent of (in)stability and frequency of the cycles, which are characterized, respectively, by the real and imaginary components of the leading eigenvalues for each of these equilibria (Otto & Day, 2007) Table 1). We find that genetic variation is depleted relative to the neutral expectation at similar rates across parameter space and does not depend significantly on the stability properties of the equilibria (Figure 2a,b). This is expected as -in the absence of migration -there is little opportunity for perturbations in allele frequency near the boundary equilibria. Although this analysis reveals that genetic diversity is insensitive to the stability properties ( Figure 2b), it provides a useful comparison to the behaviour we find in the spatially structured models where rare variants are often introduced via migration.

| Two-patch model
We introduce spatial structure by following the dynamics of coevolution in two discrete patches (demes) connected by migration. We assume that migration is symmetrical between patches with hosts and parasites dispersing to the other patch at rates and respectively. As with the one-patch model, we complement simulations of the stochastic dynamics with an analysis of the

TA B L E 1
The composite parameters and corresponding model parameters for each model, as well as the number of replicates and generation sampled for simulations. To test the consequence of changing a particular composite parameter for stability, and thus maintenance of genetic variation, we solve for the associated model parameter value that gives us particular composite parameters (e.g. those plotted in Figures 2, 4) using the function given, where the function is the model parameter written as a function of the composite parameters. We then simulate 1000 replicates for each composite parameter value and sample at the host generation indicated. Generations were chosen such that transient dynamics reached equilibrium prior to sampling.

Continuous
One-patch

Discrete
All three models a i in , nm S System 3 exhibits the same five MA equilibria as in system 2 (Appendix S1.5). The stability properties of these equilibria are, however, now characterized by three composite parameters (given in  Figure 2c).
In terms of the evolutionary parameters, coevolution maintains genetic diversity above that seen under neutrality at intermediate levels of migration coupled with strong coevolution (Figure 3b). The maintenance of diversity is a consequence, in part, of the average effect of a migrant on diversity and the frequency at which migration events occur. When host migration rate is relatively high, allele frequencies in the two patches are tightly coupled, hence on average a migrant will introduce little diversity (i.e. will not perturb the system away from a boundary equilibrium, an important force according to our search of parameter space). When migration is rare, allele frequencies between the patches are less coupled such that there are fewer opportunities to move away from the boundary equilibria (and subsequently increase diversity). Although increasing migration homogenizes the two patches, increasing selection drives divergence in local allele frequencies. We thus find the average rate of introduced diversity is highest when migration is intermediate and selection is strong (Figure 3d). The greater difference between allele frequencies under strong coevolution, relative to weaker coevolution, can be explained by considering the potential for divergence in allele frequency trajectories between synchronizing migration events.

| Island-mainland model
Because the dynamics of heterozygosity in the two-patch model are determined by the auto-correlation of allele frequencies between patches, the results of this model may not extend to systems with more complex geographical distributions. In this section, we consider an alternative geographical model, the island-mainland model, which exhibits, by definition, no spatial auto-correlation (Smith, 1970;Wright, 1943).
In this model, we track allele frequency dynamics in a single 'island' population that receives host and parasite migrants from a 'mainland' population at rates and respectively. Although the TA B L E 2 Summary of results for each of the six models considered. Co-evolutionary strength is captured by (X − Y), where X and Y are the infection rates for matching and mismatching host and parasite genotypes respectively. We consider the dynamics of genetic variation relative to neutrality.

F I G U R E 2
Maintenance of genetic diversity in continuous time depends on the stability properties of the deterministic equilibria. In the one-patch model, there is a non-linear relationship between the composite parameter i(λ in ) -which describes the cycling rate around the polymorphic, internal equilibrium -and the relative maintenance of genetic diversity, H coev − H neut (a), but this relationship is not significant (b). (c, d) In the two-patch model, the relative maintenance of genetic diversity depends on three composite parameters: internal stability, internal cycling and matching boundary stability. The most diversity is maintained relative to neutrality with greater internal equilibrium stability, stronger internal equilibrium cycling and greater matching boundary equilibria instability. (e, f) For the island-mainland model, the most diversity is maintained with strong stability at the internal equilibrium and greater instability at the boundary equilibria. For all replicates, α = 1 and δ = γ = 0.5. In the one-patch and island-mainland models, κ H = κ P = 300, whereas in the two-patch model κ H = κ P = 150. across composite parameter space, we find that genetic variation is maximized when the internal equilibrium is most stable and the cycles are most frequent (Figure 2e,f). The effect of the stability properties on diversity are, however, non-linear and hence may be best described in the biological context.
Like the two-patch case, the internal equilibrium is stabilized by sufficient migration (Table S3). Heterozygosity is initially depleted relative to neutrality; however, coevolution coupled with strong migration results in a net maintenance of diversity (Figure 3d,e).
The extent to which heterozygosity is maintained increases as the migration rate increases and the strength of selection decreases The contrasting effect of selection in the island-mainland versus the two-patch model is particularly intriguing. In the two-patch case, stronger co-evolutionary selection resulted in more diverse migrants being introduced via migration, maintaining more diversity than weaker selection. In the island-mainland case, the system spends more time near the boundary equilibria, at which two processes will occur: (1) directional selection depleting variation (resulting following fixation in the parasite) and (2) perturbations from the boundary equilibria via the introduction of diverse migrants. The strength of directional selection in the former process increases with the strength of coevolution, hence strong directional selection can lead to less diversity near these boundary equilibria. In the absence of coevolution, the system will not be pushed away from the boundary equilibria following a perturbation event, such as a migration event (the latter process). Thus, weaker, but non-zero, coevolution maintains more diversity than strong coevolution and neutrality. Below we test the robustness of these results to assumptions in life-history synchrony.

| One-patch model
As in the continuous-time case, we begin with an overview of the dynamics in a single patch. We now use a Wright-Fisher model  (Fisher, 1930;Wright, 1931) to follow the frequency of host (p H ) and parasite (p P ) type 1 genotypes at the single biallelic locus in the single patch with finite and constant host and parasite population sizes in discrete time. Following the notation of Gandon (2002), in which the analogous two-patch model was first described, in a given generation, host and parasite fitness are determined by two parameters, specificity S and virulence V. Specificity determines the probability of successful infection given contact between a non-matching host and parasite, whereas virulence gives the probability of infection induced mortality following successful infection. The continuoustime infection parameters X, Y and can be related to S and V, with For derivation of these relations, see Appendix S1.4.
The resulting recursion equations for the host and parasite allele frequencies are given by: where After simplification, we have: For a full derivation of the above expressions, see Appendix S1.3.1.
Note system 6 has the classic form of p(t + 1) = p(t)s W , where s is the selection coefficient and W is mean fitness.
System 6 exhibits the same five MA equilibria as discussed previously (Appendix S1.6; Table S4, Supplemental Materials). coevolution depletes genetic variation across parameter space. As with the continuous-time models, the extent to which heterozygosity is lost can be understood in terms of either the stability (mathematical approach, Table 1) or the two biological parameters S and V, providing complementary insights. Mathematically, the primary determinant of genetic variation is the instability of the internal equilibrium with variation decreasing with increasing instability (Figure 4a,b). The more unstable the internal equilibrium, the closer the limit cycles will be to the boundary of the state space, hence reducing the time-averaged heterozygosity. Recapitulating the results of the continuous-time, island-mainland model, increasing instability of the boundary equilibria increases heterozygosity as repulsion from the boundaries will push the limit cycles towards the centre of state space. Parallels between the continuous-and discrete-time models can also be expressed in terms of the biological parameters S and V. Increasing either specificity, S , or virulence, V, leads to stronger directional selection following pathogen fixation and hence a more rapid depletion of heterozygosity ( Figure S3).
Although Figure 4a shows the difference in heterozygosity between coevolution and neutrality at generation 125, at earlier generations, intermediate stability and non-matching boundary instability result in a net gain in diversity (Figure 4b). Early on (prior to fixation in any species), the instability of the internal equilibrium pushes the host and parasite allele frequencies towards the boundary such that the boundary instability becomes important. Since the magnitudes of matching and non-matching instability are inversely related, the values of S and V that produce sufficiently high boundary instability at both types of boundary equilibria should maintain the most diversity. This occurs when the internal stability is maximized and boundary stability minimized.
As we expect at later generations, the instability of the nonmatching boundary determines the dynamics of diversity as this quantity mediates the strength of directional selection in the host following pathogen fixation.

| Two-patch model
The two-patch model is identical to the one-patch model, with the addition of migration of hosts and parasites between two patches.
We thus now need to track the absolute fitness and allele frequency changes for hosts and parasites in both patches, with allele frequency changes resulting from both selection and migration following selection. The probabilities for hosts and parasites migration per generation are given by m H and m P . Allele frequencies in patch i at the start of generation t + 1 are: where i = 1, 2, i ≠ j and p s Hi (t) and p s Pi (t) are the alleles frequencies for hosts and parasites following selection at time t. For a full derivation of the above expressions, see Appendix S1.3.2.
The stability of the equilibria in our two-patch model is the same as in the one-patch case, with the internal equilibrium and all boundary equilibria being unstable (Table S5). Again, coevolution depletes diversity relative to neutrality and, as in the one-patch case, increasing specificity and virulence increases the rate of loss ( Figure   S4). Similar to the two-patch, continuous-time model, intermediate migration maintains more diversity than weak or strong migration ( Figure S5). The main determinants of genetic variation are (1) the stability of the internal equilibrium in and (2) the instability of the non-matching boundary equilibria nm (whose magnitude is inversely related to that of the matching boundary equilibria; Table 1). Similar to the one-patch, discrete-time model, increasing internal instability will push the limit cycles towards the boundaries of state space, reducing heterozygosity, whereas increasing boundary equilibria instability will push the system away from the boundaries, maintaining heterozygosity (Figure 4c,d) with a net result that diversity (as measured by heterozygosity) is minimized when the internal equilibrium is most unstable (i.e. large leading eigenvalue) and non-matching boundary equilibria least unstable (i.e. small leading eigenvalue greater than 1).

| Island-mainland model
Like in the continuous-time, island-mainland model, we consider allele frequency dynamics in a single 'island' population that receives host and parasite migrants from a 'mainland' population with probabilities m H and m P respectively. Again, we track the allele frequency changes in the island population only. Allele frequencies at the start of generation t + 1 are For a full derivation of the above expressions, see Appendix S1.3.3.
System 8 has only one (internal) equilibrium with the boundary equilibria appearing only in the limit of no migration (Table S6). The maintenance of diversity is determined, like the other discrete-time systems, by stability of the internal equilibrium and instability of the non-matching boundary equilibria in the limit of no migration (Table 1). We observe a net depletion of diversity relative to neutrality across parameter space, particularly with greater instability at the internal equilibrium ( Figure 4e). The effect of non-matching boundary instability is weak but aligns with the pattern seen in the other discrete-time models (Figure 4f). This weak effect is expected, as the constant input of diverse migrants maintains the population near the internal equilibrium such that dynamics at the boundary equilibria are of little importance. We observe a positive relationship between migration rate and long-term heterozygosity: strong migration maintains heterozygosity close to perfect polymorphism (H = 0.5; Figure S7). With insufficient but non-zero migration, we observe a greater decline in heterozygosity under both neutrality and coevolution, with a faster initial decline under coevolution. This initial difference between coevolution and neutrality is due to the greater instability at the internal equilibrium caused by coevolution.
With this weak, non-zero migration, both neutrality and coevolution approach new internal limit cycles.

| DISCUSS ION
Here, we explored the maintenance of diversity, as measured by expected heterozygosity, under a matching alleles model of hostparasite coevolution in finite populations with various forms of spatial structure, considering three different forms of population structure (one-patch, two-patch and island-mainland) in both discrete and continuous time. We analysed our models using complementary deterministic and stochastic methods. Exploration of equilibria stability conditions allowed us to identify which specific biological processes (e.g. migration vs. coevolution), in the deterministic limit, determined the direction and magnitude of stability. Results for the onepatch model, for which we observe a primary effect of internal equilibrium instability and a secondary effect of non-matching boundary instability. Maintenance of diversity under coevolution -relative to neutrality -is maximized with lesser internal equilibrium instability and greater boundary equilibrium instability. (c, d) Results for the two-patch model, for which we observe similar dependencies on the internal equilibrium stability and non-matching boundary equilibrium stability. (e, f) For the island-mainland model, we observe a weaker effect of the non-matching boundary equilibrium. As host heterozygosity remains closer to the internal equilibrium (i.e. host heterozygosity H = 0.5), it is expected that the internal equilibrium stability should have a stronger effect.
this effect is a consequence of recurrent introduction of diversity via migration: we only observed the maintenance of diversity (either long term or transient) relative to neutrality when the effective rate at which diversity introduced via migration was sufficiently high ( Figure 3). Our work also clarifies the specific role for coevolution in this maintenance of diversity: rather than contributing to stability at the internal equilibrium, coevolution destabilized the boundary equilibria such that perturbations by migrants with rare genotypes grow, pushing the system at least temporarily towards greater genetic diversity.
Comparing outcomes under discrete versus continuous time allowed us to parse the importance of how time is structured for the maintenance versus loss of diversity. Discrete-time models are inherently less stable than their continuous counterparts, and we see this reflected in our results. In our two-patch, discrete-time model, all equilibria were unstable and we failed to observe maintenance of diversity relative to neutrality (note the potential for transiently greater diversity, Figure 4d). In contrast, the internal equilibrium of our two-patch, continuous-time model was stabilized by migration.
This stability, coupled with the introduction of diverse migrants, resulted in parameter space in which more diversity was maintained by coevolution than neutrality (Figure 3).
In this article, we have focused on describing the maintenance of diversity in cases where the host and parasite populations are of equal size and both species have the same generation time. As noted in the Theoretical background section, differences in population size or relative generation times can impact evolutionary potential and the extent of host local (mal)adaptation (Gandon & Michalakis, 2002). As the focus on the comparison of genetic variation relative to drift emphasizes, however, neither the relative population sizes nor relative generation times can be considered in the absence of their effects on the rate of genetic drift. For example, in the absence of mutation, increasing the number of parasite generations per host generation (see top row of Figure 6 in Gandon & Michalakis, 2002) can increase the rate of drift in the parasite and lead to reduced parasite local adaptation. A deeper understanding of the local adaptation in co-evolving systems will hence require consideration of the effect of drift on evolutionary potential.
It remains to be determined the scale of population structure needed to maintain diversity in the long term under coevolution.
Our one-patch and island-mainland models reflect two extremes, with the island-mainland model functioning similarly to considering a single patch in a large meta-population, in which the average allele frequency remains near p = q = 0.5, whereas allele frequencies in individual patches may differ. Future work that considers larger, finite patch networks could help determine the minimal requirements of population structure to see the long-term maintenance of diversity.

CO N FLI C T O F I NTE R E S T S TATE M E NT
The authors declare no conflict of interest.

PE E R R E V I E W
The peer review history for this article is available at https:// www.webof scien ce.com/api/gatew ay/wos/peer-revie w/10.1111/ jeb.14207.

DATA AVA I L A B I L I T Y S TAT E M E N T
The code and data that support the findings of this study are openly available in GitHub repository "host_parasite_coevolution" at https://doi.org/10.5281/zenodo.8096853