Environmental reservoir dynamics predict global infection patterns and population impacts for the fungal disease white-nose syndrome.

Disease outbreaks and pathogen introductions can have significant effects on host populations, and the ability of pathogens to persist in the environment can exacerbate disease impacts by fueling sustained transmission, seasonal epidemics, and repeated spillover events. While theory suggests that the presence of an environmental reservoir increases the risk of host declines and threat of extinction, the influence of reservoir dynamics on transmission and population impacts remains poorly described. Here we show that the extent of the environmental reservoir explains broad patterns of host infection and the severity of disease impacts of a virulent pathogen. We examined reservoir and host infection dynamics and the resulting impacts of Pseudogymnoascus destructans, the fungal pathogen that causes white-nose syndrome, in 39 species of bats at 101 sites across the globe. Lower levels of pathogen in the environment consistently corresponded to delayed infection of hosts, fewer and less severe infections, and reduced population impacts. In contrast, an extensive and persistent environmental reservoir led to early and widespread infections and severe population declines. These results suggest that continental differences in the persistence or decay of P. destructans in the environment altered infection patterns in bats and influenced whether host populations were stable or experienced severe declines from this disease. Quantifying the impact of the environmental reservoir on disease dynamics can provide specific targets for reducing pathogen levels in the environment to prevent or control future epidemics.

Infectious diseases can drive population declines (1, 2), species extinctions (3–7), and the restructuring of ecological communities (3, 6, 8, 9). The presence of an abiotic pathogen reservoir is recognized as an important attribute of many disease systems, because it can amplify disease impacts by maintaining transmission independent of the affected hosts (10–15). Despite their importance, the effects of environmental reservoirs can be difficult to measure and, in most cases, remain unquantified (16), including for globally important diseases such as polio (17), plague (18), amphibian chytridiomycosis (19), and avian influenza (15). Although environmental sanitation practices are a hallmark of disease control (11, 20, 21), the influence of the environmental reservoir on outbreak size and population impacts remains unknown for many pathogens.

Here, we examine how disease outbreaks and population impacts are influenced by reservoir dynamics of a virulent pathogen by combining data on pathogen transmission and population growth rates from declining and stable communities across the globe. The recent introduction of the fungal pathogen Pseudogymnoascus destructans, which causes white-nose syndrome (WNS), from Eurasia to North America (22–24) has resulted in the death of millions of bats and led to >90% population declines across the range of at least 3 formerly abundant species in North America (4, 25, 26). P. destructans has been present in Eurasia for thousands of years (23), where it is known to infect bats, but has not been associated with mass mortality events (27, 28). Seasonal WNS epidemics begin each year when bats enter caves and mines to hibernate, where the fungus can survive in their absence (29–31). However, how the environmental reservoir influences population declines and whether reservoir dynamics differ between Eurasia and North America remain unknown.

We used standardized methods to quantify the P. destructans environmental reservoir on walls and ceilings in 101 hibernacula (subterranean sites where bats hibernate), as well as the corresponding infection dynamics in 39 species of bats in 8 countries on 3 continents over 8 y (21,341 total samples; ). We simultaneously collected population census data at each hibernaculum and used Bayesian hierarchical models to analyze changes in pathogen prevalence, fungal loads, and variation in population growth rates. Samples were collected from multiple time points during seasonal epidemics and data from North America included sites in multiple stages, including pre- and postpathogen invasion.

Results

Across Europe and Asia, bat populations were stable over the last 2 decades, but considerable variation was observed among sites and years (Fig. 1 and ; e.g., Myotis petax and Myotis nattereri had a mean annual population growth rate of 1.02 ± 0.19 and 1.07 ± 0.07, respectively). Similarly, in North America, prior to and during the first year of P. destructans invasion into sites, populations of Eptesicus fuscus, Myotis lucifugus, Myotis septentrionalis, and Perimyotis subflavus were stable or growing (Fig. 1 and ). However, in the 3 subsequent years following pathogen invasion, all 5 species in this region suffered severe population declines with mean annual decreases across all species ranging from 31 to 47% in years 2 to 4 (Fig. 1 and ).

