Abstract
The federally threatened southeastern beach mouse Peromyscus polionotus niveiventris occupies just a fraction of its former range along Florida’s Atlantic coast, mostly within the 72 km of continuous coastline on federal lands of the Cape Canaveral Barrier Island Complex. Although believed to be caused by loss and degradation of habitat associated with human urban expansion, no precise explanation for the reduction in southeastern beach mouse range is known. We used nine years of winter habitat occupancy survey data for a trend analysis, and a subset of four years to identify key factors influencing occupancy dynamics in coastal habitat throughout the core population. Using a maximum likelihood based model selection approach, we tested evidence for the effects of environmental and anthropogenic factors arranged into four hypotheses related to habitat structure, beach erosion, population isolation, and access to inland habitat. As predicted, sites closer to the primary dune, farther from open water, and with higher elevation had higher initial occupancy than sites with the opposite characteristics. Also as predicted, sites farther from open water, with higher elevation, more coastal strand, less ruderal vegetation, and in proximity to oak trees had lower site extinction. Contrary to predictions sites closer to infrastructure (roads, railroads, paved sites) and more isolated by narrow width or unsuitable habitat features also had lower extinction, possibly because of restriction of predator access. A Bayesian analysis showed that annual southeastern beach mouse habitat occupancy was generally high (0.63–0.97) but variable across the 72 km of coastal dune spanning most of the remaining species’ range. Simulation showed that the power of the survey was sufficient to detect a 20% decline between years. This study provides a baseline of knowledge about habitat occupancy dynamics for future monitoring and management of the southeastern beach mouse and other beach mouse subspecies.
Introduction
Wildlife of coastal lands, many already dealing with pressure from human encroachment, face increased disturbance because of climate change. A good example of this are beach mice, a group of seven subspecies of the oldfield mouse Peromyscus polionotus restricted to coastal areas in Florida and Alabama (Hall 1981). The loss or degradation of habitat because of human activities, sometimes combined with hurricanes and tropical storms, has resulted in population declines in all beach mouse subspecies. The southeastern beach mouse Peromyscus polionotus niveiventris occurs along Florida’s Atlantic coast and is federally listed as threatened under the US Endangered Species Act (US Fish and Wildlife Service 1989). Once occupying 360 km along Florida’s Atlantic coast (Stout 1992), the range of the southeastern beach mouse has been reduced to a single core population along 72 km of continuous coastline on federal lands of the Cape Canaveral Barrier Island Complex (CCBIC, Figure 1) and two small, isolated populations located at Smyrna Dunes Park (16 km to the north) and Pelican Island National Wildlife Refuge (72 km to the south). Much of the habitat occupied by the southeastern beach mouse is adjacent to eroding beaches, and accelerating sea level rise because of global climate change is threatening the remaining habitat (Foster et al. 2017).
Location of 114 sites for dynamic habitat occupancy survey of southeastern beach mouse Peromyscus polionotus niveiventris on the Cape Canaveral barrier island complex (CCBIC) during winters of the years 2010–2018. Red points are 35 sites on Cape Canaveral Space Force Station, yellow points are 16 sites on Kennedy Space Center, and green points are 63 sites on Canaveral National Seashore.
Location of 114 sites for dynamic habitat occupancy survey of southeastern beach mouse Peromyscus polionotus niveiventris on the Cape Canaveral barrier island complex (CCBIC) during winters of the years 2010–2018. Red points are 35 sites on Cape Canaveral Space Force Station, yellow points are 16 sites on Kennedy Space Center, and green points are 63 sites on Canaveral National Seashore.
The southeastern beach mouse had low levels of genetic diversity prior to recent human disturbance, and the isolated small populations at Smyrna Dunes Park and Pelican Island National Wildlife Refuge have lower genetic diversity than the large core population (Kalkvik et al. 2012). The core population on CCBIC appears to be well-mixed, which suggests that habitat fragmentation is not currently a problem (Zimmerman et al. 2015), but this could change rapidly as barriers to movement develop in the future. Most habitat within the former range is highly fragmented and in proximity to human disturbance and associated domesticated predators (i.e., cats, dogs, rodents). However, there is potential to reintroduce the southeastern beach mouse to some managed areas from which it was extirpated. Information on southeastern beach mouse habitat requirements will be useful to prioritize sites by suitability and identify management that could enhance success.
In 2010, we developed a monitoring study to learn what factors influence the dynamics of southeastern beach mouse habitat occupancy. To begin we held a workshop with biologists that had previously worked on the subspecies and together developed a set of 25 proposed mechanisms of how environmental factors influence southeastern beach mouse habitat occupancy (Table 1). To help organize our thinking, we grouped the proposed mechanisms into four sets of similar effects associated with background theory on the ecological factors that explain patterns in habitat occupancy of southeastern beach mice:
List of covariates thought to potentially explain patterns in habitat occupancy for the southeastern beach mouse Peromyscus polionotus niveiventris on the Cape Canaveral barrier island complex (CCBIC) during winters of the years 2010–2018. Covariates are grouped into four broader hypotheses (column Group: Refugia = access to inland refugia, CD = coastal dynamics, Isolation, and Habitat = habitat structure) for which they provide predictions of direction of effect (+ or −) on parameters in the dynamic habitat occupancy model (Detection, Site Initial Occupancy, Site Extinction, and Site Colonization). Covariates with Group = Detect refers to covariates thought to effect detection only. Bolded A or D identify effects that had model selection support; A agreed, and D disagreed with a priori predictions.

