key: cord-0330906-los56y94
authors: Byrne, Allison Q.; Waddle, Anthony W.; Saenz, Veronica; Ohmer, Michel; Jaeger, Jef R.; Richards-Zawacki, Corinne L.; Voyles, Jamie; Rosenblum, Erica Bree
title: Host Species is Linked to Pathogen Genotype for the Amphibian Chytrid Fungus (Batrachochytrium dendrobatidis)
date: 2021-11-24
journal: bioRxiv
DOI: 10.1101/2021.11.24.469827
sha: db37dd8fc8801a4e9dabc92beb92095fd3e11c8b
doc_id: 330906
cord_uid: los56y94

Host-pathogen specificity can arise from certain selective environments mediated by both the host and pathogen. Therefore, understanding the degree to which host species identity is correlated with pathogen genotype can help reveal historical host-pathogen dynamics. One animal disease of particular concern is chytridiomycosis, typically caused by the global panzootic lineage of the amphibian chytrid fungus (Batrachochytrium dendrobatidis, Bd), termed the Bd-GPL. This pathogen lineage has caused devastating declines in amphibian communities around the world. However, the origin of Bd-GPL and the fine-scale transmission dynamics of this lineage have remained a mystery. This is especially the case in North America where Bd-GPL is widespread, but disease outbreaks occur sporadically. Herein, we use Bd genetic data collected throughout the United States from amphibian skin swab and cultured isolate samples to investigate Bd genetic patterns. We highlight two case studies in Pennsylvania and Nevada where Bd-GPL genotypes are strongly correlated with host species identity. Specifically, in some localities bullfrogs (Rana catesbeiana) are infected with Bd-GPL lineages that are distinct from those infecting other sympatric amphibian species. Overall, we reveal a previously unknown association of Bd genotype with host species and identify the eastern United States as a Bd diversity hotspot and potential ancestral population for Bd-GPL.

Further phylogenetic analyses incorporating previously published Bd isolates shows that 223 that all Bd samples sequenced for this study belong to the Bd-GPL lineage (Figures 2,3) and that 224 the two clusters identified in our dataset correspond to the previously published Bd-GPL 225 subclades Bd-GPL1 and Bd-GPL2 (14,39). The difference between Bd-GPL1 and Bd-GPL2 226 samples is shown clearly in the combined PCA (Figure 1C The factors that may explain Bd genetic structure at each site were determined using an 245 AMOVA. We found that most of the variation in Bd genotypes (64.3% in PA, 58.4% in NV) can 246 be explained by variation between host species, rather than within an individual (34.7% in PA, 247 8.8% in NV) or between sampling sites (1.0% in PA, 32.8% in NV). Comparisons of variation 248 within individuals and variation between host species were statistically significant using 249 permutation tests to randomly permute samples between groups (999 times, p < 0.001 for each 250 test, Figure S2 ).

We 299 Bd can suppress host immune responses and some susceptible species show an over-activation of 300 the immune system that may worsen disease outcomes (41). At the same time, immunization 301 against Bd has been successful in some species, indicating the relationship between adaptive 302 immunity and Bd varies between hosts (42). Amphibian species identity is often the strongest 303 predictor of skin microbial community composition (43,44) and these microbes can be associated 304 with Bd disease outcomes (7,45,46), providing another possible mechanism for facilitating the 305 evolution of host-pathogen specificity. Finally, Bd isolates in the lab can show rapid 306 evolutionary change over short time periods (47,48). Therefore, differences between Bd 307 genotypes in the way they interact with amphibian skin microbes or immune systems, in 308 combination with Bd's evolutionary lability, could allow for the formation of the host-pathogen 309 genetic correlations that we detected.

In PA, we observed a ubiquitous pattern of association between Bd-GPL1 and bullfrogs 311 in all individual sites. Moreover, not a single bullfrog sample (N=18) in the PA region was

For each preamplification PCR reaction we used the FastStart High Fidelity PCR System 119 (Roche) at the following concentrations: 1x FastStart High Fidelity Reaction Buffer with MgCl 2 , 120 4.5mM MgCl 2 , 5% DMSO, 200µM PCR Grade Nucleotide Mix, 0.1 U/µl FastStart High Fidelity 121 Enzyme Blend. We added 1µl of cleaned DNA to each preamplification reaction and used the 122 following thermocycling profile: 95˚C for 10 min, 15 cycles of 95˚C for 15 sec

) and 124 incubated for 30 min at 30˚C, then 15 min at 80˚C. Treated products were diluted 1:5 in PCR-Juno™ LP 192.24 IFC (Fluidigm, Inc.) for library 128 preparation

