Research ArticleANIMAL ECOLOGY

Environment, host, and fungal traits predict continental-scale white-nose syndrome in bats

See allHide authors and affiliations

Science Advances  29 Jan 2016:
Vol. 2, no. 1, e1500831
DOI: 10.1126/sciadv.1500831


White-nose syndrome is a fungal disease killing bats in eastern North America, but disease is not seen in European bats and is less severe in some North American species. We show that how bats use energy during hibernation and fungal growth rates under different environmental conditions can explain how some bats are able to survive winter with infection and others are not. Our study shows how simple but nonlinear interactions between fungal growth and bat energetics result in decreased survival times at more humid hibernation sites; however, differences between species such as body size and metabolic rates determine the impact of fungal infection on bat survival, allowing European bat species to survive, whereas North American species can experience dramatic decline.

  • Emerging infectious disease
  • energetic model
  • Extinction
  • Eptesicus fuscus
  • Eptesicus serotinus
  • fungal pathogen
  • Geomyces destructans
  • hibernation
  • Myotis lucifugus
  • Myotis myotis
  • Pseudogymnoascus destructans


Central questions in disease and evolutionary ecology relate to how interactions among hosts, their parasites, and the environment affect host-parasite dynamics and why pathogenicity may differ among individuals and species. Wildlife diseases can have major impacts on human, animal, and ecosystem health. Epizootics of rinderpest affected African ecosystems for decades (1), and chytridiomycosis and ranavirus have caused worldwide decline of amphibian populations and species extinctions (2, 3). Of similar magnitude and potential for long-term impacts, white-nose syndrome (WNS) is one of the most rapidly spreading wildlife diseases ever recorded (Fig. 1) (4, 5). WNS is caused by the psychrophilic (cold-growing) fungus Pseudogymnoascus destructans (Pd; previously Geomyces destructans) (68), which grows at approximately 0° to 19.7°C (9) and invades the skin tissues of hibernating bats (10, 11). European bats infected with Pd can hibernate without obvious mortality (12); however, since the discovery of WNS in North America during 2007, WNS has been diagnosed in seven species spanning 26 U.S. states and 5 Canadian provinces (Fig. 1) and has killed millions of bats. Genetic evidence suggests that Pd was introduced to North America from Europe (7, 1315).

Fig. 1 Spread of WNS over eight winters.

The annual spread of the fungus that causes WNS, Pd, in North America. Shaded land units represent counties in the United States and voting subdistricts in Canada. Red counties show where the disease is believed to have started during the winter of 2006–2007.

Hibernating bats spend the majority of winter in deep torpor (body temperature <10°C) with intermittent arousals to euthermia (35° to 38°C) (16). Euthermic arousals at low ambient temperature consume the majority of overwinter energy reserves, and bats require specific and narrow ranges of hibernaculum temperatures to survive winter on their limited energy budgets (16, 17). Studies show that WNS mortality begins approximately 3 months into hibernation, and diseased bats typically become emaciated by late winter as arousal frequency from torpor increases (8, 18). Thus, a common hypothesis for why bats die from WNS is that Pd infection causes them to arouse with increasing frequency and deplete their required fat reserves (8, 18, 19). However, support for how this hypothesis plays out in nature is largely circumstantial, and a mechanistic understanding that would allow prediction of how the WNS epidemic will spread and which species will be most affected is lacking. Few studies explicitly integrate feedback on bat arousal from temperature-dependent fungal growth (20), despite clear patterns in how temperature influences Pd growth (9) and how Pd infection increases arousal frequency (8, 18). Furthermore, empirical evidence indicates that relative humidity (RH) influences WNS population declines (21), yet humidity-dependent fungal growth has not successfully been integrated into bat survival models. Species of bats in the northeastern United States showing greatest susceptibility to WNS are known to consistently use the wettest hibernacula (10), and conidial fungi are more likely to germinate on wet surfaces (22), suggesting that Pd may be more infectious in humid sites. Humidity has been integrated into bat hibernation models (20, 23), but only as a factor influencing evaporative water loss (EWL) in the bat hosts and not accounting for potentially nonlinear interactions among temperature, RH, and growth of the fungus in bat skin.

