Abstract

Globalized infectious diseases are causing species declines worldwide, but their source often remains elusive. We used whole-genome sequencing to solve the spatiotemporal origins of the most devastating panzootic to date, caused by the fungusBatrachochytrium dendrobatidis, a proximate driver of global amphibian declines. We traced the source ofB. dendrobatidis to the Korean peninsula, where one lineage,BdASIA-1, exhibits the genetic hallmarks of an ancestral population that seeded the panzootic. We date the emergence of this pathogen to the early 20th century, coinciding with the global expansion of commercial trade in amphibians, and we show that intercontinental transmission is ongoing. Our findings point to East Asia as a geographic hotspot forB. dendrobatidis biodiversity and the original source of these lineages that now parasitize amphibians worldwide.

Global diversity ofB. dendrobatidis

To resolve these inconsistencies, we isolatedB. dendrobatidis from all the candidate source continents and sequenced the genomes of 177 isolates to high depth, then combined our data with published genomes from three prior studies (11,12,18) to generate a globally representative panel of 234 isolates (Fig. 1A andfig. S1). This data set covers all continents from whichB. dendrobatidis has been detected to date, and spans infections of all three extant orders of Amphibia (fig. S1 and table S1). Mapped against theB. dendrobatidis reference genome JEL423, our sequencing recovered 586,005 segregating single-nucleotide polymorphisms (SNPs). Phylogenetic analysis recovered all previously detected divergent lineages (Fig. 1B andfig. S2). The previously accepted lineagesBdGPL (global),BdCAPE (African),BdCH (European), andBdBRAZIL (Brazilian) were all detected (19), but our discovery of a new hyperdiverse lineage in amphibians native to the Korean peninsula (BdASIA-1) redefined these lineages and their relationships. TheBdCH lineage, which was previously thought to be enzootic to Switzerland (11), now groups with theBdASIA-1 lineage. A second Asian-associated lineage (BdASIA-2) was recovered from invasive North American bullfrogs in Korea and is closely related to the lineage that is enzootic to the Brazilian Atlantic forest (BdBRAZIL) (20). It was not possible to infer the direction of intercontinental spread between isolates within this lineage, so it was namedBdASIA-2/BdBRAZIL. Conditional on the midpoint rooting of the phylogeny inFig. 1B, we now define the main diverged lineages asBdGPL,BdCAPE,BdASIA-1 (which includes the singleBdCH isolate), andBdASIA-2/BdBRAZIL. Previous phylogenetic relationships developed using the widely used ribosomal intragenic spacerITS-1 region do not accurately distinguishB. dendrobatidis lineages (fig. S3), and this likely explains much of the place-of-origin conflict in the literature (15–17).

Pairwise comparisons among isolates within each lineage show that the average number of segregating sites is greater forBdASIA-1 than for any other lineage by a factor of 3 (Fig. 1A andTable 1) and that nucleotide diversity (π;fig. S4) is greater by a factor of 2 to 4. Seven of our eightBdASIA-1 isolates were recently cultured from wild South Korean frogs, and the other came from the pet trade in Belgium; all eight were aclinical infections. These isolates show that the Korean peninsula is a global center ofB. dendrobatidis diversity and that East Asia may contain the ancestral population ofB. dendrobatidis, as suggested by Batailleet al. (17).

We investigated this hypothesis further using Bayesian-based haplotype clustering (21) and found the greatest haplotype sharing among isolates withinBdASIA-1 and betweenBdASIA-1 and all other lineages (fig. S5). This provides direct genetic evidence thatBdASIA-1 shares more diversity with the global population ofB. dendrobatidis than any other lineage. In an independent test of ancestry, we used OrthoMCL (22) to root aB. dendrobatidis phylogeny to its closest known relative,B. salamandrivorans, which currently threatens salamanders (23). This tree indicates that the Asian and Brazilian isolates ofB. dendrobatidis lie outside a clade comprising all other isolates (fig. S6 and table S2). To identify the signature of demographic histories across lineages, we used Tajima’sD (24). Genome scans of most lineages showed highly variable positive and negative values ofD with maximum amplitude exhibited byBdGPL (–2.6 to +6.2;Fig. 2F), indicating that these lineages (BdASIA-2/BdBRAZIL,BdCAPE, andBdGPL) have undergone episodes of population fluctuation and/or strong natural selection that are consistent with a history of spatial and host radiations. In striking contrast,BdASIA-1 shows a flat profile for Tajima’sD (Fig. 2F) indicating mutation-drift equilibrium likely reflective of pathogen endemism in this region.

