Habitat of the endangered salt marsh harvest mouse (Reithrodontomys raviventris) in San Francisco Bay

Abstract Understanding habitat associations is vital for conservation of at‐risk marsh‐endemic wildlife species, particularly those under threat from sea level rise. We modeled environmental and habitat associations of the marsh‐endemic, Federally endangered salt marsh harvest mouse (Reithrodontomys raviventris, RERA) and co‐occurrence with eight associated small mammal species from annual trap data, 1998–2014, in six estuarine marshes in North San Francisco Bay, California. Covariates included microhabitat metrics of elevation and vegetation species and cover; and landscape metrics of latitude–longitude, distance to anthropogenic features, and habitat patch size. The dominant cover was pickleweed (Salicornia pacifica) with 86% mean cover and 37 cm mean height, and bare ground with about 10% mean cover. We tested 38 variants of Bayesian network (BN) models to determine covariates that best account for presence of RERA and of all nine small mammal species. Best models had lowest complexity and highest classification accuracy. Among RERA presence models, three best BN models used covariates of latitude–longitude, distance to paved roads, and habitat patch size, with 0% error of false presence, 20% error of false nonpresence, and 20% overall error. The all‐species presence models suggested that within the pickleweed marsh environment, RERA are mostly habitat generalists. Accounting for presence of other species did not improve prediction of RERA. Habitat attributes compared between RERA and the next most frequently captured species, California vole (Microtus californicus), suggested substantial habitat overlap, with RERA habitat being somewhat higher in marsh elevation, greater in percent cover of the dominant plant species, closer to urban areas, further from agricultural areas, and, perhaps most significant, larger in continuous size of marsh patch. Findings will inform conservation management of the marsh environment for RERA by identifying best microhabitat elements, landscape attributes, and adverse interspecific interactions.

Due to these direct and indirect human impacts, many marsh wildlife species are listed as endangered or threatened pursuant to the U.S. Endangered Species Act. It is likely that many changes to ecosystems caused by shifting climate patterns have not been recognized, but habitats can respond quickly and drastically in unpredictable ways, leading to substantial uncertainty about impacts to endemic wildlife.
Given this, uncertainty defining habitat characteristics is growing increasingly important for the conservation of biodiversity, especially for at-risk species.
Many resource management agencies have established a habitat-based approach, rather than an individual species approach, for recovery of marsh endemic species (DOI, 2013;Goals Project, 2015;USFWS, 2013), making it all the more important to understand habitat relationships of at-risk species presumed to be provided by a habitat-based approach. Our study focused on the salt marsh harvest mouse (Reithrodontomys raviventris, RERA, Figure 1), federally-, state-, and IUCN-listed as endangered (CDFG, 2015; USFWS, 2013; www.iucnr edlist.org/speci es/19401/ 22385344).
RERA is endemic to the wetlands of the San Francisco Bay (SFB), California, where marsh plant communities and obligate wildlife are at risk of habitat loss from local sea level rise over the coming decades (Thorne, Buffington, Elliott-Fisk, & Takekawa, 2015;Veloz et al., 2013).
RERA currently occupies <25% of its historic range because of habitat fragmentation and loss, with >80% loss of historic wetland habitat in SFB (Statham et al., 2016). The northern subspecies (Reithrodontomys raviventris halicoetes) is found along San Pablo and Suisun Bays, and the southern subspecies (R. r. raviventris) is found in South SFB with the subspecies divide somewhere in Central Bay (Shellhammer, 1989;Statham et al., 2016;USFWS, 2013).
However, detailed habitat associations of both subspecies are not well understood. Smith, Riley, Barthman-Thompson, Woo, et al. (2018) reviewed the RERA literature to date.  recommended research priorities to further RERA recovery, including the need for data synthesis to assess regional habitat associations and species interactions. On this basis, we collated several monitoring datasets to address this priority research need.
Several independent monitoring efforts have been conducted to assess SMHM abundance with varying sampling design and effort, very low SMHM captures, and lack of marking individual captures, such that traditional occupancy models would not be appropriate.
Therefore, we used a probabilistic analysis approach with Bayesian network (BN) modeling to determine RERA correlates with environmental and habitat parameters, and to evaluate relationships with other sympatric small mammal species.
Bayesian network models link variables with probabilities calculated using Bayes' Theorem (Korb & Nicholson, 2010;Koski & Noble, 2011) and are relatively robust to zero inflated data and collinearity.
We analyzed trap data and associated vegetation, patch (elevation and inundation), and landscape covariates, including presence of other small mammal species, to determine habitat associations of RERA. Explicitly, we aimed to assess (a) microhabitat and landscape attributes determining RERA presence, and (b) if RERA presence was associated with other small mammal species. We developed all BN models from empirical field and spatial data. This modeling approach can be used to project impacts on RERA from future changes in habitat conditions, disturbances, and climate change, and to inform management of these tidal marsh environments to facilitate species recovery and conservation.

