Integrating multiple sign types to improve occupancy estimation for inconspicuous species.

Standard occupancy models enable unbiased estimation of occupancy by accounting for observation errors such as missed detections (false negatives) and, less commonly, incorrect detections (false positives). Occupancy models are fitted to data from repeated site visits in which surveyors record evidence of species presence. Use of indirect sign (e.g., scat, tracks) as evidence of presence can vastly improve survey efficiency for inconspicuous species but can also introduce additional sources of error. We developed a "multi-sign" occupancy approach to model the detection process separately for unique sign types and used this method to improve estimates of occupancy dynamics for an inconspicuous species, the American pika (Ochotona princeps). We investigated how estimates of pika occupancy and environmental drivers differed under four increasingly realistic representations of the observation process: (1) perfect detection (commonly assumed for modeling pika occupancy), (2) standard occupancy model (single observation process without possibility of false detection), (3) multi-sign with no false detections (non-false positive model), and (4) multi-sign with false detections (full model). For the multi-sign occupancy models, we modeled the detection of each sign type (fresh scat, fresh haypiles, pika calls, and pika sightings) separately as a function of climatic and environmental covariates. Estimates of occupancy processes and inferences about environmental drivers were sensitive to different detection models. Simplified representations of the detection processes generally resulted in higher occupancy estimates and higher turnover rates than the full multi-sign model. Environmental drivers also varied in their influence on occupancy models, where (e.g.) forb cover was estimated to more strongly influence occupancy in the full multi-sign model than the simpler models. As has been reported previously in other contexts, unmodeled heterogeneity in the observation process can lead to biases in occupancy processes and uncertainty in the relationships between occupancy and environmental covariates. Overall, our multi-sign approach to dynamic occupancy modeling, which accounts for spatio-temporal variation in reliability among sign types, has strong potential to generate more realistic estimates of occupancy dynamics for inconspicuous species.

Effects of traditional harvest and burning on common camas (Camassia quamash) abundance in Northern Idaho: The potential for traditional resource management in a protected area wetland.

The bulbs of common camas (Camassia quamash) were a staple food of Indigenous Peoples of western North America for millennia. Camas harvesting site productivity was encouraged through intense management. Common camas is considered a facultative wetland species, and populations have declined due to contemporary wetland drainage and land conversion. Conservation of existing habitat, as well as restoration of degraded systems, is necessary. Traditional ecological knowledge (TEK) and resource management (TRM) are often promoted as viable modes of contemporary resource management but are rarely tested or implemented. We designed a controlled experiment, informed by a born-in-the-tradition specialist, to evaluate the response of common camas populations to traditional bulb harvest, burning, and a combination of harvest and burning. We recorded camas plant counts of three life stage classes of camas plants (single-leaf seedling, multiple-leaf adult, and flowering adult) over the course of 6 years in arrays of plots and subjected to treatments or left undisturbed as control. Harvesting removed plants (>800 bulbs) and reduced aboveground counts of camas densities ( X ¯  ~ 50% of control, p < .05). Burning contributed to a reduction in single-leaf plants but had an overall positive effect ( X ¯  ~ 150% of control, p < .05) on adult camas and flowering plant abundance, and ameliorated the digging impacts. Treatment impacts tapered over the course of the study, and results indicate that a sustainable harvesting return interval of approximately 5 years may be possible when combined with fire to reduce litter and competition from pasture grasses and to accelerate the recovery of camas. Our findings support the hypotheses proposed by traditional knowledge specialists and ethnobotanists that digging and burning reduce intra- and interspecific competition and stimulate the growth of unharvested adult plants. More generally, our study supports the integration of traditional ecological knowledge into the evidence base available for protected area wetland prairie management.

NABat: A top-down, bottom-up solution to collaborative continental-scale monitoring.

Collaborative monitoring over broad scales and levels of ecological organization can inform conservation efforts necessary to address the contemporary biodiversity crisis. An important challenge to collaborative monitoring is motivating local engagement with enough buy-in from stakeholders while providing adequate top-down direction for scientific rigor, quality control, and coordination. Collaborative monitoring must reconcile this inherent tension between top-down control and bottom-up engagement. Highly mobile and cryptic taxa, such as bats, present a particularly acute challenge. Given their scale of movement, complex life histories, and rapidly expanding threats, understanding population trends of bats requires coordinated broad-scale collaborative monitoring. The North American Bat Monitoring Program (NABat) reconciles top-down, bottom-up tension with a hierarchical master sample survey design, integrated data analysis, dynamic data curation, regional monitoring hubs, and knowledge delivery through web-based infrastructure. NABat supports collaborative monitoring across spatial and organizational scales and the full annual lifecycle of bats.