Dating the emergence ofBdGPL

The broad range of previous estimates for the TMRCA ofBdGPL spanning 26,000 years (11,12) can be explained by two sources of inaccuracy: (i) unaccounted recombination and (ii) the application of unrealistic evolutionary rates. To address these, we first interrogated the 178,280-kbp mitochondrial genome (mtDNA), which has high copy number and low rates of recombination relative to the nuclear genome. To resolve the structure of the mtDNA genome, we resorted to long-read sequencing using a MinION device (Oxford Nanopore Technologies, Cambridge, UK), which allowed us to describe this molecule’s unusual configuration;B. dendrobatidis carries three linear mitochondrial segments, each having inverted repeats at the termini with conserved mitochondrial genes spread over two of the segments (fig. S7). Additionally, we sought regions of the autosomal genome with low rates of recombination to obtain an independent estimate of the TMRCA ofBdGPL.

Detection of crossover events in theB. dendrobatidis autosomal genome (18) using a subset of the isolates in this study revealed a large (1.66 Mbp) region of Supercontig_1.2 inBdGPL that exhibits several features that identified it as a recombination “coldspot”: (i) a continuous region of reduced Tajima’sD (Fig. 2F), (ii) sustained high values of population differentiation as measured by the fixation index (F_ST) when compared with all other lineages (Fig. 3A), (iii) a continuous region of reduced nucleotide diversity (π;fig. S4), and (iv) shared loss of heterozygosity (fig. S8). We expanded sampling to infer the temporal range of pathogen introductions using a broad panel of isolates with known date of isolation (n = 184, ranging from 1998 to 2016) and whole-genome RNA baiting to obtain reads from preserved amphibians that had died of chytridiomycosis. We then investigated whether our data set contained sufficient signal to perform tip-dating inferences by building phylogenetic trees using PhyML (25) (Fig. 2, A and C). We fitted root-to-tip distances to collection dates both at the whole-tree and within-lineage scales. We observed a positive and significant correlation withinBdGPL only, for both the mitochondrial and nuclear genomes, demonstrating sufficient temporal signal to perform thorough tip-dating inferences at this evolutionary scale (Fig. 2, B and D).

Tip-dating in BEAST was used to co-estimate ancestral divergence times and the rate at which mutations accumulate within theBdGPL lineage. The mean mitochondrial substitution rate was 1.01 × 10⁻⁶ substitutions per site per year [95% highest posterior density (HPD), 4.29 × 10⁻⁷ to 1.62 × 10⁻⁶]. The mean nuclear substitution rate was 7.29 × 10⁻⁷ substitutions per site per year (95% HPD, 3.41 × 10⁻⁷ to 1.14 × 10⁻⁶), which is comparable to a recent report of an evolutionary rate of 2.4 × 10⁻⁶ to 2.6 × 10⁻⁶ substitutions per site per year for another unicellular yeast,Saccharomyces cerevisiae beer strains (26). These rate estimates are faster by a factor of >300 than the rate used in a previous study (12) to obtain a TMRCA of 26,400 years forBdGPL. Accordingly, we estimate that the ancestor of the amphibian panzooticBdGPL originated between 120 and 50 years ago (Fig. 2E), with HPD estimates of 1898 (95% HPD, 1809 to 1941) and 1962 (95% HPD, 1859 to 1988) for the nuclear and mitochondrial dating analyses, respectively (Fig. 2E).

We considered an additional calibration approach for the TMRCA of the mitochondrial genome where we included informative priors on nodes around the dates for the first historical descriptions ofBdGPL detection in Australia (1978), Central America (1972), Sierra de Guadarrama (Europe) (1997), and the Pyrenees (Europe) (2000). We did not include priors for nodes where observed declines have been reported but where the lineage responsible for those declines is unknown. This mixed dating method based on tip and node calibration yielded very similar estimates [TMRCA estimates of 1975 (95% HPD, 1939 to 1989) (fig. S9)], further strengthening our confidence in a recent date of emergence forBdGPL. An expansion ofBdGPL in the 20th century coincides with the global expansion in amphibians traded for exotic pet, medical, and food purposes (27,28). Within our phylogeny, we found representatives from all lineages among traded animals (figs. S10 to S14) and identified 10 events where traded amphibians were infected with non-enzootic isolates (Fig. 4). This finding demonstrates the ongoing failure of international biosecurity despite the listing ofB. dendrobatidis by the World Organisation for Animal Health (OIE) in 2008.

