Data collection

We obtained data on the four focal bat species from state agencies in New York, Pennsylvania, West Virginia, and Tennessee. Data were from hibernation surveys completed during 1999–2011 by trained biologists as part of long-term wildlife-monitoring programs. Surveys were performed during the hibernation period from December-March in caves and mines known to serve as bat hibernacula for one or more bat species. Out of concern for negative effects of disturbing hibernating bats [12], and due to limited personnel and resources for monitoring, hibernation surveys were typically conducted once every two to three years [42]. In simple hibernacula, most accessible areas of the hibernaculum were surveyed during a single visit. In complex hibernacula, surveys were restricted to specified survey routes. Very large, complex hibernacula were divided into multiple survey routes. During surveys, bats were visually counted where they hibernated on walls and ceilings [42], [43]. In some cases, large clusters were photographed during surveys and later counted from photographs to minimize disturbance and improve accuracy of bat counts [43].

Datasets

From the full datasets provided by state agencies, which comprised data from more than 636 surveys along 163 survey routes, we selected those surveys that provided the most reliable, consistently-collected count data on hibernating bats. On the basis of survey notes, we excluded all surveys that were incomplete or inconsistent, that were focused on recording incidence of white-nose syndrome in hibernacula rather than counting hibernating bats, or that deviated from established survey routes. We also excluded routes for which only a single survey remained in the dataset. The resulting dataset included a mean of 3.95 (range: 2–10) surveys per route. Because of the potential influence of survey timing on bat counts (see Introduction), we examined the extent of within-season and across-year variation in the timing of hibernation surveys with boxplots.

We examined the dataset separately for each bat species. To reduce zero-inflation due to inclusion of unsuitable habitat in the dataset, we excluded routes in which the focal bat species was never observed. This rendered a sample of 577 surveys along 145 routes counting 982974 individual M. lucifugus, 576 surveys along 145 routes counting 68148 individual P. subflavus, 284 surveys along 62 routes counting 136386 individual M. sodalis, and 460 surveys along 109 routes counting 5206 individual M. septentrionalis. Data are provided in Appendix S1.

Global model

To evaluate temporal changes in regional abundance, we identified a global mixed-effects model a priori [63]. Route terms were assigned to random effects, to represent our sample within this population [53], [58]. Since detection probability was not explicitly estimated and metapopulation boundaries were unknown, data from hibernacula surveys represent an unknown fraction of the population rather than a complete census [43]. We therefore limited inference to relative, rather than absolute abundance [61]. To represent within-season and year-to-year variation, this model included fixed-effects terms for day-of-winter (Day), which indicated date of the survey in number of days since December 1, and year (Year). To relax assumptions of independence to accommodate hibernacula with more than one survey route, we controlled for between-site variation by nesting survey route within hibernaculum with a random grouping term [58] and modeled unmeasured differences between hibernaculum with a random-intercept term [58], [64]. Terms for the global model are represented by the following equation:(1) is the expected count for species i at time t, g is the inverse of the selected link function (in our case the inverse of the natural logarithm ln), is the mean count for a species, and are smoothing functions for Day and Year, and is a random effect for species i, at survey route k nested within location (hibernaculum) j and time t. The model assumed Poisson-distributed counts.

Model selection

We then modified the global model to create a final model for each bat species, using a step-wise reduction in fixed and random effects. We evaluated the utility of including smoothing functions (cubic regression splines) to examine temporal variation in the fixed effects in the final model for each species by comparing smoothed and unsmoothed versions of the explanatory variables Day and Year. Each candidate model was generalized using a log-link and Poisson distribution. The assumption of Poisson-distributed counts was validated through graphic comparison of results for each species using Poisson models and overdispersed models using the quasi-Poisson structure. Because little distinction was evident between these models, we proceeded with the assumption of Poisson-distributed counts, which facilitated model selection. To improve model clarity and reduce the potential for over-fitting, we used smoothing functions with a maximum basis dimension that was large enough that extrema were apparent, but small enough that curvature was simple [53]. These were maximum basis dimensions larger than half the number of time steps by year, but smaller than the total number of time steps by year. We then selected the final model as the best candidate model for each species given the data with Akaike's Information Criterion (AIC; [63], [65]). This model selection process was completed with the statistical computing language R version 2.14.0 [66].