Evidence of region-wide bat population decline from long-term monitoring and Bayesian occupancy models with empirically informed priors.

Strategic conservation efforts for cryptic species, especially bats, are hindered by limited understanding of distribution and population trends. Integrating long-term encounter surveys with multi-season occupancy models provides a solution whereby inferences about changing occupancy probabilities and latent changes in abundance can be supported. When harnessed to a Bayesian inferential paradigm, this modeling framework offers flexibility for conservation programs that need to update prior model-based understanding about at-risk species with new data. This scenario is exemplified by a bat monitoring program in the Pacific Northwestern United States in which results from 8 years of surveys from 2003 to 2010 require updating with new data from 2016 to 2018. The new data were collected after the arrival of bat white-nose syndrome and expansion of wind power generation, stressors expected to cause population declines in at least two vulnerable species, little brown bat (Myotis lucifugus) and the hoary bat (Lasiurus cinereus). We used multi-season occupancy models with empirically informed prior distributions drawn from previous occupancy results (2003-2010) to assess evidence of contemporary decline in these two species. Empirically informed priors provided the bridge across the two monitoring periods and increased precision of parameter posterior distributions, but did not alter inferences relative to use of vague priors. We found evidence of region-wide summertime decline for the hoary bat ( λ ^  = 0.86 ± 0.10) since 2010, but no evidence of decline for the little brown bat ( λ ^  = 1.1 ± 0.10). White-nose syndrome was documented in the region in 2016 and may not yet have caused regional impact to the little brown bat. However, our discovery of hoary bat decline is consistent with the hypothesis that the longer duration and greater geographic extent of the wind energy stressor (collision and barotrauma) have impacted the species. These hypotheses can be evaluated and updated over time within our framework of pre-post impact monitoring and modeling. Our approach provides the foundation for a strategic evidence-based conservation system and contributes to a growing preponderance of evidence from multiple lines of inquiry that bat species are declining.

A review of bat hibernacula across the western United States: Implications for white-nose syndrome surveillance and management.

Efforts to conserve bats in the western United States have long been impeded by a lack of information on their winter whereabouts, particularly bats in the genus Myotis. The recent arrival of white-nose syndrome in western North America has increased the urgency to characterize winter roost habitats in this region. We compiled 4,549 winter bat survey records from 2,888 unique structures across 11 western states. Myotis bats were reported from 18.5% of structures with 95% of aggregations composed of ≤10 individuals. Only 11 structures contained ≥100 Myotis individuals and 6 contained ≥500 individuals. Townsend's big-eared bat (Corynorhinus townsendii) were reported from 38% of structures, with 72% of aggregations composed of ≤10 individuals. Aggregations of ≥100 Townsend's big-eared bats were observed at 41 different caves or mines across 9 states. We used zero-inflated negative binomial regression to explore biogeographic patterns of winter roost counts. Myotis counts were greater in caves than mines, in more recent years, and in more easterly longitudes, northerly latitudes, higher elevations, and in areas with higher surface temperatures and lower precipitation. Townsend's big-eared bat counts were greater in caves, during more recent years, and in more westerly longitudes. Karst topography was associated with higher Townsend's big-eared bat counts but did not appear to influence Myotis counts. We found stable or slightly-increasing trends over time in counts for both Myotis and Townsend's big-eared bats from 82 hibernacula surveyed ≥5 winters since 1990. Highly-dispersed winter roosting of Myotis in the western USA complicates efforts to monitor population trends and impacts of disease. However, our results reveal opportunities to monitor winter population status of Townsend's big-eared bats across this region.

Improving geographically extensive acoustic survey designs for modeling species occurrence with imperfect detection and misidentification.