Access to inland refugia: because of woody structure and position farther inland, coastal oak scrub is more stable than dune and thus serves as refugia for southeastern beach mice during disturbance (Extine and Stout 1987; Swilling et al. 1998; Pries et al. 2009; Suazo et al. 2009). This hypothesis predicts that sites with or connected to oak scrub will have higher site colonization and lower site extinction rates, and thus also higher initial occupancy, whereas sites with more nearby hardwood hammock should have opposite predictions.
Isolation: if all the southeastern beach mice occupying a particular patch die or emigrate, there must then be a recolonization event before the site is re-occupied (Holler et al. 1989; Pries et al. 2009). This hypothesis predicts that more isolated (less connected) sites will have lower initial occupancy, higher site extinction, and lower site colonization than more connected sites. Although we did not initially consider it, a separate effect might be that predators have more access to sites with connections to adjacent inland habitat, leading to the opposite predictions (Wilkinson et al. 2013).
Habitat structure: southeastern beach mice prefer sites with coastal strand and oak scrub near primary dunes (Extine and Stout 1987; Bolt et al. 2019), especially higher elevation dunes (Cove et al. 2024). At a smaller scale better habitat is believed to be more open, less woody, and have more food plants (Stout 1992). Predictions are that sites with higher quality habitat will have higher initial occupancy, lower site extinction, and higher site colonization rates, whereas sites with lower quality habitat will produce the opposite effects.
Coastal dynamics: southeastern beach mice do better in sites located in coastal regions with accreting beaches than those in regions of erosion. Sites located on sections of accreting beach tend to have more open vegetation structure preferred by beach mice (Mullen et al. 2009), whereas sites on eroding beaches have denser vegetation structure. Newer beaches may also be better suited for burrowing than older sites with more mature soils. This hypothesis predicts that sites on accreting beaches will have higher initial occupancy and lower site extinction rates than sites on eroding beaches.
These hypotheses were not mutually exclusive, nor were we certain that they included all relevant hypotheses, but there was enough evidence from previous studies to warrant their consideration. In this paper we used dynamic habitat occupancy models to measure and rank the effects of environmental factors on southeastern beach mice habitat occupancy and compared these effects to predictions. A second objective was to report trends in the annual dynamics of habitat occupancy of the southeastern beach mouse in coastal locations within the range of the core population on CCBIC during 2010–2018.
Study Site
The study site consisted of 72 km of the CCBIC on the Atlantic coast of Florida (Figure 1). This region spans the north-south extent of the core population of the southeastern beach mouse. The CCBIC was formed as prograding (building seaward) barrier islands consisting of multiple dune ridges (Schmalzer et al. 1999). Land elevations range from sea level to approximately nine meters above sea level. The climate is humid subtropical with mean annual high temperatures in July (28.0°C) and low in January (15.5°C). Annual precipitation averages approximately 130 cm with 60% of precipitation occurring between June and September. Since 1996, when instrumentation records began, regional sea levels have increased by more than 11 cm and the rise has accelerated since 2009–2010 (Foster et al. 2017); this follows the pattern of other areas of the U.S. East Coast (Yi et al. 2015). The CCBIC landscape has been greatly impacted by habitat fragmentation that not only separates habitat patches but also reduces the spread of fire across the landscape (Duncan et al. 2004); both impacts reduce habitat quality for many species (Breininger et al. 1998).
Coastal habitat on CCBIC consisted of four types generally arranged in linear zones landward from the Atlantic Ocean. Dune included the dune ridges and swales within the dynamic dune topography nearest to the ocean and was open vegetation structure dominated by sea oats Uniola paniculata and other grasses including beach grass Panicum amarum. Small shrubs were common including varnish leaf Dodonaea viscosa, seacoast marsh elder Iva imbricata, and beach tea Croton punctatus. Common herbs included beach sunflower Helianthus debilis, railroad vine Ipomoea pes caprae, shoreline purslane Sesuvium portulacastrum, ground cherry Physalis walteri, and camphorweed Heterotheca subaxillaris (Schmalzer and Foster 2016). Coastal grassland was the next landward habitat in only some parts of the CCBIC, mainly south of the cape. Located between the primary dunes and predominately shrubby vegetation inland, the primary ground cover was composed of several grasses, including Muhly grass Muhlenbergia capillaris, and bluestems Andropogon spp. or Schizachyrium spp. (Schmalzer and Foster 2016).
The coastal strand, a shrub community on older dunes inland from coastal dunes, was typically the next habitat landward. Dominant vegetation included saw palmetto Serenoa repens, sea grape Cocoloba uvifera, wax myrtle Myrica cerifera, Simpson’s stopper Myrcianthes fragrans, tough buckthorn Sideroxylon tenax, and rapanea Rapanea punctata (Schmalzer and Foster 2016). The cover of this area was dense shrub with some sandy openings with sparse herbs such as beach sunflower and camphorweed. Pruning of shrubs by salt spray was common. Coastal scrub, the most inland habitat, was a shrub community where a coastal form of live oak Quercus virginiana, was dominant (Johnson and Barbour 1990). Other shrubs included saw palmetto, wax myrtle, tough buckthorn, rapanea, and snowberry Chicocoa alba (Schmalzer and Foster 2016).
Methods
Occupancy survey design
We sampled habitat occupancy February–March, 2010–2018 at 71-114 15-m radius sites, using eight tracking tubes at each site, for up to four three-night or two six-night survey periods each year (one track paper was left in the tube over the survey period). The number of temporal replicates differed with only one sample in 2010 and 2011, four samples in 2012 and 2013, and two samples for the remaining four years (Table 2). In 2012 and 2013 we sampled 71 of the 114 sites, excluding 43 sites in the central portion of Canaveral National Seashore. In 2016 and 2017 we placed tracking tubes within 3 m of the center. For analysis, we assumed that sites were closed to occupancy state changes during the sampling period within each year, and we treated the three-night surveys in 2012 and 2013 and the six-night surveys in the other years as equivalent. We aggregated detections for all tubes within sites into a single detection event. We located occupancy sampling sites within 200 m of the Atlantic Ocean in coastal dune or strand, with site locations systematically spaced approximately 600 m apart along the coastline from a randomly selected start location at the southern CCBIC boundary. At each location, we randomly selected a distance west of the beach to locate the sample site based on the width of coastal dune or strand at the site. We placed a PVC pole to mark the center of the sample site and arranged the eight tracking tube locations in cardinal compass bearings (e.g., 0°, 45°, 90°) around the center pole, each a random distance 1–15 m from the center pole. We recorded covariate data within a 1 m radius for each tube location including: 1) percent of vegetative cover; 2) coverage (absent, patchy, or continuous) of saw palmetto within 15 m; 3) presence or absence of five food plants, beach sunflower, sea oats, seaside bean Canavalia rosea, ground cherry, and varnish leaf; 4) height of the vegetation (m); and 5) composition of vegetation (herbaceous, woody, or mixed). We used these covariates at the tube level for modeling detection and aggregated them to the site level by calculating means or rate of occurrence among the eight tube locations at each site (Table 1).
Differences in sampling design between years in the number of replicate track tube surveys (temporal replicates), the number of sites, the scale of track tube deployment, and the number of track tubes deployed at each site for the habitat occupancy survey of the southeastern beach mouse Peromyscus polionotus niveiventris on Cape Canaveral barrier island complex (CCBIC) during winters of the years 2010–2018.

