DOI: 10.5061/dryad.m0cfxpp9f
title: Data from: Bat winter foraging habitat use in working forests: A multispecies spatial occupancy approach

description: Insectivorous bats in temperate zones have evolved strategies such as
 migration or hibernation to overcome challenges of reduced resource
 availability and increased energy demand during winter. In the
 southeastern United States Coastal Plain, bats are either year-round
 residents and remain active during winter or are migrants from colder
 areas seeking milder temperatures. Southeastern Coastal Plain forests also
 may represent important areas for remnant populations of species impacted
 by white-nose syndrome. Working pine (Pinus spp.) forests comprise a large
 proportion of southeastern Coastal Plain forests, yet winter bat habitat
 associations and how forest management affects bat use remain
 understudied. Hence, we used hierarchical multispecies spatial occupancy
 models to evaluate factors influencing winter bat occupancy and foraging
 habitat associations in working forests of the southeastern Coastal Plain.
 From January to March 2020–2022, we deployed Anabat Swift acoustic
 detectors and measured site- and landscape-level covariates on six working
 landscapes. We detected five species of bats and three species groups at
 93% (224/240) of sites. We observed higher species richness at sites with
 high proportions of contiguous forest and low levels of basal area. At the
 species level, occupancy patterns were influenced by site and landscape
 covariates, which had varying effects on species with distinct foraging
 strategies. Temperature was an important predictor of detectability. Our
 findings offer new insights into the ecology of bats in working forest
 landscapes during winter, where we highlight positive responses in
 occupancy with contiguous forests and lower levels of basal area, as in
 previous summer work. By providing valuable information on winter
 community composition and foraging habitat associations, we hope to guide
 management decisions for forest attributes important to these species,
 thus increasing conservation opportunities within working forests.