Fig. 1.

Population growth rates for hibernating bats across Europe, Asia, and North America. Each dot shows an estimate of annual population change at a site for a species, with the size of the point scaled by log10 population size. The red line indicates population stability, and the black points and whiskers show the model predicted posterior mean, 95% credible intervals (thin lines), and ±1 SD of the posterior mean (thick lines) for each species across all regions (). (A) Bat population growth rates for species with >2 population growth estimates combined across Eurasia. (B) Bat population growth rates in North America (Left to Right) pre-P. destructans (Pd) invasion (all years prior to pathogen arrival), during P. destructans invasion (year 1), and following the invasion of P. destructans (years 2 to 4) (Fig. 2). Color corresponds to average declines across species for each year since P. destructans arrival (individual panels).

Fig. 2.

Changes in P. destructans prevalence on bats and in the environmental reservoir from preinvasion through P. destructans establishment (years 2 to 4) in North America. Columns in each panel show the invasion and establishment of P. destructans in the environmental reservoir (A) and on bats (B) from Left to Right, over a total of 4 y of infection with all years prior to the arrival of P. destructans shown in a single column (pre-Pd invasion). Line color indicates the different type of bat species (B) or substrate type collected (A), far (>2 m from any bat), near (10 cm from the focal bat), and under (directly under or adjacent to the bat). Points indicate a prevalence estimate for an environmental sample type (A) or species (B) at a site, and the lines indicate the posterior mean for each environmental sample type () or species ().

Differences in bat infection dynamics and WNS impacts on populations were strongly linked with differences in the environmental reservoir in winter. In North America during pathogen invasion, prevalence in the environmental reservoir and on bats was low in early winter but increased to moderately high levels by the end of winter [mean (95% credible intervals); environmental prevalence: 0.4% (0,1.0) to 17.0% (14.3,20.1); bat prevalence: 2.7% (1.4,4.0) to 53.0% (49.3,56.6); Fig. 2 and ]. Despite the high fungal prevalence on bats at the end of winter, there was no evidence of WNS impacts on populations during the first year of invasion (Figs. 1 and 2 and ), likely due to the long delay (70 to 100 d) between infection and mortality (22, 24). Large population impacts became evident in subsequent years when most species (Fig. 1 and ; e.g., M. lucifugus, M. septentrionalis, Myotis sodalis, and P. subflavus) became infected immediately upon returning to contaminated hibernacula where the environmental reservoir had established during previous outbreaks [Fig. 2; early winter bat and environmental prevalence years 2 to 4: 75.9% (72.6,77.2) and 41.1% (38.6,43.6); ].

Environmental contamination of hibernacula in North America steadily increased each winter after the initial pathogen invasion (Fig. 2 and ), and contamination levels remained stable or slightly increased over most summers following invasion [; average change over summer in prevalence (monthly slope): 0.03 (0.01,0.05)]. Pathogen dynamics across Europe and Asia were similar in the environment and on bats over the years sampled (Fig. 3, country panels and ). Prevalence in the environment at the start of winter was much lower in Eurasia than in established sites in North America (years 2 to 4), despite P. destructans likely being present in Eurasian hibernacula for thousands of years (23) [Eurasia: 19.0% (17.5,20.6); North America: 41.1% (38.6,43.6); Fig. 3, country panels and Fig. 4 and ]. Prevalence of P. destructans in the environment subsequently increased over the winter in Eurasia (Fig. 3, country panels and Fig. 4 and ), as was observed in North America (Figs. 2 and 4). However, in contrast to North America (Fig. 4 and ), P. destructans decreased in the environment in Eurasia over summer from 28.3% (26.6,29.9) to 19.0% (17.5,20.6), on average [Fig. 4 and ; monthly slope coefficients: China: −0.17 (−0.21,−0.13); Hungary: −0.19 (−0.26,−0.14); United Kingdom: −0.21 (−0.45,−0.04)]. Data from fixed marked locations in the environment sampled over time in the United States and China provided further support that P. destructans is stable and persistent in the environment in the United States, while both the presence and quantity of the fungus decayed in the environment over the summer in China ().

