Spatial spread of white-nose syndrome in North America, 2006-2018

White-nose syndrome has caused massive mortality in multiple bat species and spread across much of North America, making it one of the most destructive wildlife diseases on record. This has also resulted in it being one of the most well-documented wildlife disease outbreaks, making it possible to look for changes in the pattern of spatial spread over time. We fit a series of spatial interaction models to the United States county-level observations of the pathogenic fungus, Pseudogymnoascus destructans, that causes white-nose syndrome. Models included the distance between caves, cave abundance, measures of winter length and winter onset, and species richness of all bats and hibernating bats only. We found that the best supported models included all of these factors, but that the particular structure and most informative covariates changed over the course of the outbreak, with winter length displacing winter onset as the most informative measure of winter conditions, and evidence for the effects total species richness and hibernation varying from year to year. We also found that weather had detectable effects on spread. While the effect sizes for cave abundance and species richness were relatively stable over the length of the outbreak, distance became less important as time went on. These findings indicate that although models produced early in the outbreak captured important and consistent aspects of the spatial spread of white-nose syndrome, there were also changes over time in the factors associated with spread, suggesting that forecasts may be improved by iterative model refinement.


Introduction
The outbreak of white-nose syndrome (WNS) in North American bats is driving massive declines in multiple hibernating bat species (Frick et al., 2010(Frick et al., , 2015Langwig et al., 2012, Thogmartin et al., 2012 and continues to spread across the continent (Lorch et al., 2016;Hoyt et al., 2021). The etiological agent of WNS is the psychrophilic fungus Pseudogynmoascus destructans (Lorch et al., 2011;Verant et al., 2012;Warnecke et al., 2012), which is known to infect at least 12 species of North American cave-dwelling bats Frick et al., 2015Frick et al., , 2016Hoyt et al., 2021), and has caused declines of up to 99% in abundance of Myotis septrionalis, M. lucifugus, M. sodalis and Perimyotis subflavus populations (Frick et al., 2010(Frick et al., , 2015Langwig et al., 2012;Thogmartin et al., 2012). Since its earliest detection in Schoharie County, NY in 2006, the pathogen has been detected in hibernacula, places including caves and mines where bats hibernate, in 39 states (Hoyt et al., 2021). This combination of host mortality and spatial extent makes WNS one of the most damaging wildlife diseases known, along with chytridiomycosis in frogs (Rohr et al., 2008) and facial tumor disease in Tasmanian devils (McCallum et al., 2007).
The well-documented spread of WNS since 2006 offers a valuable opportunity to consider how our understanding of the spatial dynamics of this disease has changed over time. A comparison of multiple models fit to the pattern of spatial spread during the first 6 years of WNS expansion in the United States found that a gravity model based on distance between counties and the number of hibernacula in each county better explained the observed pattern of spatial spread than models based solely on distance (Maher et al., 2012). The best-supported model also included a covariate for the average winter length in each county whereby longer winters were associated with a greater probability of occurrence of P. destructans (Maher et al., 2012). In addition to determining the role of distance, hibernacula abundance, and climate in the observed spread of WNS, this model was used to create forecasts for the future spread of WNS, summarized as the median number of new counties infected each year and the distance from origin over time (Maher et al., 2012). These forecasts also captured the order in which counties recorded infections in many cases (USFWS; www.fws.gov/whitenosesyndrome/maps.html, 2019).
Pseudogynmoascus destructans has continued to spread steadily across North America, providing an ever-lengthening time series of newly infected counties. (It appears that once a county has been infected, it never becomes "uninfected" again.) Studies of the disease dynamics within and between hibernacula have highlighted the role of environmental conditions within hibernacula (Wilder et al., 2011;Langwig et al., 2012;Lilley et al., 2018), finding that longer time spent hibernating leads to higher fungal loads (Langwig et al., 2020) and mortality (Lorch et al., 2011, Warnecke et al., 2012. Ongoing research also aims to determine whether autumn swarming or winter movements are more important to transmission between hibernacula (Kramer et al., 2019;Langwig et al., 2015Langwig et al., , 2020. In this context, the previous finding that winter duration increases spread rate raises the question of whether the average length of winter (as tested by Maher et al. in 2012) matters because of regional differences in bat behavior or whether the specific conditions of each year drive spread. If the latter is the case, differences in spread rate could correlate with the number of hibernacula visited during autumn swarming or outcomes of mid-winter movements. The larger number of winter observations, compared to when the first model was created, offers additional information for answering this question. This raises the question whether the long distance spread to Washington was an event consistent with the existing understanding of spatial dynamics of spread. Or did the short period of disease progression to that point prevent effective forecasts of future dynamics? While the dispersal of P.
destructans to Washington has been suggested to be human-mediated (Hoyt et al., 2021), humanaided dispersal may have been a factor in earlier spread as well, and therefore implicitly recognized by the model.
Here, we revisit the spatial spread of WNS to assess how well earlier models matched subsequent spatial dynamics and factors influencing disease transmission. We extend the models previously developed to focus on three specific questions: 1) Do year-to-year weather conditions have a detectable relationship to spatial spread? 2) Was bat species richness, which was not well supported in the earlier modeling approach, a more important factor in subsequent spread? and 3) Were the effects of various factors in the original model stable during subsequent spread in space and time?

Data
Time of infection by county was obtained from the U.S. Fish and Wildlife Service (USFWS; www.fws.gov/whitenosesyndrome/maps.html, 2019), recording the first year in which either WNS was observed, the infectious agent P. destructans was discovered (e.g. in a soil sample), or both. Given that fall and winter are the crucial periods for transmission (Langwig et al., 2020;Hoyt et al., 2021) and detection (Warnecke et al., 2012;Langwig et al., 2015), the epidemic year was assumed to begin June 1 st and end May 31 st .
Number of caves in each county in the contiguous United States was provided by Maher et al. (2012) and the Euclidean distance between county centroids from the same reference (Maher et al., 2012). County-level bat species richness was estimated from NatureServe (www.natureserve.com) for the 46 bat species occurring in the contiguous United States at any point in the year, as provided by Maher et al. (2012). Each species was classified as hibernating or non-hibernating based on USGS information (USGS, 2020).
We quantified average and annual winter conditions using the minimum daily ftp://ftp.ncdc.noaa.gov/pub/data/ghcn/). Length of winter was calculated at each weather station as the number of days with temperatures below 10°C between June 1 st and May 31 st of the following calendar year. A smooth representation of the length of winter for the contiguous USA (NAD83 projection, resolution of 0.124 degrees) was obtained using anisotropic ordinary kriging performed with ArcGIS (ver 10.7, www.esri.com). The mean winter length was then calculated for each county using the package 'raster' (ver 2.9-23) in R (Hijmans 2019 citations). This was the same procedure used in Maher et al. (2012). We used the average length of winter for each county over the full time period as a measure of the winter climate of each county.
In addition to winter length, we also considered the number of days from June 1 until the onset of winter as another indicator of winter conditions, where onset of winter was defined as the first two consecutive days with a minimum daily temperature reading below 10°C for each operational weather station. The number of days until the start of winter for each county was interpolated as above and rounded to the nearest whole number.

Model
Following Maher et al. (2012), gravity models were fit to the observed time of infection to estimate the effects of weather, distance between counties, cave density, and bat species richness on the spatial spread of the epidemic. All models were of the form where " #$ is the probability that county i does not become infected from previously infected county j and is a function describing the inverse of transmission intensity from j to i. The basic version of f is: where #$ is the distance between county i and county j, # is the number of caves in county i and $ is the number of caves in county j. The fit parameter 0 is the background infection rate, and 2 and 3 adjust the distance and product of cave densities respectively. Alternative models were constructed by adding additional terms to the function f (Supplementary Table 1). We fit a total of 19 models with each winter condition measure crossed with the two species richness estimates.
A change from Maher et al. (2012) is the inclusion of time-varying environmental variables.
Previous models used fixed values for all variables representing each county, including the average length of winter in each county over the period of epidemic spread. Here, we also considered models that replaced the average winter conditions with the annual conditions corresponding to each year, as we hypothesized that specific annual conditions might affect either bat behavior or infection intensity and thus influence pathogen transmission (Verant et al., 2012;Langwig et al., 2015;Langwig et al., 2020). We also considered models where the probability of transmission depended on conditions in the prior winter, with the rationale that the observation of WNS in a bat colony (symptomatic disease) may reflect arrival of the fungus (the pathogen) in the previous winter.
Unknown parameters were estimated using maximum likelihood; model selection was performed using Akaike's Information Criterion (AIC). To consider the role of accumulating information and the influence of unique spread events on model fit, we calculated the AIC value annually for each model. We considered models with consistently low AIC values to be the most robust, but recognize that fluctuations in model performance from year-to-year provide additional insight into the interaction between epidemic pattern and covariates, and guard against model selection being dominated by the conditions of a single, possibly anomalous year at which evaluation occurred. The difference in AIC between the model with the lowest AIC in that year and each of the other models in that year (ΔAIC) was calculated for years 2008-2018. Similarly, we looked at the estimated parameter values for years 2009-2018, focusing on the models that had AIC = 0 at some point during the epidemic (Figure 3). We excluded the first two years of estimates because these early estimates were less reliable and exhibited larger changes than for later years (Supplementary Table 2).

Results
Annual measures of winter severity in U.S. counties were more strongly correlated with average winter severity (Pearson's r = 0.89) than winter start (Pearson's r = 0.88, Figure S1). Estimated values of parameters associated with background transmission ( 4 ), winter severity ( 5 ), distance ( 6 ), winter severity ( 7 ), and bat species richness ( 8 )varied over the course of the outbreak (Figure 3). Early in the outbreak the models including winter length had a higher estimate of background infection rate, suggesting the other covariates in these models were not capturing as much of the pattern of spread, but in later years the estimates of background transmission were similar across models. (Figure 3a). As expected, winter length was positively correlated with the probability WNS was observed. Similarly, winter start had a negative relationship (Figure 3b), with the magnitude of the effect of winter conditions declining over time (moving towards zero) for both measures. The estimated effect of distance on spread declined by almost an order of magnitude over time in all models (Figure 3c). Cave density had a more stable influence on spread patterns (Figure 3d). Finally, the estimated species richness parameter was stable and similar whether considering total bat richness or hibernating bat richness (Figure 3e).

Discussion
The well-resolved spatial dynamics of WNS allowed us to distinguish aspects of the dynamics that remained consistent as spread continued across the continent from those that changed as the outbreak proceeded. Winter conditions remained an important factor in spatial dynamics throughout the outbreak, but the most informative measure of winter conditions changed over the course of the epidemic, which could be due to differing drivers in different regions or correlations with unmeasured causal mechanisms related to winter severity . However, our analysis indicates that year-to-year variation in winter conditions is playing a measurable role in epidemic dynamics despite strong correlations with average climate in each county. We found that geographic patterns in species richness are related to the probability of transmission between counties, and detection of this relationship has become more reliable as the outbreak has proceeded. In addition we found that the effect of distance on transmission declined from estimates by Maher et al. (2012), and this decline was gradual, rather than an abrupt response to the 2016 detection of P. destructans in Washington State (Lorch et al., 2016).
The long distance transmission event to Washington State was of particular interest as this either represented the realization of an event (albeit one of low probability) encompassed by the model of Maher et al. (2012), or an unforeseen possibility demonstrating that models fit early in the outbreak missed important aspects influencing spread. The suggestion that this spread was human-mediated due to its extreme distance from the nearest known infected hibernacula (Hoyt et al., 2021) could be consistent with either explanation. The inclusion of species richness in most of the supported models in recent years represents one possible signal of the spread to the western United States. The community of bats in the western U.S. has different species and higher richness than the region where WNS emerged, and these areas are separated by a relatively species-poor area through the middle of the country (Maher et al., 2012). At the same time we did not observe an abrupt change in the parameter controlling the influence of distance in the model after this long distance event, instead accumulating data was already indicating distance was less of a barrier than found early in the epidemic.
It is important to note how the geographic pattern in species richness may correlate with additional factors we were unable to consider here. The sign of the estimated parameter means that the model estimated probability of transmission went down as species richness increased.
The continental pattern of spread has been from an area of low richness in the Northeast, along the Appalachian Mountains, then towards the Ozark Plateau, where more species ranges overlap.
West of the Ozark Plateau, species richness declines into the Great Plains, and then more species occur west of the 100th meridian. As such, the model may reflect a bottleneck of suitable hibernacula that corresponded to a relative increase in richness. Further, the geography of species richness reflects a combination of ecological and evolutionary processes (Miller-Butterworth et al., 2014) and ignores population heterogeneity (e.g. Wilder et al., 2015). Thus, it is possible that the specific mechanism associated with this covariate has nothing to do with the ecological community and it would be premature to conclude that counties with greater species richness in the western U.S. necessarily have low risk of pathogen transmission.
The influence of winter conditions on spatial spread of WNS is of interest not only for understanding relative risk of different bat populations (Maher et al., 2012;Wilder et al., 2011;Lilley et al., 2018), but also because of the importance of environmental conditions in disease severity (Langwig et al., 2012;Hayman et al., 2016) and hypothesized pathways of transmission between hibernacula (Langwig et al., 2020). Research suggests that bat movements during the winter may be a driver of transmission between hibernacula (Langwig et al., 2020) and ambient temperature affects bat activity during winter (Bernard & McCracken, 2017). Our findings suggest that aspects of winter severity related to winter length and the timing of winter onset both improve model performance, with winter onset appearing more informative early in the epidemic, while length of winter is favored over the full period up to the present. We also found evidence that the conditions in particular winters drive transmission risk and that this is detectable despite the strong correlation with average conditions. However, although the previous winter is more valuable for predicting the risk of transmission when looking at the full period to the present, the tendency for the best fit model to vary between years suggests we should not overestimate the strength of this inference. A signal of the previous winter could indicate detection occurring the year after fungal arrival in a hibernacula (Langwig et al., 2015).
Improving this aspect of the model may depend on improved understanding of which factors lead to the onset of hibernation or movement among hibernacula.
Examining how models of the spread of WNS have changed over time provided a richer understanding of spatial spread in this ongoing outbreak. At all time points, distance, cave density, and winter conditions best explained the spread, suggesting that even models fit to the exponential phase of an epizootic can provide information sufficient for assessing risk and designing interventions. At the same time, the models best able to explain the observations changed significantly as the outbreak proceeded, as did the relative strength of multiple factors in those models. Hence, the effort of updating the model provides additional insight into WNS dynamics and will facilitate improved forecasts of future spread. This aspect of disease modeling is likely important to any long term outbreak event and highlights the value of revisiting even high performing models.