description: We conducted our study on six working forest
 landscapes across six states (Florida, Georgia, Louisiana, Mississippi,
 North Carolina, and South Carolina) of the southeastern United States
 Coastal Plain during 2020–2022 (Fig. 1). All study areas consisted primarily of planted loblolly pine (<em>P. taeda</em>) stands interspersed with streamside management zones (predominantly mature hardwood trees), roads, and wildlife openings, with other non-forest areas accounting for the remaining land area. We selected study areas &gt;3 000 ha and comprised primarily of upland planted pine with &lt;15% in forested wetlands. Management activities were typical of commercial forestry operations in the region, including clear-cutting at 20–35 years, mechanical and/or chemical site preparation, and planting 182–283 pine trees ha<sup>−1</sup> (Gresham, 2002). Competing vegetation was temporarily suppressed through herbicide applications, prescribed fire, or mechanically, with most stands being thinned at least once. On each study area, we created a 900 × 900 m grid and used ArcGIS Pro 2.8.0 (ESRI, Redlands, California, USA) to randomly select grid intersections as sampling points. The grid spacing was selected to ensure that the distance between sampling points encompassed a core area that constituted much of an individual bat's foraging movements (Morris, Miller, & Conner, 2011; Bender <em>et al</em>., 2015). We surveyed 40 sampling points randomly selected from the grid on each study area to ensure enough samples to adequately represent variation in stand age, stand size, and management history. We sampled all points at each study area within a 1-month period. We defined January–March as the winter sampling season as mean nightly temperatures are lowest (typically &lt;10°C) during this time throughout most of the Coastal Plain region (NOAA Climate.gov. https://www.climate.gov/maps-data/data-snapshots/averagetemp-monthly-1981-2010-cmb-0000-02-00?theme=Temperature). At each sampling point, we deployed Anabat Swift acoustic detectors with omnidirectional ultrasonic microphones US-OV2 and US-OV3 (Titley Electronics, Ballina, New South Wales, Australia; Appendix S1: Table S1) for three consecutive nights, recording from 30 min before sunset to 30 min after sunrise (Reichert <em>et al</em>., 2018). If rain occurred during the sampling period, we left detectors out for additional nights to ensure three nights of rain-free sampling. We placed detectors on poles with microphones 3 m above the forest floor pointed in the direction of the least vegetation clutter (Weller & Zabel, 2002). We coupled each detector with a temperature logger (HOBO Pendant G Acceleration Data Logger, Onset Computer Corp., Pocasset, Massachusetts, USA) programmed to record hourly temperature. We used auto ID software and subsequent visual vetting to identify calls to species, as recommended by the North American Bat Monitoring Program (NABat; Reichert <em>et al</em>., 2018). We first used Kaleidoscope Pro 5.4.1 software (Wildlife Acoustics Inc., Maynard, Massachusetts, USA) to filter noise files. We selected default filter setting parameters for bat analysis specifying a signal of interest between 8 and 120 kHz, 2 to 500 ms, and at least 2 pulses per sequence. We used the batch function in Kaleidoscope Pro to split each sequence to a maximum duration of 10 s for standardization. We selected the auto classifier of Kaleidoscope Pro with a balanced sensitivity level for classification to assist the visual vetting. Subsequently, we manually analyzed all remaining files using call structure, frequency of minimum and maximum energy, characteristic frequency, duration, inter-pulse interval, and slope (O'Farrell & Gannon, 1999; Szewczak <em>et al</em>., 2011). We grouped bat passes into species groups for <em>Lasiurus borealis</em>/<em>L. seminolus</em>, <em>Eptesicus fuscus</em>/<em>Lasionycteris noctivagans</em>, and <em>Myotis austroriparius</em>/<em>M. septentrionalis</em> due to overlap in acoustic call characteristics between these species (Grider <em>et al</em>., 2016; Johnson & Chambers, 2017; Kunberger & Long, 2022). We measured three components of vegetation structure at each sampling point (Appendix S1: Table S2). First, we used a convex spherical densiometer (Forestry Suppliers Inc., Jackson, Mississippi, USA) to measure percent canopy openness, which can be managed via planting density, by averaging measurements taken at the acoustic point and four additional locations in each cardinal direction 5 m from the point. Second, we characterized vegetation clutter, which can relate to forest management through mechanical, chemical, and prescribed burning practices, using methods based on Nudds (1977) and modified by Bender <em>et al</em>. (2015). To do so, we estimated average percent coverage of a 1 m<sup>2</sup> panel raised 4.5 m above the ground and 5 m from the acoustic point in each cardinal direction and in the direction the microphone was oriented. Third, we used a 10-factor prism (Husch, Beers, & Kershaw, 2003) centered at the acoustic detector point to estimate basal area (m<sup>2</sup> ha<sup>−1</sup>) of overstory trees, which again can relate to planting density, thinning, and other forest management activities. We used ArcGIS Pro and Fragstats v4.2 (McGarigal <em>et al</em>., 2015) to calculate landscape metrics from landowner-provided and publicly available data (Appendix S1: Table S2). Although variables at this scale cannot be managed directly, they may be important for managers to consider when implementing landscape-scale planning. We measured proportions of forest and wetland cover types and determined total edge (m) as landscape composition metrics within a 450-m-radius circular buffer around sampling points. The 450-m buffer area represented the area that did not overlap with the buffers of neighboring sampling points. We defined edge as the boundary between any two of six cover types reclassified from the National Land Cover Database (Dewitz & U.S. Geological Survey, 2021). We grouped forest stands into growth stages (hereafter, stand age; 0–3 [early establishment], 4–7 [closing canopy], 8–13 [closed canopy, pre-thinned], 14–20 years [mid-rotation thinned], or 21+ years old [mature forest, semi-closed canopy; including streamside management zones/bottomland hardwood forests]) as it can relate to forest management activities (e.g., thinning, final harvest) and is easily interpreted by forest managers (Marshall <em>et al</em>., 2022). Lastly, we measured distance (m) from sampling points to roads and permanent water using the Near tool in ArcGIS Pro. We implemented the hierarchical multispecies spatial occupancy model developed by Doser <em>et al</em>. (2022). The hierarchical model, which consists of an ecological process model and an observation sub-model, accounts for residual species correlation in a joint species distribution model framework while considering imperfect detection. The model quantifies the probability of occupancy for each species by accounting for factors influencing detection (MacKenzie <em>et al</em>., 2018). This hierarchical approach, in which species-specific effects are treated as random effects arising from a common community-level distribution, allows for inference of management effects on individual species and overall communities (Zipkin <em>et al</em>., 2010). The ecological process model is <em>z</em><sub><em>i</em>,<em>j</em></sub>, the true state of presence or absence of species <em>i</em> at sites <em>j</em>. Similar to Tikhonov <em>et al</em>. (2020), this model uses a spatial factor model along with Nearest Neighbor Gaussian Processes (NNGP; Datta <em>et al</em>., 2016) to ensure computational efficiency of species assemblages at different spatial locations. The observational sub-model (detection sub-model hereafter) separately models imperfect detection from the latent ecological process (see Doser, Finley, & Banerjee, 2022 for the modeling framework). Occupancy covariates included a combination of site- (basal area, canopy openness, and vegetation clutter) and landscape-level (total forest, total wetland, total edge, distance to freshwater, distance to roads, and stand age). We expected the influence of covariates on bat species to differ depending on their foraging strategy (Appendix S1: Table S2). Detection covariates included basal area, temperature at sunset, vegetation clutter, and year. We standardized all continuous covariates for both ecological and survey processes to a mean of 0 and a standard deviation equal to 1 (Zipkin, DeWan, & Royle, 2009; Kéry & Royle, 2015). We tested for correlation among continuous predictor variables using Pearson's correlation coefficient to ensure that highly correlated (<em>r</em> ≥ |0.7|) variables were not included in the same model. We fit our models using Polya-Gamma data augmentation (Polson, Scott, & Windle, 2013) for computational efficiency in R version 4.4.1 (R Core Team, 2020) via package spOccupancy (function sfMsPGOcc; Doser, Finley, & Banerjee, 2022). Accommodating sources of spatial dependence among observations is key to obtaining valid inferences about species occupancy (Doser, Finley, & Banerjee, 2022), thus we fit a spatial factor model to control for spatial correlations and residual spatial variation in species occurrence. We implemented spatial models using three replicate Markov chain Monte Carlo (MCMC) iterations to generate 10 000 samples from the posterior distribution of each model after discarding a “burn-in” of 5 000 samples, with a thinning rate of 50. We selected an exponential covariance to model spatial dependence structure among observations (Banerjee, Carlin, & Gelfand, 2014). We estimated model parameters and community summaries, setting default vague prior hyperparameter values: <em>hypermeans</em> to 0 and <em>hypervariances</em> to 2.72 (Banerjee, Carlin, & Gelfand, 2014) in Normal priors, and scale and shape parameters to 0.1 (Lunn <em>et al</em>., 2013) in inverse-Gamma priors. To control spatial autocorrelation, the spatial decay <em>phi</em> for each latent factor followed a uniform <em>Unif</em> (0, 10) distribution. We determined model convergence of Markov chains using R-hat statistic values (&lt;1.1) for all parameters within the models (Brooks & Gelman, 1998). We used the Widely Applicable Information Criterion (WAIC; Watanabe, 2010) to compare our set of models and shortlist the best-performing models, with models with a ΔWAIC &lt; 2 being biologically plausible and relevant. To evaluate detection covariates, we constructed models of single and all possible additive combinations of variables and compared them by including an occupancy sub-model with only the spatial structure, and no covariates. Temperature at sunset was the top-ranked detection model (Appendix S1: Table S3) and was subsequently included as the only covariate in the detection sub-model. We then developed 25 spatial models that included single and additive combinations of covariates, along with null and global models, in the occupancy sub-models and temperature at sunset in the detection sub-model (Appendix S1: Table S4). We calculated posterior mean and standard deviation of the model coefficients with 95% Bayesian credible intervals (BCI). Parameter estimates of covariates with BCI that did not cross 0 were considered important predictors of species occupancy, as this was reflective of a consistent relationship within model iterations. However, we also considered covariates as biologically meaningful if estimated 75% BCIs did not overlap zero, although the 95% BCIs overlapped zero (Cumming & Finch, 2005; Nakagawa & Cuthill, 2007; Tilker <em>et al</em>., 2020). We computed Bayesian <em>P</em>-values with Freeman-Tukey statistic to assess model fit, where a model with a good fit to the data had a value near 0.5, while values &lt;0.1 or &gt;0.9 suggested poor model fit (Gelman, Meng, & Stern, 1996; Hobbs & Hooten, 2015).