Sequence Data Analysis 133 We pre-processed all sequencing data as described in (23) and generated consensus

These trees allow us to explore the relationship of the Bd 140 collected for this study to previously-published Bd samples representing all known Bd lineages 141 (N=13; 13,25). First, we trimmed our consensus sequence dataset to eliminate loci that had more 142 than 50% missing data for either the NV or PA datasets

We then collapsed all branches in all trees with <10 bootstrap 147 support and used Astral-III to generate a consensus tree

To further explore the genetic variation of Bd collected at our sites, we generated a 150 dataset of single nucleotide polymorphisms (SNPs). First, we aligned all raw sample reads (PA 151 and NV) and reads from 26 previously published Bd-GPL reference samples (13,25) using bwa the 50% missing data cutoff described above. Variants were then filtered using vcftools 156 (v.4.2) (32) to only include variants with a minor allele frequency > 0.01, quality > 30, less than 157 10% missing data, and minimum depth of 5. The final variant set

We 168 then ran the DAPC analysis specifying the optimal number of clusters. To determine the optimal 169 number of PCs to use for this analysis we ran the adegenet function optim.a.score which 170 calculates the alpha score as (P t -P r ), where P t is probability of reassignment to the true cluster 171 and P r is the probability of reassignment to randomly permuted clusters

Samples with posterior probability of inclusion in one cluster < 0.99 were labeled "unassigned

We plotted the first two PCs which together explain 68.2% of the variation in our data. We then 176 computed a PCA separately for each locality to show within locality diversity across species and 177 sites

To evaluate the relative contribution of sampling site, host species, and individual 179 variation, we ran an analysis of molecular variance (AMOVA) (36) using our SNP data. This 180 analysis was conducted separately for the PA and NV datasets and determined the relative 181 influence of within individual variation

We ran this analysis using the poppr package in R (37) and only used SNPs 183 with <5% missing data. This resulted in 40 high quality SNPS for PA and 126 high quality SNPs 184 for NV. We then tested for statistical significance using the "randtest

To evaluate the frequency of mixed infections on individuals, we searched our set of 187 amplicon sequences for those that had diagnostic alleles that discriminate between the Bd-GPL1

We identified amplicons that had at least 1 SNP that distinguished 190 GPL1 samples from GPL2 (as defined by both the DAPC and phylogenetic analyses described 191 above). The GPL1 Bd isolates in this reference set were CJB4

These samples were all sequenced from pure 195 cultures grown in the lab

From these samples we identified 24 loci that 197 occur on five separate Bd chromosomes that have specific Bd-GPL1 or Bd-GPL2 alleles. For 198 the diagnostic alleles on chromosome 1, heterozygotes (those having both alleles present) are 199 characteristic of Bd-GPL2 isolates (see Figure S3). We then manually scored each sample from 200 PA and NV at each as having either the Bd-GPL1

In the combined PCA, we see a clear division of Bd genotypes into 207 two clusters, with the presence of intermediate "unassigned" samples (Figure 1C). By evaluating 208 each site separately, we see that PA harbors much more genetic variability within samples 209 (accounting for all "unassigned" Bd genotypes, Figure 1). We confirmed that "unassigned" 210 samples were not simply the product of missing data by showing that each genotype group