Tracking tube design and assembly followed Loggins et al. (2010). However, in 2015 we modified the tube design, and inserted a PVC reducer into the 90° elbow joint to decrease the opening from 2 inches (5.08 cm) to 1 inch (2.54 cm) to reduce access to the tubes and tracking papers by larger non-target species including raccoons Procyon lotor, eastern spotted skunks Spilogale putorius, and cotton rats Sigmodon hispidus (Oddy et al. 2018). We deployed tracking tubes and measured footprints as described in Stolen et al. (2014). We followed American Society of Mammalogists guidelines (Sikes and Gannon 2011), had our procedures approved by a NASA Institutional Animal Care and Use Committee, and conducted our work under US Fish and Wildlife Service permit number TE089075-3.
Habitat analysis
We measured the extent of habitat surrounding the sample sites from a land cover map (SJRWMD Land Use and Land Cover GIS Layer 2004) with refined categorization, an increased extent, and delineation of new landscape features identified with high-resolution imagery acquired during December 2003, and site-specific ground knowledge. We calculated the composition of features within sites using an unsupervised classification of a 15 m radius area around each sample site where the source data were one-foot resolution (0.305 m) color infrared aerial imagery taken in December 2008. For each site we analyzed the proportion of unvegetated land within the 15 m radius (Table 1), and the distances to: 1) nearest hard infrastructure (e.g., train track, paved road), 2) nearest water barrier (e.g., ditch, estuary, large pond), and 3) nearest oak scrub. We also recorded the amount of each of six land cover classes within a 200 m radius: 1) oak scrub, 2) hammock (hardwood and cabbage palm), 3) coastal strand, 4) marsh (fresh and salt), 5) ruderal (mowed grass), and 6) infrastructure (paved regions). We also included as a covariate the combined area of non-wetland vegetation. Based on a digital terrain model derived from LIDAR data collected during 2006 for Volusia County (O’Neill 2007) and 2007 for Brevard County (Maune 2009) we recorded the elevation of the center of the sample site and the distance and direction (landward or seaward) to the primary dune. To develop features for coastline change, we mapped the location of the interface between the primary dune and open beach and lines representing the apparent high-water mark, using on-screen digitizing over rectified aerial photography at a display scale of 1:4000, using the best available imagery. Dune interface change compared 2005 and 2010. Shoreline change compared 1994 to each of 1999, 2000, 2003, and 2007, and then used the average change for these years, to smooth out random fluctuations. Three variables were developed for each feature (Table 1): 1) the difference in m between years, 2) a three-level factor (bottom 0.1 quantile, middle 0.1-0.9 quantile, and upper 0.9 quantile), and 3) a three-level factor (bottom 0.25 quantile, middle 0.25-0.75 quantile, and upper 0.75 quantile). We obtained illumination due to the moon using ephemeris data (Solar System Dynamics. Horizons system 2017) as follows: we first determined the fraction of moon illuminated and above the horizon during each 15-minute interval between sunset and sunrise during each survey interval, and then summed these amounts to calculate a total moon illumination during the session. We obtained daily weather data from the NOAA global historical climatology network-daily database (Menne et al. 2012), and processed it to quantify the total precipitation, and the maximum and minimum temperature during each survey period. Finally, we developed several categorical landscape variables based on the 2003 land cover map, aerial imagery, and field visits each year: 1) habitat type (dune, coastal grassland, coastal strand, and coastal scrub), 2) region dividing the area into three or four parts, 3) width of habitat west of the point (less than or greater than 150 m), 4) connection to habitat west of point with levels: no barrier, barrier between adjacent habitat, or no adjacent habitat present, and 5) back habitat availability (whether or not the site was adjacent to suitable inland habitat). For the width and back habitat variables we included two separate versions to include or exclude hardwood hammock in the determination.
Occupancy modeling to identify ecological predictors of habitat occupancy
We used dynamic habitat occupancy models (Mackenzie et al. 2003) to evaluate evidence for the effects of the covariates relating to the broader hypotheses groups. We used annual habitat occupancy surveys for 2012–2015, the four years in which we measured the ecological covariates at sites. Occupancy models were fit using maximum likelihood methods implemented in the unmarked package (Fiske and Chandler 2011; Kellner et al. 2023) in R version 4.3.1 (R Core Team 2023). Differences in the level of sampling between years, because of changes in sampling design or accidents, was handled by the model estimation software by not including them in the calculation of the model likelihoods. We evaluated evidence for occupancy models using information-theoretic model selection based on the adjusted Akaike Information Criterion, AICc (Burnham and Anderson 2002). For dynamic occupancy models we began by testing all univariate models on the detection parameter with no covariates for the other parameters. We judged univariate detection models as supported if they were included in the set of models with combined Akaike weight of 0.95 (Anderson 2008). We dropped variables identified at this stage if their effect was estimated imprecisely, as indicated by a 95% confidence interval that included zero (Arnold 2010). We retained all the supported detection covariates for inclusion in the next round. We next examined covariate models for the occupancy parameter, because the first colonization and extinction events are conditioned on the initial occupancy state. We repeated the procedure used for selecting the covariate structure of the detection parameter for covariates on the occupancy parameter, including all models deemed important in the previous round for detection. This resulted in a set of models for occupancy which were then used in the next stage. We repeated the covariate selection procedure for the extinction parameter and finally the colonization parameter, reasoning that we had the least information to model colonization given the high level of initial occupancy in our study.
Following model selection, we evaluated support for our four hypothesis groups (Inland refugia, Coastal geomorphology, Isolation, and Habitat quality) based on the degree to which the predicted effects of individual covariates agreed. To assess the strength of covariate effects we calculated model-averaged parameter estimates for the dynamic occupancy models using a balanced set of additive models with the best supported covariates for each state parameter (initial occupancy, extinction, and colonization), limiting the models to those having a single covariate for each parameter (Burnham and Anderson 2002). We scaled the continuous covariates by subtracting the mean and dividing by two standard deviations; this puts the regression coefficients for numerical and categorical covariates on similar scales, allowing comparison of effect sizes (Gelman 2008). Code and data to conduct the analysis is provided as supplemental material (Text S1, Text S2, and Text S6).
We conducted a power analysis based on simulation methods by (Guillera-Arroita and Lahoz-Monfort 2012) designed to determine the probability of detecting a change in the proportion of sites occupied between two surveys. We simulated 10,000 data sets based on the monitoring survey design and parameter estimates in typical years described in this paper (114 sites, two temporal replicates, detection probability = 0.8, occupancy probability = 0.8). We set alpha = 0.05 and beta = 0.08 and the change in occupancy to 0.2.
Trends in annual habitat occupancy on CCBIC
Where pijk = the probability of a detection at site i in year j and survey k. Covariates effecting each parameter (Ψi, φi, γi, pijk) were modeled using a logit link and included all those supported in the analysis of 2012–2015 data. In addition to the selected sequence covariate on detection, we included a scale covariate to allow detection to differ during 2016 and 2017 to reflect changes in the number and scale of distances between tubes in those years.
To evaluate evidence for covariate effects we used an indicator variable model selection approach (Kuo and Mallick 1998). We included latent binary indicator variables for the beta parameters for each of the covariates. These indicator variables were specified as independent binary variables with Bernoulli (0.5) priors. Each binary indicator variable was multiplied by the beta parameter in the linear formula, which thereby determined if the covariate was included or excluded during a particular MCMC iteration. To ensure that the resulting beta parameters were not overly prone to omission because of the MCMC getting stuck within unlikely regions of the parameter space, we used a spike-and-slab prior approach (Mitchell and Beauchamp 1988). We first fit the full model with all covariates included to estimate the mean and standard deviation of each beta parameter when included. We then used normal prior distributions with the estimated means and standard deviations for the spikes, and diffuse Normal (0, 1000) priors for the slabs, with the indicator variables determining which prior was included in the likelihood. This had the effect of keeping the MCMC chains in reasonable regions of parameter space when particular variables were excluded from the model. To evaluate model support we tabulated the number of times each linear combination of variables were included in an MCMC iteration, with the resulting proportions providing relative model support (Royle et al. 2013). Because the parameter estimates for covariate effects are influenced by the priors when the variable is not included, we used parameter estimates from the full model for inferences about effect sizes. This issue does not affect the estimates of annual habitat occupancy (or any parameters not subject to model selection) so inference for those parameters were based on the indicator variable version of the model (Royle et al. 2013).
To account for unmodeled effects associated with site conditions we included a group level (random) effect of site in the linear models for rate of initial occupancy. To account for unmodeled effects associated with yearly conditions we included group level (random) effects of year in the linear models for the rates of extinction and colonization. Because the mechanisms for extinction and colonization may overlap to some degree, we allowed these random effects to be correlated (Gelman and Hill 2006). We also included an “individual level random effect” with a level for each replicate survey, to help account for non-independence between detections (Harrison 2014). Overdispersion is a concern in models of detections of individuals when there is a possibility of non-independence between individuals. All random effects were given a normal distribution with mean 0 and precision = 0.5 (SD = 1.4) as a prior distribution; such weakly informative priors can speed up convergence but do not influence inference if random effects are present (Gelman et al. 2008).
We implemented the hierarchical Bayesian models using JAGS 4.3.0 (Plummer 2003) and R (R Core Team 2023) with the package jagsUI (Kellner 2021). For each analysis, we ran an adaptation phase of up to 10,000 iterations decided automatically by JAGS, and then discarded samples as burn-in until Gelman-Rubin convergence diagnostic (R-hat) were less than 1.01 and the number of effective samples was estimated to be greater than 4,000 for all parameters, except for individual levels of the random effects. We did not thin the MCMC posteriors (Link and Eaton 2012). We tested model goodness of fit using a Bayesian p-value approach that compared the fit of the observed data to the equivalent fit of data simulated under the model with estimated parameters (Kéry and Royle 2020). Code and data to conduct the analysis is provided as supplemental material (Text S3, Text S4, and Text S5).
Results
Habitat analysis
The spatial distribution of the sample site habitat covariates was largely determined by landscape-scale patterns, as shown in Figure 2. Starting from the north, the first 16 points were within narrow (<150 m) coastal strand bordered by cabbage palm and oak hammock to the east. Moving south, the next 44 points were on a narrow barrier island (<0.5 km) with intermittent impounded salt marsh to the east. The next 11 points were also bounded by wetlands to the east, although the barrier island here was wider and there was potential connection to oak scrub inland along a few remnant dune ridges. This region had the highest rate of dune erosion. The next 29 sites occurred within a region of wide coastal strand with no adjacent wetlands. To the east in this region the barrier island was broad with extensive oak habitat. The southernmost 14 points occurred along an area of accreting coast with broad coastal strand with many open sandy areas. Although there was some spatial correlation within individual covariates, the covariates each showed different patterns between points, and no pairs were strongly correlated.
Spatial distribution of values for the best supported covariates in dynamic habitat occupancy models for southeastern beach mouse Peromyscus polionotus niveiventris on the Cape Canaveral barrier island complex (CCBIC) during winters of the years 2010–2018.
Spatial distribution of values for the best supported covariates in dynamic habitat occupancy models for southeastern beach mouse Peromyscus polionotus niveiventris on the Cape Canaveral barrier island complex (CCBIC) during winters of the years 2010–2018.
Occupancy modeling to identify ecological predictors of habitat occupancy
Over the four years the proportion of sites with detections (naïve occupancy rate) was 0.83–0.93 and all but three of the sites had detections in at least one year. The goodness of fit test of the full model indicated no lack of fit (Probability ([Chi-square statistic of bootstrap > Chi-square statistic of data] = 0.72). The model building procedure identified Sequence and Session as the only variables affecting detection (Table 3); detection increased between each sampling session with the largest increase between the first and second sessions (Figure 3A). Because both variables produced nearly identical estimates for detectability, in subsequent stages of modeling we retained only the linear effect of Sequence, which required fewer degrees of freedom. We investigated including a post hoc quadratic Sequence effect, but the quadratic term was not significant, so we did not retain this model. Ten covariates on initial site occupancy were identified as having support based on our model selection criteria (Table 3). In agreement with our a priori predictions (Table 1), initial site occupancy decreased with increasing distance from the dune (Figure 3B), increased with increasing distance from open water (Figure 3C), and increased as elevation of the center of the site increased (Figure 3D). Seven other variables were also selected but all had large variance estimates and were dropped from further consideration, including the aerial extent (hectares within 200 m of the center point) of infrastructure, the distance to oak scrub, extent of hammock, presence/absence of hammock, extent of marsh, distance to infrastructure, and presence of habitat westward of the site.
Predicted effects on probability of detection and initial site occupancy for the best supported covariates in dynamic habitat occupancy models for southeastern beach mouse Peromyscus polionotus niveiventris on the Cape Canaveral barrier island complex (CCBIC) during winters of the years 2010–2018. A) predicted detection as a function of survey sequence, B) predicted initial occupancy rate as a function of distance of the site from the primary dune, C) predicted initial occupancy rate as a function of distance of the site from open water > 0.2 ha in extent, D) predicted initial occupancy rate as a function the elevation of the center of the survey site. Graphs give the predicted mean (point and solid line) and 95% confidence interval (bars and dotted lines).
Predicted effects on probability of detection and initial site occupancy for the best supported covariates in dynamic habitat occupancy models for southeastern beach mouse Peromyscus polionotus niveiventris on the Cape Canaveral barrier island complex (CCBIC) during winters of the years 2010–2018. A) predicted detection as a function of survey sequence, B) predicted initial occupancy rate as a function of distance of the site from the primary dune, C) predicted initial occupancy rate as a function of distance of the site from open water > 0.2 ha in extent, D) predicted initial occupancy rate as a function the elevation of the center of the survey site. Graphs give the predicted mean (point and solid line) and 95% confidence interval (bars and dotted lines).
Model selection results for dynamic habitat occupancy models for all supported covariates for southeastern beach mouse Peromyscus polionotus niveiventris on Cape Canaveral barrier island complex (CCBIC) during winters of the years 2010–2018. Ψ() = covariate model on initial occupancy, γ() = covariate model on site colonization, ε() = covariate model on site extinction, p() = covariate model on detection, K = number of parameters estimated in model, ΔAICc = difference between adjusted Akiakie’s information criterion of model and the best supported model within group, Wi = model weight, ΣWi = sum of model weights for given model and all models with more support within group.