description: # Data from: Bat winter foraging habitat use in working forests: A
 multispecies spatial occupancy approach
 [https://doi.org/10.5061/dryad.m0cfxpp9f](https://doi.org/10.5061/dryad.m0cfxpp9f) ## Description of the data and file structure ## Overview and Objectives Working forest owners and managers are increasingly committed to conserving biodiversity, as evidenced by voluntary enrollment in sustainable forestry certification programs that include biodiversity principles. Given the geographic scale and economic and social importance of privately owned working forests, understanding how biodiversity can be conserved in managed landscapes is imperative. However, limited data on foraging ecology and selection of foraging areas by bats in working forest landscapes, especially outside the growing season, hinders our ability to evaluate management decisions. Thus, investigations into winter habitat use are needed to ensure that management actions provide suitable habitat conditions year-round. Hence, we used a multispecies spatial occupancy modeling approach that explicitly accounts for imperfect detection, spatial autocorrelation, and species correlations to examine winter bat associations on working forest lands across the southeastern United States Coastal Plain. To do so, we examined the influence of site- and landscape-level habitat characteristics on species richness to study the effect on winter bat community composition and determined species-specific winter foraging habitat occupancy at site- and landscape-levels. ## Data The dataset includes: * **Code for Models**: Scripts used for the Spatial Factor Multi-Species Occupancy Models. More information about sp0occupancy package and the sfMsPGOc function can be found [https://doserlab.com/files/spoccupancy-web/reference/sfmspgocc](https://doserlab.com/files/spoccupancy-web/reference/sfmspgocc) * **Code for Figures**: Scripts used to generate figures included in this study. * **Data**: Excel spreadsheet including temperature and vegetation structure data for each site.  Please note that data collected on or derived from privately owned properties are sensitive and cannot be publicly shared. For more detailed information, feel free to contact me. ## Contact For inquiries about the data or experimental methods, please reach out to: * **Dr. Santiago Perea**: [santip1320@gmail.com](mailto:santip1320@gmail.com) * **Prof. Steven B. Castleberry**: [scastle@uga.edu](mailto:scastle@uga.edu) ### Files and variables #### File: sfMSOM\_winterbats\_figures.R **Description:** Scripts used to generate figures for the study. #### File: sfMSOM\_winterbats.R **Description:** Script used for the Spatial Factor Multi-Species Occupancy Models #### File: sfMSOM\_data\_available.excel **Description:** Available data; the rest of data collected on or derived from data collected on privately owned properties are sensitive and cannot be provided publicly. The dataset includes presence/absence data for each bat species detected (EPFU: *Eptesicus fuscus*, LABO/LASE: *Lasiurus borealis* / *Lasiurus seminolus*, LACI: *Lasiurus cinereus*, MYO: *Myotis austroriparius* / *Myotis septentrionalis*, NYHU: *Nycticeius humeralis*, PESU: *Perimyotis subflavus*, TABR: *Tadarida brasiliensis*) across three survey nights per site, three years and spanning working forest lands of six states: Florida, Georgia, Louisiana, Mississippi, North Carolina, and South Carolina. To account for factors influencing detectability (*p*), the dataset includes temperature at sunset (°C) for each survey night as well as the mean sunset temperature (°C). Temperature data were recorded using HOBO Pendant G Acceleration Data Loggers (Onset Computer Corp., Pocasset, Massachusetts, USA), programmed to record hourly temperature at each sampling site. For occupancy (*ψ*), the dataset includes three key components of vegetation structure at each sampling site: * **Canopy Openness** – Estimated using a convex spherical densiometer (Forestry Suppliers Inc., Jackson, Mississippi, USA) by averaging measurements taken at the acoustic point and four additional locations (5 m away in each cardinal direction). * **Vegetation Clutter** – Assessed using a modified version of Nudds (1977) as described by Bender et al. (2015). Clutter was estimated by averaging percent coverage of a 1 m² panel positioned 4.5 m above the ground and 5 m from the acoustic point in each cardinal direction, as well as in the direction the microphone was oriented. * \**Basal Area **– Measured using a 10-factor prism (Husch, Beers, & Kershaw, 2003) centered at the acoustic detector location, providing an estimate of basal area (m² ha⁻¹). **References** Bender, M.J., Castleberry, S.B., Miller, D.A. & Wigley, T.B. (2015). Site occupancy of foraging bats on landscapes of managed pine forest. *For. Ecol. Manage.* 336, 1–10. Husch, B., Beers, T.W. & Kershaw, J.A. (2003). *Forest mensuration*, 4th edn. New York, NY: John Wiley. Nudds, T.D. (1977). Quantifying the vegetation structure of wildlife cover. *Wildl. Soc. Bull.* 5, 113–117. ## Code/software We implemented the hierarchical multispecies spatial occupancy model developed by Doser *et al*. (2022). We fit our models using Polya-Gamma data augmentation for computational efficiency in R version 4.4.1 (R Core Team, 2020) via package spOccupancy (function sfMsPGOcc; Doser, Finley, & Banerjee, 2022). All figures were created using ggplot2 in R.  **References** Doser, J.W., Finley, A.O., Kéry, M. & Zipkin, E.F. (2022). spOccupancy: an R package for single-species, multi-species, and integrated spatial occupancy models. *Methods Ecol. Evol.* 13, 1670–1678. R Core Team. (2020). *R: a language and environment for statistical computing*. Vienna: R Foundation for Statistical Computing. Doser, J.W., Finley, A.O. & Banerjee, S. (2022). Joint species distribution models with imperfect detection for high-dimensional spatial data. arXiv, 2204.02707