Fig. 3.

Global host and environmental reservoir dynamics of P. destructans over the winter (country panels). Solid lines indicate the mean of the posterior distribution for the prevalence of P. destructans for each species in a country sampled. Points represent prevalence estimates for a species or environmental sample at a site and the size of the point indicates the sample size. Red points on the map indicate the location of sample collection, and black lines link regions with corresponding graphs. Circles in the graphs show pathogen prevalence for a species or the environmental reservoir at a site, and lines show the mean of the posterior distribution for each species or environmental type in each region (). Dashed lines indicate predictions beyond the sampling dates. Prevalence of P. destructans in the environment is shown as the red line in each panel. Four-letter species codes correspond to the first 2 letters of the genus and species names for bats (). The data from North America show the changes in prevalence for the first year of P. destructans invasion at a site and all subsequent years following the invasion of the fungus (years 2 to 4 combined from Fig. 2) for comparison of dynamics across all regions.

Fig. 4.

Changes in P. destructans prevalence in the environmental reservoir over the summer (A) and on bats over the winter (B). Solid lines indicate the mean of the posterior distribution for the prevalence of P. destructans for each country sampled, and transparent ribbons indicate ±1 SD of the posterior distribution. Points represent prevalence estimates for a species or environmental sample at a site, and the size of the point indicates the sample size. (A) Dynamics of P. destructans in the environmental reservoir over the summer. Data for all years over the summer are combined into “US Established” and shown as individual summer periods (1 to 3) in . (B) The change in prevalence on bats over the winter period. Data for the United States are split into 2 groups with separate lines (invasion and established 2 to 4 y postpathogen invasion) and the species are combined for each region ().

Across Eurasia, prevalence on bats in early winter was consistently low for all species [bat prevalence: 21.1% (19.0,23.2); Fig. 3, country panels and Fig. 4] but increased over the hibernation period [late winter bat prevalence: 62.0% (60.0,64.0); Fig. 3, country panels and Fig. 4 and ]. Importantly, the change in prevalence of P. destructans on most Eurasian bat species across all years was similar to bats in North America during the first year of invasion (Fig. 4 and ), when environmental reservoir levels were equally low (Fig. 4 and ) and population growth rates were stable (Fig. 1). For example, increases in P. destructans prevalence over winter on a North American species, M. lucifugus, during the first year of pathogen invasion was similar to 2/3 of the species examined across Eurasia (e.g., M. petax populations in China and Mongolia and M. nattereri in Hungary and the United Kingdom; ). The fungal loads on M. lucifugus during the first year of invasion were also similar to 90% of the Myotis spp. examined across Eurasia (e.g., M. petax, Myotis daubentonii, Myotis gracilis, Myotis blythii, etc.; ), but differed from all species in North America, including other M. lucifugus during the established disease phase (years 2 to 4) (). Although winter length varied among some regions (), differences in the length of hibernation (combined effect of days below a mean daily temperature and latitude) could not explain the observed pathogen prevalence patterns ().