Predictive GAMMs

We produced separate GAMMs from the final model for each species [54], [67] in R [66]. Sample R codes are provided in Appendix S2. Final model terms were fit using penalized quasi-likelihood methods in the R software package mgcv [68]. Library mgcv computes effective degrees of freedom for smoothed terms from the trace of the GAMM influence matrix, for computing AIC values [62]. We examined the within-season response of relative abundance to smoothed Day graphically. To examine variation across years, we extracted parameters from the GAMMs to populate predictive models, then fixed the value for Day at its median for each species, calculating trajectories as if they had all been sampled on the same day-of-winter. We compared trajectories with and without correction for variation in survey date. We calculated expected values and approximate confidence intervals using the R function predict [69]. Confidence intervals were estimated at plus and minus twice the standard error [62]. Because the interpretation of confidence intervals at multiple time periods within a repeated measures context is subject to debate ([52]; P. de Valpine, T.E. Ingersoll & W. Rainey, unpublished), we used the R function anova [66] to calculate p-values for the smoothed Day and smoothed Year terms, testing the null hypothesis of unchanging relative abundance over time.

Because our GAMMs estimate relative, rather than absolute abundance [61], we sought to avoid the perception that estimates of relative abundance were informative of absolute abundance. We therefore normalized estimates for after calculating trajectories [52], providing a common scale of relative abundance for all species. To avoid selection of an arbitrary baseline year from which to normalize counts and measure population changes over time, we calculated the relative abundance of a species by dividing predicted values by the maximum expected value for that species. Thus, our normalization procedure set the maximum abundance estimate for a species during the study period equal to 1.0. Because the GAMMs produced a complex series of additive terms, predicted relative abundances and confidence intervals were rendered graphically for ease of interpretation.

To evaluate the influence of bias from within-season survey date on estimates of long-term population trajectories, we compared corrected trajectories (which accounted for variable survey date) and naïve uncorrected trajectories (where survey date was not included in the model). Corrected trajectories were from the final models for each species, and models for uncorrected trajectories excluded the fixed effect for Day. In one species (P. subflavus), Day was not selected in the final model (see Results), so comparison was between the final model and the best alternate model that included Day.

Survey timing

The timing of hibernacula surveys varied substantially within the hibernation season in every year, and exhibited a systematic shift towards later dates (Fig. 1). Half of all surveys in each year were conducted during a 3–6 week period from late January to early March, but some surveys were conducted as early as mid-December and as late as the end of March. Within-season variability in survey date was highest in 2008. The median survey date was early February in most years, with a trend toward later date with time starting in 2006 (Fig. 1). In three years, median survey date was particularly late; ∼1 week later than other years in 2010, and ∼2 weeks later in both 2002 and 2011.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Figure 1. Timing of hibernation surveys across years.

Box plots showing date of hibernacula surveys during 1999–2011.

https://doi.org/10.1371/journal.pone.0065907.g001

Model selection

AIC model selection (Table 1) resulted in the selection of the fixed effects terms Day and Year for M. lucifugus, M. septentrionalis and M. sodalis. Smoothed Day was selected for these three species, smoothed Year was selected for M. lucifugus and linear Year was selected for M. septentrionalis and M. sodalis. For P. subflavus, linear Year alone was selected. Random intercept terms for Route were selected for all species, and Route within Location was selected for M. lucifugus.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Table 1. Model selection.

https://doi.org/10.1371/journal.pone.0065907.t001

Predictive GAMMs

Relative abundance varied non-linearly with survey date within the hibernation season, in all species except P. subflavus (smoothed Day terms, M. lucifugus, F_6.48 = 10.03, p<0.001; M. sodalis, F_5.90 = 21.88, p<0.001; M. septentrionalis, F_4.34 = 2.93, p = 0.033; Fig. 2; Appendix S3), suggesting that, if not accounted for, systematic heterogeneity in survey timing would bias relative abundance estimates in three of the four species. When comparing long-term trajectories corrected for survey timing (blue traces in Fig. 3) versus naïve uncorrected trajectories (red traces in Fig. 3), we found that uncorrected models underestimated relative abundance for M. lucifugus in most years (Fig. 3A). Declines in M. sodalis were underestimated in uncorrected models (Fig. 3C). Declines in M. septentrionalis were overestimated in uncorrected models (Fig. 3D). As expected, due to the insignificant effect of survey day for this species, trajectories from corrected versus uncorrected models for P. subflavus were nearly indistinguishable (Fig. 3B).

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Figure 2. Within-season temporal variation in bat counts.