Figure 1: Map of sampling localities for Bd in southern Nevada (A) and northwestern

Symbols and ellipses colored as in maps above. 221 312 infected with Bd-GPL2, or included in the "unassigned" group. The finding that host species was region and has coevolved with native amphibians, including bullfrogs. We can see additional 316 evidence for this hypothesis at individual sites; at PA Site 4 bullfrogs and the closely related 317 green frog (Rana clamitans) were both preferentially infected with Bd-GPL1, while the more 318 distantly related northern leopard frog (Rana pipiens) was infected with Bd-GPL2 (49)

An alternative possibility that could produce patterns of host species-Bd genotype

In NV the association 327 between Bd genotypes and host species is not as consistent as in PA, indicating the coexistence 328 of different Bd genotypes in the NV region may be a more recent phenomenon. For example, at 329 NV Site 4 Bd genotypes were strictly partitioned by species, with Bd-GPL1 only present on 330 bullfrogs; however, at NV Site 1 all individuals sampled across species were infected with the 331 same Bd-GPL2 genotype

Supporting this hypothesis, one of Bd-GPL? 343 The implied coexistence of Bd and amphibians at evolutionary time scales in PA could 344 indicate that the ancestral population of Bd-GPL existed in this part of the world. This possibility 345 has been previously proposed, following reports of higher diversity in Bd samples collected from 346 bullfrog populations (53,54), and evidence that bullfrogs may have facilitated Bd spread from 347 east to west in the US (55)

In other 350 regions, such as Panama and Mexico, genetic diversity of Bd-GPL is either spatially structured 351 across the region (Mexico) or genetically homogenous (Panama), the latter indicating a more 352 recent introduction and spread (15,16)

US (within the native range of the bullfrog) may have been the ancestral population for Bd-GPL

What can the "unassigned" samples tell us about Bd dynamics

These samples, all of 360 which were from swabs collected in PA, show high levels of heterozygosity across the entire 361 panel of diagnostic alleles. There are two possible explanations for the occurrence of a highly 362 heterozygous, intermediate Bd genotype. The first is that these individuals were coinfected with 363 both Bd-GPL1 and Bd-GPL2. The second is that they were infected with a hybrid lineage 364 resulting from recombination between Bd-GPL1 and Bd-GPL2. Indeed, if we look at the 365 reference Bd isolates in Figure S3 we see that some Bd-GPL2 pure isolates are heterozygous at 366 certain alleles

Our data show not only the coexistence of multiple Bd-GPL lineages at a single site, but 370 also highlights areas where Bd lineage mixing on individuals is rare (in NV, and PA Sites 3 and 371 4), or common (PA Site 1). The most common frog species caught at PA Site 1 was the green 372 frog (R. clamitans), in fact this species accounts for 71% of all

At this point it is unclear whether there is something unique about PA Site 1, or R. clamitans, 374 that may allow for such consistent lineage mixing or the potential presence of hybrid lineages

To our knowledge no chytridiomycosis outbreaks have been reported in R. clamitans, and 376 individuals from the same site in PA showed no increase in stress hormones when 377 experimentally infected with Bd (58). Probably, this species is a tolerant carrier of Bd, like the 378 closely related bullfrog. This host-pathogen relationship is worth further exploration as the 379 potential for

This process has been to exploit hosts (61). Parasexual reproduction often results in aneuploidy, which is a 388 common phenomenon in Bd (25,62). Furthermore, elevated chromosomal copy number has been 389 linked to increased virulence in Bd (62), and hybrids of Bd have been more deadly than parent 390 lineages when paired with certain host species (63). Indeed, it was originally proposed that Bd-391 GPL could have originated from recombination between closely related Bd lineages (12)

