Susceptibility of amphibians to chytridiomycosis is associated with MHC class II conformation

The pathogenic chytrid fungus Batrachochytrium dendrobatidis (Bd) can cause precipitous population declines in its amphibian hosts. Responses of individuals to infection vary greatly with the capacity of their immune system to respond to the pathogen. We used a combination of comparative and experimental approaches to identify major histocompatibility complex class II (MHC-II) alleles encoding molecules that foster the survival of Bd-infected amphibians. We found that Bd-resistant amphibians across four continents share common amino acids in three binding pockets of the MHC-II antigen-binding groove. Moreover, strong signals of selection acting on these specific sites were evident among all species co-existing with the pathogen. In the laboratory, we experimentally inoculated Australian tree frogs with Bd to test how each binding pocket conformation influences disease resistance. Only the conformation of MHC-II pocket 9 of surviving subjects matched those of Bd-resistant species. This MHC-II conformation thus may determine amphibian resistance to Bd, although other MHC-II binding pockets also may contribute to resistance. Rescuing amphibian biodiversity will depend on our understanding of amphibian immune defence mechanisms against Bd. The identification of adaptive genetic markers for Bd resistance represents an important step forward towards that goal.


Introduction
The emerging infectious skin disease chytridiomycosis, caused by the chytrid fungus Batrachochytrium dendrobatidis (denoted Bd) [1,2], is a primary driver of global amphibian population declines and species extinctions [3,4]. Bd colonizes the keratinized epithelial cells of its host's skin [2]. In susceptible amphibians, acute infections can disrupt membrane potentials, leading to paralysis and heart failure [5]. Bd has been detected in more than one-third of all amphibian species surveyed around the world [3,4]. Disease progression, and whether infection culminates in morbidity or mortality, depends on many factors, especially the capacity of the host immune system to respond to the infection [6]. Thus, even in the midst of epizootics, variation exists among individuals in their vulnerability to the disease. In amphibian communities, some species thrive even as others around them perish. The major histocompatibility complex (MHC) encodes receptors at the cell surface that induce and regulate acquired immune responses against pathogens in all jawed vertebrates [7]. Resistance of amphibians to bacterial and viral diseases that cause population die-offs has been demonstrated to be conferred by particular MHC alleles [8,9]. A recent study identified possible MHC allele-specific disease resistance to chytridiomycosis [10]. These findings suggest that features of MHC molecules may improve their capacity to bind Bd antigens and induce adaptive immune responses. If the results have some generality, alleles that encode these molecules should be strongly selected in infected populations.
Many variable amino acid residue positions associated with disease susceptibility in vertebrates are situated in exon 2 of the MHC class II (MHC-II) B gene [11,12], which encodes the b1 segment of the antigen-binding groove that presents epitopes to T cells [7]. The stability of the antigen-MHC complex depends mainly on deep pockets within the binding groove that interact directly with antigen residues [12]. Amino acid changes in this segment result in important structural modifications of the binding groove that can change the affinity of the MHC molecule for particular antigens [11,12]. The MHC-II b1 domain of several amphibian species has been sequenced [10,[13][14][15][16] and residue positions important for antigen binding may be under positive selection [10,15]. Analyses of Bd resistance among amphibians as a function of their MHC-II b1 domain offer a means to examine the evolution of resistance to Bd and to identify specific MHC conformations that confer Bd resistance.
Here, we use a combination of comparative and experimental approaches to demonstrate selection for MHC conformations that foster the survival of Bd-infected amphibians across four continents. We show that resistant amphibians share amino acids within the MHC-II binding groove. We assess how specific sites in this region affect the conformation and binding affinities of MHC molecules that are important in conferring Bd resistance. We then test our predictions by conducting laboratory challenge experiments on a threatened Australian tree frog species, comparing survivorship of MHC genotypes in populations that have evolved with Bd to those that remain naive to the pathogen.

Material and methods (a) Sample collection
In South Korea, we obtained 30 toe-clip samples from wild-caught Asiatic toads (Bufo gargarizans) in Geumsan (3688.249 0 N, 127822.876 0 E) and Jeonju (35847.241 0 N, 12788.348 0 E), and 10 samples from wild-caught oriental fire-bellied toads (Bombina orientalis) in Chuncheon (37858.664 0 N, 127836.146 0 E) and Chiaksan (37823.676 0 N, 128803.221 0 E) (electronic supplementary material, table S1). These two species historically have been infected by endemic Bd strains in Korea and show no evidence of morbidity nor mortality attributable to chytridiomycosis despite extensive fieldwork having been undertaken [17]. In Australia, we obtained toe-clip samples from 30 wild alpine tree frogs (Litoria verreauxii alpina) in each of three similar but geographically separate populations in Kosciuszko National Park, New South Wales. Two of the three populations were long-exposed to the pathogen (site A, Kiandra, 35852.335 0 S, 148829.994 0 E; site B, Ogilvies Creek, 3682.175 0 S, 148819.327 0 E) and a third never had been exposed to Bd (site C, Grey Mare Range, 36819.010 0 S, 148815.567 0 E).
(b) MHC-II b1 genotyping and structure comparison DNA was extracted from toe clips using a salting-out extraction method with ammonium acetate [18]. We used existing primers for B. gargarizans and newly designed primers (see detailed methods in electronic supplementary material) for L. v. alpina and B. orientalis to genotype samples at the b1 domain of one MHC-II locus. Full details of primer sequences and PCR protocols are given in electronic supplementary material, table S2. MHC-II b1 amplicons were purified and cloned using the RBC A&T cloning kit and accompanying HIT-DH5 a competent cells (RC001 and RH617, RBC Bioscience, Taipei, Taiwan) following the manufacturer's protocol. Between 8 and 20 clones were amplified by PCR (electronic supplementary material, table S2) and sequenced using M13 primers. Allelic identity was confirmed when the corresponding DNA sequence was found in at least two clones. Only one or two alleles were recovered from any of our samples, confirming that our primers targeted one single MHC-II locus in all species.
Sequences obtained were aligned with MHC-II b1 sequences of 17 amphibian species from around the world that vary in susceptibility to Bd [19][20][21][22][23][24][25] (length of sequences: 37-81 amino acids; complete list in electronic supplementary material, figure S1; GenBank accession numbers KJ679288-KJ679331) using CLUSTAL W [26]. Sequences then were translated into amino acids using BIOEDIT [27]. We examined the amino acid composition at 15 codon positions known to affect the properties of the P4 (seven sites), P6 (four sites) and P9 (four sites) pockets of the MHC-II peptide-binding groove in humans [11,12]. Codon b9 of pocket P9 was not evaluated because it was missing from more than half of our dataset. The most frequent amino acid across MHC-II b1 sequences of Bd-resistant amphibians was recorded for each of the 15 positions.
We hypothesized that these most frequent amino acid compositions represent pocket conformations associated with Bd resistance. Allelic frequencies within populations could not be considered because such data are not available for most amphibian species. Nevertheless, if Bd induced strong selection for particular MHC-II conformations, we would expect to find a trend. Failure to find a relationship, however, would not allow us to draw any conclusion. For each MHC-II binding pocket, the proportion of b1 sequences containing all the amino acids most frequent in resistant amphibians was compared with that characteristic of susceptible species. Statistical significance was assessed by two-sample tests of equality of proportions, correcting for multiple tests with the Benjamini and Hochberg's procedure [28], using the statistical platform R v. 3.0.2 [29].

(c) Microsatellite genotyping and data analysis
A panel of nine microsatellite markers, including two markers isolated from the Litoria ewingii complex [30] and seven from a L. v. alpina partial genomic DNA library (see the electronic supplementary material), was organized into two fluorescently labelled multiplexes, and multiplex PCR was performed with Qiagen multiplex PCR master mix following the method described in electronic supplementary material, table S3 [31]. Details of microsatellite primers, GenBank accession numbers and PCR protocols are given in electronic supplementary material, table S3. Allele sizes were assigned using an ABI3730 DNA analyser at NICEM (Seoul National University, South Korea) and PEAK SCANNER v. 1.0 (Applied Biosystems, Carlsbad, CA, USA). Conformity to Hardy-Weinberg equilibrium and linkage disequilibrium was determined with GENEPOP v. 4.0 [32]. Heterozygosity values and frequency of null alleles were estimated with CERVUS [33].
Frequency of null alleles was high (F . 0.10) for four markers (Livea-GT1, Livea-AG3, Le2 and Le4) in the three populations genotyped. Therefore, we statistically adjusted allele frequencies assuming a single new allele size by generating a new genotyping dataset corrected for null alleles in FREENA [34]. Subsequent rspb.royalsocietypublishing.org Proc. R. Soc. B 282: 20143127 analyses were run both with the five-marker dataset and the full nine-marker dataset corrected for null alleles.

(d) Experimental infection
Fifteen clutches of L. v. alpina were collected from the three populations described above. Frogs were reared in Bd-free quarantine conditions to adulthood. Frogs were confirmed to be Bd-negative by qPCR [35] prior to the commencement of the exposure experiment. Two hundred adults were randomly chosen from each clutch and population to be assigned to treatments that were inoculated with 750 000 infective Bd zoospores (strain AbercrombieNP-L.booroolongensis-09-LB-P7). The experimental design, involving the utilization of frogs from each population and clutch together with details of the blind randomized block design used for allocation of treatment groups, is outlined in electronic supplementary material, table S4. As a control, another 56 frogs were sham-infected with culture medium in dilute salt solution. Numbers of frogs available for treatments were subject to actual clutch sizes and natural attrition during growth and development.
Frogs were maintained in a quarantine room at temperatures between 18 and 208C under a 12 L : 12 D regimen. Subjects were housed individually in small plastic tubs on an angled rack with drainage holes, a loose pebble floor and a gauze-covered top. They were fed twice weekly with vitamin-dusted (calcium and herptivite alternately) crickets, and tubs were cleaned daily by flushing with fresh filtered water.
Subjects were monitored daily for clinical signs of infection (dullness, lethargy, peripheral erythema and increased skin shedding). Bd infection intensity data were collected weekly by swabbing all subjects, each with fresh gloves and a new sterile dry swab (MW100; Medical Wire and Equipment, Corsham, UK). Infection intensity was estimated as Bd zoospore equivalents (ZSE) in each swab sample using a TaqMan real-time qPCR assay [35] and standards with known Bd quantity. Individual swabs were analysed in triplicate with an internal positive control. Once the last surviving subjects effectively cleared themselves of Bd infection (0-10 ZSE) over four successive weeks, the experiment was terminated.
Frogs showing marked clinical signs of chytridiomycosis were euthanized throughout the experiment using tricaine methanesulfonate (MS-222). Subjects to be sacrificed first were swabbed to quantify Bd infection intensity by qPCR. Foot or toe-clip samples for MHC analysis were collected post-mortem soon after sacrificing subjects and were placed into 90% ethanol. We genotyped the MHC-II b1 domain in 100 L. v. alpina that we had infected with Bd. This selection included the only six individuals that survived the experiment and 94 individuals randomly chosen and encompassing all the clutches collected from the three populations (8-10 individuals per clutch; electronic supplementary material, table S4).

(e) Survival statistical analysis
We used Cox proportional hazard models [36] to identify variables that affected survival of L. v. alpina during the experiment. Parameters included in the initial model were site of origin, clutch, maximum Bd infection load, mass at the start of the experiment, change in mass over the course of the experiment, MHC heterozygosity, allelic divergence (measured as p-distance) and pocket residue composition (table 1). We also examined two-way interaction terms between maximum infection load and the other variables, and between masses and the other variables sequentially. We used likelihood ratio tests to calculate the predictive power of each variable. Independent variables that did not show a significant association with length of survival were excluded from the model. The significance of each variable's effect on hazard risk was evaluated with Wald z statistics. Analyses were conducted with the 'survival' package in R v. 3.0.2.  [37] implemented in LOSITAN [38]. Markers with F ST values significantly higher than the simulated distribution (posterior probability more than 0.95) were considered to be under positive selection. Analyses were done on the full dataset including the three populations, and on pairs of populations, to identify population-specific differences in selection pressure. Balancing or directional selection affecting MHC-II in wild L. v. alpina populations was also assessed using two tests based on deviation of allele frequencies from neutral expectation (Watterson's homozygosity test [39] and Slatkin's exact test [40]) conducted in ARLEQUIN v. 3.5 [41]. Selection pressures on the MHC-II b1 domain were assessed by estimating the ratio of non-synonymous to synonymous nucleotide substitutions using multiple methods available on the Datamonkey website [42]. First, we confirmed the absence of recombination in our sequences using single breakpoint recombination (SBR) and genetic algorithms for recombination detection (GARD) [43]. To estimate the ratio of non-synonymous to synonymous nucleotide substitutions, we used the partitioning approach for robust inference of selection (PARRIS) [44] (v ¼ dN/dS; positive selection: v . 1). Analyses were conducted across MHC-II b1 sequence alignments to identify the selection model best fitting our data. Positively selected sites were detected using single-likelihood ancestor counting (SLAC), fixed effects likelihood (FEL), random-effects likelihood (REL) analysis [45] and mixed-effects model of evolution (MEME) [46]. We also tested whether specific biochemical properties were driving substitutions at the peptide-binding sites using the property-informed model of evolution (PRIME) method [42], which includes the Holm-Bonferroni procedure to control for familywise false positives within sites.

(a) Comparison of MHC-II b1 in worldwide amphibians
The MHC-II b1 alleles of amphibian species least affected by Bd infection [21][22][23][24][25] consistently presented the same amino acids, or amino acids with similar chemical properties, at all 15 pocket residues (figure 1 and table 1; electronic supplementary material, figure S1). By contrast, these pocket-specific amino acid compositions were less frequent in species susceptible to Bd [19,20]   Electrically charged (e.g. negative Asp/Glu at 28b or 66b; positive Arg/Lys at 71b) or hydrophobic (e.g. Tyr/Phe at 26b, 30b or 37b) amino acids were found for 12 of the 15 codon positions for amphibians with low susceptibility to Bd (table 1). Amino acids with these properties can modulate anchor residue specificity by changing the charge characteristics of the pockets [12]. Other common amino acids, Valb11 and Serb13, increase the space available within the pockets for large anchor residues owing to their very small side chain [12]. Codon b56 encodes proline, the only such amino acid within the peptide-binding region of this MHC-II b locus (table 1; electronic supplementary material, figure S1). Proline is an amino acid with a unique conformational rigidity strongly affecting protein secondary structure such as alpha helices [48].

(c) Experimental infection of Litoria verreauxii alpina
After inoculating frogs with Bd, subjects from long-exposed site A survived significantly longer than those from the other long-exposed site B and naive site C (Kaplan-Meier p ¼ 0.001; figure 2; electronic supplementary material, table S6). Many subjects were heavily infected and died within five weeks (figure 2). Infection loads of some individuals from sites A and C stabilized after five weeks but increased afterwards, leading to morbidity and mortality after eight weeks. Six of 200 frogs (five from site A, one from site C) demonstrated greatly reduced loads after three weeks and survived until the end of the experiment (figure 2). All non-exposed, control frogs survived until the end of the experiment.
We genotyped the b1 domain of one MHC-II locus in 84 individuals, including the six survivors, that we had inoculated with Bd. Twenty-two MHC-II b1 alleles were recovered, with eight alleles identified among the surviving individuals (Livea-1, 2, 3b, 5a, 5b, 11, 13 and 14; figure 3; electronic supplementary material, table S6 and figure S1). These eight alleles had identical residues at five codon positions associated with the P9 pocket, corresponding to the composition identified in resistant amphibians worldwide (table 1 and figure 1; electronic supplementary material, S1). Subjects that died had significantly lower frequencies of alleles with the specific P9 pocket composition than those that survived (x 2 1 ¼ 7:18, p ¼ 0.004; table 1 and figure 3). None of the eight alleles had the P4 pocket composition identified in the worldwide MHC-II b1 alignment, and only two alleles had the targeted P6 pocket composition (table 1). We identified another 15 alleles in individuals that did not survive infection. Six of these alleles presented the P9 pocket composition identified in Bd-resistant amphibians (figures 1 and 3; electronic supplementary material, S1).
We used Cox proportional hazard models to determine whether P6 and P9 pocket residue composition, MHC heterozygosity and other variables affected the survivorship of subjects during the experiment. The model that best fitted the data included P6 and P9 pocket residue compositions, the clutch and site of origin, and the interaction between P9 composition and the maximum Bd infection load (  The best-fitting selection model allowed for site-specific positive selection (v ¼ 3.14, LRT ¼ 6.58, p ¼ 0.04). Along with other peptide-binding residues, residues b37 and b57 of the P9 pocket were identified as being under positive selection using REL analysis ( posterior probability more than 0.99 for both residues). Residue b57 also was identified as being under positive selection using MEME ( p ¼ 0.002; electronic supplementary material, table S5). If MHC alleles with the specific P9 pocket conformation are advantageous to L. v. alpina against Bd infection, we should also observe signs of directional selection pressure on this locus in wild L. v. alpina exposed to this pathogen. We sampled 30 individuals from each of the three sites used as source populations for the infection experiment. We genotyped these individuals for the MHC-II b1 domain and nine microsatellite markers. Contrary to our predictions, the frequency of P9-specific alleles in exposed site A (37.5%) was lower than that of other alleles (x 2 1 ¼ 5:04, p ¼ 0.50), whereas P9-specific alleles were more frequent in exposed site B (93.8%, x 2 1 ¼ 70:04, p , 0.0001) and in naive site C (97.5%, x 2 1 ¼ 84:38, p , 0.0001). Levels of MHC heterozygosity were high in all populations (60-65%). Genetic variation among populations (measured as F ST ) was significantly higher at the MHC-II locus than at the microsatellite loci among the populations (F ST ¼ 0.147, posterior probability more than 0.99; electronic supplementary material, table S7), suggesting that directional selection is affecting the MHC class II locus across populations [49]. However, Watterson's F-test and Slatkin's exact p-test did not detect any significant deviation from neutrality in MHC allele frequencies of the three populations

Discussion
We have shown that chytridiomycosis appears to select for specific properties of the MHC class II molecule P9-binding pocket. This is apparent both in our experimental infection study of L. v. alpina and our surveys of wild populations of Bd-resistant Bufo and Bombina species in Korea. While our comparative study also points to specific P4 and P6 pocket conformations that correlate with Bd resistance worldwide, these conformations did not improve survival of infected L. v. alpina, nor are they abundant in Korean Bufo and Bombina populations. Possibly, the importance of the various antigen-binding pockets in triggering an efficient immune response varies with hosts, infecting Bd strains and environmental factors.
The protective P9 pocket conformation includes an aromatic b37, an acidic Aspb57, a Prob56 and a hydrophobic b60 residue. We could not fully assess the importance of b9, another P9 peptide-binding residue, because the position is missing from sequences of many of the species that we examined. Residue Prob56, although not directly binding to anchor residues of epitopes, may be of especial importance for efficient Bd antigen binding by influencing how the peptide-binding residue b57 is positioned within the P9 pocket [48].
Our experimental study on L. v. alpina demonstrates adaptive immune-system-mediated recovery during infection solely in individuals bearing two MHC alleles with the advantageous P9 conformation. MHC alleles are co-dominantly expressed [50], so two doses of MHC molecules binding preferentially to Bd peptides may be required to offset the pathogen's mechanism of attack on the adaptive immune system [51]. Alternatively, MHC alleles with the protective P9 conformation may interact more efficiently with the specific repertoire of T-cell receptors or T-cell subsets that respond to Bd antigens [52]. Many individuals with the advantageous P9 pocket conformation nonetheless died. The protective effect of MHC class II b1 molecules may be modulated by other factors regulating immune system function. The high virulence of the Bd strain to the frogs, indicated by massive infection loads borne by some subjects (figure 2; electronic supplementary material, table S6), may have resulted in effective inhibition of adaptive immunity [51].
Results from wild L. v. alpina populations show only weak evidence of selection for specific alleles of the MHC class II locus in infected populations. Inhibition of immune response may prevent selection of advantageous MHC alleles in the wild. However, the differential survival of subjects from the two historically infected sites suggests that selection for resistance to Bd varies among sites. Genetic variation and environmental factors may influence how strongly resistance alleles are selected in infected populations. L iv e a -1 L iv e a -2 L iv e a -3 a L iv e a -3 b L iv e a -4 L iv e a -5 a L iv e a -5 b L iv e a -5 c L iv e a -6 L iv e a -7 a L iv e a -8 a L iv e a -8 b L iv e a -1 0 L iv e a -1 1 L iv e a -1 2 L iv e a -1 3 L iv e a -1 4 L iv e a -1 5 L iv e a -1 6 L iv e a -1 7 L iv e a -1 8 L iv e a -2 1 P 9 ( 2 /2 ) P 9 ( 1 /2 ) P 9 ( 0 /2 ) P 6 ( 2 /2 ) P 6 ( 1 /2 ) P 6 ( 0 /2 ) (a) ( b) frequency proportion of survivors Overall disease resistance may correspond to levels of MHC heterozygosity [10,49]. The variation at the P9 pocket residues that we observed in amphibian populations long exposed to Bd is consistent with this view, although we did not observe any evidence of heterozygote advantage in our infection experiment. Amphibian populations with high frequencies of resistance alleles together with substantial MHC heterozygosity should be those most likely to recover from infection by Bd. Yet resistance alleles may have pleiotropic deleterious effects, including decreased growth, development or survival [8]. Heterozygosity may mitigate these effects in some conditions. It also may afford enhanced resistance against different Bd strains, as well as other fungal, bacterial and viral pathogens [8,9]. However, heterozygotes are not always more resistant than homozygotes to single infectious pathogens [53] such as Bd.
Further work is needed to assess the prevalence of the three pocket conformations in other Bd-resistant species. Also, our work suggests that the frequency of these conformations should increase in susceptible species as populations recover from chytridiomycosis epizootics. These MHC markers may make possible the identification of those amphibian populations most susceptible to Bd. Knowing which species are most at risk should facilitate the prioritization of conservation efforts and allow more accurate projections as to how Bd will affect complex host assemblages [54]. For those species now dependent on ex situ management for their very survival, selective breeding for Bd resistance may make possible successful re-introductions even into those areas in which Bd is enzootic. The molecular bases of resistance also may be critical to the development of immunization strategies. But much remains to be learned before these actions become practicable.
For most species included in our study, we have compared MHC-II b1 sequences obtained from one specific locus. Yet the MHC complex of many amphibians consists of at least two class II loci [14,15]. The structure and role of those additional loci in conferring disease resistance need to be further studied. In addition, the capacity of Bd to produce factors inhibiting amphibian adaptive immune response represents a major issue for any mitigation strategy based on adaptive immunity [51] that needs to be further elucidated.
Recent studies demonstrate that innate and adaptive immune defence systems can confer resistance against Bd [6,55,56]. Our study demonstrates for the first time that susceptibility to chytridiomycosis is associated with selection for specific immunogenetic traits across a wide variety of amphibians. Around the world, some amphibian species appear to have evolved immunity to chytridiomycosis through a common mechanism. Rescuing amphibian biodiversity will depend on our understanding of amphibian innate and adaptive immune defence mechanisms against Bd. The identification of adaptive markers for Bd resistance is an important step forward towards that goal.