Relative abundance and approximate 95% confidence intervals during December-March for (A) M. lucifugus, (B) P. subflavus, (C) M. sodalis, and (D) M. septentrionalis. Relative abundance was set equal to 1.0 at the maximum expected value.

https://doi.org/10.1371/journal.pone.0065907.g002

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Figure 3. Long-term population trajectories.

Expected relative abundance and approximate 95% confidence intervals during 1999–2011 for (A) M. lucifugus, (B) P. subflavus, (C) M. sodalis, and (D) M. septentrionalis. Relative abundance was set equal to 1.0 at the maximum expected value. Two trajectories are shown for each bat species: the trajectory with abundance estimates corrected for survey date of bat counts (in blue), and the uncorrected trajectory (red).

https://doi.org/10.1371/journal.pone.0065907.g003

Corrected estimates of regional relative abundance demonstrated marked declines across the study period for each species (smoothed Year terms, M. lucifugus, F_4.912 = 42.69, p<0.001; P. subflavus, F₁ = 30.02, p<0.001; M. sodalis, F₁ = 18.68, p<0.001; M. septentrionalis, F₁ = 3.90, p = 0.049; blue traces in Fig. 3; Appendix S3, S4). Some growth in the regional population of M. lucifugus was evident prior to 2004 (Fig. 3A; Appendix S4). Dynamics of the final models were remarkably similar among the other three species, exhibiting slow, steady declines in all cases (Figs. 3B, 3C, 3D; Appendix S4). Cumulative declines in regional relative abundance by 2011 from peak levels were 71% (with 95% confidence interval of ±11%) in M. lucifugus, 34% (±38%) in P. subflavus, 30% (±26%) in M. sodalis, and 31% (±18%) in M. septentrionalis (Appendix S4).

Population trajectories

Our results clearly show that within our four-state study area of the eastern United States, the regional populations of M. sodalis and M. septentrionalis were in decline (Fig. 3C, 3D; Appendices S3, S4) and the regional population of M. lucifugus was in sharp decline (Fig. 3A; Appendices S3, S4). Our data are strongly suggestive that the regional population of P. subflavus was also in decline (Fig. 3B; Year term in Appendix S3), though our estimates of decline for this species did not exceed 95% confidence intervals (Fig. 3B; Appendix S4). It is unclear why M. lucifugus would be declining faster than the other species, but each of these species varies in the location and environmental conditions of preferred foraging and hibernation sites [70], and these factors may be important in determining susceptibility to disturbance [12] and population-level responses to climate change [18] and disease [71]. As the most common hibernating species in this region, and one that clusters in large aggregations [70], M. lucifugus may also be particularly susceptible to disturbance [12] and bat-to-bat transmission of diseases such as white-nose syndrome [71], [72]. It also is the recipient of less conservation attention than the other species that aggregates in large clusters in the region, M. sodalis, which is federally-protected [45].

Population dynamics varied by species. We estimated that the relative abundance of M. lucifugus fluctuated slightly and gradually during the beginning of the study period (although still well within confidence intervals during this period), reached a peak in the study area during 2004, then declined severely and persistently thereafter (Fig. 3A; Appendix S4). Trajectories for all other species suggested a persistent declining trend throughout the study period (Fig. 3B, 3C, 3D; Appendix S3, S4). These different model dynamics may have occurred because the greater amount of data on M. lucifugus enabled us to determine finer-scale temporal variation in regional population trajectories. Alternatively, these dynamics may represent ecological differences among the species in response to environmental conditions or to ongoing threats. The pre-2004 increase in M. lucifugus could suggest a period of favorable environmental conditions against a background of overall decline, bat recovery after past declines, colonization of new winter habitat (perhaps associated with climate change), or changing levels of threat over time.

