Two distinct host-specialized fungal species cause white-nose disease in bats

Abstract

The emergence of infectious diseases, particularly those caused by fungal pathogens, poses serious threats to public health, wildlife and ecosystem stability¹. Host-fungus interactions and environmental factors have been extensively examined^2-4. However, the role of genetic variability in pathogens is often less well-studied, even for diseases such as white-nose in bats, which has caused one of the highest disease-driven death tolls documented in nonhuman mammals⁵. Previous research on white-nose disease has primarily focused on variations in disease outcomes attributed to host traits or environmental conditions^6-8, but has neglected pathogen variability. Here we leverage an extensive reference collection of 5,479 fungal isolates from 27 countries to reveal that the widespread causative agent is not a single species but two sympatric cryptic species, each exhibiting host specialization. Our findings provide evidence of recombination in each species, but significant genetic differentiation across their genomes, including differences in genome organization. Both species contain geographically differentiated populations, which enabled us to identify the species introduced to North America and trace its source population to a region in Ukraine. In light of our discovery of the existence of two cryptic species of the causative agent of white-nose disease, our research underscores the need to integrate the study of pathogen variability into comprehensive disease surveillance, management and prevention strategies. This holistic approach is crucial for enhancing our understanding of diseases and implementing effective measures to prevent their spread.

Fig. 1. Multilocus microsatellite typing reveals two clades of P.destructans.

a, Sampling locations in Eurasia (North American sites not shown; n = 255). b, Phylogenetic tree, based on the D_A distance, representing the relationship between all 1,866 unique multilocus genotypes (North American isolates are indicated with an asterisk) originating from the 5,479 isolates from 264 sites based on 18 microsatellite loci (for the D_A distance, see the section ‘Analyses of multilocus genotypes’ in the Methods; the monophyly of Pd-1 and Pd-2 was consistently recovered when jackknifing loci). For better visualization of both clades, the section containing Pd-2 was magnified. c, Principal component analysis (bottom) of isolates (Pd-1 was subsampled to ensure even sampling between clades and to maximize geographical coverage, which resulted in 234 isolates for Pd-1 and 92 isolates of Pd-2; Supplementary Table 1). Density (top) of principal component 1 (PC1) coefficients after random subsampling of Pd-1 to obtain the same number of sampling sites as for Pd-2 (51 sites, replicated 100,000 times). Dotted lines represent the absolute edges of the distribution (see the section ‘Analyses of multilocus genotypes’ in the Methods for details).

Fig. 2. Genomic differentiation between clades.

a, Sampling locations of the 18 isolates used for phylogenetic and sequence divergence analyses. b, Phylogenetic tree of 664 BUSCO genes with 1,000 bootstraps and site concordance factors for nodes of interest (as percentages). The branch to the outgroup Gd267 has been shortened for visualization purposes. c, Boxplot of the pairwise genetic distances between isolates (each isolate is compared with the 17 other isolates) for the 664 BUSCO genes, partitioned between intra-clade and inter-clade distance. The darker line indicates the median, and the lower and upper hinges represent the first and third quartiles, respectively. The whiskers extend to 1.5 times the interquartile range, with any data points beyond this range marked as outliers. d, Genomic differentiation between pools of 69 and 63 individuals from clades Pd-1 and Pd-2, respectively, estimated using F_ST across a window size of 200 SNPs, using Gd293 (Pd-1 clade) as the reference genome, with its 18 contigs successively coloured. e, Heatmap of phylogenomic profiling depicting clade-specific synteny network clusters, with isolates in rows and clusters in columns (n = 1,365; profile for run 13, Supplementary Table 17). Grey denotes the absence of a gene, whereas other colours indicate the number of gene copies (see the key).

Fig. 3. Strong population differentiation in Pd-1.

a, Estimation of effective migration (m) surfaces based on 2,261 isolates from Pd-1 in Europe (all sites excluding Russia and the United States after clone correction, n = 225 sites). For visualization, results from eight independent runs (each with 8 million iterations and between 100 and 450 demes) were combined. Different shades of colour represent variable levels of high (blue) or low (brown) effective migration rates. Sampling locations are represented by red dots. b, Distribution of the distance between the true and assigned site of each Eurasian isolate of Pd-1 (2,191 isolates; dataset limited to 20 isolates per site) for the observed and randomized datasets of DAPCs (bin width, 100 km).