Eleven covariates on the site extinction parameter sub-model were identified as having support based on our model selection criteria (Table 3); two of these (area of hammock and year fixed effect) had low precision for parameter estimates and were thus not considered further. The covariates retained were distance from water, point elevation, distance from infrastructure, area of coastal strand, barrier class, habitat width, area of vegetated, area of disturbed ruderal, and presence of oak. In agreement with a priori predictions (Table 1) site extinction decreased with increasing distance from open water (Figure 4A) and decreased with increasing elevation (Figure 4B). Contrary to our prediction, site extinction increased with increasing distance from human infrastructure (Figure 4C). As predicted, sites with more coastal strand and more vegetation had lower rates of annual site extinction, and sites with more disturbed ruderal cover had higher extinction (Figures 4D, 4G, and 4H). Contrary to our predictions, sites with barriers had lower extinction than did either sites with no barrier to adjacent habitat or sites that lacked adjacent habitat (Figure 4E). Also contrary to our prediction, sites with wide habitat had higher extinction than did sites with narrow habitat (Figure 4F). As predicted, site extinction was lower for sites with oak within 200 m (Figure 4I). Only the effect of distance from human infrastructure on site colonization was identified as having support based on our model selection criteria (Table 3). This effect agreed with our a priori prediction but was imprecisely estimated.
Predicted effects on probability of site extinction for the best supported covariates in dynamic habitat occupancy models for southeastern beach mouse Peromyscus polionotus niveiventris on the Cape Canaveral barrier island complex (CCBIC) during winters of the years 2010–2018. A) distance of the site from open water > 0.2 ha in extent, B) elevation of the center of the survey site, C) distance of the site from paved surfaces, D) amount of coastal strand within 200 m of the site, E) isolation of the site from surrounding habitat, F) width of the barrier island at the site, G) amount of vegetated within 200 m of the site, H) amount of ruderal within 200 m of the site, I) presence or absence of oak adjacent to the site.
Predicted effects on probability of site extinction for the best supported covariates in dynamic habitat occupancy models for southeastern beach mouse Peromyscus polionotus niveiventris on the Cape Canaveral barrier island complex (CCBIC) during winters of the years 2010–2018. A) distance of the site from open water > 0.2 ha in extent, B) elevation of the center of the survey site, C) distance of the site from paved surfaces, D) amount of coastal strand within 200 m of the site, E) isolation of the site from surrounding habitat, F) width of the barrier island at the site, G) amount of vegetated within 200 m of the site, H) amount of ruderal within 200 m of the site, I) presence or absence of oak adjacent to the site.
Model-averaged parameter estimates for the selected covariates on each of the four parameters are shown in Figure 5. Model-averaged estimates for covariate effects on initial site occupancy agreed with predictions from individual covariate models for distance from dune and distance from water, but the model-averaged predicted effect of point elevation was weaker. The model average covariate effects on site extinction for the distance to infrastructure and area of vegetated were weaker and showed more uncertainty than did the effects in the individual covariate models. Other model average covariate effects on site extinction were like those from the individual covariate models. The model average effect of distance to infrastructure on site colonization was weak and showed much uncertainty. The power analysis predicted a power of 0.68 to detect a change in occupancy of 0.2 between subsequent years.
Model-averaged estimated covariate effects on detection, occupancy, extinction, and colonization for the 95% confidence set of dynamic habitat occupancy models for southeastern beach mouse Peromyscus polionotus niveiventris on the Cape Canaveral barrier island complex (CCBIC) during winters of the years 2010–2018.
Model-averaged estimated covariate effects on detection, occupancy, extinction, and colonization for the 95% confidence set of dynamic habitat occupancy models for southeastern beach mouse Peromyscus polionotus niveiventris on the Cape Canaveral barrier island complex (CCBIC) during winters of the years 2010–2018.
Trends in annual habitat occupancy on CCBIC
Over the 9 years of occupancy surveys the proportion of sites with detections (naïve occupancy rate) was 0.29–0.93 and all sites had detections in at least one year. The mean number of years of detections among sites was 5.4 (SD = 1.8). One site had detection in only one year, whereas four sites had detections in every year. Support for including covariates for the time series data agreed well with the 4-year analysis. The variables most often included by the indicator variable method were distance to dune and distance to open water on initial occupancy, and point elevation and the barrier state on site extinction (Table 4). The model with only these variables was chosen in 32% of the MCMC iterations, and these variables were included in 82% of the iterations. The Bayesian model did not show evidence of a lack of fit based on the goodness of fit check (Bayesian p-value value for the closed portion of the model was 0.51 and for the open portion 0.56).
Summary statistics for the posterior parameter estimates from the Bayesian dynamic habitat occupancy model using indicator variables to allow selective inclusion of covariates for southeastern beach mouse Peromyscus polionotus niveiventris on Cape Canaveral barrier island complex (CCBIC) during winters of the years 2010–2018. Three variables were included in the model for initial occupancy (psi), ten in the model for site persistence (phi) and one for site colonization (gamma). For each of these effects an indicator variable determined if the effect was included during each iteration of the Bayesian model. The detection model included the covariates sequence (linear effect of survey order) and scale (to adjust for changes in sampling scale) with no model selection. Correlated random effects for year were included in the persistence and colonization models, a site level random effect was included in the initial occupancy model, and a replicate level random effect was included in the detection model to account for potential overdispersion. Values give the estimated mean and standard deviation (SD) and the lower and upper quantiles for a 95% credible interval.