Our estimates of decline were less severe than previous estimates [22], [32]. This difference is likely because, unlike previous studies, our study analyzed data from hibernacula from across the region, including not only hibernacula infected with white-nose syndrome, but also hibernacula where white-nose syndrome had not been detected. Additionally, our estimates utilized the first multi-species, regional estimates of long-term trajectories in abundance of bat populations in North America to account for several factors believed to bias estimates of abundance from wildlife count data: non-linearity with time (Figs. 2, 3), non-Gaussian error distributions, variation among sites, and correlated errors [48], [49], [50]. In addition, although we were not able to explicitly estimate heterogeneous detection probability with these data, we did correct for an important source of detection heterogeneity [48], namely variation in sampling date (Figs. 1, 2, 3). Nonetheless, our results indicate considerable cause for concern for the regional populations of all four hibernating bat species studied. Although substantial uncertainty remains in estimates for some species, mean estimates of regional declines from their peak abundances were equal to or greater than 30% by 2011 in all species studied (Fig. 3; Appendix S4). The regional population of M. lucifugus is of particular concern, as mean estimates of decline from its peak abundance reached 71% by 2011 (Fig. 3A; Appendix S4). If trajectories continue along their current paths (Fig. 3), substantial further declines can be expected for each species in the near future [22]. Because of the critical ecological roles of bats in directly limiting populations of nocturnal arthropods through predation [3], [73], [74], indirectly limiting herbivory by predated species of arthropods in forests [73], [75] and agricultural systems [5], [76], and providing the nutrient subsidies [77] that are essential to supporting diverse and specialized cave fauna [78], the regional declines in the abundances of hibernating bat species that we observed could have important cascading effects in ecological communities.

Within-season variation in abundance

Unlike previous studies, our methods corrected long-term population trajectories for within-season variation. Such correction is needed to reduce bias, especially given the high variability in dates of hibernacula surveys, and recent systematic shifts to later date-of-winter for hibernacula surveys near the end of the survey period (Fig. 1). In P. subflavus, survey date had almost no effect on within-season estimates of relative abundance (Table 1B; Fig. 2B), and correction did not influence long-term estimates of relative abundance (Fig. 3B). However, survey date had an important effect on within-season estimates of relative abundance in the other three species (Tables 1A, 1C, 1D; Figs. 2A, 2C, 2D). As a result, our correction for survey date substantially changed estimates of rates of decline and cumulative amounts of decline in M. sodalis (Fig. 3C) and M. septentrionalis (Fig. 3D), and changed estimates of patterns and timing of population changes (but not cumulative amounts of decline over the study period) in M. lucifugus (Fig. 3A).

The influence of survey date on within-season and across-year estimates of relative abundance in these bat species may arise as a result of changes in the size of colonies within hibernacula over the course of the hibernation season. Hibernation imposes severe physiological challenges on bats, including progressive energy and evaporative water loss over the hibernation period [79], [80], and these challenges are exacerbated by disturbance in hibernacula [12] and white-nose syndrome [81]. As a result, progressive mortality is expected, and should lead to a gradual reduction of the colony (and declines in regional relative abundance) over the hibernation period [24]. This could explain the within-season declines we observed throughout the hibernation period in M. sodalis (Fig. 2C), and the within-season declines late in the hibernation period in M. septentrionalis (Fig. 2D). In addition, bat species vary in the dates and duration of hibernation [55], frequency of between-hibernacula movements [57], and propensity to emerge from hibernation during intervals of warm weather during winter [57]. Variation among bat species in hibernation timing would cause a portion of the population to be present outside of the hibernacula, and thus be undetectable during hibernation surveys. Such variation could result in the non-linearities in within-season patterns of relative abundance that we observed in three of the four species in this study (Tables 1A, 1C, 1D; Figs. 2A, 2C, 2D).

The bat species' differing roosting habits during hibernation may also affect detectability, which could also help explain the influence of within-season variation in survey date on bat counts. In particular, the small size, distinctive forearms, and style of roosting away from large aggregations make P. subflavus easy to identify, and its habit of roosting as singles or in very small clusters hanging on walls or ceilings make them relatively easy to count throughout the hibernation season [82]. In contrast, the roosting behavior of the other three species may vary throughout the hibernation season in ways that affect detectability and thus bat counts. M. septentrionalis also roosts as singles or small groups, but tends to seek out crevices in walls or ceilings. Because differences in temperature or environmental factors may affect how deep they move into crevices [83], this species may vary in detectability across the hibernation season, particularly in hibernacula with high ceilings or deep crevices. M. lucifugus and M. sodalis, in contrast, often hibernate in large aggregations. Although they typically form large clusters (10's to 1000's) that are easy to locate, these clusters may vary considerably in density of packing throughout the hibernation season [82], [84], [85], perhaps also due to temperature or other environmental factors [83]. These species are also typically found on high ceilings, which makes estimating such within-season changes in cluster density difficult unless high-resolution photography is used [51].