We model the growth dynamics of Pd and energetic requirements of WNS-affected hibernating bats under a range of environmental conditions. Populations of the little brown bat (Myotis lucifugus) in the northeastern United States and Canada have been more affected by WNS than any other species (5). We compare model outcomes from M. lucifugus to another species less affected by WNS in North America, the big brown bat (Eptesicus fuscus), and to two apparently unaffected European species, the serotine (Eptesicus serotinus) and greater mouse-eared bats (Myotis myotis). We model Pd growth as a function of body temperature and RH and incorporate this into an energetic model across a range of ambient temperatures. We predict mortality times, based on critical depletion of stored body fat, for specific combinations of environmental conditions. We then use climate data from within the distributions of each species to predict survival times of each species compared to spatially varying winter durations within its range, with and without fungal infection.


Our modified energetic model (16) is generalizable across species that prevent their body temperatures from dropping below (defend) different temperatures (Ttor-min) during overwinter hibernation. Our generalization ensures that if the temperature is at the optimal torpor temperature (Ttor-min), then energy use is minimized [minimum torpor metabolic rate (TMRmin)] so that time in torpor (ttor) is maximized, whereas if the temperature is away from Ttor-min, ttor is reduced such that it is equal for given TMR above and below Ttor-min (fig. S1). Our results demonstrate that bats are energetically constrained to narrow ranges of cold temperatures for surviving long winters (for example, 6 months) even in the absence of fungal infection (fig. S2). We additionally model fungal growth dynamics and incorporate these dynamics into the energetic model by assuming that Pd infection affects bats by increasing the consumption of stored fat with increasing frequency of arousal from deep torpor.

Higher arousal frequency increases overwinter mortality in the environmental conditions within which Pd grows well, especially warmer and wetter hibernacula (Fig. 2). The model predicts that hibernating M. lucifugus infected with Pd can make energy reserves last 6 months at both ambient temperatures (Ta) between 1° and 6°C and <98% RH. However, bats that defend different minimum temperatures (Ttor-min) have different optimal torpor temperatures, and this and other host traits affect the duration of winters they can survive (Fig. 2).

Fig. 2 Predicted times to deplete overwinter energy reserves for bats infected with Pd.

(A to D) Model predictions in months (surface colors and contours) are shown for a range of RH percentages (88 to 99%) and ambient hibernacula temperatures (0° to 19.4°C) at which Pd grows for M. lucifugus (A), E. fuscus (B), M. myotis (C), and E. serotinus (D). The arrow in (A) shows 3.3 months at 7°C and 97% RH for M. lucifugus.

Support for our model predictions includes independent experimental observations, in which M. lucifugus experimentally infected with European and North American isolates of Pd died within 70 to 120 days at 7°C and 97% RH (8). Under those conditions, our model predicts that M. lucifugus can survive 100 days, matching the experimental results (Fig. 2). The model predicts survival to 185 days (6 months) under these same conditions without Pd infection (fig. S2). Additional evidence of model validity includes predicted surface areas infected (fig. S3) and maps of predicted overwinter survival across the distributions of affected species (Fig. 3) that qualitatively match observed disease patterns, as well as observations of greater mortality in hibernacula with warmer temperatures (21).

Fig. 3 Comparison of winter fat depletion in North American bats.

(A to F) Difference between winter duration and predicted time to deplete overwinter energy reserves for M. lucifugus (A to C) and E. fuscus (D to F) within their distributions. Differences in months are shown before (A and D) and after (B and E) the arrival of the fungus Pd, the cause of WNS. Blue indicates that bats are predicted to have more than enough energy reserves to survive a typical winter (+ values), with white (no difference, all energy reserves used) and red (− values) indicating that bats are unable to survive winters with enough energy reserves to survive through hibernation. Dark gray indicates that energy reserves are greater than an 8-month upper limit in (D). Distributions of the difference between winter duration and the model’s predicted overwinter survival times are in the histograms (C and F), with median values (dashed lines) for pre-Pd (blue) and post-Pd (pink) infection and zero difference (black line) shown.