Species that roosted in areas with lower levels of P. destructans contamination in the environment in early winter had lower infection prevalence and fungal loads at the end of winter across all 3 continents [Fig. 5; late winter bat prevalence intercept: 0.45 (−0.43,1.39); early winter environmental prevalence slope: 1.91 (1.18,2.66); Fig. 5; late winter bat load intercept: −3.26 (−3.66,−2.81); early winter environmental prevalence slope: 1.33 (1.14,1.53); ]. Similarly, across these regions, bat populations roosting in areas with higher levels of P. destructans in the environmental reservoir during early winter had higher fungal loads (Fig. 5 and ) and lower annual population growth rates, suggesting that higher levels of P. destructans in the environment resulted in higher population impacts from WNS [Fig. 5; population growth rate (λ) intercept: 0.18 (−0.14,0.51); environmental prevalence slope: −1.47 (−2.01,−0.90); ]. There was support for this effect (decreasing population growth rate as environmental prevalence increased) in both North American (intercept: 0.02 [−0.24,0.30]; slope: −1.18 [−1.69,−0.70]) and Eurasian bat populations (intercept: 0.43 [−0.01,0.86]; slope: −2.13 [−3.46,−0.77] [country-specific results: China population intercept: 0.57 {−0.18,1.38}; slope: −2.01 {−3.74,−0.23}]), and the slopes were similar among these regions (). For example, 1 population of M. lucifugus that declined sharply in the United States (λ: 0.18 or population decline of 82%) roosted in areas with the highest levels of pathogen contamination (early environmental prevalence: 90 ± 0.13%) and had the highest fungal loads at the end of winter (late winter fungal loads: −0.07 log10 ng DNA; Fig. 5 and ). In contrast, a population of M. petax in China roosted in areas with lower pathogen contamination (early environmental prevalence: 33 ± 0.20%) and had low pathogen loads (late winter fungal loads: −2.95 log10 ng DNA) and nearly stable population growth rates (λ: 0.93).

Fig. 5.

Relationship between late winter P. destructans prevalence (A), fungal loads (B), and annual population growth rates (C) of bats and environmental reservoir prevalence in early winter. Black lines show the posterior mean and the gray ribbons show 95% credible intervals. The size of the points in each panel represents the sample size from bats (A and B) and the environment (C) for a species or the substrate sample at a site. The shape indicates the genus of the bat species for each point. (A and B) The relationship between late winter bat prevalence and early winter prevalence of P. destructans in the environment (within 10 cm of bats) [A; intercept: 0.45 (−0.43,1.39); slope: 1.91 (1.18,2.66)] and late winter fungal loads on bats and early winter prevalence of P. destructans in the environment [B; intercept: −3.26 (−3.66,−2.81); slope: 1.33 (1.14,1.53)]. Each point shows a prevalence estimate or mean fungal load for a species at a site and the contamination of P. destructans 10 cm around that species at a site. (C) The relationship between population growth rates for a bat species at a site and pathogen prevalence in the environment 10 cm around that species at a site in early winter [intercept: 0.18 (−0.14,0.51); slope: −1.47 (−2.01,−0.90)]. Each dot shows the prevalence of P. destructans in the environment 10 cm from around individuals of that species at a site in early winter (x axis) and the y axis shows the population growth rate (N/N) from the previous late winter, N, to the following late winter, N. The horizontal red line indicates population stability. bR2 are Bayesian R2 values, calculated as the variance of the predicted values divided by the variance of predicted values plus the variance of the errors.

Discussion

White-nose syndrome has caused declines in bat populations between 50 and 90% for species in North America, while many bat species across Eurasia are stable (25, 32). Our results suggest that differences in environmental reservoir dynamics account for the large contrast in population growth rates of bats in both North America and Eurasia. Data from both native regions where bat populations have coexisted with P. destructans for thousands of years, and invading disease regions where populations have suffered extensive mortality, reveal that a crucial difference between stable and declining bat populations is the extent of the environmental reservoir in early winter. The seasonal decay of P. destructans in the environment over the summer in Eurasia (Fig. 4 and ) leads to a reduced reservoir of P. destructans in early winter, delayed infection, a shorter period of pathogen growth on bats, and lower fungal loads at the end of winter (Figs. 4 and 5 and ). The long period (70 to 100 d) between infection and mortality for WNS (22, 24, 33–36) means that while most bats eventually become infected by the end of winter across Eurasia, they survive until spring, when they can emerge from hibernation and clear infection (29). In contrast, bats in WNS-established areas of North America become rapidly infected when they return to an extensive and heavily contaminated environmental reservoir in early winter (Fig. 4 and ) and, given these longer periods of infection, have higher fungal burdens () and greatly increased mortality (Fig. 1).

