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 et al., 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 <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 et al., 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 et al., 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 et al., 2011). We grouped bat passes into species groups for Lasiurus borealis/L. seminolus, Eptesicus fuscus/Lasionycteris noctivagans, and Myotis austroriparius/M. septentrionalis due to overlap in acoustic call characteristics between these species (Grider et al., 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 et al. (2015). To do so, we estimated average percent coverage of a 1 m² 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² ha⁻¹) 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 et al., 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 et al., 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 et al. (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 et al., 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 et al., 2010). The ecological process model is z_i,j, the true state of presence or absence of species i at sites j. Similar to Tikhonov et al. (2020), this model uses a spatial factor model along with Nearest Neighbor Gaussian Processes (NNGP; Datta et al., 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 (r ≥ |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: hypermeans to 0 and hypervariances to 2.72 (Banerjee, Carlin, & Gelfand, 2014) in Normal priors, and scale and shape parameters to 0.1 (Lunn et al., 2013) in inverse-Gamma priors. To control spatial autocorrelation, the spatial decay phi for each latent factor followed a uniform Unif (0, 10) distribution. We determined model convergence of Markov chains using R-hat statistic values (<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 < 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 et al., 2020). We computed Bayesian P-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 <0.1 or >0.9 suggested poor model fit (Gelman, Meng, & Stern, 1996; Hobbs & Hooten, 2015).