The modeled mechanism of feedback between increased fungal growth, decreased torpor time, and bat energetics driven by temperature and RH seems reasonably supported for M. lucifugus, but the robustness of this mechanism is further supported across multiple species. For example, our model predicts that E. fuscus can hibernate for >6 months in the absence of Pd infection (fig. S2) and over a larger range of environmental conditions than M. lucifugus. This trend is reflected in the extended survival times for E. fuscus when Pd infection is incorporated into the model (Fig. 2). We find that European M. myotis is limited to a smaller environmental parameter space than the other species analyzed (Fig. 2). However, when we compare predicted potential hibernation duration to estimated winter duration (hereafter “survival capacity”) of each of the four species across their distributional ranges using climate data (Figs. 3 and 4), results indicate that although M. myotis has a lower median survival capacity (+3.2 months) than E. serotinus (+4.3 months), the median survival capacities of both infected European species are higher than those predicted for both infected North American species (M. lucifugus, −1.2 months; E. fuscus, +2.2 months). Modeled survival capacities of M. lucifugus and E. fuscus without Pd infection were +1.5 and +4.9 months, respectively (Fig. 3), again indicating substantially higher overwinter survival in the absence of the fungus. Bootstrapped Harrell-Davis quantiles of predicted survival capacity times are statistically higher for E. fuscus than M. lucifugus (+0.6 month versus −3.2 months, P < 0.001; +2.1 months versus −1.0 month, P < 0.001; +3.3 months versus +0.9 month, P < 0.001).

Fig. 4 Comparison of winter fat depletion in European bats.

(A to D) Difference between winter duration and predicted time to deplete overwinter energy reserves for M. myotis (A and C) and E. serotinus (B and D) within their distributions. Results in months are shown with the fungus Pd, the cause of WNS in American bats. Blue indicates that bats are predicted to have more than enough energy reserves to survive a typical winter (+ values), with white (no difference, all energy reserves used) and red (− values) indicating that bats are unable to survive winters with enough energy reserves to survive through hibernation. Distributions of the difference between winter duration and model’s predicted overwinter survival times are in the histograms (C and D), with median values (dashed red lines) and zero difference (black line) shown.

We assumed that European bat parameters were measured from populations with Pd infection. Given this assumption, our model suggests that behavioral and/or physiological traits may have evolved or been preadapted in the European species to increase survival with Pd infection, whereas in M. lucifugus and E. fuscus, such traits have not evolved. Predictions suggest that E. fuscus may be better suited to survive WNS infection than M. lucifugus. Our findings suggest that environmental conditions and basic host traits alone may explain much of the variability in disease outcomes among species of bats infected by Pd in North America and Europe.

Our results show how incorporating temperature- and humidity-dependent fungal growth into models to predict bat hibernation energetics, and thus overwinter survival, is critical because these key features of fungal ecology are inextricably linked to bat mortality from WNS. Sensitivity analysis assessing that the proportion of examined parameter space survival is predicted over after 6-month winters (16) reveals that survival is significantly affected by numerous parameters, but that our model is most responsive to fungal growth changes related to RH (Fig. 5).

Fig. 5 Sensitivity analysis of model parameters.

Sensitivity analysis results using partial rank correlation coefficients (PRCCs) for the output variable predicted mortality after 6 months for a range of RH percentages (88 to 99%) and ambient hibernacula temperatures (0° to 19.4°C) for M. lucifugus. Positive PRCCs indicate that increasing a parameter increases mortality. Parameter definitions and PRCC values are in Table 2. Bat parameters are light gray, temperature-related fungal growth parameters are mid-gray, and humidity-related fungal growth parameters are dark gray. Significance at α = 0.05 is demarcated by the dashed line.

The importance of RH to Pd growth is plausible, because conidial fungi such as Pd are more likely to germinate and degrade nutrient substrates in the presence of high moisture levels in their environments (22). These observations and results suggest that individuals, colonies, or species of bats that use microhabitats with lower RH will be less susceptible to WNS. In most caves, RH reaches 100% far from entrances but can vary throughout, and RH is affected by Ta, airflow, and atmospheric pressure (24). Typically, bats use hibernation sites with 90 to 100% RH (25); however, three North American species that seem less severely affected by WNS—Myotis sodalis, Myotis leibii, and E. fuscus—tend to select drier areas within hibernacula, whereas the three species most dramatically affected by WNS—M. lucifugus, Myotis septentrionalis, and Perimyotis subflavus—consistently roost in the most humid locations within hibernacula and are regularly observed with condensation on their fur (10, 21). Changes in RH are predicted to affect bat survival independent of Pd infection [for example, Cryan et al. (10) and Willis et al. (23)]. We have not incorporated RH into the bat energetic model to allow for a parsimonious modeling approach; however, we predict that increased RH will likely improve the survival of M. lucifugus in the absence of Pd by reducing EWL (26, 27), whereas our model predicts that increasing humidity in hibernacula with Pd generally decreases bat survival through WNS, a result supported by field data (21). The single most affected species, M. lucifugus, is in positive water balance at 2° and 4°C only at ≥99% RH (25), which is the ideal growing conditions for Pd. Monitoring RH and understanding the interactions between water vapor pressure and EWL (26, 27) in bat hibernacula may be a key element of future disease surveillance and research efforts, because this environmental variable might explain a considerable amount of the variability in Pd pathogenicity observed among species, hibernation environments, and continents.