We explored whether environmental conditions and host colony size could explain the seasonal decay or persistence of P. destructans in the environment. Comparisons of hibernaculum temperatures between North America and Eurasia showed that only 3 of the 7 regions varied from the United States (Mongolia and the United Kingdom were lower and Israel was higher; ) and were unable to account for the continental differences in pathogen decay that were observed (). In addition, the higher prevalence of P. destructans in the environment in North American established sites was not driven by higher densities of bats, as we found little difference between colony size among regions () and found no support for a relationship between colony size and pathogen contamination in the environment at the start of the next winter ( and results in the figure legend).

In Eurasia, the observed pathogen reduction over summer in the environment may involve species that compete with or consume P. destructans. Many invasive plants and animals have fewer natural enemies in their introduced range, which allows them to proliferate unconstrained by predation and competition (37, 38), and this may have occurred with P. destructans in North America. Our results indicate that if environmental prevalence could be reduced below ∼20% each summer in North America by active control measures or natural competition and predation, it could result in stable bat populations on average (Fig. 5), similar to those observed in Eurasia.

Understanding environmentally mediated transmission and its contribution to disease outbreaks is critical for reducing epidemic potential (21, 39). We have shown that the extent of the environmental reservoir not only influences pathogen prevalence and loads on bat hosts but corresponds to the degree of population declines. Control efforts for pathogens with environmental reservoirs often aim to reduce indirect transmission by managing pathogen levels in the environment (11, 21, 40) but the relative importance of the environmental reservoir is rarely quantified, leading to a lack of defined targets. Characterizing the environmental reservoir across space and time will provide more accurate and targeted control efforts to break chains of transmission and prevent disease outbreaks.

Materials and Methods

Data Collection.

We quantified bat population growth rates and the presence of P. destructans on bats and in the environment at 101 sites in 8 countries on 3 continents across the Northern Hemisphere. We visited sites twice per winter over a total of 8 y (2012 to 2019) to measure changes in infection prevalence in the United States (8 winters; 2012 to 2019), the United Kingdom (2 winters; 2016 to 2018), Hungary (2 winters; 2016 to 2018), Japan (1.5 winters; 2017 to 2018), Mongolia (1.5 winters; 2016 to 2017), and China (3.5 winters; 2015 to 2018). In addition, we collected samples from a single time point in late winter from Israel (2016) and Georgia (2017).

At each site, we estimated population size by counting all individuals of each species. Hibernating bats primarily roost on cave and mine surfaces, making near-complete counts of individuals possible. We collected epidermal and environmental swabs during early and late winter. We targeted winter sampling because WNS epidemics only occur during winter when bats enter caves and mines to hibernate (29). Bats become infected with the fungus in the fall when they return to hibernacula where the pathogen can persist in the absence of bats (29–31), and mortality from WNS occurs in mid to late winter, ∼70 to 100 d after initial infection (22, 24, 29, 34). If bats survive until spring, they generally clear infections when they emerge from hibernation, and remain uninfected until the following fall (29, 41, 42).

Sample Collection.