Hybridization between recontacting lineages ofB. dendrobatidis

To determine the extent to which the four main lineages ofB. dendrobatidis have undergone recent genetic exchange, we used the site-by-site–based approach implemented in STRUCTURE (29). Although most isolates could be assigned unambiguously to one of the four main lineages, we identified three hybrid genotypes (Fig. 3B), including one previously reported hybrid (isolate CLFT024/2) (20), and discovered two newly identified hybrids ofBdGPL andBdCAPE in South Africa. Furthermore,BdCH (isolate 0739) appears to be a chimera of multiple lineages that may represent unsampled genomic diversity residing in East Asia, rather than true hybridization. These hybrid genomes demonstrate thatB. dendrobatidis is continuing to exchange haplotypes among lineages when they interact after continental invasions, generating novel genomic diversity. We analyzed isolate clustering using principal components analysis on a filtered subset of 3900 SNPs in linkage equilibrium, revealing an overall population structure that is consistent with our phylogenetic analyses (Fig. 3C). In addition, the putatively identified hybrid isolates ofB. dendrobatidis were shown to fall between main lineage clusters (Fig. 3C), further strengthening our hypothesis of haplotype exchange occurring during secondary contact between lineages.

Associations among lineage, virulence, and declines

Genotypic diversification of pathogens is commonly associated with diversification of traits associated with host exploitation (30) and is most commonly measured as the ability to infect a host and to cause disease post-infection. We tested for variation of these two phenotypic traits across fourB. dendrobatidis lineages by exposing larval and postmetamorphic common toads (Bufo bufo). Larvae are highly susceptible to infection but do not die before metamorphosis, in contrast to postmetamorphic juveniles, which are susceptible to infection and fatal chytridiomycosis (31). In tadpoles, bothBdGPL andBdASIA-1 were significantly more infectious thanBdCAPE andBdCH (fig. S15 and tables S3 and S4). In metamorphs,BdGPL was significantly more infectious than the other treatments, relative to the control group, and was significantly more lethal in experimental challenge than the geographically more restrictedBdCAPE,BdASIA-1, andBdCH (Fig. 2G). We further tested for differences in virulence among lineages by using our global data set to examine whether chytridiomycosis was nonrandomly associated withB. dendrobatidis lineage. We detected a significant difference (P < 0.001) in the proportion of isolates associated with chytridiomycosis among the three parental lineages (BdASIA-1 andBdASIA-2/BdBRAZIL were grouped because of low sample sizes), and post hoc tests indicated significant excess in virulence in bothBdGPL andBdCAPE lineages relative to the combinedBdASIA-1 andBdASIA-2/BdBRAZIL (allP < 0.05). However, we did not detect a significant difference betweenBdGPL andBdCAPE (fig. S16 and table S5). These data suggest that althoughBdGPL is highly virulent, population-level outcomes are also context-dependent (32); under some conditions, other lineages can also be responsible for lethal amphibian disease and population declines (33).

Historical and contemporary implications of panzootic chytridiomycosis

Our results point to endemism ofB. dendrobatidis in Asia, out of which multiple panzootic lineages have emerged. These emergent diasporas include the virulent and highly transmissibleBdGPL, which spread during the early 20th century via a yet unknown route to infect close to 700 amphibian species out of ~1300 thus far tested (34). With more than 7800 amphibian species currently described, the number of affected species is likely to rise. The international trade in amphibians has undoubtedly contributed directly to vectoring this pathogen worldwide (Fig. 4) (35,36), and within our phylogeny we identified many highly supported (≥90% bootstrap support) clades on short branches that linked isolates collected from wild amphibian populations across different continents (Fig. 4 andfigs. S10 to S14). However, the role of globalized trade in passively contributing to the spread of this disease cannot be ruled out. It is likely no coincidence that our estimated dates for the emergence ofBdGPL span the globalization “big bang”—the rapid proliferation in intercontinental trade, capital, and technology that started in the 1820s (37). The recent invasion of Madagascar by Asian common toads hidden within mining equipment (38) demonstrates the capacity for amphibians to escape detection at borders and exemplifies how the unintended anthropogenic dispersal of amphibians has also likely contributed to the worldwide spread of pathogenic chytrids.