Sensitivity analyses revealed that host traits had smaller effects on survival than changes in humidity-dependent growth parameters (Table 2 and Fig. 5). However, torpor is already constrained thermally (fig. S2) (16), and host traits and species distributions interact to explain species-specific differences observed in survival times (Figs. 2 and 3). Decreasing time in torpor and increasing euthermic times also increased mortality in the sensitivity analysis as expected. With respect to host traits, survival was most sensitive to changes in body mass (BM), with increasing mass decreasing mortality (possibly due to increased fat available and decreased thermal conductance), followed by the lower temperature at and above which bats are thermally neutral (Ttor-min) and minimal resting metabolic rate (RMR) occurs, Tlc. However, there was little natural variation in Tlc [coefficient of variation (CV), 0.01] among species (Table 1), and greater variation in other parameters that significantly affected the model results (Fig. 5) [for example, TMRmin (CV, 1.21) and Ttor-min (CV, 0.66)] suggests that these may be higher priority for study because they might provide greater predictive ability regarding differences in survival between species. Attempts to predict the necessary host trait parameters (for example, Tlc and TMRmin) from more frequently measured parameters (for example, BM and Ttor-min) were unsuccessful (r2 = 0.01 to 0.84), highlighting a more general need to study bats further.

Table 1 Winter energy expenditure model and parameters related to fungal growth.

Data are given for little brown (M. lucifugus), big brown (E. fuscus), greater mouse eared (M. myotis), and serotine (E. serotinus) bats with the calculated CV for each parameter that was varied.

View this table:
Table 2 Sensitivity analysis results for the E. fuscus–Pd model.

Positive PRCCs indicate that increasing a parameter increases mortality. Default parameter values are in Table 1 and were sampled from uniform distributions from the default values to a minimum of 10% lower than default values.

View this table:

Our model results predict that those species more likely to survive Pd in North America will have a combination of physical and behavioral traits, including larger body sizes and hibernation in drier sites and/or colder sites. The results suggest that within populations, there will be strong selection pressure for these traits.

Our results indicate that we can predict bat mortality using a model that incorporates bat traits, fungal growth components, and environmental conditions. We recognize, however, that other factors can play a role in the disease process. Increased metabolic rates may occur through alternative mechanisms in the absence of increased arousal frequency (28). In addition to energy expenditure, other mechanisms that may lead to decreased survival of WNS-affected bats include altered physiological processes during winter, such as fungal damage to wing membranes potentially disrupting blood circulation, water and electrolyte balance, or immune function (10, 2932). All or some of these processes may influence arousals from torpor, but we were able to explain and predict bat mortality using a parsimonious host energetic model, and some of these additional aspects may be captured phenomenologically in our model, because Pd growth is what influences arousal (8, 18). The overwhelming effect of increased arousals depleting energy reserves is persuasive given that one arousal bout of M. lucifugus hibernating at 5°C consumes the same amount of fat energy as 67 days spent in torpor (17).

Previous studies have modeled the temporal dynamics of WNS as it spreads (33). We assume that because Pd is a saprobic fungus, it can persist in the underground environments inhabited by overwintering bats. Therefore, we did not model transmission dynamics within the bat population itself. Host density and social behaviors likely influence the transmission dynamics of WNS (21) and could be incorporated into modeling frameworks along with host-specific fungal growth as more data become available (34). Our mechanistic model may provide a framework for potentially integrating species distribution modeling (16), climate change, and the impacts of WNS (for example, fig. S4). However, such efforts may require a better understanding of how to scale macroscale environmental changes, such as climate change, to microscale environmental change within caves—our analyses suggest that simple scaling (for example, fig. S4) may not be representative (35, 36).

The results of our model paint a bleak picture for the American bat species M. lucifugus, as predicted increased arousal frequency across most of their distribution suggests that they will struggle to survive (Fig. 3). In conclusion, our model results allow us to understand the interactions between host and pathogen at the individual level, and their interaction with the environment, as well as allow us to scale these up to regional levels to understand and predict widespread species survival or decline. These results are important steps toward integrating knowledge gained through infection traits, host biology, and climate to increase understanding of WNS disease outcomes in bats, and they demonstrate the value of modeling dynamics of the “epidemiologic triad” (environment, pathogen, and host) to predict species-specific differences in infection outcome.