In general, the effect size and direction of covariate effects estimated by the full Bayesian model (all covariates included) for the time series data agreed with model average estimates for the focal years (Figure 6; compare with Figure 5 and Table 5). The estimates from the Bayesian model had better precision for the barrier state effect on site extinction and the effect of distance to infrastructure on site colonization. For several covariates on site extinction with high uncertainty in both data sets (distance to infrastructure, areas of vegetated and ruderal) the results differed in direction of effect (but all still overlapped zero). For the effect of habitat width on site extinction the effect was predicted to be much lower by the Bayesian model for the time series data. The annual estimated site occupancy ranged from 0.65 to 0.98 (Figure 7), and the occupancy rate appeared more similar between consecutive years than more distant years, suggesting fluctuations over several years.
Posterior credible intervals for effects of covariates on the probability of initial site occupancy (psi), site extinction (phi) and colonization (gamma) for southeastern beach mouse Peromyscus polionotus niveiventris on Cape Canaveral barrier island complex (CCBIC) during winters of the years 2010–2018. The central dot shows the location of the mean, the thicker bar the 80% credible intervals, and the thinner bar the 95% credible intervals.
Posterior credible intervals for effects of covariates on the probability of initial site occupancy (psi), site extinction (phi) and colonization (gamma) for southeastern beach mouse Peromyscus polionotus niveiventris on Cape Canaveral barrier island complex (CCBIC) during winters of the years 2010–2018. The central dot shows the location of the mean, the thicker bar the 80% credible intervals, and the thinner bar the 95% credible intervals.
Annual probability of site occupancy for southeastern beach mouse Peromyscus polionotus niveiventris on Cape Canaveral barrier island complex (CCBIC) during winters of the years 2010–2018, predicted with the Bayesian hierarchical model. The bars show the 95% credible interval of the predictions. The red dots show the observed proportion of sites with detections. For 2012 and 2013 when not all sites were sampled, occupancy was inferred from covariates and model parameters, resulting in an observed occupancy above the predicted lower credible limits in one of those years (2012).
Annual probability of site occupancy for southeastern beach mouse Peromyscus polionotus niveiventris on Cape Canaveral barrier island complex (CCBIC) during winters of the years 2010–2018, predicted with the Bayesian hierarchical model. The bars show the 95% credible interval of the predictions. The red dots show the observed proportion of sites with detections. For 2012 and 2013 when not all sites were sampled, occupancy was inferred from covariates and model parameters, resulting in an observed occupancy above the predicted lower credible limits in one of those years (2012).
Summary statistics for the posterior parameter estimates from the Bayesian dynamic habitat occupancy model including all covariates for southeastern beach mouse Peromyscus polionotus niveiventris on Cape Canaveral barrier island complex (CCBIC) during winters of the years 2010–2018. Three variables were included in the model for initial occupancy (psi), ten in the model for site persistence (phi) and one for site colonization (gamma). The detection model included the covariates sequence (linear effect of survey order) and scale (to adjust for changes in sampling scale) with no model selection. Correlated random effects for year were included in the persistence and colonization models, a site level random effect was included in the initial occupancy model, and a replicate level random effect was included in the detection model to account for potential overdispersion. Values give the estimated mean and standard deviation (SD) and the lower and upper quantiles for a 95% credible interval.