Samples were analyzed using a previously described swabbing protocol (27, 29) and were standardized across all regions. We collected a temperature measurement for each bat sampled using a Fluke 62 MAX IR thermometer (Fluke) taken directly adjacent to the bat. To sample for P. destructans on bats, we dipped polyester swabs (Puritan) in sterile water and swabbed the wing and muzzle 5 times back and forth (43). We collected environmental pathogen samples by swabbing a section of substrate that was similar to the length of a bat’s forearm (36 to 40 mm) using the same methods, but without dipping the swab in sterile water. Three types of environmental samples were collected: under or directly adjacent to hibernating bats (“under”), ∼10 cm from a hibernating bat (“near”), and >2 m from any hibernating bat but in locations where bats might roost (“far”). In addition to these environmental samples, we also collected samples from fixed marked locations over time at a subset of sites in the United States and China (4 sites in China and 5 sites in the United States). At sites in China and the United States, we installed between 15 and 20 reflective markers with unique identification numbers on the substrate in each site. Samples were collected radiating out from the marker (36 to 40 mm), consecutively in a clockwise fashion to avoid removing P. destructans in an area that would be sampled at a later time. These stations were installed in areas where bats would normally roost and were sampled during consecutive late and early winter sampling visits. All samples were stored in salt preservation buffer (RNAlater; Thermo Fisher Scientific) immediately following collection to stabilize and preserve DNA.

Sample Testing.

We extracted DNA from samples using a modified Qiagen DNeasy Blood & Tissue Kit (29) and tested for the presence of P. destructans using qPCR (44). All samples were run in duplicate, with 13 negative controls and quantification standards on each plate. All quantification standards were within a consistent range and all negative controls had no fungal detection.

Permits.

All research was conducted under Institutional Animal Care and Use Committee protocols: Virginia Polytechnic Institute: 17-180; University of California, Santa Cruz: Kilpm1705 and FrickW1106; Wisconsin Endangered/Threatened Species Permit 882-886; US Fish and Wildlife Service Threatened & Endangered Species Permit TE64081B-1; Natural England Project Licenses 2016-26874-SCI-SCI and 2017-32554-SCI-SCI in the United Kingdom (England); Israel Nature and Parks Authority (2016/41242); Wildlife Science and Conservation Center of Mongolia by the Ministry of Environment and Tourism of Mongolia; National Inspectorate for Environment and Nature Hungary (14/2138-7/2011; PE-KTFO/4384-24/2018); Japan Ministry of Environment (1603013, 1703084) and prefectural governments (Shizuoka: 27-21-1 and 28-15-1; Kochi: 1; Iwate: 1642 and 2494).

Analyses.

We present analytical methods in identical order as they appear in the main text and figures. In addition, an appendix for analyses, figures, and all results is provided (), which indicates where the results can be found (e.g., supplemental table, figure legend, or article text), whether the results are displayed in a figure (and if so, the figure number), the parameters of the model, the data that were used, and the main findings.