Fungal growth model

We assume that Pd infection occurs at the beginning of the hibernation period, and we modeled Pd growth as a function of bat body temperature, Tb, and ambient RH through a series of torpor-arousal bouts during winter hibernation.

Two independent studies report the growth of Pd over time for different temperatures in culture (9, 37). We used one that cultured the fungus across more temperatures and reported the thermal performance function results to obtain the temperature-dependent function shape (9). To obtain temperature-dependent Pd growth rates, we first simulated the expected growth (Gweek) of Pd over 5 weeks for a range of temperatures (9). We simulated the growth using the thermal performance equation (Eq. 1) (9)Embedded Image(1)to obtain the predicted temperature-dependent Pd growth at each week for 5 weeks using reported estimates of (α1 to α4) (fig. S5A) (9).

We used the results of these simulations to calculate the relationship between Pd growth rates and temperature through time. The increase in colony area was a linear function of time (fig. S5B), so we used linear model coefficients to estimate the rate of change in Pd colony area per week. Calculated maximum Pd growth from the two temperatures in Chaturvedi et al.’s study (37) was fourfold greater than that in Verant et al.’s study (9) (fig. S5C), and we increased the maximum Pd growth under optimal conditions to match Chaturvedi et al.’s study. To use a single function for Pd temperature-dependent growth in the model, we then fit a nonlinear function to the different temperature-dependent linear model coefficients to model temperature-dependent growth (GT) with the functional formEmbedded Image(2)where T is temperature, Tmin is the minimum temperature at which the fungus grows (fixed at 0°C) (9, 37), β1 and β2 are shape parameters, and Tmax is the maximum temperature at which Pd grows (fig. S4D). The model was fit and parameters were estimated by nonlinear (weighted) least-squares estimates [nls function in R (38)] and converted to an hourly rate to match the parameter units for hibernating bats.

RH rates >80% are usually necessary for fungal growth, though this ranges from 76 to 96% depending on temperature, substrate, and species (3942), and this requirement is true for Geomyces, close relatives of Pseudogymnoascus, even in extreme conditions (43). In the absence of experimental data from which to estimate humidity-dependent parameters, we estimated the effect of RH on Pd growth by assuming the rate of bat colony declines from WNS at different RH directly correlates to the effect of RH on Pd growth. M. sodalis is the only species for which data were reported, and the most dramatic declines occurred at high RH (21). We use the data from M. sodalis and assume that Pd colonizes and grows in the same way across all species. We extracted the data for our analyses using the digitize function (44) in R. To model the effect of RH on Pd growth, we fit a saturating Michaelis-Menton (MM; Eq. 3) function to the data using nls to obtain the Pd growth rate and humidity relationshipEmbedded Image(3)where μ1 and μ2 are the MM shape parameters. We set maximum growth to occur at maximum RH, which corresponded with greatest population declines. We did not constrain the MM function to go through zero at lower levels of RH because one colony of M. sodalis still experienced a decline at the lower RH levels we modeled. We fit the MM function to the log10 (λ) (21) (λ, the finite rate of population change), which is similar to the instantaneous rate of population change (r, the natural log of λ) (fig. S5E). To integrate temperature and humidity relationships with growth, we assumed that the greatest growth in culture media (37) represented maximum growth at saturation (100% RH) and scaled our estimated Pd growth for each temperature as a function of RH (fig. S4E).

Thus, to incorporate the effects of Tb and RH on fungal growth, we modeled hourly Pd growth as a function of Tb (GT, Eq. 2), scaled by the effects of humidity (GRH, Eq. 3), for the duration of winter (twinter) asEmbedded Image(4)

Hibernation energetic model

To model the energetic requirements of a hibernating bat with WNS in a way that allowed us to incorporate a reduction in the time in torpor (ttor) due to increased arousal frequency, we adapted a bat energetic model (16).

