Challenging a host–pathogen paradigm: Susceptibility to chytridiomycosis is decoupled from genetic erosion

Abstract The putatively positive association between host genetic diversity and the ability to defend against pathogens has long attracted the attention of evolutionary biologists. Chytridiomycosis, a disease caused by the chytrid fungus Batrachochytrium dendrobatidis (Bd), has emerged in recent decades as a cause of dramatic declines and extinctions across the amphibian clade. Bd susceptibility can vary widely across populations of the same species, but the relationship between standing genetic diversity and susceptibility has remained notably underexplored so far. Here, we focus on a putatively Bd‐naive system of two mainland and two island populations of the common toad (Bufo bufo) at the edge of the species’ range and use controlled infection experiments and dd‐RAD sequencing of >10 000 SNPs across 95 individuals to characterize the role of host population identity, genetic variation and individual body mass in mediating host response to the pathogen. We found strong genetic differentiation between populations and marked variation in their susceptibility to Bd. This variation was not, however, governed by isolation‐mediated genetic erosion, and individual heterozygosity was even found to be negatively correlated with survival. Individual survival during infection experiments was strongly positively related to body mass, which itself was unrelated to population of origin or heterozygosity. Our findings underscore the general importance of context‐dependency when assessing the role of host genetic variation for the ability of defence against pathogens.