We examined the population growth rate (λ) between pairs of late winter counts (N) separated by T years, [λ = (N/N)1/], of bat species across Europe, Asia, and North America. Population counts were collected at each site during sampling for P. destructans, and some sites were also counted in previous years as part of ongoing bat population monitoring. We included counts from the years 2000 to 2019 for all sites where data were available. We fit a Bayesian hierarchical model using population growth rate as our response variable with a gamma distribution and log link. Species was included as a population-level effect (equivalent to a fixed effect in a frequentist framework) and site as a group effect (equivalent to a random effect in a frequentist framework). For North America, we included an interaction for the stage of pathogen invasion (e.g., number of years with P. destructans) to examine differences in population trends before and after disease invasion (Fig. 1). We fit all models (unless otherwise noted) using the No-U-Turn Sampler (NUTS), an extension of Hamiltonian Markov chain Monte Carlo. All Bayesian models were created in the Stan computational framework (https://mc-stan.org/) accessed with the “brms” package in program R (45). To improve convergence and avoid overfitting, we specified weakly informative priors (a normal distribution with a mean of 0 and SD of 10). Models were run with a total of 4 chains for 2,000 iterations each, with a burn-in period of 1,000 iterations per chain resulting in 4,000 posterior samples, which, given the more efficient NUTS, was sufficient to achieve adequate mixing and convergence. All values were less than or equal to 1.01, indicating model convergence.

We analyzed changes in P. destructans prevalence on bats and the environment across North America and Eurasia by fitting separate Bayesian hierarchical models with binomial distributions and logit links to each region. For all analyses examining pathogen prevalence, the detection of the fungus, P. destructans, in each sample was our Bernoulli (0/1) response variable. In North America, each sample type (substrate in Fig. 2 or bat in Fig. 2) was analyzed separately, and we included a 3-way interaction between bat species or environmental sample type (samples collected under, 10 cm from, or >2 m from bats), a scaled winter date of sample collection, and years since P. destructans detection. We treated year as categorical for bat prevalence models and as continuous for substrate prevalence models, because bats clear infections between years whereas fungal dynamics were continuous in the environment. All models included site as a group-level effect.

We separated the data from North America over the summer to examine the individual slopes and intercept for each year following the initial invasion of the fungus into sites. We analyzed the change in P. destructans prevalence in the environment over the summer for each summer since P. destructans arrival using the same methods as above, but we included a 2-way interaction with Julian day and summer since invasion with site and substrate sample type as a group-level effect (). In addition to the change in environmental prevalence over the summer, we also examined the change in fungal loads over the summer for different countries (United States, United Kingdom, Hungary, and China). We included fungal loads (log10 ng of DNA) as our response variable and a 2-way interaction between season (late winter year t and early winter year t + 1) and country with sites as a group effect ().

To further explore P. destructans dynamics in the environmental reservoir over the summer, we analyzed the change of P. destructans prevalence and fungal loads at fixed marked locations in the environment in both the United States and China. Our response variable was the detection of P. destructans in a sample (0|1) and we included a 2-way interaction between sampling date and country with sampling station ID nested within site as a group-level effect. For the changes in fungal loads, we ran the identical model structure described above with a Gaussian distribution and fungal loads as our response variable ( and results in the figure legend).

To examine changes in P. destructans prevalence in bats and the environment over the winter for each country across Eurasia, we ran the same models described above with a 3-way interaction between winter sampling date, species, and winter year (the winter when the data were collected, e.g., 2015/2016, 2016/2017, etc. and equivalent to the years since P. destructans detection in models described above) with site as a group-level effect. We included substrate in the analysis as a “species” and combined under and near substrate sample types. For countries where there were fewer than 2 y of data, we excluded winter year from the analysis (Mongolia and Japan). We included winter year as a categorical population-level effect for China, Hungary, and the United Kingdom so that variation in species between winter periods could be visualized. For countries where data were collected during a single time period during the winter (Israel and Georgia), we excluded sampling date and winter year from the analyses (Fig. 3, country panels and ).

To compare patterns of prevalence in the environment over the summer among continents, we used a model with Julian date interacting with country of sample collection, and site and environmental sample type as group effects (Fig. 4 and ). In addition, we analyzed the change in P. destructans prevalence on bats over the winter using the same methods as above, but we included a 2-way interaction with winter day and country with site and species as group-level effects. North America was broken down into 2 different stages, first year of P. destructans invasion or all years following invasion (Fig. 4 and ).

We analyzed the difference in fungal loads on bats during late winter by first fitting a model with fungal loads on bats as our response variable, and sampling date interacting with species sampled in a country as a population effect and with site as a group effect. The US data were grouped by the first year of pathogen invasion (year 1) and after P. destructans establishment (years 2 to 4). The model with date was fit to account for differences in sampling date among the regions (). We then extracted the posterior mean, SD, and 95% credible intervals for each species in a country on March 1, shown in . The pairwise comparison for each species combination was then calculated to determine differences among species at the end of winter (). An estimate of fungal loads could not be obtained for species with no positive detections, and positives at only a single time point or fewer than 9 samples over winter were dropped from the analysis.

To examine the differences in hibernation length, we retrieved data from the National Oceanic and Atmospheric Administration’s National Climatic Data Center. Mean minimum daily temperatures from the swab sampling study period (2013 to 2018) were retrieved from multiple climate stations within 200 km of each site in each country where samples were collected. We fit a Bayesian generalized additive model to the data for each country and year, with date predicting mean daily temperature to fit a smoothed model of mean minimum daily temperatures. We then calculated the days below 0 °C and 5 °C for each country and year based on the model fit to quantify approximate hibernation length for bats in the different regions sampled (). To determine the relationship between climate (a proxy for days in hibernation) and prevalence of P. destructans in early winter on bats, we fit a model with the P. destructans detection as our response variable and additive effects of latitude and days under 5 °C for a region as our predictor variables using a Gaussian distribution and logit link, and we included site as a group-level effect ( and results in the figure legend).

We examined the relationship between early winter prevalence of P. destructans in the environmental reservoir 10 cm from a species of bat at a site and late winter prevalence or fungal loads. Early winter P. destructans prevalence and fungal loads for a species were calculated using samples collected before December 31 (bats across most regions typically begin hibernation between October and December) and late winter bat loads and prevalence were calculated using samples collected after February 1. In addition, we examined the relationship between early winter prevalence of P. destructans in the environmental reservoir near a species of bat at a site and early winter bat prevalence. We modeled these relationships by fitting a Bayesian hierarchical model with prevalence or loads of P. destructans on bats in late winter as our response variable and environmental reservoir prevalence as the predictor variable using a binomial or Gaussian distribution, respectively. We fit a Bayesian hierarchical model with prevalence and loads of P. destructans on bats in late winter (Fig. 5 ) and in early winter ( and results in the figure legend) as our response variables and early winter environmental reservoir prevalence as the predictor variable using a binomial distribution and logit link, and we included site, species, and country as group-level effects. We also examined the relationship between bat population growth rate (as calculated above) and early winter prevalence in the environment (combining under and near substrate samples) for a given species at a site. We fit a Bayesian hierarchical model with a gamma distribution and log link, and included site and species as group-level effects (Fig. 5 and results in the figure legend). In addition, we ran the same set of 3 models described above, but with fungal loads in the environment during early winter as our predictor and against the same response variables in . The same model described above examining the relationship between bat population growth rate (as calculated above) and early winter prevalence in the environment was repeated but with a population-level effect for the continent (North America or Eurasia) to examine differences in this relationship between regions ( and results in the figure legend).

We examined differences in hibernaculum temperatures in relation to the change in P. destructans over summer at all sites. First, we estimated differences in winter roosting temperatures among countries, and fit a hierarchical model using the temperature measurement from each bat at each site as our response variable and country as a population-level effect, with site and species as group-level effects to determine whether hibernation temperatures differed among countries (). Then, to quantitatively assess the relationship between roosting temperatures within a site and the decay or persistence of P. destructans over the summer in that site, we fit a Bayesian hierarchical model with the change in P. destructans prevalence over the summer in samples taken >2 m from a bat [PNovember − PMarch/(1 − PMarch), where P is prevalence] as our response variable and mean roosting temperature at that site as the predictor variable using a Gaussian distribution and logit link, and we included site as a group-level effect ( and results in the figure legend).

To determine whether there were differences in colony sizes of bats among regions, we fit a hierarchical model with total colony size of a site as our response variable and country as our population-level effect. Data for the United States were broken down into 2 subcategories (prepathogen, establishment and established) and site was included as a group effect (). We also examined whether the total size of a bat colony predicted early substrate prevalence by fitting a model using early substrate prevalence 10 cm around a species at t + 1 as our response variable and colony size for a species at a site in late winter (t) as our population effect with site as a group effect ( and results in the figure legend). All statistical tests were carried out using R 3.3.2.

Data Availability.

Raw data points are contained within the figures and data files when possible. The data and code, not accessible in the manuscript, have been made available at GitHub, https://github.com/hoytjosephr/global-wns.