The energy required per torpor bout (Ebout) is a function of both euthermic energy costs (Eeu) and duration (teu) and torpor energy cost (Etor) and duration (ttor,Pd), plus the additional energy usage required to arouse from torpor (Ear) for the arousal duration (tar)Embedded Image(5)Energy costs during euthermic periods (Eeu) are RMRs plus energy cost due to euthermic thermal conductance (Ceu) at hibernacula temperatures (Ta) below the lower critical temperatures (Tlc, Eq. 6). The Tlc is the temperature at and above which bats are thermally neutral and minimal RMR occurs, and below which additional cost occurs (16) [note that in the parameter ranges we explore, Tlc (32°C) is always >Ta]. ThusEmbedded Image(6)

During torpor, TMR is substantially lower than RMR, and metabolic processes during torpor (Etor) decrease exponentially as a function of the hibernacula temperature (the Q10 effect) until a temperature at which the bats defend their body temperature (Ttor-min) is reached. At this point, Etor is analogous to Eeu and depends on TMR and the thermal conductance at torpid temperatures (Ct) below Ttor-min. We model the energy use per unit time in torpor, Etor, asEmbedded Image(7)Embedded Image(8)where TMRmin is the minimum metabolic rate in torpor that occurs at temperature Ttor-min. Q10 is a quadratic in temperature (Table 1).

Experiments have recently demonstrated that EWL is a key driver in determining arousal frequency in hibernating bats (26, 27), but we used Ta to drive ttor phenomenologically because Ta and RH are highly correlated and only below 2°C (and perhaps as low as −2°C) is EWL a better predictor of torpor duration (25). Although EWL is thought to ultimately drive arousal during torpor, and so ttor, Pd does not grow below freezing, and above freezing, the duration of ttor can be predicted equally well using Ta [R2 = 0.53 for Ta and R2 = 0.5 for EWL (25, 45)]. Thus, we model ttor,Pd (torpor time with Pd infection) as a function of Ta. Because Pd infection decreases the duration of bat torpor (ttor) (8, 18), we model the reduction in ttor due to Pd growth (ttor,Pd) as ttor/Pd (Eqs. 9 and 10). We further explored this relationship using sensitivity analysis, as described below.

Given Etor, we suppose that ttor is such that the total energy expended over time is constant, that is, if the temperature Ta is at the minimum temperature that the bat defends (Ttor-min), then the energy use is minimized so that the duration of ttor is maximized, whereas if the temperature is away from this optimal time, energy use is higher and thus ttor declines. In the original model, the authors include a parameter, k, an analytical constant that attempts to yield equal values of ttor for a given TMR above and below Ttor-min. However, this works for single values below Ttor-min only (fig. S1). To generalize the model and remove this scaling parameter, which would need to be reestimated for every point below Ttor-min, we let Etorttor = W for some constant W. Thus, ttor = W/Etor, and with Pd, ttor,Pd may be given asEmbedded Image(9)Embedded Image(10)where ttor-max is the maximum time in torpor that occurs at Ttor-min (that is, ttor-max = W/Etor) (fig. S1). Removing k means no constant need be estimated and the model is free to be used for any species for which the parameters in Eqs. 9 and 10 are available. Removing Pd from Eqs. 9 and 10 gives ttor with no Pd infection.

The energy costs of arousal (Ear) are modeled as a function of Tb and tissue-specific heat capacity (S) of hibernating bats (16), and we model Ear asEmbedded Image(11)where the change in Tb during arousal from torpid body temperature (Ttor) to euthermic temperatures (Teu) isEmbedded Image(12)Embedded Image(13)

To model Pd growth as a function of Tb, we modeled Tb as torpid body temperature (Ttor) (Eq. 13), unless the bat arouses to euthermic temperatures (Teu) (Eq. 14); thusEmbedded Image(14)Embedded Image(15)

The total energy requirement for the winter, Ewinter, is modeled as a function of the duration of winter (twinter), the expected duration of each torpor bout including arousal time (ttor,Pd + teu + tar), and the energy required for each torpor bout (Ebout)Embedded Image(16)

To convert energy requirements, Ewinter, to g fat required for winter hibernation, we used g fat = Ewinter × 20.1/(39.3 × 100) (46). To calculate species-specific survival times for a hibernating bat, we then ran the parameterized model for a species for a range of Ta and RH to estimate the length of winter that animals could survive given their prehibernation fat reserves. We assumed that the proportion of body fat available for hibernation was a general feature of bats (17, 4750), which report means of 25 to 30% BM as fat deposits for hibernating bats, and so we used 30% as a conservative metric and assume that this trait scales with body size.