We observed that the association of Bd-GPL1 with bullfrogs is 400 consistent in NV where bullfrogs have been more recently introduced, and that the patterns in 401 this region are consistent with more recent introductions of novel Bd lineages. We also found 402 evidence of either mixed Bd infections or a potential hybrid lineage that was common at a single 403 site in PA and on green frogs across all PA sites. Given the high genetic diversity and evidence 404 of evolutionary relationships between Bd-GPL lineages and particular host species

North America may be the ancestral origin of Bd-GPL. Overall, this study presents evidence that

the US, reveals more about the 407 relationship of the bullfrog to Bd, and points to a possible cradle of Bd-GPL genetic diversity in

Smith for their assistance with 411 funding acquisition and management for the Nevada study. We thank the following people for 412 their assistance with field and laboratory work in Nevada: Rebeca Rivera

We thank the following 415 people for assistance with field and laboratory work in Pennsylvania: Laura Brannelly

Reproducibility Statement

All data used for this study and R code used for the analyses are available on github

Red Queen Meets Santa Rosalia: Arms Races and the Evolution of Host 426 Specialization in Organisms with Parasitic Lifestyles

Parasite and host assemblages: embracing the 429 reality will improve our knowledge of parasite transmission and virulence

Evolutionary implications of host-pathogen specificity: fitness 433 consequences of pathogen virulence traits

Genomic 435 surveillance reveals multiple introductions of SARS-CoV-2 into Northern California

Batrachochytrium Dendrobatidis gen. et sp. nov., a 439 Chytrid Pathogenic to Amphibians

Riding the Wave: Reconciling the 442 Roles of Disease and Climate Change in Amphibian Declines

PREDICTED DISEASE 445 SUSCEPTIBILITY IN A PANAMANIAN AMPHIBIAN ASSEMBLAGE BASED ON 446 SKIN PEPTIDE DEFENSES

Amphibian declines in Ecuador: overview and first report of 449 chytridiomycosis from

Investigating the Population-Level 451 effects of Chytridiomycosis: An Emerging Infectious Disease of Amphibians

Chytridiomycosis causes amphibian mortality associated with population declines in the 456 rain forests of Australia and Central America

Recent Asian 459 origin of chytrid fungi causing global amphibian declines

Multiple 463 emergences of genetically diverse amphibian-infecting chytrids include a globalized 464 hypervirulent recombinant lineage

Cryptic 467 diversity of a widespread global pathogen reveals expanded threats to amphibian 468 conservation

Novel, 471 panzootic and hybrid genotypes of amphibian chytridiomycosis associated with the 472 bullfrog trade

475 Divergent regional evolutionary histories of a devastating global amphibian pathogen

Early presence of 478 Batrachochytrium dendrobatidis in Mexico with a contemporary dominance of the global 479 panzootic lineage

Minimising 481 exposure of amphibians to pathogens during field studies

Batrachochytrium dendrobatidis and the Decline and Survival of the Relict Leopard Frog

Rapid quantitative detection of 487 chytridiomycosis ( Batrachochytrium dendrobatidis ) in amphibian samples using real-488 time Taqman PCR assay

Population-level 491 resistance to chytridiomycosis is life-stage dependent in an imperiled anuran

Systematic approach 494 to isolating Batrachochytrium dendrobatidis

Inheritance of DNA methylation in Coprinus cinereus

Unlocking the 499 story in the swab: A new genotyping assay for the amphibian chytrid fungus 500

Batrachochytrium dendrobatidis

503 Batrachochytrium salamandrivorans sp. nov. causes lethal chytridiomycosis in 504 amphibians

Complex 507 history of the amphibian-killing chytrid fungus revealed with genome resequencing data

MUSCLE: multiple sequence alignment with high accuracy and high 512 throughput

Geneious 515 Basic: An integrated and extendable desktop software platform for the organization and 516 analysis of sequence data

ASTRAL-III: polynomial time species tree 522 reconstruction from partially resolved gene trees

Fast and accurate short read alignment with Burrows-Wheeler transform

Haplotype-based variant detection from short-read sequencing. 528 arXiv Prepr arXiv12073907

