Evolutionary rescue in a host–pathogen system results in coexistence not clearance

Abstract The evolutionary rescue of host populations may prevent extinction from novel pathogens. However, the conditions that facilitate rapid evolution of hosts, in particular the population variation in host susceptibility, and the effects of host evolution in response to pathogens on population outcomes remain largely unknown. We constructed an individual‐based model to determine the relationships between genetic variation in host susceptibility and population persistence in an amphibian‐fungal pathogen (Batrachochytrium dendrobatidis) system. We found that host populations can rapidly evolve reduced susceptibility to a novel pathogen and that this rapid evolution led to a 71‐fold increase in the likelihood of host–pathogen coexistence. However, the increased rates of coexistence came at a cost to host populations; fewer populations cleared infection, population sizes were depressed, and neutral genetic diversity was lost. Larger adult host population sizes and greater adaptive genetic variation prior to the onset of pathogen introduction led to substantially reduced rates of extinction, suggesting that populations with these characteristics should be prioritized for conservation when species are threatened by novel infectious diseases.


| INTRODUCTION
Infectious diseases are rapidly becoming a principal threat to biodiversity (Daszak, 2000;Smith, Sax, & Lafferty, 2006). Myriad factors such as global climate change and habitat fragmentation are intensifying disease outbreaks in natural systems, and the outcomes can be severe (Aguirre & Tabor, 2008;Harvell et al., 2002). Treatments for infected populations, if available, are often cost prohibitive, logistically intractable or result in substantial side effects to nontarget species (e.g., Bosch et al., 2015;Pederson & Fenton, 2015). One effective preventative measure could be to maintain genetic diversity within natural populations such that infected host populations can quickly respond to the selective pressures imposed by disease outbreaks. However, the genetic, demographic and environmental conditions that allow host populations to rapidly evolve in response to a pathogen remain largely unknown. Understanding the evolutionary processes that affect disease dynamics on ecologically relevant timescales is essential for identifying which populations are at the greatest risk for diseasedriven declines and for predicting the long-term persistence of host populations.
From a mechanistic perspective, host populations can respond to infectious disease by evolving resistance or tolerance to the pathogen (reviewed in Altizer, Harvell, & Friedle, 2003;Penczykowski, Forde, & Duffy, 2011). Recent empirical work has revealed that evolutionary responses can happen very quickly, sometimes within only a handful of generations (Christie, Marine, Fox, French, & Blouin, 2016;Hendry, 2016;Stockwell, Hendry, & Kinnison, 2003). When an environmental change impacts a population, a rapid response to selection can prevent population extinction in a process known as "evolutionary rescue" (Carlson, Cunningham, & Westley, 2014). Although much of evolutionary rescue research has focused on abiotic environmental changes (e.g., climate change, pollution; reviewed in Gonzalez, Ronce, Ferriere, & Hochberg, 2012;Carlson et al., 2014), biotic perturbations, such as the introduction of a novel pathogen, likely have similar outcomes.
However, the introduction of a pathogen is different from abiotic stress because the effects of pathogens on host populations are often density-dependent (e.g., for many pathogens with direct transmission). Therefore, the selective pressures from a pathogen will change depending on the density, abundance and susceptibility of the host population (e.g., Burdon & Chilvers, 1982;Duffy et al., 2012). Furthermore, evolutionary rescue could result in two possible outcomes in disease systems, host-pathogen coexistence or clearance (pathogen extinction), and it is unclear under what conditions and how frequently each of these two outcomes will occur. Thus, there is a need to understand evolutionary rescue in host-pathogen systems (Bonte, Hovestadt, & Poetke, 2009;Gandon, Hochberg, Holt, & Day, 2012), particularly when populations are confronted with novel or introduced pathogens.
One parameter that is critical for improving our understanding of evolutionary rescue in host-pathogen systems is the amount and type of genetic variation present within host populations. Some disease models predict that increased host genetic variation has little effect on the spread of pathogens within populations (Nath, Woolliams, & Bishop, 2008;Springbett, MacKenzie, Woolliams, & Bishop, 2003;Yates, Antia, & Regoes, 2006), while others show that high genetic variation can directly reduce disease spread (Lively, 2010). However, the majority of experimental studies investigating the effects of genetic variation in host susceptibility have used single generations and have not allowed host evolution to occur (e.g., Zhu et al., 2000;Hughes & Boomsma, 2004; but see Altermatt & Ebert, 2008). Thus, while host genetic variation can play a significant role in short-term disease dynamics (i.e., single host generations), much remains unknown regarding the relationship between variation in host susceptibility and the ability of host populations to evolve over multiple generations. For example, it is unknown how the within-population distribution of genetically based host susceptibilities interacts to influence the likelihood of evolutionary rescue. Additionally, rapid host evolution may fail to rescue populations from going extinct in some scenarios (i.e., there are limits to evolutionary rescue; Bell, 2013;Osmond & Mazancourt, 2013;Stewart et al., 2017). Thus, more studies are necessary to understand the conditions that allow for evolutionary rescue in response to pathogens.
One pathogen that has had large negative effects on host populations is the fungus, Batrachochytrium dendrobatidis (Bd), which infects the keratinized structures of amphibians and can result in high rates of mortality (Garner et al., 2009;Searle, & Gervasi, et al., 2011).
Global amphibian population declines and extinctions have been directly linked to the introduction of this pathogen (Skerratt et al., 2007). However, some amphibian species and populations are seemingly unaffected by Bd, suffering no detectable population declines in its presence (e.g., Daszak et al., 2004Daszak et al., , 2005Reeder, Pessier, & Vredenburg, 2012). Variation in disease dynamics after Bd introduction may be driven by a number of variables including rapid evolution of amphibian hosts in response to Bd. There is evidence that Bd epidemics can impose directional selection for MHC genes incurring resistance (Bataille et al., 2015;Savage & Zamudio, 2011, 2016 and individuals from populations with a history of exposure to Bd can experience reduced pathogen loads compared to individuals from naïve populations (Knapp et al., 2016). One recent study in toads found that estimates of narrow-sense heritability for Bd-load were fairly substantial (h 2 = 0.22; Palomar, Bosch, & Cano, 2016), suggesting that there is heritable variation for susceptibility in natural populations and that these populations may be able to rapidly respond to the strong selection imposed by Bd. Critically, however, it remains unknown whether the evolution of amphibian hosts in response to Bd can lead to evolutionary rescue. Several previous models of amphibian-Bd dynamics have been developed (Briggs, Knapp, & Vredenburg, 2010;Briggs, Vredenburg, Knapp, & Rachowicz, 2005;Converse et al., 2016;Drawert, Griesemer, Petzold, & Briggs, 2017;Grogan et al., 2016;Wilber, Langwig, Kilpatrick, McCallum, & Briggs, 2016), but none focuses on evolution of the amphibian host. The goals of our study were to concentrate on the amphibian-Bd system to (i) identify how variation in host susceptibility (range and distribution of host susceptibilities) and host genetic diversity (number of host genotypes) affects the ability of populations to rapidly evolve in response to a pathogen and (ii) quantify how this host evolution affects population recovery (i.e., evolutionary rescue).

| MATERIALS AND METHODS
We constructed a spatially explicit, forward-time, individual-based model of an amphibian population to determine the relationships between variation in host susceptibility, host evolution and population persistence after the introduction of Bd. We chose an individual-based model because we could vary the susceptibility of each individual within a population (DeAngelis & Grimm, 2014), and manipulate a range of population and disease parameters (e.g., adult carrying capacity, transmission rates). Many of our response variables are at the population level because population-level responses are often what is important from a conservation standpoint. However, selection in the model only occurs at the individual level. Our model followed lifehistory characteristics consistent with a majority of temperate anuran species (Wells, 1977(Wells, , 2007 Table S1).

| Model overview
Our model consisted of a series of sequential steps (Figure 1), where each simulation began with the creation of two spatially discrete habitats: (i) an aquatic breeding habitat where adults congregated to breed and where tadpoles were born, and (ii) a terrestrial overwintering habitat, where individuals were solitary. We began with unrelated adults at a population size equal to their carrying capacity (K adults ). Multilocus diploid genotypes for each individual were created in accordance with Hardy-Weinberg Equilibrium and consisted of 100 neutral bi-allelic loci (i.e., single nucleotide polymorphisms; SNPs) with an initial minor allele frequency of 0.2. We also assigned a Bd-susceptibility value for each individual, referred to as genotype-specific mortality (GSM), which had values ranging from 0 to 1, indicating the probability that a genotype dies each year from Bd if infected. We parameterized the model based on existing literature (Table S1).

| Reproduction
The first step in the model was reproduction, which occurred by randomly pairing all adults and creating a maximum of 100 offspring for each pair. As most offspring did not survive to adulthood, this procedure created high variance in reproductive success among pairs, which is commonly observed in natural populations (Christie, Marine, French, Waples, & Blouin, 2012;Hedrick, 2005). Varying the maximum number of offspring per pair from 10 to 1,000 resulted in little change to the number of offspring from each family that survive to metamorphosis (Fig. S1). After offspring were created, multilocus genotypes were assigned to each individual in accordance with Mendelian inheritance (i.e., the diploid offspring inherited one randomly selected allele from each parent). To assign GSM values for each offspring in a family, we sampled a random deviate from a normal distribution where the mean of the distribution was equal to the mid-parent GSM value and the variance equal to half of the genetic variance of the initial population (e.g., Dunlop, Baskett, Heino, & Dieckmann, 2009;. Variation in phenotype driven by the environment would increase this variance.

| Tadpole mortality
After reproduction, the offspring (tadpoles) moved randomly throughout the pond and density-dependent mortality of tadpoles occurred based on tadpole carrying capacity (K tads ). This mortality was random with respect to GSM, neutral genotype, parentage and location. High mortality at the tadpole stage is consistent with most temperate amphibians which have high fecundity and type III survivorship (e.g., Anderson, Hassinge, & Dalrymple, 1971;Calef, 1973).

| Transmission
The model ran for 49 years in the absence of Bd, and at year 50, we introduced a single Bd-infected adult. Because we employed an individual-based model, keeping track of each zoospore shed by each host (up to 68 zoospores per minute; Reeder et al., 2012) was computationally intractable. Thus, we modelled transmission as a densitydependent process by setting the maximum distance that a zoospore could travel to reach a neighbouring host (in a 100 × 100 m environment). We simulated a season of transmission where transmission F I G U R E 1 Schematic representation of the steps in our model. The population is initiated in the aquatic environment with adult population size at carrying capacity. Adults form pairs randomly and susceptibility to Bd has a heritable component that is passed on to the offspring. Reproduction is followed by density-dependent mortality of tadpoles (based on tadpole carrying capacity) and transmission of Bd. The population then transitions to the terrestrial environment with tadpoles completing metamorphosis (becoming metamorphs: "metas.") and adults leaving the pond. Densitydependent adult mortality occurs (based on adult carrying capacity) followed by Bdassociated mortality (based on individual genotype-specific mortality; GSM), after which the population returns to the aquatic environment to begin a new year. A single infected individual is introduced in year 50; in years 1-49 both the transmission and Bd mortality steps are omitted

| Adult mortality
After the transmission step, tadpoles became metamorphs (newly metamorphosed individuals), both adults and metamorphs moved to the terrestrial environment, and density-dependent mortality of the adults and metamorphs occurred based on adult carrying capacity (K adults ). The number of adults was recorded at the beginning of each year in the model (immediately before reproduction) and mortality on metamorphs and adults occurred at the end of the year (after adult emigration) to keep the combined adult and metamorph population near K adults in the absence of any Bd-related mortality. Variation in density-dependent mortality was introduced each year by sampling a random deviate from a normal distribution (μ = N−N t+1 ; SD = 7). All tadpoles completed metamorphosis in one season, which is common for the majority of temperate anurans (Wells, 2007; data from Bancroft et al., 2011).

| Bd mortality
Once an individual became infected, its probability of dying from in-

| Population outcomes
For each simulation, the population was monitored until it either went extinct or 100 years elapsed after Bd introduction (150 years total) at which point the population outcome was recorded. Populations were classified as "extinct" if there were no remaining amphibians, "cleared" if only uninfected amphibians remained (i.e., the population cleared infection), or "coexisting" if infected individuals were present in the population (i.e., both the host and pathogen were present). All model simulations and analyses were run in R version 3.2 (R Core Team 2016).

| Fixed within-population susceptibility
We first modelled scenarios where all adults and offspring within a population had the same GSM, such that host evolution could not occur. This lack of variation is likely present in some populations (e.g., those with genetic constraints or no additive genetic variation) and allowed us to systematically examine the interactions between GSM, population demography and transmission on population outcome. We varied GSM from 0 to 1 (by increments of 0.1), K adults as 100, 200 or 300 individuals, and transmission from low to high (0.2, 0.3 and 0.4 m).
K tads was set to 1,000 individuals. For each unique combination of parameters (Table S2), we performed 100 replicate model simulations.

| Mixed within-population susceptibility
We  Table S2). We performed 100 replicates for each combination of parameters.

| Exploring parameter space
We next modelled 7,098 combinations of parameter values to determine the effect of each variable in contributing to population outcome. We varied transmission from low to very high (13 levels),  (Table S2) was replicated 80 times.
To illustrate the main effects, we plotted each variable against the population outcome without holding the other variables constant.
This procedure, which can be considered a type of sensitivity analysis, reflects a summary of each variable across a large exploration of the total parameter space. We did not perform any significance tests on our data (e.g., calculate p-values) as such analyses can be biased by the large number of replicates (White, Rassweiler, Samhouri, Stier, & White, 2014).

| Fixed within-population susceptibility
The highest proportion of populations that cleared infection occurred at high levels of GSM, coexistence was highest at low levels of GSM, and extinction was greatest at intermediate levels of GSM ( Figure 2).
Additionally, there were substantially more population extinctions as Bd transmission increased; higher transmission increased both the frequency of population extinction at a given GSM value and the range of GSM values at which extinctions occurred (cf. Figure 2a-c). Similar relationships between GSM and population outcome were found when K adults was set to 100 or 300 (Fig. S3) or when we increased the number of infected individuals initially introduced into the population (Fig. S4).
These population outcomes were observed because at very low levels of GSM, there was little to no Bd-related mortality and adult population sizes remained close to K adults (Figure 3, Fig. S5, S6), leading to high rates of coexistence. At intermediate GSM levels, enough individuals survived each year to maintain the pathogen in the population, but mortality was high enough that adult population sizes were greatly reduced, increasing the chances of extinction. When GSM was near one, many infected hosts died before transmitting Bd, leading to the population clearing infection so quickly that at most, a small, temporary reduction in adult population size was observed (e.g., Figure 3f).
In all scenarios, populations that cleared Bd eventually returned to a population abundance near K adults (Figure 3).
The presence of Bd also caused the greatest reduction in neutral genetic diversity at intermediate levels of GSM (Figure 4, Fig. S7).

| Mixed within-population susceptibility
Host populations with wider range in GSM values experienced a greater evolutionary response, leading to a greater proportion of host populations coexisting with Bd ( Figure 5). The uniform distribution for GSM resulted in the greatest evolutionary response ( Figure 5a,b), followed by the log-normal distribution (Fig 5g,h) and then the normal distribution (Figure 5d

| Exploring parameter space
coexisted with Bd, we found that higher initial genetic diversity led to a smaller reduction in adult abundance (Fig. S8).
With greater transmission, the proportion of populations that cleared Bd decreased rapidly and the proportion of extinctions increased in a sigmoidal fashion (Figure 6b). Rates of coexistence peaked at a transmission value of 1.0 m (i.e., when a Bd zoospore could travel up to 1.0 metres to infect a new host). At very low levels of transmission, the pathogen was unable to persist. As transmission increased slightly (e.g., from 0.25 to 1.0 m), Bd was able to persist in populations, but rates of infection were low enough that extinction was uncommon, allowing coexistence to occur. However, as transmission continued to increase (over 1.0 m), many individuals became infected, increasing rates of extinction, reducing rates of clearance and slightly reducing rates of coexistence.
K adults had a much larger effect on population outcome than K tads .
Doubling K adults led to a 14.1% increase in the number of populations coexisting with Bd, and a 2.6% and 9.7% decrease in the number experiencing extinction and clearance, respectively (Figure 6c). K tads had relatively little effect on population outcome ( Figure 6d); for each additional tadpole, there was a 5.1 × 10 −6 and 4.1 × 10 −6 increase in coexistence and clearance, respectively, and a 9.2 × 10 −6 decrease in the proportion of extinctions.

| DISCUSSION
A rapid evolutionary response to selection is one mechanism by which host populations can persist after the introduction of a novel pathogen. Our analyses demonstrate that increased variation in amphibian host susceptibility allowed populations to respond faster and to a greater extent to Bd-induced mortality (Figures 5 and 6a). This evolutionary response greatly decreased rates of host extinction and Genetic diversity (heterozygosity expected) Response of host populations to Bd depending upon their within-population composition of host susceptibility (GSM; genotypespecific mortality). We explored three possible distributions labelled on the x-axis with a "U," "N" or "L," to indicate a uniform (a-c), normal (d-f) or log-normal (g-i) distribution, respectively. The first column (a,d,g) represents the average distribution of GSM within populations prior to the introduction of Bd (year 49; inset depicts a representative histogram of one population). For each distribution, six different within-population variances in GSM were tested. The second column (b,e,h) illustrates the average distribution of GSM 100 years after the introduction of Bd or immediately before host population extinction. The third column (c,f,i) illustrates the population outcomes for each population distribution. Host populations with a wider range of susceptibility values can respond to selection and have a much greater rate of coexistence than populations that do not have sufficient adaptive genetic variation to respond to selection Genotype-specific mortality Genotype-specific mortality G enotype-specific mortality Genotype-specific mortality GSM distributions within an amphibian population will directly influence the chances of evolutionary rescue after the introduction of a pathogen ( Figure 5).
Given the importance of host susceptibility (GSM) in our study, it is imperative that we gather more empirical estimates for this parameter.
Susceptibility is a phenotypic trait which is determined by both genetic and environmental contributions. Thus, the first step for comparing GSM across populations is to develop a standardized assay to measure this phenotype, although disentangling the genetic and environmental contributions to GSM is challenging in systems without genetic clones. While quantitative genetic breeding designs may circumvent some of these challenges, identifying panels of genetic markers that correlate with susceptibility may be a more practical approach (Bataille et al., 2015;Hirschhorn & Daly, 2005;Savage & Zamudio, 2016). With enough individuals, this procedure (or conceptually similar alternatives) will allow for the characterization of the within-population variation in host susceptibility. Regardless of the chosen approach, there remains a need to empirically quantify the variation in susceptibility within and among amphibian populations and species in order to predict which populations will be able to respond to Bd.
When we prevented host evolution from occurring in our model (i.e., when all individuals in a population had the same GSM value), extinction was highest at intermediate levels of GSM (Figure 2). In these scenarios with fixed GSM values throughout the population, GSM behaves similarly to pathogen virulence, which is also predicted to have the largest effects on populations at intermediate levels (Anderson, 1979(Anderson, , 1982 (Bell & Gonzalez, 2011). Thus, if there is no (or limited) genetic variation for susceptibility within a host population (e.g., in small host populations), then populations with very high or low levels of susceptibility will be the least likely to go extinct.
Because we focused on evolution of the host, none of our model simulations allowed for evolution of the pathogen. While there is some evidence that pathogen evolution could occur in this system (Piovia-Scott et al., 2015;Refsnider, Poorten, Langhammer, Burrowes, & Rosenblum, 2015), most amphibian communities contain single or very few strains of Bd, such that evolution of the pathogen is expected to be limited (e.g., Jenkinson et al., 2016;Morgan et al., 2007). However, if Bd can evolve reduced virulence in a declining host population, we would expect rates of coexistence to increase.
Clearing infection did not require host evolution and, in fact, rapid host evolution promoted coexistence and not clearance (Figures 5 and   6a). Although pathogen clearance from a host population occurred in many of our simulations, clearance of Bd has only been documented in the field after substantial intervention (Bosch et al., 2015).
which requires long-term monitoring for successful detection. For example, some populations in our study did not clear infection until over 50 years after Bd introduction (Figure 3c,d). Thus, our results highlight the need for long-term studies to identify reservoirs of Bd and to understand why Bd clearance has been so rarely documented in natural populations.
Many of the simulations that led to populations clearing Bd also led to reduced neutral genetic diversity. Reduced neutral genetic diversity can affect long-term population persistence by reducing the ability of populations to respond to future perturbations (e.g., introduction of a different pathogen strain or species). Because neutral genetic diversity and analogous measures (e.g., effective population size) are often used as a proxy for future adaptive potential (Waples, 1990;Reed & Frankham, 2003; but see Mittell, Nakagawa, & Hadfield, 2015), these measures may be useful as a coarse overview of the number of adaptive genotypes present within a population. Additionally, if the introduction of Bd is documented within a population, a rapid reduction in genetic diversity may suggest high rates of Bd-associated mortality.
In addition to variation in host susceptibility, transmission was an important factor in our model simulations, in part because our range of transmission values was not bounded above zero. In contrast, many of our other parameters were constrained by definition (e.g., GSM can only have values from 0-1) or design (e.g., we chose K adults values that represent natural population sizes: Table S1). In natural systems, variation in transmission could occur due to environmental conditions or Bd strain. Our results show that if Bd transmission is low, the chances of host extinction are also low, regardless of the mean and distribution of susceptibilities (Figures 2 and 6b). Additionally, changing the transmission function of our model to include environmental transmission could lead to different results for host evolution and population outcome. In general, allowing for survival of a pathogen in the environment is likely to increase rates of extinction and reduce the chances of pathogen clearance, because transmission can remain high even if host density becomes extremely low. In the Bd system, a previous study using a mathematical model found that the longer Bd can survive in water, the greater the chances of host extinction (Mitchell, Churcher, Garner, & Fisher, 2008 Managing amphibian populations to have large adult population sizes and high genetic diversity may also help reduce the chances of population extinction after the introduction of Bd. These are the same conditions that are predicted to increase the chances of evolutionary rescue in previous studies (reviewed in Carlson et al., 2014). In our model, rates of extinction were lowest when populations consisted of a wider range of susceptibility values ( Figures 5   and 6a), although these populations coexisted with Bd and therefore maintained population sizes below carrying capacity (Figure 3). We also found that increased genetic diversity benefited populations even when the population outcome was the same; populations with greater numbers of host genotypes were able to maintain larger population sizes when coexisting with Bd (Fig. S8). Thus, we recommend implementing policies to maintain genetic diversity of threatened populations to reduce their extinction risk from Bd. However, managing populations for potential invasion by a novel pathogen can result in conflicting trade-offs; small populations can be less likely to be colonized by an invasive pathogen than large populations due to the smaller number of potential hosts and reduced capacity for transmission (Burdon, Ericson, & Muller, 1995;Stapp, Antolin, & Ball, 2004). Additionally, higher gene flow into a population will generally increase genetic diversity, but also increased risk of pathogen introduction. Thus, management policies must consider trade-offs between the risk of Bd introduction versus the risk of extinction once Bd has already arrived in a population. In the ab- However, it is important to note that in our model, hosts varied in their response to infection (i.e., pathogen-induced mortality; GSM), but not in their likelihood of becoming infected. In systems where host defences primarily involve reducing the chances that an individual becomes infected, evolutionary rescue may lead to an increase in clearance and not coexistence. Thus, we predict that many of the qualitative results we found in this study would also be found in other systems, but the commonality of our findings to systems that lack variation in pathogen-induced mortality or that have different modes of pathogen transmission (e.g., a pathogen that can survive long periods of time in the environment) requires further validation.
Because individual-based models are able to incorporate a multitude of different host, pathogen and environmental scenarios, they continue to be a powerful tool for the study and management of infectious disease in systems with a wide range of characteristics.

| Conclusion
In summary, our analyses demonstrate that possible withinpopulation variation in the susceptibility of amphibian hosts to Bd could play a large role in evolutionary rescue. Rapid evolutionary responses have been increasingly documented in diverse taxa such as plants (Franks, Sim, & Weis, 2007), fishes (Christie et al., 2016;Lescak et al., 2015), lizards (Campbell-Staton et al., 2017) and birds (Grant & Grant, 2006), suggesting that many host species have the potential to rapidly respond to pathogen-mediated selection. However, here we show that in host-pathogen systems, evolutionary rescue can lead to increased rates of coexistence and decreased rates of extinction. Such evolutionary rescue, however, is not without trade-offs; decreased population sizes and lower rates of clearing infection are the costs associated with this process. Nevertheless, for a virulent pathogen that has extirpated large numbers of populations, the drawbacks associated with coexistence greatly outweigh increased risk of extinction. Thus, the successful conservation and management of amphibians in the wake of a Bd epidemic should focus on maintaining genetic diversity, bolstering adult population sizes, and minimizing disease transmission to prevent future extinctions and extirpations of threatened and vulnerable amphibians.