Threats contributing to bat population declines

It is not possible to definitively attribute population declines to specific causes solely from count data, but the timing of changes in population trajectories may be suggestive of key threats associated with these changes [52]. White-nose syndrome has often been assumed to be the key cause of regional population declines in hibernating bat species in the region we studied, due to repeated observations of mass mortality events associated with the disease [22], [32], [34], [86]. Our population trajectories, indicating persistent declines in each species in recent years (Fig. 3; Appendix S3, S4) are consistent with the interpretation that white-nose syndrome is a severe and pervasive threat that has led to extensive declines in regional bat abundance.

Our results (Fig. 3; Appendix S4), however, indicate changes in regional bat abundances throughout the study period, including declines in all species that began prior to the wide proliferation of white-nose syndrome in the region beginning in 2008–2010, and prior even to the first detection of white-nose syndrome in 2006 [22], [34]. In M. lucifugus, we observed declines that began in 2004, two years prior to the first recorded instance of the disease in the region [34]. AIC-based inference in which linear models were selected suggests that declines in all other species persisted across the study period, and so preceded discovery of WNS by at least seven years. Furthermore, insufficient evidence existed to support a conclusion that rates of decline increased in the other species following 2006 (Tables 1B, 1C, 1D). Mass mortality events due to white-nose syndrome are associated with large numbers of bat carcasses in hibernacula that should have been evident during hibernacula surveys [22], so it is unlikely that the epizootic went undetected for several years prior to 2006. Rather, additional persistent threats likely contributed substantially to recent declines.

Which other threats could have caused recent bat population declines remains unclear. Bat mortalities from in-flight collisions with turbines have increased with expanding wind energy development [20], [21] since 2000 [87], but these have been primarily associated with fatalities of migratory, foliage- and tree-roosting bats rather than hibernating bats [20]. Recent climate changes could also have played a role [17], [18], but climate changes have been gradual in the region for at least thirty years [18], including during periods in which bat populations were increasing (Fig. 3A). Other long-standing threats, such as changes to critical roosting [10], [11], [12] or foraging habitat [9], [13], [14], in-flight collisions with vehicles or buildings [19], or effects of pesticides and other chemical compounds [15], [16], [26] may also have influenced regional bat abundances, though it is unclear whether any single threat could have produced the persistent declines we observed (Fig. 3). The temporally-variable trajectory of M. lucifugus in particular may reflect synergistic or interacting effects of multiple threats [52], perhaps in conjunction with bat population responses to variable environmental conditions [17]. Further study is needed to clarify the relative importance of each threat for bat populations.

Implications of regional declines in bat abundances for bat conservation status

M. lucifugus, P. subflavus, and M. septentrionalis are not currently federally-protected in the United States, whereas M. sodalis is listed as Endangered [88], [89], the most protected designation under the federal Endangered Species Act. Several proposals are currently being evaluated for listing the three currently-unprotected species at the state or federal level within the United States (e.g., [33], [90], [91]). We do not evaluate the merits of these proposals here. However, our results indicating multi-year, continuing declines suggest a need for enhanced conservation efforts for all four bat species, and for hibernating bat species more generally.

The IUCN Red List thresholds for assigning Vulnerable or Endangered conservation status are, respectively, 30% or 50% declines in population size, where causes of decline may not have ceased or may not be understood [92]. According to the IUCN Red List, M. sodalis is currently Endangered with a decreasing population trend [30], and M. lucifugus, P. subflavus, and M. septentrionalis all currently have an IUCN conservation status of Least Concern with a stable population trend [30]. Our estimates (Fig. 3; Appendix S3, S4), however, provide new evidence in our study area of regional-scale declines in these species, and especially in M. lucifugus. Specifically, in our study area, M. lucifugus exceeded the 30% (Vulnerable) threshold of population decline by 2008 and the 50% (Endangered) threshold by 2009. In both cases, the upper confidence limit on M. lucifugus abundance estimates exceeded these thresholds in the following year, indicating substantial certainty in these results. The 30% threshold of population decline was exceeded in our study area by 2010 in P. subflavus, and by 2011 in M. sodalis and M. septentrionalis, though upper confidence limits have not yet exceeded the 30% threshold in these three species. Note, however, that the declines in M. sodalis are on top of persistent historical declines that occurred in the species up through 2003 [30], [31], [89].