The variant 530 call format and VCFtools

vcfr: a package to manipulate and visualize variant call format 533 data in R

Discriminant analysis of principal components: a new 536 method for the analysis of genetically structured populations

adegenet: a R package for the multivariate analysis of genetic markers

Analysis of molecular variance inferred from metric 541 distances among DNA haplotypes: application to human mitochondrial DNA restriction 542 data

Poppr: an R package for genetic analysis of 544 populations with clonal, partially clonal, and/or sexual reproduction

Fighting a 558 Losing Battle: Vigorous Immune Response Countered by Pathogen Suppression of Host 559 Defenses in the Chytridiomycosis-Susceptible Frog Atelopus zeteki. G3 (Bethesda)

Amphibian 563 resistance to chytridiomycosis increases following low-virulence chytrid fungal infection 564 or drug-mediated clearance

The 566 amphibian skin-associated microbiome across species, space and life history stages

Co-habiting amphibian 570 species harbor unique skin bacterial communities in wild populations

Proportion of individuals with anti-573

Batrachochytrium dendrobatidis skin bacteria is associated with population persistence in 574 the frog Rana muscosa

Greater Species Richness of Bacterial Skin Symbionts Better Suppresses the Amphibian 578 Fungal Pathogen Batrachochytrium Dendrobatidis

A Fungal 581 Pathogen of Amphibians, Batrachochytrium dendrobatidis, Attenuates in Pathogenicity 582 with In Vitro Passages

Genomic 585 Correlates of Virulence Attenuation in the Deadly Amphibian Chytrid Fungus

Spatiotemporal Diversification of the True Frogs (Genus Rana): A Historical 590 Framework for a Widely Studied Group of Model Organisms

Phylogeny and geography predict pathogen community 593 similarity in wild primates and humans

Amphibians and Reptiles in Nevada

Opening the file drawer: 598 Unexpected insights from a chytrid infection experiment

Rapid 601 Global Expansion of the Fungal Disease Chytridiomycosis into Declining and Healthy 602 Amphibian Populations

Amphibian 605 chytridiomycosis in Japan: distribution, haplotypes and possible route of entry into Japan

Introduced bullfrog facilitates pathogen 609 invasion in the western United States

A century of 612 Batrachochytrium dendrobatidis in Illinois amphibians (1888-1989)

Prevalence of Batrachochytrium 616 dendrobatidis in 120 archived specimens of Lithobates catesbeianus (American bullfrog) 617 collected in California

Relationships between 619 glucocorticoids and infection with Batrachochytrium dendrobatidis in three amphibian 620 species

Amphibian-killing chytrid in Brazil comprises both locally 624 endemic and globally expanding populations

Hybridization Facilitates

Sex and Virulence of Human Pathogenic Fungi

Chromosomal 633 copy number variation, selection and uneven rates of recombination reveal cryptic genome 634 diversity linked to pathogenicity

636 Hybrids of amphibian chytrid show high virulence in native hosts

Chytrid fungi and global amphibian declines

A) PCA as shown in Figure 1C colored by proportion missing data. (B) BIC chart for 644 DAPC analysis showing K = 2 clusters as the optimal number of clusters (or "elbow" in the 645 chart). (C) Density plot showing that unassigned

Results from AMOVA testing the relative importance of variation within samples 649 (top), variation between host species (middle), and variation between sites (bottom) in explaining

Bd genetic variance for northwestern Pennsylvania (left) and southern Nevada (right)

Each histogram shows the results from 653 randomly permuting the sample matrices 999 times using the rand.test function in the R package 654 ade4. The black dot indicates the value obtained in this study with a red asterisk indicating 655 statistical significance

A) Astral phylogenetic tree and (B) allele plot for 26 pure-isolate Bd-GPL reference 658 samples used to identify the 24 Bd-GPL1 and Bd-GPL2 diagnostic alleles shown in Figures 6 659 and 7. Tree displayed as in Figures 3 and 4