Basic parameters for four species—M. lucifugus, E. fuscus, M. myotis, and E. serotinus—are in Table 1. We note that E. serotinus has recently been subdivided into multiple species (51); however, we treat it as a single species here given that its parameters are reported in the literature from when it was recognized as a single species and are from the current distribution (51). To obtain the necessary host trait parameters typically not reported for other species (for example, Tlc and Ttor-min) from more frequently measured parameters (for example, BM) reported in the literature, we used linear regression for prediction.

Model validation

We validate our model by comparing the results of our model predictions for M. lucifugus to experimental study results (the only species for which such data were available), comparing the predictions for four species against observations on two continents (see Predicting survival times within species distributions) and by modeling the percent surface area (SA) of the bat infected. We calculated the percent SA of the bat covered by the predicted Pd infection by calculating the total cm2 Pd grew to ensure the predicted percent SA infected with Pd was plausible (that is, ≥0 and ≤100%). To estimate the SA in m2 for bats, we estimated the wing (w) SA using (52)Embedded Image(17)

Total SA was then estimated given that wings are approximately 85% of a bat’s total SA (53). We model a single introduction of Pd and do not account for increased SA enabling increased infection; thus, the modeled percentage of skin SA infected by Pd ranged from 0 to 40% excluding wing SA and from 0 to 6% including wing SA for M. lucifugus (fig. S3), generally corresponding to observed proportions of affected wing area [see figures and histological sections within previous studies (4, 5, 10, 18)].

Sensitivity analysis

To assess the response to variation in parameters, we defined our response as the proportional change in g fat required to overwinter over the complete environmental parameter space for M. lucifugus. We converted model-derived g fat required to survive hibernation into probabilities, using the mean and SD (47), such that if >~34% BM was required, then survival (S) was zero, and if <~26% BM was required, then survival was assured (S = 1) (fig. S5F). We assigned intermediate values based on the cumulative probability of bats being those weights (47). To compare the changes in g fat required to overwinter given the change in temperature-dependent growth and the overall effect of Pd on ttor,Pd, we then compared the proportion of environmental parameter space for the combined ranges of hibernacula Ta and percent RH for which survival was certain (S = 1).

For our sensitivity analyses, we present results for survival predicted after a 6-month winter. We chose 6 months because hibernation durations up to 6.4 months are reported in Canada (16), suggesting that these are upper limits, but model results for 5 months are qualitatively similar. We present the results as the change in survival over the whole environmental parameter space relative to our default environmental values for changes to host- and fungal-related parameters (Table 1).

For the sensitivity analysis, we used a multiparameter sensitivity analysis proposed by Blower and Dowlatabadi (54). We constructed 100 random parameter sets using stratified random samples from uniform distributions extracted with Latin hypercube sampling [lhs package in R (55)] spanning a range of potential values from the default (often maximum values) to 10% lower than default values (see Table 1). Note that increasing the MM function parameter μ2 in the fungal growth model (Eq. 3) led to asymptotic behavior, and so for comparison values for all, parameters were sampled between the default and 10% lower values only for constancy. PRCCs between each parameter and model output determined the relative importance of each parameter. In a K + 1 by K + 1 symmetric matrix, C, where K is the number of parameters, and 1 to N is the rank of each column defined by the set (r1i, r2i, … rki, Ri), where i is the run number and τ is the average rank [(1 + N)/2], the matrix C is defined with elements Cij, such thatEmbedded Image(19)Ri replaces rij and ris for the Cj, K+1 elements, and the inverse of C becomes bij for matrix B. The PRCC between the ith parameter and yth variable is thenEmbedded Image(20)

In the PRCC analysis, a positive PRCC indicates that as the parameter increases, the output variable of interest increases, and a negative PRCC indicates that the variable decreases with increasing parameter values (54, 56). Significance (tiy) of PRCCiy is determined by a Student’s t distribution with N − 2 degrees of freedom; thusEmbedded Image(21)

Because some parameters were not statistically independent, we compared the PRCC results to univariate analyses and results were qualitatively similar, so only PRCC results are presented. We checked for nonlinear, nonmonotonic relationships among the output and the LHS parameter inputs (fig. S6) (57).

Predicting survival times within species distributions