The implications of our regional-scale results for the conservation status of each species as a whole depend on whether the declines we observed in our four-state study area are pervasive across each species' geographic range. The species M. sodalis [93] and M. septentrionalis [94] and the subspecies M. lucifugus lucifugus [95] and P. subflavus subflavus [96] occur principally in eastern and northern North America, where many of the threats we identified above are widespread and increasing [7], [9]. This includes white-nose syndrome, which has already spread from a single known hibernaculum to 22 eastern, southeastern, and Midwestern states and five eastern Canadian provinces in only seven years [97], and has the potential to rapidly spread across the remaining entire geographic ranges of each taxon studied here [93], [94], [95], [96]. Thus our results likely have strong implications for conservation status at the species or subspecies level.

Sampling approach

The hibernacula surveys we used in this study and others conducted by state wildlife-monitoring programs are the most reliable and consistent datasets currently available for long-term, regional studies of North American bats. Sampling effort across the region each year is extensive, and the long-term nature of these data provides an important historical record of population trends in these bat species. Despite the great utility of this dataset, however, it exhibited several potential sampling problems that are typical of long-term regional monitoring programs [41]. It also exhibited potential sampling problems that likely result from the cryptic nature of bats and the difficulty of monitoring them.

In this paper, we addressed several of these potential problems that would otherwise have biased our results. However, we were unable to account for other potential problems, in part due to limitations of the dataset. For instance, the survey routes in the state monitoring programs that provided data for this study were established in known hibernacula where hibernating bats could be reliably observed, and as a result easily-surveyed and consistently-used hibernacula may be overrepresented in the data. The extent and net impact of such potential sampling bias (towards over or underestimation) from the use of these particular sites for monitoring cannot be determined from our data. In addition, due to concerns about disturbance of bats during the hibernation period [12], surveys were not repeated on the same survey route in the same season (i.e., at most one survey was conducted per route per year). This aspect of the dataset limited our ability to examine interactions between temporal variables and to improve confidence estimates via modeling of detection probabilities. Also, our ability to identify key threats associated with population changes was limited by a lack of data collected concurrently on factors that may be correlated with bat counts. This limited us to using counts as an index of abundance, an approach that could affect interpretation of abundance changes if detection probability for a species changed over time. Finally, although we used a 13-year dataset, further historical data would have been beneficial to fully characterize population dynamics prior to the onset of recent declines and to more precisely identify the baseline from which to calculate cumulative declines.

Such potential sampling problems could be addressed through changes in the selection of survey routes and improvements in the survey methods employed in the monitoring programs. Questions about representativeness could be addressed through comprehensive efforts to document all hibernacula followed by implementation of random sampling of hibernacula within the entire set of hibernacula in the region. Improved survey methods could also be developed that use limited-disturbance survey methods (e.g., thermal imaging video recording) that would allow some within-season repeated measures to improve estimation of bat detectability [61]. Improved survey methods could include the use of double-observer methods for bat counts during hibernation surveys that distinguish between observer error and natural variability [49], [98], [99]. Survey efforts could also be expanded to include the recording of standardized detection covariates (e.g., hours of effort, length of survey route, and distance from observer to bat clusters) during hibernation surveys [49], [61]. This information could increase the accuracy and precision of estimates of bat abundance, and facilitate efforts to distinguish the relative importance of multiple threats driving declines. Continued monitoring efforts are also needed to improve the estimation of population dynamics over time and directly and accurately document further changes in conservation status. Other steps that could increase the accuracy and consistency of hibernation surveys are discussed in detail elsewhere (e.g., [41]).

Modeling approach