Acoustic recording units (ARUs) enable geographically extensive surveys of sensitive and elusive species. However, a hidden cost of using ARU data for modeling species occupancy is that prohibitive amounts of human verification may be required to correct species identifications made from automated software. Bat acoustic studies exemplify this challenge because large volumes of echolocation calls could be recorded and automatically classified to species. The standard occupancy model requires aggregating verified recordings to construct confirmed detection/non-detection datasets. The multistep data processing workflow is not necessarily transparent nor consistent among studies. We share a workflow diagramming strategy that could provide coherency among practitioners. A false-positive occupancy model is explored that accounts for misclassification errors and enables potential reduction in the number of confirmed detections. Simulations informed by real data were used to evaluate how much confirmation effort could be reduced without sacrificing site occupancy and detection error estimator bias and precision. We found even under a 50% reduction in total confirmation effort, estimator properties were reasonable for our assumed survey design, species-specific parameter values, and desired precision. For transferability, a fully documented r package, OCacoustic, for implementing a false-positive occupancy model is provided. Practitioners can apply OCacoustic to optimize their own study design (required sample sizes, number of visits, and confirmation scenarios) for properly implementing a false-positive occupancy model with bat or other wildlife acoustic data. Additionally, our work highlights the importance of clearly defining research objectives and data processing strategies at the outset to align the study design with desired statistical inferences.

Occupancy modeling species-environment relationships with non-ignorable survey designs.

Statistical models supporting inferences about species occurrence patterns in relation to environmental gradients are fundamental to ecology and conservation biology. A common implicit assumption is that the sampling design is ignorable and does not need to be formally accounted for in analyses. The analyst assumes data are representative of the desired population and statistical modeling proceeds. However, if data sets from probability and non-probability surveys are combined or unequal selection probabilities are used, the design may be non-ignorable. We outline the use of pseudo-maximum likelihood estimation for site-occupancy models to account for such non-ignorable survey designs. This estimation method accounts for the survey design by properly weighting the pseudo-likelihood equation. In our empirical example, legacy and newer randomly selected locations were surveyed for bats to bridge a historic statewide effort with an ongoing nationwide program. We provide a worked example using bat acoustic detection/non-detection data and show how analysts can diagnose whether their design is ignorable. Using simulations we assessed whether our approach is viable for modeling data sets composed of sites contributed outside of a probability design. Pseudo-maximum likelihood estimates differed from the usual maximum likelihood occupancy estimates for some bat species. Using simulations we show the maximum likelihood estimator of species-environment relationships with non-ignorable sampling designs was biased, whereas the pseudo-likelihood estimator was design unbiased. However, in our simulation study the designs composed of a large proportion of legacy or non-probability sites resulted in estimation issues for standard errors. These issues were likely a result of highly variable weights confounded by small sample sizes (5% or 10% sampling intensity and four revisits). Aggregating data sets from multiple sources logically supports larger sample sizes and potentially increases spatial extents for statistical inferences. Our results suggest that ignoring the mechanism for how locations were selected for data collection (e.g., the sampling design) could result in erroneous model-based conclusions. Therefore, in order to ensure robust and defensible recommendations for evidence-based conservation decision-making, the survey design information in addition to the data themselves must be available for analysts. Details for constructing the weights used in estimation and code for implementation are provided.

Variation in subsurface thermal characteristics of microrefuges used by range core and peripheral populations of the American pika (Ochotona princeps).