The hyperdiverse hotspot identified in Korea likely represents a fraction of theBatrachochytrium genetic diversity in Asia, and further sampling across this region is urgently needed because the substantial global trade in Asian amphibians (39) presents a risk of seeding future outbreak lineages. Unique ribosomal DNA haplotypes ofB. dendrobatidis have been detected in native amphibian species in India (40,41), Japan (16), and China (42). Although caution should be observed when drawing conclusions about lineages based on short sequence alignments (fig. S3), other endemic lineages probably remain undetected within Asia. It is noteworthy that the northern European countryside is witnessing the emergence ofB. salamandrivorans, which also has its origin in Asia. The emergence ofB. salamandrivorans is linked to the amphibian pet trade (43), and the broad expansion of virulence factors that are found in the genomes of these two pathogens is testament to the evolutionary innovation that has occurred in these AsianBatrachochytrium fungi (23).

Our findings show that the global trade in amphibians continues to be associated with the translocation of chytrid lineages with panzootic potential. Ultimately, our work confirms that panzootics of emerging fungal diseases in amphibians are caused by ancient patterns of pathogen phylogeography being redrawn as largely unrestricted global trade moves pathogens into new regions, infecting new hosts and igniting disease outbreaks. Within this context, the continued strengthening of transcontinental biosecurity is critical to the survival of amphibian species in the wild (44).

Supplementary Material

Acknowledgments

DNA sequencing was carried out in the NBAF GenePool genomics facility at the University of Edinburgh, and we thank the GenePool staff for their assistance. This work used the computing resources of the UK MEDical BIOinformatics partnership - aggregation, integration, visualization and analysis of large, complex data (UK MED-BIO) which is supported by the Medical Research Council (grant number MR/L01632X/1). We thank J. Rhodes for the provision of flow cells and reagents for MinION sequencing, the staff at Oxford Nanopore Technologies for admission to the MinION Early Access Programme, and three anonymous reviewers for constructive comments and suggestions.

Funding: S.J.O., T.W.J.G., L.Br., A.Lo., A.A.C., D.S.S., E.A.C., C.M., J.B., D.M.A., F.C., and M.C.F. were supported through NERC (standard grant NE/K014455/1). S.J.O. acknowledges a Microsoft Azure for Research Sponsorship (subscription ID: ab7cd695-49cf-4a83-910a-ef71603e708b). T.W.J.G., A.Lo., A.A.C., D.S.S., E.A.C., C.M., J.B., D.M.A., F.C., and M.C.F. were also supported by the EU BiodivERsA scheme (RACE, funded through NERC directed grant NE/G002193/1 and ANR-08-Biodiversa-002-03) and NERC (standard grant NE/K012509/1). M.C.F., E.A.C., and C.M. acknowledge the Nouragues Travel Grant Program 2014. R.A.F. was supported by an MIT/Wellcome Trust Fellowship. T.W.J.G. was supported by the People’s Trust for Endangered Species and the Morris Animal Foundation (D12ZO-002). J.M.G.S. and M.C.F. were supported by the Leverhulme Trust (RPG-2014-273) and the Morris Animal Foundation (D16ZO-022). F.B. was supported by the ERC (grant ERC 260801–Big_Idea). D.M.A. was funded by Wellcome Trust grant 099202. J.V. was supported by the Hungarian Scientific Research Fund (OTKA K77841) and Bolyai János Research Scholarship, Hungarian Academy of Sciences (BO/00579/14/8). D.J.G. was supported by the Conservation Leadership Programme (grant 0134010) with additional assistance from F. Gebresenbet, R. Kassahun, and S. P. Loader. C.S.-A. was supported by Fondecyt Nº11140902 and 1181758. T.M.D.-B. was supported by the Royal Geographical Society and the Royal Zoological Society of Scotland with assistance from M. Hirschfeld and the Budongo Conservation Field Station. B.W. was supported by the National Research Foundation of Korea (2015R1D1A1A01057282). L.F.T. was supported by FAPESP (#2016/25358-3) and CNPq (#300896/2016-6). L.Be., L.F.S., and R.J.W. were supported by the Australian Research Council (FT100100375, DP120100811). A.A.C. was supported by a Royal Society Wolfson Research Merit award. J.H., A.La., and S.M. were funded by the Swedish Research Council Formas (grant no. 2013-1389-26445-20). C.W. was funded by the National Research Foundation, South Africa. T.Y.J. and T.S.J. acknowledge NSF grant DEB-1601259. W.E.H. was funded by the NSERC Strategic and Discovery grant programs.

Footnotes

References

Associated Data

This section collects any data citations, data availability statements, or supplementary materials included in this article.