Fig. 4. Assignment of the source of the North American introduction of P.destructans using Bayesian inference.

a, Bayesian inference (SPASIBA analysis; trained on 5,162 Pd-1 isolates) was performed independently for each of the 33 North American isolates across the continuous landscape (>5 million square kilometres). We then calculated the log-likelihood of the assignment across all 33 isolates (Supplementary Table 5), depicted on the map by a colour scale from blue (lowest log-likelihood) to red (highest log-likelihood). Green dots represent sites from which Pd-1 samples were collected, and the two yellow dots, in Podillia, Ukraine, indicate the inferred origin of the introduced ancestor of the 33 North American isolates (Supplementary Table 16). The zoomed-in area shown in b is represented here by the black rectangle. b, Details of the region with the most probable assignment. This area has a 6 × 10⁹ greater likelihood of assignment than any pixel occurring outside this region (Supplementary Table 5). The central black-bordered pixel has the highest likelihood of assignment and contains the site to which the DAPC assigned all 33 North American isolates. The remaining black-bordered pixels have a likelihood within an order of magnitude of the central black-bordered pixel. Scale bar, 50 km.

Extended Data Fig. 1. Percentages of swab samples collected from the 6 most frequently sampled bat species or combination of bat species.

(Myotis myotis/M. blythii, Myotis nattereri/M. crypticus/M. escalerai, M. mystacinus, Myotis daubentonii, Myotis dasycneme, M. brandtii) and all other species or combination of species combined (“Other”) per clade (Eurasian sites only). Morphologically cryptic/highly similar species were treated together due to the difficulty of reliable species identification during winter hibernation when bats are not handled to minimise disturbance. Samples from substrates other than bats (n = 267) or without bat species information (n = 1) were not included in this figure, resulting in data from 1,388 and 92 swabs for Pd-1 and Pd-2, respectively. Note that 17 swabs out of 1,463 harboured isolates from both clades and are thus used to calculate percentages in both graphs. See the section ‘Statistical analyses’ in the Methods for statistical analyses formally testing the relationship between clade identity (Pd-1 or Pd-2) and environmental factors, including bat species.

Extended Data Fig. 2. Colouration of culture medium by growth of Pd-1 and Pd-2.

Left: Difference in colour (darkness) of culture medium after 7 weeks using pixel density as a proxy. Photos were taken of 45 and 34 isolates of Pd-1 and Pd-2 respectively and analysed in R (see section ‘Analysis of culture darkness’ in the Methods and Table S8). To calculate the difference, the median pixel density at 8 weeks was subtracted from the pixel density at 1 week whereby a positive value indicates an increase in darkness. The black line indicates the median, while the lower and upper hinges represent the first and third quartiles, respectively. The whiskers extend to 1.5 times the interquartile range. Right: Examples of analysed photos (original) for Pd-1 and Pd-2. Clades Pd-1 and Pd-2 differ significantly in terms of colouration of the agar medium (Week 1 – Week 8; two-sided t-test: t = −11.58, df = 60.06, p < 0.001).

Extended Data Fig. 3. Temperature and absolute humidity recorded in Eurasian hibernacula.

Results are shown as raw data (dots) and as violin plots (coloured shading) showing the probability density estimate of the variables per clade. Temperature data were obtained from 152 and 26 sites in which clades Pd-1 and Pd-2 were sampled while absolute humidity was recorded in 91 and 22 sites per clade. There was no significant difference between the clades, either for temperature (two-sided t-test, t = 1.65, df = 35.52, p = 0.11) or absolute humidity (two-sided t-test, t = −0.56, df = 37.96, p = 0.58).

Extended Data Fig. 4. Bi-weekly size of cultures belonging to clades Pd-1 and Pd-2.

Measurement of culture sizes for 45 and 34 isolates of Pd-1 and Pd-2 respectively recorded for a growth period of 7 weeks (after which point growth ceases) at 15 °C. The size was measured from photos using R software (for more information see section ‘Analysis of growth’ in the Methods). The black line indicates the median, while the lower and upper hinges represent the first and third quartiles, respectively. The whiskers extend to 1.5 times the interquartile range. There was no significant difference in culture size of isolates belonging to Pd-1 compared to Pd-2 (two-sided t-test, p ranging from 0.07 to 0.97 for each week).