Finally, to predict the survival times of hibernating bat species across their ranges, we obtained climate data from the U.S. National Oceanic and Atmospheric Administration Earth System Research Laboratory site using the RNCEP R package (58) and species distributions from the International Union for Conservation of Nature ( We used these species distributions in the absence of good winter hibernacula data (16). We used mean maximum winter (January to April, 2001–2002) RH, as a proxy for cave RH, and mean annual surface temperatures, which predict cave temperatures (24). Because of the revised phylogeny of E. serotinus and thus geographic range, we created a new coverage file based on the new distribution (fig. S7) (51). To ensure that Pd remains positive as mean annual surface temperatures exceed the maximum or minimum required for Pd growth, we fixed ttor beyond Tmax and Tmin. We estimated the difference between predicted survival time given environmental conditions within a species range and average winter length, determined by the period of the year where mean nightly temperature from 12 a.m. to 6 a.m. was below freezing. This method of winter duration was determined a priori (fig. S8). The differences between winter duration and predicted survival times were then plotted within species distributions (Figs. 3 and 4).

To assess the significance of the changes in the distributions of the winter durations between North American species, we used a Harrell-Davis distribution free estimator with a percentile bootstrap (59) to test for differences between quantiles of the different distributions of survival times.

We coded all analyses in R (38). For data manipulation, interpolation, and plotting, we used the RNCEP (58), lattice (60), MASS (61), plyr (62), tgp (63), akima (64), maptools (65), raster (66), sp (67), adehabitat (68), maps (69), rgdal (70), shape (71), maptools (72), ggplot2 (73), rgeos (74), ggmap (75), mapdata (76), mapproj (77), grid (38), gridExtra (78), PBSmapping (79), data.table (80), deSolve (81), and WRS (82) packages in R. Code and data files are available at


Supplementary material for this article is available at

Fig. S1. We generalized the energetic model by revising the time in torpor (ttor) such that energy use is equal for values of ttor for a given minimum TMR above and below a lower minimum set-point temperature that is defended (Ttor-min).

Fig. S2. Predicted times to deplete overwinter energy reserves for bats without infection by the fungus Pd.

Fig. S3. Predicted percent of the SA of hibernating M. lucifugus infected by the fungus Pd based on a single introduction of Pd.

Fig. S4. Histograms of the difference between the predicted winter duration and the modeled hibernation durations within the two North American (A and B) and two European (C and D) bat species’ distributions under different climate change scenarios.

Fig. S5. Summary of model fitting and parameter estimation.

Fig. S6. Parameter input values plotted against the output values used in our sensitivity analysis.

Fig. S7. The distribution of serotine bats (E. serotinus) (51) used in this study.

Fig. S8. The predicted winter duration in months within the two North American and two European bat species distributions that underlie the model predictions.

This is an open-access article distributed under the terms of the Creative Commons Attribution-NonCommercial license, which permits use, distribution, and reproduction in any medium, so long as the resultant use is not for commercial advantage and provided the original work is properly cited.


Acknowledgments: We thank the Webb and Pulliam laboratories, especially A. Luis, C. Leach, and C. Pearson, for useful discussions. M. Verant (USGS National Wildlife Health Center) offered helpful input on earlier versions of this manuscript, and J. Juste (Estación Biológica de Doñana, Spain) confirmed the distribution of serotine bats. Any use of trade, firm, or product names is for descriptive purposes only and does not imply endorsement by the U.S. Government. Although software programs have been used by the U.S. Geological Survey (USGS), no warranty, expressed or implied, is made by the USGS or the U.S. Government as to the accuracy and functioning of the programs and related program materials nor shall the fact of distribution constitute any such warranty, and no responsibility is assumed by the USGS in connection therewith. Funding: D.T.S.H. acknowledges funding from a David H. Smith Fellowship in Conservation Research from the Cedar Tree Foundation and Society of Conservation Biology. D.T.S.H., J.R.C.P., and C.T.W. have support from the Research and Policy for Infectious Disease Dynamics (RAPIDD) program of the Science and Technology Directorate (U.S. Department of Homeland Security) and the Fogarty International Center of the U.S. NIH. Author contributions: D.T.S.H., J.R.C.P., P.M.C., and C.T.W. designed the study; D.T.S.H., J.R.C.P., J.C.M., and C.T.W. developed the model, code, and statistical analyses; D.T.S.H. analyzed the data and led the writing with contributions from all authors. Competing interests: The authors declare that they have no competing interests. Data and materials availability: All data needed to evaluate the conclusions in the paper are present in the paper and/or the Supplementary Materials. Additional data related to this paper may be requested from the authors. Code and data files are available at
View Abstract

Stay Connected to Science Advances

Navigate This Article