Discussion
We used both frequentist and Bayesian methods to model the effects of covariates on habitat occupancy of the southeastern beach mouse, using the strengths of each approach as needed (Kéry and Kellner 2024). Using maximum likelihood-based model selection and multi-model inference for the four years for which we had an expanded collection of data, we identified several covariates that explained patterns in habitat occupancy of the southeastern beach mouse and found support for three of the four ecological explanations considered. This principled statistical modeling approach based on a set of a priori hypotheses about effects helped us avoid being misled by spurious relationships (Anderson 2008). We then included additional years of monitoring data within a more flexible Bayesian analysis that accommodated changes in study design across nine years and including the best supported covariate effects. This analysis provided more details about how the key variables explained patterns in habitat occupancy dynamics.
The dynamic habitat occupancy model includes three parameters (initial occupancy, colonization, and extinction) which can be modeled with covariate sub-models. The initial occupancy parameter informs about what factors affect occupancy in the first year and the colonization and extinction parameters account for changes in site occupancy between sampling occasions (years for our study). These sub-models reveal how factors influence occupancy dynamics as environmental conditions change and populations respond to disturbance. By modeling dynamic occupancy as a function of habitat variables thought to be related to our a priori hypotheses, we quantified evidence of effects on occupancy parameters for variables related to hypotheses within our four groups (habitat structure, inland refugia, isolation, and coastal geomorphology).
Habitat structure
It is generally believed that beach mice are strongly associated with primary dune habitat (Extine and Stout 1987), but recent attention has been on the value of adjacent inland habitat, especially oak scrub. We included this hypothesis group to try to identify important features of dune habitat. When applied carefully, initial site occupancy can be used to infer habitat suitability (Yirka et al. 2018; Northrup et al. 2022). Three covariates explained patterns in the initial occupancy of the southeastern beach mouse and agreed with our a priori predictions about the effects of habitat. We found an increase in initial site occupancy with decreasing distance of sites from the primary dune, which makes sense given that dune is recognized as highly suitable for beach mice (Extine and Stout 1987). We also found that greater elevation of sites was associated with higher initial occupancy, which agrees with reports that large primary dune structures improve suitability (Stoddard et al. 2019). However, when applied to the nine-year time series data the effect of elevation was much weaker than was found during the variable selection phase with just four years of data. The association of higher elevation with increased site occupancy therefore warrants further attention. The decrease in site extinction rate as site elevation increased was as predicted and supports the value of higher dunes. We also found lower site extinction rates for sites with more coastal strand, more non-wetland vegetation, and less disturbed ruderal land cover as predicted, which can be interpreted as benefits of increasing availability of habitat.
Access to inland refugia
The only effect we found related to the hypothesis that adjacent inland habitat may allow beach mice to persist after habitat disturbance was a negative effect on extinction for sites with the presence of adjacent oak habitat, as predicted. A reason why there was not more evidence for this hypothesis could be that there were no direct hurricane impacts on the CCBIC during the study. Although there were a few hurricanes that passed nearby, and several winter storms, we did not investigate storm impact at a temporal or spatial scale fine enough to detect the impact on site occupancy.
Isolation
The finding that sites closer to open water had lower initial occupancy can be interpreted in two ways. Open water with surrounding wetland fringes may expose beach mice to novel predators such as aquatic vertebrates. Also, open water features on the CCBIC generally occur as physical barriers (long wide ditches, large borrow ponds) that could isolate sites. The decrease in site extinction rate as distance to water increased was as predicted and can be interpreted similarly to their effects on initial occupancy. The three covariates that had effects on site extinction that disagreed with our a priori predictions (width, barrier state, distance to infrastructure) were surprising. We predicted that the sites located in areas with more extensive habitat (wider) or with more connection to other habitat would have lower rates of extinction because of demographic rescue effects (Lehtinen 2023). However, we did not consider that increased connectivity might also increase predation pressure. Since we observed similar strength of effects for both data sets it is not easy to dismiss these findings as spurious. The implications for management could be significant, so understanding these effects warrants further attention, perhaps using camera traps to directly measure predator occurrence (Hooker et al. 2024). The effect of distance to infrastructure on site extinction rate is less clear because the prediction for the focal data had a predicted effect contrary to, whereas the time series data agreed with, our a priori prediction. The effect of distance to infrastructure on site colonization was as predicted for both data sets, although the estimates were very uncertain. Perhaps the best we can say is that the spatial arrangement of infrastructure that may serve as a barrier to beach mice deserves some attention in habitat management.
Coastal dynamics
There was no evidence for this hypothesis among the variables we tested. This may be because the scale of the investigation was too limited spatially and temporally to detect differences in the important factors. Unfortunately, since the CCBIC is the last remaining part of the SEBM range we may not be able to directly test this hypothesis. There was also a potential confounding between the geomorphology hypothesis and the spatial distribution of habitat. It is worth noting that the SEBM evolved on a high-energy dune system and may be adapted to changes in beach morphology over time, at least historically when the range was not restricted.
A key element of successful habitat occupancy studies is to correctly measure factors effecting detectability (Royle 2006). In our survey, detectability of beach mice by tracking tubes was high so it was not surprising that none of the numerous environmental, microhabitat, or anthropogenic variables we tested were useful in modeling detectability. Detectability increased between the sampling sessions (three or six nights of tube exposure) within the primary periods (Years). We believe this could be explained by the movements of southeastern beach mice within their home ranges. Since southeastern beach mouse home ranges are probably many times larger than the area of the track tube detector array (1.5 ha versus 0.07 ha), individual mice simply might need time to travel near enough to be detected. If true, this change in availability between surveys alters the definition of detection probability estimated in our model to a combined (confounded) availability and detection during the survey (Guillera-Arroita et al. 2017). The definition of site occupancy is also altered to one of site use rather than continuous occupancy. This seems appropriate for a mobile animal that moves around a home range many orders of magnitude larger than itself. Because our study design included multiple levels of replication (spatial and temporal) this could be more thoroughly investigated in the future.
Habitat occupancy by beach mice was high over the entire core range and appeared stable over the nine years of surveys. This was not initially expected given the extensive range reduction that occurred over the last several decades (Stout 1992). However, while occupancy is appropriate for understanding species-habitat relationships, it does not provide precise information about abundance. Estimation of abundance of beach mice is labor intensive, limiting the extent of spatial variation sampled, and methods such as trapping can be hazardous (Suazo et al. 2005). Monitoring beach mouse habitat occupancy is an important tool for conservation and may also provide information about changes in abundance (Wilkinson et al. 2012), although more information is needed relating beach mouse occupancy with abundance. Power analysis simulations showed that our long-term habitat occupancy survey could detect a decline in beach mouse habitat occupancy rate for the core population on CCBIC of 20% or greater between any two years. This is important because detecting problems early allows more time for managers to identify and mitigate the threats to populations (Caughley and Gunn 1996). Understanding the causes of the extensive reduction in range for the subspecies, and rapid extirpations from many sites, is important for recovery planning. The population is now much more concentrated geographically than it was before the recent reduction in range, and so is more susceptible to a single catastrophic event such as a hurricane (Oli et al. 2001). If the causes of recent local extirpations could be determined, some of the potential suitable habitat remaining within the historic range of the population might be suited to reintroduction attempts.
Supplemental Material
Please note: The Journal of Fish and Wildlife Management is not responsible for the content or functionality of any supplemental material. Queries should be directed to the corresponding author for the article.
Text S1. ColextFitter_fn.txt. An R script with accessory function to fit a set of occupancy models using the unmarked package.
Available: https://doi.org/10.3996/JFWM-24-023.S1 (2.58 KB TXT)
Text S2. Dynamic Occupancy Model Selection.txt. An R script with code to conduct dynamic occupancy model selection for Southeastern beach mouse on the Cape Canaveral barrier island complex (CCBIC) during winters of the years 2010–2014 using unmarked package.
Available: https://doi.org/10.3996/JFWM-24-023.S2 (24.6 KB TXT)
Text S3. Hierarchical Bayesian model selection.txt An R script with code to conduct dynamic occupancy model selection for Southeastern beach mouse on the Cape Canaveral barrier island complex (CCBIC) during winters of the years 2010–2018 using a hierarchical Bayesian model implemented in JAGS and R.
Available: https://doi.org/10.3996/JFWM-24-023.S3 (30.6 KB TXT)
Text S4. Hierarchical Bayesian Dynamic Occupancy Model for Southeastern beach mouse.txt. An R script with code to fit a dynamic occupancy for Southeastern beach mouse on the Cape Canaveral barrier island complex (CCBIC) during winters of the years 2010–2018 using a hierarchical Bayesian model implemented in JAGS and R.
Available: https://doi.org/10.3996/JFWM-24-023.S4 (22.6 KB TXT)
Text S5. sebm_Bayesian_data.txt. An R data file with data for the hierarchical Bayesian analysis of the dynamic occupancy model for Southeastern beach mouse on the Cape Canaveral barrier island complex (CCBIC) during winters of the years 2010–2018.
Available: https://doi.org/10.3996/JFWM-24-023.S5 (25.3 KB TXT)
Text S6. sebm_data.txt. An R data file with data for the unmarked analysis of dynamic occupancy model for Southeastern beach mouse on the Cape Canaveral barrier island complex (CCBIC) during winters of the years 2010–2018.
Available: https://doi.org/10.3996/JFWM-24-023.S6 (29.7 KB TXT)
Acknowledgments
This project was funded directly and by in-kind services by the NASA Environmental and Medical Contract (NEMCON Contract # 80KSC020D0023), United States Air Force, 45th Space Wing, Natural Assets Office at Cape Canaveral Space Force Station, United States Fish and Wildlife Service at Merritt Island National Wildlife Refuge and National Park Service at Canaveral National Seashore. We thank the Environmental Monitoring Branch of NASA’s John F. Kennedy Space Center (Lynn Phillips, Don Dankert, James Brooks, and Jeff Collins). We also thank the personnel at Merritt Island National Wildlife Refuge, Cape Canaveral Space Force Station, and Canaveral National Seashore for their support. We thank Carlton Hall, Karen Holloway-Adkins, B. Bolt, J. Provancha, and T. Kozusko for their helpful comments. We thank the anonymous Associate Editor and two anonymous peer reviewers for many helpful suggestions that greatly improved the manuscript.
Any use of trade, product, website, or firm names in this publication is for descriptive purposes only and does not imply endorsement by the U.S. Government.
References
Author notes
The findings and conclusions in this article are those of the author(s) and do not necessarily represent the views of the US Fish and Wildlife Service.