| Study area
The SFB ecosystem in the central coast of California has a Mediterranean climate with warm dry summers and cool rainy winters. North SFB (38°08′N, 122°24′W), comprised of San Pablo and Suisun Bays, receives most of its freshwater from the Sacramento and San Joaquin Rivers in the form of summer snowmelt from the Sierra Nevada mountains. This region has a mixed semidiurnal tide with a mean tide range of 1.17-1.63 m and a mean spring tide range of 1.57-2.09 m (Parker, Callaway, Schile, Vasey, & Herbert, 2012). Most tidal marshes are found above mean tide level (MTL, defined as the arithmetic mean of mean high water and mean low water; Parker et al., 2012) whereas marsh plain elevations are closer to or above mean high water (MHW, defined as the average of all the high water heights observed over 19 years; Takekawa et al., 2013) (https ://shore line.noaa.gov/gloss ary.html).
The marsh plain is dominated by pickleweed (Salicornia pacifica, syn. Sacrocornia pacifica, formerly Salicornia virginica; Calflora, 2018), and the lower marsh area is dominated by Spartina foliosa, a native of the region.
Examples of local marsh inhabitants are salt marsh common yellowthroat (Geothlypis trichas sinuosa), San Pablo song sparrow (Melospiza melodia samuelis), and San Pablo vole (Microtus californicus sanpabloensis), state-listed species of special concern; California black rail (Laterallus jamaicensis coturniculus), a state threatened species; and California Ridgway's rail (Rallus obsoletus, formerly clapper rail, Rallus longirostris obsoletus). However, RERA is the only mammal species in the world that is entirely endemic to coastal marshes  and is found only in tidal marshes in SFB ( Figure 1).
It currently occupies <25% of its historic range because of habitat fragmentation and loss, with >80% loss of historic wetland habitat in SFB (Statham et al., 2016). The northern subspecies (R. r. halicoetes) is found along San Pablo Bay, and the southern subspecies (R. r. raviventris) is found in South SFB with the subspecies divide somewhere in Central Bay (Shellhammer, 1989;Statham et al., 2016;USFWS, 2013).
Sherman live traps (H.B. Sherman Traps, Inc.) were placed in three trap layout patterns (Table 2): grid, transect, and random placement. Grids varied in size, covering arrays of 5 × 5, 5 × 10, 7 × 7, and 2 × 25 traps, depending on local conditions and marsh area. Where local marsh conditions narrowed, not providing adequate area for a grid, one to three transects were used, each consisting of ten traps per transect. Random placements of traps were used only in the Tubbs Island Setback site because tidal marsh habitat and the upland transition zone were narrow with variable width such that the random placement allowed better coverage by which to associate RERA captures to habitat. No one site used all three trap layout patterns, and all sites used the grid or transect patterns except for Tubbs Setback (Table 2).
Traps were spaced at 10-m intervals following Jones, McShea, Conroy and Kunz (1996). This spacing was also consistent with findings from Bias and Morrison (1999) who reported that RERA moved short distances (mean 11.9 m) between consecutive 2-hr telemetry observations, and that RERA home ranges averaged just 2,133 m 2 .
Traps were initially deployed twice a year in spring and summer at Tolay Creek to capture a range of reproductive conditions. At other sites, traps were deployed once a year in late summer or early fall to maximize the capture of juveniles who would have been born in the spring of the same year. Although individual RERA movements patterns were generally highest in summer (June) and lowest in early winter (November), movements of pregnant females were low more consistently throughout the seasons (Bias & Morrison, 1999). Therefore trap distance was kept at a 10-m interval regardless of season.
Each trapping session consisted of three consecutive nights (with limited exploratory sampling conducted for four nights at Tolay Creek to increase detection of RERA). Traps were baited with a mix of crushed walnuts, birdseed, and mealworm (for insectivorous shrews) and opened in the evenings, checked each morning at sunrise, and closed during the day.
Polyester batting was placed within each trap to keep small mammals warm. Wooden shingles were placed on top of each trap to protect captured animals from exposure. Species identification, sex, age, mass (mg), reproductive condition, body length, tail length, and presence of wounds or parasites were recorded for all individuals. Reproductive condition in males was characterized by presence and development of the testes. Reproductive condition in females was characterized by the presence and development of mammaries and whether the animal was pregnant. Animals captured and identified to the genus Reithrodontomys also included records of tail width 20 mm from the base of the tail, hind foot length, ear length, venter coloration of tail and belly, bicoloration of tail, and behavior (e.g., aggressiveness). Individuals were marked by clipping fur with small scissors to identify recaptures.

| Microsite and patch covariates
We recorded microsite and patch conditions in terms of vegetation, topography, and inundation patterns at trap locations. Each trap location was defined by a set of covariates including trap data, location, plant species, elevation, distance to natural or anthropogenic features, and marsh patch size (Appendix S1: Table S1; see Metadata S1 for the list and definitions of covariates), as follows.

| Vegetation
We visually assessed percent cover of each vascular plant species within a 0.25 m 2 vegetation quadrat centered on each trap and measured the mean and maximum height (to nearest cm) of each species. In wetland restoration sites, habitat at each trap was described using the closest and the second-closest vegetation quadrats for Tolay Creek , Tubbs Setback (Woo, Takekawa, Rowan, Gardiner, & Block, 2007), Guadalcanal (Woo, Takekawa, Gardiner, Dembosz, & Bishop, 2009), and Benicia-Martinez (Woo, Bishop, & Takekawa, 2010). For our covariate data set, we identified from the vegetation plots the plant species or other cover categories having the three highest individual cover percentages, maximum heights, and average heights; we used all these covariates in the initial modeling, discussed below. Other cover categories recorded and use in initial modeling included presence and percent cover of green algae, brown algae, bare ground, litter, dead plants standing, and dead plants not standing.
Total plant cover in a quadrat could exceed 100% due to vegetation layering. We followed Jepson Flora Project (2018) for vascular plant nomenclature. See Figure 3 for examples of vegetation and environmental conditions.

| Topography
We measured topographic variation at each site by surveying marsh surface elevation along transects perpendicular to the major tidal sediment source, with survey points taken every 12.5 m, and transects separated by 50 m using a Real Time Kinematics (RTK) GPS (Leica Geosystems Inc.). We used the Geoid09 model to calculate orthometric heights from ellipsoid measurements (m, NAVD88; North American Vertical Datum of 1988) and projected all points to NAD83 UTM Zone 10 using Leica GeoOffice v7.0.1 (Leica Geosystems Inc.).
We then combined the elevation survey data to create a digital elevation model (DEM) at each site in ArcGIS 10.2.1 Spatial Analyst (ESRI 2011) with exponential ordinary kriging methods (5 × 5 m cell size) after adjusting model parameters to minimize the root-mean-square error (RMS). Using the DEM data, we were then able to assign a specific elevation to each trap location.

| Inundation
To determine inundation patterns and calculate site-specific tidal datums, we used a mix of available local data sources. We deployed water level data loggers (Model 3001; Solinst Canada Ltd.) at Corte Madera, Tolay Creek and Tubbs Island Setback Marshes (Takekawa, Thorne, Buffington, & Freeman, 2014;Takekawa et al., 2013). For other sites (Benicia-Martinez, Fagan, Guadalcanal), we used data from the nearest NOAA tidal station. Sites with locally collected data had one or two loggers that were placed at the mouth and upper reaches of second-order tidal channels to capture high tides and determine seasonal inundation patterns. Water level readings were collected every 6 min starting on the date of deployment and continued for at least one year to calculate tidal datums for each trapping location.
We used data from the lowest elevation logger at each site to develop local hydrographs and inundation rates, and we surveyed loggers with RTK GPS at the time of deployment and at each data download to correct for any vertical movement and to reference water levels to NAVD88 elevations. All raw water level data were corrected with a local time series of barometric pressure.
For Solinst loggers, we deployed independent barometric loggers (Model 3001; Solinst Canada Ltd.). We used water level data to estimate local tidal datums for all sites using procedures outlined in the NOAA Tidal Datums Handbook (NOAA 2003), and calculated only local mean high water (MHW) and mean higher high water (MHHW) because the loggers were positioned in the intertidal and therefore could not be used to compute lower datums. We then assigned MHW and MHHW values to each trap location based on the nearest local water monitoring location and its associated local tidal datums ( Figure 2).

| Correlation among variables
We first conducted bivariate Pearson correlations among the covariates to determine which variables may be highly correlated so as to exclude them from further modeling. Although the BN models we used are more tolerant of multicollinearity among covariates than are frequentist multivariate models (Pawson, Marcot, & Woodberry, 2017), it is still useful to narrow the set of covariates to reduce unnecessary model complexity.
For each of these two response variable outcomes, we developed 19 variants of BN models, thus 38 total models, to represent different combinations of covariates. We grouped covariates into logical sets representing trap sampling design, location (latitude-longitude), vegetation and cover, elevation, distance to natural or anthropogenic features, and marsh patch size, and then developed BN models using variables drawn from one or more of those sets (Table 3).

| BN model accuracy, complexity, selection, and sensitivity
We identified the best-performing models as those with low complexity and high accuracy (low overall prediction error. and low Type I error of false positive). The overall aim of denoting model calibration accuracy and complexity was to provide information by which to choose best models, which would be those with highest classification accuracy (lowest error) and lowest model complexity, akin to use of Akaike information criterion (AIC) metrics in frequentist statistical modeling (Akaike, 1973) that is strictly not applicable in BN modeling. Selection of best Bayesian models can be a somewhat subjective process (Hooten & Hobbs, 2015 We include spherical payoff because it provides complementary information to, and can be poorly correlated with, overall confusion error rate (Marcot, 2012). Lastly, we conducted sensitivity tests of the best models to determine the degree to which the predicted occupancy probabilities were sensitive to uncertainty about each covariate. We used sensitivity analysis procedures and metrics detailed in Marcot (2012) entailing calculations of entropy reduction which depict incremental responses of an outcome variable given incremental changes in each covariate. Results are useful for comparing relative sensitivity among covariates within a model, and relative sensitivity of a given covariate among models.

| Habitat relationships
We  with Tolay Creek, also furthest from urban areas; and Fagan was furthest from the bay and from levees ( Figure S1.2). Tolay Creek had the largest average habitat patch size, and Guadalcanal had the largest extended habitat patch size not impeded by barriers or channels >3 m wide ( Figure S1.3).

| Environmental and vegetation covariates
Of the 78 Pearson correlations r among 13 continuous environmental covariates (elevation, distance, and patch size variables), nine were statistically significant (p < .05) with |r| > 0.75 (Table S1.3). From these results, we eliminated three covariates (MHHW, distance to closest urban area and distance to closest agriculture), leaving the remaining 10 to include in the BN models in various combinations with low correlation. The 12 vegetation covariates were not part of the correlation analysis due to sparsity of these data (low numbers of plots in which all variables were present), but we included all of them in the BN models.

| Small mammals
Field trapping resulted in capture of nine species of small mammals (  Figure 2).

| Bayesian network models
We developed 38 variants of BN models consisting of 19 different combinations of covariates, each with two variations of trap-result response variables: RERA presence, models 1-19; and presence of each small mammal species, models 20-38 (Table 3; Figure S1.4).
Among the RERA presence models ( Among the all-species presence models (Table 7, Table S1.5), model complexity also varied widely among the model variants.   For RERA presence, the best-performing models were models 2, 6, and 7 (Table 6), and for all-species presence, these were models 25 and 26 (  Table S1.1 for covariate definitions).
Models 1 and 20, based on just study site location, sampling time, and sampling method, performed fairly well. However, they were more complex and incurred greater Type I error rates than the best-performing models. Models 2 and 21, based on just latitude and longitude, fared well. This indicted that, for the general study region, site location is a fair predictor of RERA presence, although these models provided no information on habitat or environmental conditions. However, models 3, 4, 22, and 23, based solely on vegetation, fared very poorly with high Type I errors. Clearly, something other than, or in addition to, vegetation-viz., elevation, distance to roads, and marsh patch size-influenced RERA presence. The remaining models were generally increasingly complex with various mixes of covariate sets, but improvements to calibration error rates were minor at best and not offset by the greatly increased levels of model complexity.
Cross-validation results (Table S1.6) were essentially the same as the calibration error analyses (Tables 6 and 7), conveying that the best-performing BN models were not overfit. Sensitivity analyses of best-performing models (Table S1.7) suggested the following regarding presence of RERA. RERA presence was more sensitive to longitude (easting) than to latitude (northing). Also, presence of all-species, compared with presence of RERA, was slightly more sensitive to distance to roads and much more sensitive to patch size.
None of the 19 BN models of RERA presence provided a classification confusion error rate of <19%. Four models resulted in no

| Habitat differences of salt marsh harvest mouse and California vole
Reithrodontomys raviventris and MICA generally spanned a similar range of site attributes ( Figure S1.4). However, some differences in distributions were suggested by results of unpaired two-sample t tests of covariate values at locations where RERA and MICA were both trapped (Table 8; Table S1.8). At sites where both species were trapped, values of ten covariates, denoting microsite conditions, differed significantly between the species. Attributes of trap sites with RERA captures differed from those with MICA captures in the following ways. RERA capture sites had significantly greater percent cover and a greater maximum height of the most dominant plant species (pickleweed). RERA sites also had a higher marsh elevation, higher MHW and MHHW levels (these two variables were significantly positively correlated; see Table S1.2), and higher elevation of trap location relative to MHW. RERA sites also had greater distance to agriculture, closer distance to nearest urban area, and greater marsh patch area not impended by barriers or channels >3m wide.

| D ISCUSS I ON
Our study found that most of the RERA captures were associated with sites having pickleweed as the dominant ground-cover category.
Fewer RERA sites had bare ground as the dominant ground-cover category, although about half of the RERA captures were associated with bare ground as the secondary ground-cover category.
Reithrodontomys raviventris have long been associated with dense wetland cover with pickleweed as preferred habitat (Shellhammer, 1989;Shellhammer et al., 1982). Shellhammer et al. (1982) Basson (2009) reported finding the species to be randomly distributed with no association with pickleweed salinity, pickleweed height, distance to levees, distance to dry or filled water bodies, percent cover of vegetation, or sympatric rodents. Differential use of pickleweed sites by sex and season, as well as movement across various conditions (Geissel et al., 1988), might explain the findings by Basson (2009), as well as those by Botti, Warencyia, and Becker (1986) of several RERA in areas lacking pickleweed in the Suisun Marsh region of SFB.
We found that RERA presence correlates positively with marsh elevation, MHW level, and elevation of the trap location compared to MHW, particularly as compared with distributions of sympatric MICAs. This is consistent with other findings that RERA remained within the marsh in tall vegetation rather than swim to upland areas for high tide refugia sites (Smith, Barthman-Thompson, Gould, & Mabry, 2014).
We found evidence of some degree of tolerance of RERA with levees and water (<200 m distance), and also with agricultural areas and paved roads (<2 km distance). This was consistent with Bias and Morrison's (1999) radio-telemetry study of RERA at Mare Island, Solano County, California (Figure 2), who concluded that the species crossed roads, levees, and canals. However, RERA in our study did not apparently select for such conditions, as they were also associated with larger marsh patch sizes not interrupted by levees and roads.
We captured nine species of small mammals and found no spe- We found that RERA presence, as compared with MICA presence, is characterized as being higher in marsh elevation, higher in MHW level, and higher in elevation compared with MHW. Also, as compared with MICA, RERA occurred in sites with greater height and percent cover of pickleweed, that were closer to urban areas but further from agricultural areas (these two variables were strongly negatively correlated; see Table S1.3), and that were larger in marsh patch size. Similarly, in a study of Spartina marshes along the southwest Atlantic coast of the United States, Canepuccia, Pascual, Biondi, and Iribarne (2015) found that occurrence of small mammals was related to vegetation cover and diversity, and that small mammal species composition varied by landscape context and configuration.
A Type II error may be more egregious if models predict RERA nonpresence than when this State-and Federally listed endangered species is actually present. That is, false absences (Type II error) of an endangered species have a greater potential negative consequence for conserving population viability than do false presences (Type I error). Predicting absence when the species is present could result in no specific habitat management being directed to such locations, potentially resulting habitat degradation or loss for the species. In such cases, models cannot substitute for on-the-ground surveys. This type of habitat association model is better suited to assess RERA presence than to conclude absence. We note, too, that we have not analyzed any temporal trends or differences in habitat associations over the course of the trapping seasons and years, assuming that the availability of environmental conditions in each trapping location remained constant over time. We also note that studies are lacking on identifying specific threats to population viability of the species, other than inferences made on threats of habitat changes, which aligns with the recommended research priorities to further RERA recovery, including range-wide population estimates, demographics, and dynamics .
Additional threats to small mammals in tidal marshes and potentially to RERA conservation may yet appear or confound recovery, such as uptake of environmental toxins as documented by Clark, Foerster, Marn, and Hothem (1992). BN models for habitat associations can be informative for tidal marsh restorations that may be able to accommodate higher marsh elevations in relation to inundation and larger marsh patches in their designs. Ultimately, management and conservation of RERA, particularly under threats of climate change, sea level rise, and anthropogenic alteration of habitats, may need to take an ecosystem-level approach (e.g., Stagg et al., 2016) that includes tracking increased inundation and flooding of pickleweed habitat from rising sea level (Field, Bayard, et al., 2016;Kirwan et al., 2010;Rosencranz et al., 2018).
The current study could not compare RERA occurrence with attributes of random-site locations, so resource selection functions of the mouse for particular habitat attributes could not be determined.
This could be a valuable topic for future studies, to better understand RERA habitat selection for guiding site management, especially considering landscape level changes associated with sea level rise and climate change. Further studies could also address existing or emerging threats to the species.

ACK N OWLED G M ENTS
We thank John Y. Takekawa