Some disadvantages of our modeling approach were also apparent. Our results suggest that, despite our use of data from 636 surveys from 163 survey routes across four states, regional population estimates were subject to substantial uncertainties (Fig. 3). In addition, although temporal smoothing assisted with identifying long-term regional patterns in highly-variable count data, the smoothing could have had the side effect of obscuring single-year changes in population trajectories. Other factors than the observation heterogeneities we accounted for in this study, such as within- and among-year climatic variation, could also have affected bat detectability and thus estimates of population abundance. In addition, our modeling approach would need to be modified before being applied to datasets that exhibit significant deviations from the Poisson assumption such as overdispersion.

Despite these caveats, our modeling approach provided several important benefits. First, we were able to provide estimates of not only relative abundance but also uncertainty (Fig. 3; Appendix S3, S4). Both of these elements are essential to effective management decisions, but few studies of bat populations have presented estimates of uncertainty [40]. Our confidence intervals are wide in some places (Fig. 3; Appendix S4), reflecting uncertainty in point estimates due to the high spatial variability and irruptive temporal variability (P. de Valpine, T.E. Ingersoll & W. Rainey, unpublished) present in bat counts at individual hibernacula. Methods for confidence interval estimation and interpretation for time series data in GAMMs are incompletely developed, and should be the focus of future study (P. de Valpine, T.E. Ingersoll & W. Rainey, unpublished). Nonetheless, our results demonstrated strong evidence of change at the regional scale in the abundance of each species over time (Fig. 3; Year or smoothed Year terms, Appendix S3). Thus, our modeling approach reduced bias and improved interpretability of long-term regional population changes (Fig. 3), despite high short-term and local variability.

Second, our modeling approach provided estimates over time – not just before/after snapshots – while accommodating non-linearities in population abundances over time. This aspect of the modeling approach is important because there is no reason to assume, as is routine [49], that changes in population abundance over time will be linear [53]. Rather, non-linearity is especially likely to be observed in long-term studies of populations that may be responding to the combined influences of multiple environmental factors and potentially complex suites of threats. Our use of a modeling approach that accommodated non-linearity enabled us to identify time-dependent changes in the population trajectories and clarify the timing and extent of increases and declines in population abundance, allowing inference of changes in trend.

Third, trajectory estimates over time accounted for non-independent repeated-measures along sampling routes. The selection of ‘Route within location’ for M. lucifugus (Table 1A) but not the other species (Tables 1B, 1C, 1D) likely reflects both environmental preferences for roost sites and sample size differences. M. lucifugus appeared more commonly in more complex hibernacula that were surveyed with several routes. In addition, the greater amount of data for M. lucifugus (colony sizes as well as numbers of surveys and routes in which it was detected) enabled route-specific differences to be detected for this species in these more complex hibernacula that were surveyed with several routes. More generally, the use of the random grouping term for location enabled us to relax independence assumptions, allowing us to integrate the survey data from diverse sites – including data from the single routes that are sufficient to survey simple hibernacula and data from the multiple, non-independent routes that are needed to adequately survey complex hibernacula – in the same analysis. It also improved estimation of temporal changes in abundance by accounting for some of the detection heterogeneity associated with survey route.

Fourth, our modeling approach controlled for observation heterogeneities. Controlling for changes in within-year variation in survey date was particularly important, due to the potential of such changes to systematically bias results over time. For example, we observed a shift to later survey dates beginning near the end of the study period (Fig. 1). According to survey notes associated with hibernacula surveys, this shift was an intentional effort by managers of the monitoring programs to facilitate the detection of white-nose syndrome (which is easier to observe later in the hibernation season; T.E. Ingersoll, unpublished) after its discovery in 2007. However, such variation in survey dates (Fig. 2) could combine with systematic, within-season changes in bat abundances [24] and movements within and among hibernacula [55], [56], [57] to affect long-term, regional estimates of relative abundance (Fig. 3). Specifically, our models suggested that, if not accounted for in models, changing survey dates would have biased abundance estimates for M. lucifugus, M. sodalis, and M. septentrionalis (Fig. 3A, 3C, 3D; Appendix S3). Our focus on addressing observation heterogeneities explicitly addressed this potential source of bias and thereby improved the accuracy of our results.

Thus, our modeling approach proved useful in clarifying long-term, regional trajectories of populations of hibernating bat species on the basis of count data, by providing estimates of uncertainty, accommodating non-linearities, and accounting for observation heterogeneities. Such factors are often important in studies of population abundances [40], [48], [49], [50], [54], and thus our approach has wide potential applicability to other studies of wildlife populations.