in which a potentially lethal pathogen is novel, the neutral genetic diversity of host populations is expected to be positively associated with survival because of the tendency for heterogenous host genotypes to limit pathogen spread (e.g. Ekroth et al., 2019;Gibson & Nguyen, 2020; see also Ostfeld & Keesing, 2012) and the protection conferred by higher genome-wide heterozygosity at the individual level (e.g. Coltman et al., 1999;Pearman & Garner, 2005;Spielman et al., 2004). Due to previous parasite-mediated selection and adaptation, systems with a long history of infection are, however, governed by different effects of host gene pools on disease outbreaks than systems that have not been previously severely affected, even leading to situations where selection can lead to decreased functional genetic diversity paralleled with increased disease resistance (e.g. Leeds et al., 2010;Savage & Zamudio, 2016).
Plants and animals are generally unevenly distributed in space.
The extent to which they exist as partitioned populations depends on their natural histories as well as the distribution of suitable habitats and broadly increases from the core of a species' range towards its edge (e.g. Sexton et al., 2009). A particularly high level of natural fragmentation pertains to populations that are distributed across islands, resulting in high levels of local erosion of genetic diversity combined with high degrees of distinctiveness as a function of population size and island connectivity (e.g. Whittaker & Fernández-Palacios, 2007). Particularly when compared to mainland populations, a given taxon that occupies a series of adjacent islands thus represents an amenable system to study the effects of hierarchical population structure on phenotypic evolution and ecological processes, including the extent to which populations respond to pathogens (Grenfell & Harwood, 1997;Millins et al., 2018;Parratt et al., 2016;Whiteman et al., 2006). Due to low vagility and life history attributes such as water dependency, amphibians are particularly dispersal limited and often reside in genetically distinct, spatially structured populations (e.g. Hancock & Hedrick, 2018). With marked population declines in all biogeographic regions, they also represent a particularly prominent example of the rapid and severe loss of biodiversity currently underway (Houlahan et al., 2000;Stuart et al., 2004;Wake & Vredenburg, 2008). A main driver for these declines is chytridiomycosis, an emerging infectious disease caused by the fungal pathogen Batrachochytrium dendrobatidis (hereafter Bd). Chytridiomycosis is regarded as the most devastating vertebrate disease ever recorded and implicated in the declines or extinctions of hundreds of species (O'Hanlon et al., 2018;Scheele et al., 2019).
Batrachochytrium dendrobatidis is however not universally destructive, with the severity of disease outbreaks depending on Bd strain, amphibian taxonomy and a multitude of environmental factors (for a review see Fisher & Garner, 2020). The impact of infection also varies widely among individuals of the same species, and controlled laboratory experiments have revealed a generally positive association between resistance to Bd and body size (Bielby et al., 2015;Carey et al., 2006;Garner et al., 2009;Meurling et al., 2021;Searle et al., 2011). Marked intraspecific variation in Bd susceptibility further occurs between distinct populations, suggesting that innate genetic factors are acting (Tobler & Schmidt, 2011;Luquet et al., 2012;Bataille et al., 2015;Bradley et al., 2015;see also Palomar et al., 2016). After an outbreak has occurred, field surveys linking Bd prevalence with host neutral genetic variability have yielded mixed findings, with positive (Addis et al., 2015;Horner et al., 2017), negative (Savage et al., 2015) and neutral (Wagner et al., 2017) associations all reported. Given Bd's impact as a pathogen when it invades an area for the first time, the role played by levels of standing genetic variation prior to exposure in predicting the immediate responses of host populations is, therefore, a critical question that has, however, so far remained underexplored.
Despite recent declines, the common toad (Bufo bufo) remains one of the most widespread amphibians in northern Europe, where fragmented population structures on islands and adjacent coastlines have been documented (O'Brien et al., 2021;Roth & Jehle, 2016;Seppä & Laurila, 1999;Tuncay et al., 2018). As a result of its intermediate susceptibility to Bd (e.g. Kärvemo et al., 2019), B. bufo has also served as a model species to unravel key aspects of chytridiomycosis, including the discovery of the highly virulent BdGPL lineage (Farrer et al., 2011;Fisher et al., 2009) and the protective role played by body mass and development to reduce Bd susceptibility (Bielby et al., 2015;Garner et al., 2009Garner et al., , 2011Meurling et al., 2021).
In the present study, we sample B. bufo individuals from mainland and inshore islands at the edge of the species' range and characterize these populations using a high-density panel of de novo assembled SNP markers. Using controlled infection experiments, we investigate whether population-level genetic diversity or individual multi-locus heterozygosity affect the capacity to survive exposure to Bd, taking body mass as a known phenotypic determinant of Bd susceptibility into account.

| Sample collection and animal husbandry
Our study was based on four B. bufo populations in western Scotland ( Figure 1), which were presumed to be Bd-naïve as a national survey revealed that the closest site with a positive record was over 280 km south of the study area (Cunningham & Minting, 2008). Two populations (MAK, MAT) were situated on the mainland, one population (SKB) was situated on the large (1656 km 2 ) island of Skye separated from the mainland by about 300 m, and one population (CRO) was located on the small island of Crowlin (approximately 1 km 2 in area and 1.5 km from the mainland), isolated by sea at least since the last period of glacial activity ended approximately 9500 years ago (Lambeck, 1993). Geographic distances between the populations ranged from approximately 3 km (MAT-CRO) to approximately 12 km (MAT-SKB, Figure 1).
Bufo bufo is an explosive synchronous breeder, and all animals used in the experiment were equivalent in age, being young-of-theyear collected as spawn over a three-day period in March 2016.
Nineteen clutches (egg strings, 3-6 per population) were sampled by collecting approximately 10 cm sections from each of them. Eggs were transported to the Institute of Zoology, Zoological Society of London, where they were hatched in UK Home Office-approved outdoor animal facilities and housed separately in 9 L polypropylene boxes (Really Useful Products LTD) containing dechlorinated, aged tap water and an air-driven sponge filter system with cover and vegetation. Partial (approximately 20%) water changes, as well as feeding using crushed Tetra Tabimin ® tablets, were carried out every two days. After metamorphosis, individuals were housed by clutch in 360 × 210 × 160 mm terrestrial tanks (Faunarium, ExoTerra ® ) lined with moist paper towelling and containing a cover object and water container. All individuals were fed hatchling crickets ad libitum and were subject to daily health checks, and biosecurity and sanitation were maintained throughout animal raising and experimental procedures.

| Bd exposure trials
Experimental procedures took place in a temperature-controlled room at 18°C with a 12:12-h light:dark cycle. Ninety days after the start of metamorphosis, toadlets were weighed to the nearest 0.001 g and transferred to individual housing in 0.7 L polypropylene boxes (Really Useful Products LTD) lined with a moist paper towel and containing a cover object. Following a 2-week period of acclimatization, 202 toadlets were randomly allocated into one of two experimental treatments: (I) exposure to an active Bd culture or (II) a control group exposed to culture media alone. Treatments broadly followed Garner et al. (2009Garner et al. ( , 2011. Each toadlet was placed into an individual Petri dish filled with 30 ml of aged tap water and 450 μl of culture media for four hours. In order to ensure a level of exposure sufficient to cause mortality, nine such treatments were carried out across 21 days. All individuals in the exposed group received the same dosage of the Bd isolate UK CORN'12 3.1, part of the hypervirulent global panzootic lineage BdGPL (Farrer et al., 2011). Dosages were calculated by counting live spores using a haemocytometer and varied between 40 500 and 382 500 zoospores depending on session. The experiment proceeded for 50 days after the first exposure (i.e. a 21-day period of inoculations followed by a further 29 days of daily monitoring of mortality). If a toadlet reached a humane end point (i.e. was unresponsive or was unable to support its own weight), it was euthanised by licenced personnel in accordance with the Animal (Scientific Procedures) Act 1986 using the non-schedule 1 method of immersion in buffered tricaine methanesulphonate (MS222) followed by fixation in 70% ethanol. All animals surviving to the end of the experiment were euthanised and fixed as described above.

| Genotyping and Bd detection
Toe clips from the hind foot of all experimental individuals were sampled for Bd screening after death, using the qPCR assay protocol described in Boyle et al. (2004). DNA was extracted using Prepman™ Ultra (Applied Biosystems), and samples were analysed in duplicate using a TaqMan ® qPCR assay (Applied Biosystems) with four dilution standards representing 100, 10, 1 and 0.1 Bd zoospore genomic equivalents (BdGE), and a negative control. In the case where only one replicate amplified, that sample was reanalysed. This occurred with three samples, which were negative in subsequent analyses.
DNA for genotyping was extracted from hind leg muscle of toadlets using a Qiagen DNEasy extraction kit following the man- Once obtained from Floragenex, a de novo catalogue of SNP loci was constructed using StackS 2.0 (Catchen et al., 2013), given that the B. bufo genome (Streicher & Darwin Tree of Life Consortium, 2021) was not yet available at the time of the analyses. The raw sequences were filtered and demultiplexed using the process_radtags pipeline. Reads with an uncalled base were discarded, as were reads containing a 15 bp window in which the average quality dropped below a phred score of 10 (i.e. a 90% probability of being correct).
Barcodes and RAD-tags containing one mismatch to an expected sequence were retained. To identify an optimal set of parameter values in StackS, we followed procedures described in Paris et al. (2017) and Rochette and Catchen (2017). Based on these considerations,

| Data analyses
Observed (H o ) and expected (H e ) heterozygosities for each population were calculated using the R package adegenet (Jombart, 2008).
Deviations from Hardy-Weinberg equilibrium and pairwise F ST values were calculated using the software genepop 4.4 (Rousset, 2008; input files were generated from vcf files using pgdSpider 2.1.1.5, Lischer & Excoffier, 2012). Individual heterozygosities were calcu- equalling the number of samples used. We used a chain length of 10 000 000 generations sampling every 1000 trees and visualized the resulting tree using denSitree (Bouckaert, 2010) as implemented in beaSt. An absolute dating of nodes derived from Snapp is difficult to achieve without detailed biogeographic information and can be further biased due to gene flow, non-dependence of loci as well as selection (e.g. Barratt et al., 2018;Stange et al., 2018). We, therefore, refrained from a time calibration of the obtained tree.
To evaluate the factors influencing survival in response to Bd exposure, Cox Proportional Hazard (CPH) analyses were performed incorporating the explanatory variables population of origin, individual mean heterozygosity and pre-exposure body mass (cg). Each variable was tested to confirm adherence to model assumptions.
These tests were carried out using the Survival package version 2.38 (Therneau, 2015) in R 3.5.0 (R Core Team, 2018). To further illustrate the impact of mass on survival, individuals were divided into quartiles based on mass, and Kaplan-Meier survival curves were plotted for these groupings as well as for populations.

| Genetic structure of populations
The final filtered SNP data set contained a total of 10 884 biallelic loci. As expected, the mainland populations were characterized by

| Individual-and population-level response to experimental Bd exposure
Within 50 days of the initial exposure, 72 of 101 (71%) toadlets exposed to Bd died, whereas all 101 control animals survived. All experimental mortalities and 35% (10/29) of exposed animals surviving to the end of the experiment were recorded as infected, whereas no Bd was detected in the unexposed control group. Infection prevalence for exposed individuals at point of death was 70% for MAT, 80% for CRO and 86% for MAK and SKB. Survival was 18% (4/22) for SKB, 19% (7/37) for MAK, 27% (4/15) for CRO and 52% (14/27) for MAT (Figure 3a, with median survival times of 28, 28, 42 and >50 days, respectively) and significantly differed between populations (logrank test = 12.6, d.f. = 3, p = 0.007). When compared to the population in which survival was highest (MAT), population membership significantly determined mortality when controlling for the effect of body mass (CPH coefficients ranging from 0.926 to 1.046, see Table 1); the membership of clutch, on the contrary, did not show a significant effect (at however low sample sizes per clutch, data not shown).
Individual mass of exposed individuals did not vary between populations (ANOVA, F(3, 97) = 1.211, p = 0.31), but had a highly significant effect on survival during trials, with an additional 1 cg of mass leading to an approximately 6% decrease in mortality risk per day when controlling for population (CPH coefficient in combined model = −0.063, p < 0.001, Table 1). This was also reflected in the significant difference in survival times between mass quartiles (logrank test = 39.4, d.f. = 3, p < 0.001. Figure 3b). SNP heterozygosity was not associated with mass (either overall or within population, Spearman rank correlation p-values >0.10 in all cases), nor was it found to be associated with survival when evaluated independently in a CPH (see also Figure 4). In the combined model (Table 1), however, SNP multi-locus heterozygosity was found to have a significant negative relationship with survival when accounting for population and mass.

| DISCUSS ION
The present study combined controlled infection trials with highthroughput genotyping to investigate whether host genome-wide genetic variation is linked to disease susceptibility in the chytridamphibian system. Our results confirm that island populations were subject to elevated levels of genetic erosion compared with populations on the mainland and demonstrate that the study populations responded differentially to pathogen exposure. We also replicate the observation that larger body size protects against the negative effects of Bd exposure. However, contrary to the fundamental paradigm that host genetic diversity impedes pathogen success (Coltman et al., 1999;Pearman & Garner, 2005;Spielman et al., 2004), the observed variation in survival was not driven by isolationmediated population genetic erosion and even showed a negative association with individual heterozygosity. Our study underscores the importance of context-dependency when assessing the inter- Our finding of a pronounced spatial genetic structure aligns with other studies on northern European B. bufo populations (e.g. Brede & Beebee, 2004;Thörn et al., 2021), with the high-density panel of SNPs unequivocally demonstrating a reduced level of heterozygosity in island populations compared with the mainland, an effect which was less apparent in previous studies based on fewer loci (Roth & Jehle, 2016;Seppä & Laurila, 1999). The early putative split of the population on the small island (CRO) confirms evidence from microsatellites and mtDNA that supports its natural origin likely through F I G U R E 2 Genetic structure of the four study populations. (a) SNP-based phylogenomic tree as inferred by the Snapp algorithm; node labels denote relative node ages, with bars showing 95% highest posterior density intervals. (b) Scatterplot of the first two principal components of DAPC analysis. Each individual is coloured according to its population of origin rafting events during the melting of glaciers as opposed to a recent human introduction (O'Brien et al., 2021).
The ability to carry the costs incurred by an infection intuitively depends on the phenotype. While a positive relationship between body mass and resistance or resilience to disease is not universally present across all host-parasite systems (for a meta-analysis see Sánchez et al., 2018), our finding of increased survival in larger individuals is consistent with previous Bd exposure experiments conducted with B. bufo (Bielby et al., 2015;Garner et al., 2009Garner et al., , 2011Meurling et al., 2021) and a range of other anuran species (Bradley et al., 2015;Carey et al., 2006;Kriger et al., 2006;Searle et al., 2011;Tobler & Schmidt, 2011). Bd infection is confined to the epidermis, and smaller individuals exhibit higher relative rates of energetically costly skin sloughing and ion loss than their larger counterparts (Wu et al., 2018). Despite an initially smaller potential area of exposure, they also may generally be more vulnerable because both surface area and metabolic rate scale allometrically with mass (Klein et al., 2016;White et al., 2006). Our SNP data revealed that the observed individual differences in mass are not mediated by an effect of heterozygosity on attained size. This finding contrasts a previous study which demonstrated a link between genetic variation and developmental rates in B. bufo populations suffering from urbanization (Hitchings & Beebee, 1998; see also Rowe et al., 1999;Lesbarrères et al., 2007 for similar findings in other anuran species), suggesting F I G U R E 3 Kaplan-Meier survival curves comparing the survival of (a) four Bd challenged populations and controls, and (b) mass quartile groupings of Bufo bufo individuals experimentally exposed to Bd. All control animals survived (grey lines)  (Coles et al., 2019;Reading, 2007), coupled with a projected expansion of the range of Bd in temperate zones of the Northern Hemisphere associated with a warming climate (Xie et al., 2016) suggesting that the role and even the isolated island population CRO. Our unexpected finding of a negative association between survival and heterozygosity could be caused by advantageous additive genetic variants, whose homozygote state is reflected in low genome-wide levels of genetic variation. While a GWAS approach proved unsuited due to samplesize limitations (See Table S1), future studies could also more specif-

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

AUTH O R ' S CO NTR I B UTI O N S
DO'B and JH provided fieldwork expertise; DS, RJ, DO'B and JH conducted the fieldwork; DS, CS and LMB performed the animal husbandry work; DS performed the experimental infections; DS, TWJG, XAH and RJ conceived the study, TWJG, XAH and RJ coordinated the study; DS and RJ analysed the data and wrote the manuscript, with input from all authors.

F I G U R E 4
Relationship between individual heterozygosity and disease trial performance (days survived) of Bufo bufo individuals exposed to Bd. CRO and SKB are island populations, and MAK and MAT are mainland populations

PE E R R E V I E W
The peer review history for this article is available at https://publo ns.com/publo n/10.1111/jeb.13987.

DATA AVA I L A B I L I T Y S TAT E M E N T
The data that support the findings of this study are openly available