Study design
This study began in 2013, before P. destructans invaded Michigan and Wisconsin, and continued through 2020, after all Michigan and Wisconsin hibernacula were invaded by the pathogen. Throughout this period, we visited 22 hibernacula twice per winter during bat hibernation to quantify bat colony sizes, individual bat roosting temperatures, fungal loads (to which we added a constant 0.0001 before transforming to the log10 scale), and recapture probabilities (Fig. 2). We sampled bats during early hibernation in November, when more than 95% of swarm activity was expected to be finished51, and again during late hibernation in March, when less than 1% of spring emergence activity was expected to have begun51,52. For each sampling event, we counted all bats of all species within the site. We focused our analyses on the little brown bat (Myotis lucifugus). This species suffered severe declines due to WNS24,53, and was the only species abundant enough to provide sufficient sample sizes for our analyses. We divided bat counts by sections within the hibernaculum (i.e., “rooms”) that potentially varied in microclimate, and we used HOBO Pro v2 data loggers to continuously record temperatures every 3 h in these sections (Supplementary Figs. 2 and 7).
After counting all bats, we haphazardly sampled 20–25 individual little brown bats stratified across sections and roughly in proportion to the number of bats in each section (Fig. 2). For each sampled bat, we used a Fluke 62 MAX IR laser thermometer to quantify the temperature of the substrate directly adjacent to the bat (<1 cm; the “roosting temperature”). We then used a previously developed protocol to collect a standardized swab of each bat’s forearm and muzzle, which we stored in RNAlater until we could quantify fungal loads using qPCR54. Finally, we banded all swabbed bats that could be safely removed from their roosts (mean = 73%, median = 100% of all swabbed bats) with an aluminum lipped-band (2.9 mm; Porzana Ltd., Icklesham, E. Sussex, U.K.), so that they could be re-sighted and individually-identified during subsequent visits (Fig. 2).
All sites and sampling procedures were covered by the appropriate state and federal permits, Virginia Tech IACUC protocol #17-180, and University of California, Santa Cruz IACUC protocol Kilpm1705. We followed field hygiene and decontamination protocols in accordance with United States Fish and Wildlife Service.
Over winter changes in fungal loads on bats
For infected, banded bats that were sampled during early hibernation and recaptured during late hibernation, we compared estimated fungal loads between the two time points to quantify how changes in fungal loads were affected by early hibernation roosting temperatures. In particular, we calculated fungal load change rates (λ) on individual bats using our estimates of fungal loads during early hibernation (N0), fungal loads during late hibernation (Nt), and the number of weeks between collection of the early and late samples (t), where λ = ln((Nt/N0)^(1/t)). We excluded bats that were uninfected in November because if they had become infected between November and March, we could not have determined how long they had been infected (t in the λ equation) (Fig. 2).
Though the relationship between changes in fungal loads and early roosting temperatures on recaptured bats was linear over most of the range of early roosting temperatures that we observed for banded bats, fungal growth rates should be hump-shaped over large temperature ranges, with a minimum threshold temperature, an optimum temperature (Topt), and a maximum threshold temperature (Tmax). One such nonlinear temperature-dependent growth model is the Logan-10 function55:
$$lambda left( T right) = aleft( {frac{1}{{1 + ce^{ – pT}}} – e^{ – frac{{T_{max } – T}}{{Delta T}}}} right)$$
(1)
where T is temperature, ∆T is the width of the upper boundary layer (Tmax – Topt), and a, p, and c are shape parameters (Eq. 1). We fit this curve to our field-collected fungal load change data.
This fungal load analysis did not include bats that were not recaptured, and thus the best-fitting temperature-dependent growth curve based on recaptured bats might be shallower than the curve for the whole population. We accounted for the potential missing bat bias using two methods that addressed the limited data from bats that roosted at relatively warm temperatures and bats that began hibernation with relatively high fungal loads. First, we expected that we would recapture relatively few bats that roosted near the optimum temperature for fungal growth, which could lead to model convergence issues. Therefore, we fit the temperature-dependent fungal growth curve in a Bayesian framework, using uniform priors for all parameters (except for the maximum temperature for fungal growth, Tmax) that were weakly informed by data from a previously-published laboratory fungus growth experiment (see below)17. Second, we expected that we would recapture relatively few bats that had high fungal loads during early hibernation, and that any recaptured bats with high early loads would have had lower changes in fungal loads that allowed them to survive their initially high loads. Therefore, we allowed the intercept of the temperature-dependent fungal growth curve to vary with early fungal load, where the curve for the lowest early fungal loads would be the least likely to be affected by disease-induced mortality.
We acquired reasonable ranges for Topt, Tmax, and the shape parameters using previously published fungal growth data from a laboratory experiment17, where P. destructans was grown on Sabouraud dextrose agar plates for five weeks at each of nine average growth chamber temperatures: 0.8, 2.0, 4.6, 7.2, 11.9, 16.0, 17.6, 19.0, and 21.4 °C. The 7.2 °C treatment was repeated in three separate trials, which we combined for our analyses. We assumed that all temperatures had the same fungal carrying capacity (K) and selected a carrying capacity that minimized the mean squared error of model fits across all temperatures: K = 0.85 cm2. Using this value and assuming that starting colony sizes (N0) were the sizes observed at Week 1 during the laboratory experiment, we fit the logistic growth model (K/(1 + ((K/N0) − 1)*exp(−r*Day)) to the fungal colony size data from Weeks 2–5 (28 days of growth) to estimate intrinsic fungal growth rates (r) at each temperature. We then fit the temperature-dependent growth curve (Eq. 1) to the nine r estimates using uninformative uniform priors for a, p, and c and less flexible normal priors for Topt (mean = 14, SD = 5) and Tmax (mean = 21, SD = 1) in a Bayesian framework, which provided us with point estimates and 95% credible intervals for a, p, c, Topt, and Tmax. These models were run using three MCMC chains of 60000 iterations each, with a burn in of 30000 iterations. Model convergence was assessed using visual inspection of MCMC runs and by confirming that (hat R) values were less than 1.01. All analyses were run with package ‘R2jags’ in R56.
We used the lab-derived 95% credible intervals as uniform priors for all parameters (except Tmax) in a second Bayesian analysis that fit the Logan-10 curve to the change in fungal load data from recaptured bats using three MCMC chains of 160000 (burn in = 80000). We also used an uninformative uniform prior for the load-dependent intercept term [b~unif(−0.5,0.5)]. Instead of estimating the thermal maximum, we set Tmax in the Logan-10 curve for the field data to a constant (the estimated thermal maximum from the laboratory study), because no bats roosted close to Tmax. To confirm that these priors only weakly informed the model and that any reasonable curve shape would be possible, we ran a simulation to visualize the curve shapes that would be considered during the model fitting process. In particular, we plotted Logan-10 growth curves for all combinations of the minimum, maximum, and median of the ranges for each parameter (Supplementary Fig. 1).
We used early roosting temperatures in our model because we expected that large changes in fungal loads during early hibernation would be more detrimental to bat survival than large changes in fungal loads during late hibernation57. We calculated how much variability in fungal load change (log10λ) the best-fitting load-dependent Logan-10 curve explained (R2) by running a linear regression of observed versus predicted values. We also confirmed that there were significant correlations between log10λ and early roosting temperatures and between log10λ and early fungal loads using Spearman’s rank correlation tests, which were performed in a frequentist framework without using priors from the laboratory experiment. We then repeated all analyses with late roosting temperatures instead of early temperatures and confirmed that late roosting temperatures always explained less variation (i.e., lower R2), if any, in log10λ than early roosting temperatures.
Bat recapture rates
Bats experiencing WNS symptoms often leave their hibernacula, where they are likely to die from exposure or starvation on the landscape58. In our study, these would be bats banded during early hibernation that left the hibernaculum before our March survey, which occurred weeks before average spring emergence dates51. Factors other than disease-induced mortality could have resulted in a failure to recapture bats (e.g., emigration/immigration, handling disturbance, non-disease related mortality), which would have reduced the likelihood of finding a significant effect of our predictors. However, only 5% of bats that were banded/sighted during early hibernation and missed during our late hibernation surveys were later re-sighted in a different year, and always in the same hibernaculum. Our model yielded qualitatively identical results whether we included these bats are recaptured or not in a given year. Some bats might also have died after our late hibernation survey, but before spring emergence, which would cause us to underestimate overall mortality. Therefore, we assume that recapture success was a conservative proxy for overwinter survival22,59,60.
We used a mixed effects logistic regression with a logit link to quantify how early hibernation roosting temperatures and fungal loads affected the probability that bats infected during early hibernation were (=1) or were not (=0) re-sighted in late hibernation the same year. To quantify the predictive capabilities of this model, we performed 5-fold cross-validation on 1000 random divisions of our dataset, and calculated the average Area Under the Curve (AUC) for the resulting 1000 Receiver Operating Characteristic (ROC) curves, which can be interpreted as the percent of the test data that was successfully predicted by the training models.
Across invasion years, there was a trend towards increased recapture probabilities with time since WNS arrival (Fig. 5a), but including this covariate in the logistic recapture model with early hibernation temperatures and fungal loads resulted in very high variance inflation factors (VIF = 33.5). To determine whether invasion year explained any additional variation in recapture success, as might be expected if bats were evolving physiological resistance or tolerance to P. destructans, we ran a partial regression comparing the Pearson residuals for recapture success given early hibernation temperatures and fungal loads (logistic GLMM described above) to the Pearson residuals for invasion year given early hibernation temperatures and fungal loads (Poisson regression).
Bat distribution shifts
By 2020, there were 12 hibernacula that had harbored P. destructans for at least three years (>2 years is defined here as the post-invasion period) and that we had sampled during both the pre-invasion period (up to and including year 0, when WNS was first detected) and invasion period (years 1–2 of invasion). From pre-invasion to post-invasion, bat temperature distributions might have changed across these 12 hibernacula via two mechanisms: (1) bat preferences could have remained the same while populations declined in hibernacula and/or microsites with unfavorable temperatures, causing distribution shifts via selection; and/or (2) bat preferences could have changed towards avoidance of hibernacula and/or microsites with high temperature-mediated disease mortality, causing distribution shifts via selection or learning. In both cases, we would expect the average regional bat roosting temperatures to change over time. In contrast, we would only expect the average bat roosting temperature in any given site (or section) to change over time if the site provided a wide range of temperatures to select from, which was not always the case in the hibernacula we studied. Therefore, we quantified how the regional mean early hibernation roosting temperatures for all counted bats changed across all 12 of these hibernacula from pre-invasion to post-invasion (2013–2020, up to four years after the fungus was first detected) using a Gaussian regression, without including a random site effect. We visualized this change using density plots with 1 °C smoothing bin widths.
Though we stratified bat roosting temperature sampling by section, roughly in proportion to the total number of bats present in a section, we only sampled up to 25 little brown bats per site, and we could not sample exactly in proportion to the total bats present. Therefore, calculating the mean temperatures based on only sampled bats might be misleading. Instead, we used the observed temperatures from our sampled bats to estimate the roosting temperatures for bats that were counted but not sampled. In particular, we randomly generated an additional set of temperatures for all unmeasured, counted bats based on the mean and standard deviation of the roosting temperatures of sampled bats in the same section during the same survey. We then used both known (sampled bats) and estimated (counted bats) temperatures in our Gaussian regression to determine how average temperatures changed from pre-invasion to post-invasion. We also confirmed that qualitative results were the same if we only used the temperatures from sampled bats, without estimation. By considering how bats’ early hibernation temperature distributions changed throughout invasion, regardless of infection status, we were quantifying bat roosting preferences within and between hibernacula, rather than behavioral responses to infection within a season (i.e., “behavioral fever”).
Finally, we used the number of bats counted per hibernaculum during early hibernation (November) and late hibernation (March) in all hibernacula that still had bats in post-invasion years to calculate two demographic parameters: over summer population growth rates (reflecting immigration and recruitment to sites between March and November) and over winter mortality rates (population growth rates from November to March). We used a Gaussian linear model to determine whether the log10 population growth rates were correlated with average early hibernation roosting temperatures for sampled bats in each hibernaculum, the demographic season (over summer versus over winter), and the interaction between temperature and demographic season. A significant interaction term would indicate that temperature-correlated bat immigration and recruitment rates were decoupled from temperature-mediated population declines.
All analyses described above were performed in R version 3.5.161 with packages ‘ggplot2’62, ‘lme4’63, ‘loo’64, ‘pROC’65, and ‘R2jags56. Visual inspection of residuals and predictions plots confirmed acceptable model fits for all regression models.
Reporting summary
Further information on research design is available in the Nature Research Reporting Summary linked to this article.
Source: Ecology - nature.com