Extended Data Fig. 5. Density of multi-locus F_ST (a), π for Pd-1 (b) and π for Pd-2 (c) when mapping the 18 individually tagged isolates (11 Pd-1 and 7 Pd-2) on Gd293 reference genome (blue) or on Gd45 reference genome (orange).

To better visualize π values, values greater than 0.01 were omitted (representing less than 0.25% of values). 95% highest density intervals for multi-locus F_ST are 0.50—0.95 and 0.51—0.95 when using ref Gd293 and ref Gd45 respectively. For π, the means for Pd-1 are 4.2 × 10⁻⁴ and 7.1 × 10⁻⁴, and the 95% highest density intervals for Pd-1 are 6.5 × 10⁻⁵ — 8.3 × 10⁻⁴ and 0 — 0.9 × 10⁻⁴ when using ref Gd293 and ref Gd45 respectively. For Pd-2, the means are 4.6 × 10⁻⁴ and 7.3 × 10⁻⁴, and the 95% highest density intervals for π are 2.1 × 10⁻⁴ — 1.5 × 10⁻³ and 1.9 × 10⁻⁴ — 1.5 × 10⁻³ when using reference genome Gd293 and reference genome Gd45 respectively. See section ‘Diversity and differentiation calculation’ in the Methods, and Table S10 for further information on methodology and isolates, respectively.

Extended Data Fig. 6. Boxplot of the pairwise distance between isolates for 11 full genomes, partitioned between intra- and inter-clade distance.

Intra-clade Pd-1 & Pd-2 divergence are coloured in blue (Pd-1) and orange (Pd-2) while inter-clade divergence is coloured in green (as per Fig. 2c). The black line represents the median while the lower and upper hinges correspond to the first and third quartiles and the whiskers extend to 1.5 times the interquartile range with data points beyond this range shown as outlier points.

Extended Data Fig. 7. Density of multi-locus F_ST across the genome when mapping the Pool-seq data on Gd293 reference genome (blue) or on Gd45 reference genome (orange).

Note that the distributions are extremely similar independently of the reference genome. Indeed, the median multi-locus F_ST values are 0.88 (95% highest density interval [hdi], 0.50—0.99) and 0.88 (95% hdi, 0.49—0.97) when using reference genome Gd293 and reference genome Gd45 respectively.

Extended Data Fig. 8. Genomic location of within-clade recombination breakpoints.

Based on the four-gamete test using SNPs from clade Pd-1 (a, c; n = 11 isolates) and Pd-2 (b, d, n = 7 isolates) when using Gd293 as reference genome (a, b), or Gd45 (c, d). When using Gd293 reference genome: 264 out of 331 and 130 out of 215 windows with at least two SNPs show evidence of recombination in Pd-1 and Pd-2 respectively. When using Gd45 reference genome: 261 out of 309 and 117 out of 210 windows with at least two SNPs show evidence of recombination in Pd-1 and Pd-2 respectively. The Φw test of recombination significantly rejected clonality (p = 0.0) in all four instances (within each of the two Pd clades, whether considering Gd45 or Gd293 as reference genome). Contiguous regions alternate between blue and red at break points estimated by the four-gamete test. The population recombination rate (r = 2 Ne r; see section ‘Analyses of recombination’ in the Methods) was estimated at 1.0 × 10⁻⁴ and 4.6 × 10⁻⁵ for Pd-1 and Pd-2, respectively. The recombination rate was lower in Pd-2 than in Pd-1, confirming the result obtained with the FGT test.

Extended Data Fig. 9. Population differentiation in Pd-2.

a, Estimation of effective migration surfaces based on 107 isolates from Pd-2 in Europe (all sites excluding Russia after clone correction). For visualization, results from eight independent runs (each with 8 million iterations and between 100 and 450 demes), were combined. Different shades of a colour represent variable levels of high (blue) or low (brown) effective migration rates. Sampling locations are represented by black dots. b, Distribution of distance between true and assigned site of each isolate (n = 279) for the observed and randomized datasets (binwidth = 100 km). The mean distance of assignment was 42.88 km with a median of 0 km. In the Null-DAPC with randomization of sites before assignment, the mean and median were 913.68 km and 774.05 km, respectively.