Microrefuges provide microclimates decoupled from inhospitable regional climate regimes that enable range-peripheral populations to persist and are important to cold-adapted species in an era of accelerated climate change. However, identifying and describing the thermal characteristics of microrefuge habitats is challenging, particularly for mobile organisms in cryptic, patchy habitats. We examined variation in subsurface thermal conditions of microrefuge habitats among different rock substrate types used by the American pika (Ochotona princeps), a climate-sensitive, rock-dwelling Lagomorph. We compared subsurface temperatures in talus and lava substrates in pika survey sites in two US national park units; one park study area on the range periphery and the other in the range core. We deployed paired sensors to examine within-site temperature variation. We hypothesized that subsurface temperatures within occupied sites and structurally complex substrates would be cooler in summer and warmer in winter than unoccupied and less complex sites. Although within-site variability was high, with correlations between paired sensors as low as 47%, we found compelling evidence that pikas occupy microrefuge habitats where subsurface conditions provide more thermal stability than in unoccupied microhabitats. The percentage of days in which microhabitat temperatures were between -2.5 and 25.5°C was significantly higher in occupied sites. Interestingly, thermal conditions were substantially more stable (p < .05) in the lava substrate type identified to be preferentially used by pikas (pahoehoe vs. a'a) in a previous study. Our study and others suggest that thermal stability appears to be the defining characteristic of subsurface microrefuges used by American pikas and is a likely explanation for enigmatic population persistence at the range periphery. Our study exemplifies an integrated approach for studying complex microhabitat conditions, paired with site use surveys and contextualized with information about gene flow provided by complementary studies.

Replicated landscape genetic and network analyses reveal wide variation in functional connectivity for American pikas.

Landscape connectivity is essential for maintaining viable populations, particularly for species restricted to fragmented habitats or naturally arrayed in metapopulations and facing rapid climate change. The importance of assessing both structural connectivity (physical distribution of favorable habitat patches) and functional connectivity (how species move among habitat patches) for managing such species is well understood. However, the degree to which functional connectivity for a species varies among landscapes, and the resulting implications for conservation, have rarely been assessed. We used a landscape genetics approach to evaluate resistance to gene flow and, thus, to determine how landscape and climate-related variables influence gene flow for American pikas (Ochotona princeps) in eight federally managed sites in the western United States. We used empirically derived, individual-based landscape resistance models in conjunction with predictive occupancy models to generate patch-based network models describing functional landscape connectivity. Metareplication across landscapes enabled identification of limiting factors for dispersal that would not otherwise have been apparent. Despite the cool microclimates characteristic of pika habitat, south-facing aspects consistently represented higher resistance to movement, supporting the previous hypothesis that exposure to relatively high temperatures may limit dispersal in American pikas. We found that other barriers to dispersal included areas with a high degree of topographic relief, such as cliffs and ravines, as well as streams and distances greater than 1-4 km depending on the site. Using the empirically derived network models of habitat patch connectivity, we identified habitat patches that were likely disproportionately important for maintaining functional connectivity, areas in which habitat appeared fragmented, and locations that could be targeted for management actions to improve functional connectivity. We concluded that climate change, besides influencing patch occupancy as predicted by other studies, may alter landscape resistance for pikas, thereby influencing functional connectivity through multiple pathways simultaneously. Spatial autocorrelation among genotypes varied across study sites and was largest where habitat was most dispersed, suggesting that dispersal distances increased with habitat fragmentation, up to a point. This study demonstrates how landscape features linked to climate can affect functional connectivity for species with naturally fragmented distributions, and reinforces the importance of replicating studies across landscapes.

A goodness-of-fit test for occupancy models with correlated within-season revisits.

Occupancy modeling is important for exploring species distribution patterns and for conservation monitoring. Within this framework, explicit attention is given to species detection probabilities estimated from replicate surveys to sample units. A central assumption is that replicate surveys are independent Bernoulli trials, but this assumption becomes untenable when ecologists serially deploy remote cameras and acoustic recording devices over days and weeks to survey rare and elusive animals. Proposed solutions involve modifying the detection-level component of the model (e.g., first-order Markov covariate). Evaluating whether a model sufficiently accounts for correlation is imperative, but clear guidance for practitioners is lacking. Currently, an omnibus goodness-of-fit test using a chi-square discrepancy measure on unique detection histories is available for occupancy models (MacKenzie and Bailey, Journal of Agricultural, Biological, and Environmental Statistics, 9, 2004, 300; hereafter, MacKenzie-Bailey test). We propose a join count summary measure adapted from spatial statistics to directly assess correlation after fitting a model. We motivate our work with a dataset of multinight bat call recordings from a pilot study for the North American Bat Monitoring Program. We found in simulations that our join count test was more reliable than the MacKenzie-Bailey test for detecting inadequacy of a model that assumed independence, particularly when serial correlation was low to moderate. A model that included a Markov-structured detection-level covariate produced unbiased occupancy estimates except in the presence of strong serial correlation and a revisit design consisting only of temporal replicates. When applied to two common bat species, our approach illustrates that sophisticated models do not guarantee adequate fit to real data, underscoring the importance of model assessment. Our join count test provides a widely applicable goodness-of-fit test and specifically evaluates occupancy model lack of fit related to correlation among detections within a sample unit. Our diagnostic tool is available for practitioners that serially deploy survey equipment as a way to achieve cost savings.

Habitat availability and gene flow influence diverging local population trajectories under scenarios of climate change: a place-based approach.

Ecological niche theory holds that species distributions are shaped by a large and complex suite of interacting factors. Species distribution models (SDMs) are increasingly used to describe species' niches and predict the effects of future environmental change, including climate change. Currently, SDMs often fail to capture the complexity of species' niches, resulting in predictions that are generally limited to climate-occupancy interactions. Here, we explore the potential impact of climate change on the American pika using a replicated place-based approach that incorporates climate, gene flow, habitat configuration, and microhabitat complexity into SDMs. Using contemporary presence-absence data from occupancy surveys, genetic data to infer connectivity between habitat patches, and 21 environmental niche variables, we built separate SDMs for pika populations inhabiting eight US National Park Service units representing the habitat and climatic breadth of the species across the western United States. We then predicted occurrence probability under current (1981-2010) and three future time periods (out to 2100). Occurrence probabilities and the relative importance of predictor variables varied widely among study areas, revealing important local-scale differences in the realized niche of the American pika. This variation resulted in diverse and - in some cases - highly divergent future potential occupancy patterns for pikas, ranging from complete extirpation in some study areas to stable occupancy patterns in others. Habitat composition and connectivity, which are rarely incorporated in SDM projections, were influential in predicting pika occupancy in all study areas and frequently outranked climate variables. Our findings illustrate the importance of a place-based approach to species distribution modeling that includes fine-scale factors when assessing current and future climate impacts on species' distributions, especially when predictions are intended to manage and conserve species of concern within individual protected areas.

Woodland dynamics at the northern range periphery: a challenge for protected area management in a changing world.

Managers of protected natural areas increasingly are confronted with novel ecological conditions and conflicting objectives to preserve the past while fostering resilience for an uncertain future. This dilemma may be pronounced at range peripheries where rates of change are accelerated and ongoing invasions often are perceived as threats to local ecosystems. We provide an example from City of Rocks National Reserve (CIRO) in southern Idaho, positioned at the northern range periphery of pinyon-juniper (P-J) woodland. Reserve managers are concerned about P-J woodland encroachment into adjacent sagebrush steppe, but the rates and biophysical variability of encroachment are not well documented and management options are not well understood. We quantified the rate and extent of woodland change between 1950 and 2009 based on a random sample of aerial photo interpretation plots distributed across biophysical gradients. Our study revealed that woodland cover remained at approximately 20% of the study area over the 59-year period. In the absence of disturbance, P-J woodlands exhibited the highest rate of increase among vegetation types at 0.37% yr(-1). Overall, late-successional P-J stands increased in area by over 100% through the process of densification (infilling). However, wildfires during the period resulted in a net decrease of woody evergreen vegetation, particularly among early and mid-successional P-J stands. Elevated wildfire risk associated with expanding novel annual grasslands and drought is likely to continue to be a fundamental driver of change in CIRO woodlands. Because P-J woodlands contribute to regional biodiversity and may contract at trailing edges with global warming, CIRO may become important to P-J woodland conservation in the future. Our study provides a widely applicable toolset for assessing woodland ecotone dynamics that can help managers reconcile the competing demands to maintain historical fidelity and contribute meaningfully to the U.S. protected area network in a future with novel, no-analog ecosystems.

The idiosyncrasies of place: geographic variation in the climate-distribution relationships of the American pika.

Although climate acts as a fundamental constraint on the distribution of organisms, understanding how this relationship between climate and distribution varies over a species' range is critical for addressing the potential impacts of accelerated climate change on biodiversity. Bioclimatic niche models provide compelling evidence that many species will experience range shifts under scenarios of global change, yet these broad, macroecological perspectives lack specificity at local scales, where unique combinations of environment, biota, and history conspire against generalizations. We explored how these idiosyncrasies of place affect the climate-distribution relationship of the American pika (Ochotona princeps) by replicating intensive field surveys across bioclimatic gradients in eight U.S. national parks. At macroecological scales, the importance of climate as a constraint on pika distribution appears unequivocal; forecasts suggest that the species' range will contract sharply in coming decades. However, the species persists outside of its modeled bioclimatic envelope in many locations, fueling uncertainty and debate over its conservation status. Using a Bayesian hierarchical approach, we modeled variation in local patterns of pika distribution along topographic position, vegetation cover, elevation, temperature, and precipitation gradients in each park landscape. We also accounted for annual turnover in site occupancy probabilities. Topographic position and vegetation cover influenced occurrence in all parks. After accounting for these factors, pika occurrence varied widely among parks along bioclimatic gradients. Precipitation by itself was not a particularly influential predictor. However, measures of heat stress appeared most influential in the driest parks, suggesting an interaction between the strength of climate effects and the position of parks along precipitation gradients. The combination of high elevation, cold temperatures, and high precipitation lowered occurrence probabilities in some parks, suggesting an upper elevational limit for pikas in some environments. Our results demonstrate that the idiosyncrasies of place influence both the nature and strength of the climate-distribution relationship for the American pika. Fine-grained, but geographically extensive, studies replicated across multiple landscapes offer insights important to assessing the impacts of climate change that otherwise may be masked at macroecological scales. The hierarchical approach to modeling provides a coherent conceptual and technical framework for gaining these insights.

Assessing the status and trend of bat populations across broad geographic regions with dynamic distribution models.

Bats face unprecedented threats from habitat loss, climate change, disease, and wind power development, and populations of many species are in decline. A better ability to quantify bat population status and trend is urgently needed in order to develop effective conservation strategies. We used a Bayesian autoregressive approach to develop dynamic distribution models for Myotis lucifugus, the little brown bat, across a large portion of northwestern USA, using a four-year detection history matrix obtained from a regional monitoring program. This widespread and abundant species has experienced precipitous local population declines in northeastern USA resulting from the novel disease white-nose syndrome, and is facing likely range-wide declines. Our models were temporally dynamic and accounted for imperfect detection. Drawing on species-energy theory, we included measures of net primary productivity (NPP) and forest cover in models, predicting that M. lucifugus occurrence probabilities would covary positively along those gradients. Despite its common status, M. lucifugus was only detected during -50% of the surveys in occupied sample units. The overall naive estimate for the proportion of the study region occupied by the species was 0.69, but after accounting for imperfect detection, this increased to -0.90. Our models provide evidence of an association between NPP and forest cover and M. lucifugus distribution, with implications for the projected effects of accelerated climate change in the region, which include net aridification as snowpack and stream flows decline. Annual turnover, the probability that an occupied sample unit was a newly occupied one, was estimated to be low (-0.04-0.14), resulting in flat trend estimated with relatively high precision (SD = 0.04). We mapped the variation in predicted occurrence probabilities and corresponding prediction uncertainty along the productivity gradient. Our results provide a much needed baseline against which future anticipated declines in M. lucifugus occurrence can be measured. The dynamic distribution modeling approach has broad applicability to regional bat monitoring efforts now underway in several countries and we suggest ways to improve and expand our grid-based monitoring program to gain robust insights into bat population status and trend across large portions of North America.

Estimating temporal trend in the presence of spatial complexity: a Bayesian hierarchical model for a wetland plant population undergoing restoration.

Monitoring programs that evaluate restoration and inform adaptive management are important for addressing environmental degradation. These efforts may be well served by spatially explicit hierarchical approaches to modeling because of unavoidable spatial structure inherited from past land use patterns and other factors. We developed bayesian hierarchical models to estimate trends from annual density counts observed in a spatially structured wetland forb (Camassia quamash [camas]) population following the cessation of grazing and mowing on the study area, and in a separate reference population of camas. The restoration site was bisected by roads and drainage ditches, resulting in distinct subpopulations ("zones") with different land use histories. We modeled this spatial structure by fitting zone-specific intercepts and slopes. We allowed spatial covariance parameters in the model to vary by zone, as in stratified kriging, accommodating anisotropy and improving computation and biological interpretation. Trend estimates provided evidence of a positive effect of passive restoration, and the strength of evidence was influenced by the amount of spatial structure in the model. Allowing trends to vary among zones and accounting for topographic heterogeneity increased precision of trend estimates. Accounting for spatial autocorrelation shifted parameter coefficients in ways that varied among zones depending on strength of statistical shrinkage, autocorrelation and topographic heterogeneity--a phenomenon not widely described. Spatially explicit estimates of trend from hierarchical models will generally be more useful to land managers than pooled regional estimates and provide more realistic assessments of uncertainty. The ability to grapple with historical contingency is an appealing benefit of this approach.