Study area
Bakersfield is located in Kern County in the southern San Joaquin Valley in central California, USA. The city currently has an area of 392 km2 and current human population is ca. 390,000. Average elevation is 124 m, with little topographic variation. Climate is characterized by hot, dry summers and cool winters with infrequent precipitation in the form of rain. Average high and low temperatures are 13.7 C and 3.9 C in December and 36.2 C and 21.4 C in July. Mean annual precipitation is 164 mm (National Oceanic and Atmospheric Administration 2020). Bakersfield is bounded by occupied SJKF habitat to the northeast and southwest (Cypher et al. 2013) with irrigated agriculture bordering the city elsewhere.
Habitat attributes
To assess habitat attribute preferences by SJKF in Bakersfield, we established a sampling grid over the city consisting of 357 1-km2 grid cells that contained at least 50% human development. We used satellite imagery to quantify the occurrence of 13 urban landscape attributes in each of the cells. The attributes were those identified in Deatherage et al. (2021) as being important for urban SJKF and are described in Table 1. To estimate land use composition and index feature abundance, we superimposed a point grid (10×10 grid, 100 points total) over each cell in Google Earth Pro. We used Google Earth Pro imagery dated 26 April 2018 at an eye altitude of 300-m above ground level to characterize grid points. We characterized each point by the attribute that best described the location of the point (i.e., the terrestrial land use type on which the majority of the point was located or closest to), and recorded if the point fell on a mature tree. If a point appeared to fall equally on two different landscape types, we split the proportion of the point equally between the attributes (0.5:0.5). The proportion of points in each land use in a given cell was considered to approximate the proportion of that cell consisting of that land use. Water bodies generally were small and were recorded as present or absent for each cell.
We tested for pairwise correlations between landscape attributes using Spearman’s Rank tests and adjusted the resulting P values using the method proposed by Legendre and Legendre (1998) to account for the inflated risk of a type I statistical error when running multiple tests. Correlations were completed in Minitab 19 statistical software (Minitab LLC, State College, PA) and an α-level of 0.05 was used to determine statistical significance.
Grid use by SJKF
We conducted surveys for SJKF in a sample of grid cells using camera stations. We used the randomization function in Excel 2010 (Microsoft, Redmond WA) to randomly choose 120 grid cells to potentially sample from the 357 available cells. Within each selected grid cell, we then identified locations (1) that were accessible to kit foxes and (2) where the risk of camera theft was low (i.e., locations with restricted public access or where a camera could be placed in a cryptic location). Kit foxes are able to access most locations within the urban landscape, but generally avoid high-density residential areas due to the presence of fences, walls, dogs, and high levels of disturbance (Frost 2005; Cypher 2010). Consequently, most camera stations were placed in locations such as school campuses, city or county storm water drainage basins, municipal facilities, churches, golf courses, private businesses (with owner permission), and undeveloped parcels.
Within each sampled grid cell, we employed an automated camera station design and methodology developed specifically to survey for kit foxes and other sympatric carnivores (Westall and Cypher 2017). We used Cuddeback Digital Black Flash IR cameras (Model 1255, Non Typical Inc., Green Bay, WI) that employ a “black flash” infrared LED flash that creates almost no light visible to humans and that take high-resolution images (20 megapixels). We secured the cameras to 1.2-m U-posts using zip-ties. At some locations where cameras might be discovered by the public, we placed the cameras in protective cases (“CuddeSafe” Model 3327, Non Typical Inc., Green Bay, WI) that were then secured with a cable lock to fences, trees, or other immobile structures. To attract foxes, we placed several drops of a scent lure (Carman’s Canine Call Lure, New Milford, PA) in front of the camera and on surrounding vegetation. We also staked a 163-ml can of cat food to the ground approximately 2 m in front of each camera using 30-cm nails. The cat food cans were perforated to allow scent to void but limit access to the food. The staked cat food cans functioned as a further attractant for foxes and also caused them to remain in the camera’s field of view for an extended period as foxes attempted to access the food in the cans.
The average diameter of a kit fox home range in Bakersfield is 1 km, based on an estimated mean home range size of 0.78 km2 (CSUS Endangered Species Recovery Program, unpublished data). Camera station locations typically were separated by at least 1 km to reduce the probability of a given kit fox being detected at more than one station. Because 97.1% of kit foxes are typically detected at camera stations within six nights (Westall and Cypher 2017), cameras were deployed at each location for seven nights. Images were then downloaded from each camera and examined for visits by kit foxes each of the seven days. Because kit foxes are primarily nocturnal, each new day started at 1200. The sampled grid cells with and without kit fox detections were recorded. This survey was conducted each summer from 2015 to 2019. However, only the 2015 results were used for the habitat analyses because an epidemic of sarcoptic mange (Cypher et al. 2017) reduced the abundance of SJKF in each subsequent year (see population estimation results below), and these latter results were less reflective of actual grid cell use by the foxes.
Suitable habitat
The quantity of suitable habitat for SJKF in Bakersfield was estimated using occupancy modeling based on the 2015 camera station survey results in conjunction with habitat attributes measured in each of the 1-km2 cells. We used kit fox detection histories in single-species, single-season occupancy models in program PRESENCE (Hines 2006) to produce probability estimates of kit fox occupancy (ψ) in relation to landscape covariates while accounting for detection uncertainty (MacKenzie et al. 2018). We assumed no unmodeled heterogeneity in our data and that occupancy state in each cell did not change over the seven days of the survey. To develop our candidate models, we fit models first with no occupancy covariates (null model), followed by a global model including all covariates considered most likely to influence kit fox occupancy, then with each covariate individually, and finally with pairwise additive combinations of covariates following three a priori modeling categories representing general landscapes that may influence kit fox occupancy: campuses or campus-like landscapes, roads, and commercial or commercial-like landscapes. We allowed detection parameters to vary for each survey day. Covariates that were strongly correlated at rs> ±0.60 were not included in the same multi-covariate model (Burnham and Anderson 2002).
We calculated a Pearson χ² goodness of fit statistic and an overdispersion factor, , in PRESENCE to evaluate the fit of our global model with 1,000 bootstraps (Burnham and Anderson 2002). We assessed goodness of fit significance at an α-level of 0.05 and considered > 1 as overdispersion in models (Burnham and Anderson 2002). We used Akaike’s Information Criterion (AIC) and beta (β) estimates to evaluate models and used quasi-corrected AIC (QAIC) to adjust for poor model fit and overdispersion (Burnham and Anderson 2002). We considered sample size/maximum number of model parameters (K) < 40 as small, and used QAICc to further adjust for small sample sizes (Burnham and Anderson 2002). We derived a model confidence set based on QAICc weights (w) of the most supported predictors of kit fox occupancy and calculated model-averaged estimates of occupancy and detection parameters as well as 95% confidence intervals across multiple competing models in Excel 2019 (Burnham and Anderson 2002).
We calculated a model-averaged estimate of kit fox occupancy in each 1-km2 cell in PRESENCE. We used the Natural Breaks function in ArcGIS Desktop 10.6 (ESRI, West Redlands, CA) to classify kit fox occupancy probability values into three categories of habitat suitability: high, medium, and low. This function uses the Jenks Natural Breaks algorithm to identify breaks in the data that best groups similar values together and maximizes the differences between classes (de Smith et al. 2021). We summed the number of cells in each classification to determine the total area of each habitat category within the study area. We mapped these results to display the distribution of the habitat categories across the city. Finally, we created a more fine-scale representation of SJKF habitat suitability in Bakersfield. We used the Polygon to Raster conversion tool in ArcGIS to create a raster with a cell size 100 m2 (10x10 m). This cell size approximates the area that was characterized by each grid point when measuring the habitat attributes. The ArcGIS Spatial Analyst Focal Statistics tool to recalculate the SJKF occupancy probability for each 100-m2 cell. For each cell, we calculated its value as the mean of cell values in a 1,000x1,000-m rectangular window around the cell. In this final map, each suitability category was divided into two to produce a representation that better depicts the uncertainty in precise boundaries between the suitability categories.
Carrying capacity and population estimation
We estimated a conceptual carrying capacity for SJKF in Bakersfield. We calculated an estimate two ways using the suitable habitat quantities derived from the occupancy analysis. For one estimate, we divided the quantity of highly suitable habitat by 0.78 km2, which is the average home range size for kit foxes in Bakersfield (CSUS Endangered Species Recovery Program, unpublished data). This provided an estimate of the number of home ranges that potentially could be present in Bakersfield. Each social group occupying a home range minimally consists of an adult male and an adult female (Cypher 2010). Therefore, we multiplied the number of home ranges by two to produce an estimate of adult carrying capacity.
For the second estimate, we conducted the calculation above. In addition, we divided the quantity of medium quality habitat by the average home range size to estimate the number of home ranges in this habitat. However, because this was lower quality habitat, we assumed that the number of foxes might be lower (e.g., either foxes require more area than the average home range size, not all areas were occupied, not all areas might have a pair of adult foxes). Thus, we multiplied the number of home ranges by one. The resulting number was added to that derived for the high quality habitat to produce a potential carrying capacity of adult foxes.
To assess the plausibility of the carrying capacity estimates, we estimated the number of SJKF in Bakersfield each year from 2015 to 2019 using the systematic camera station survey described above. We conservatively estimated the number of individuals based on the presence of multiple animals in a single image or distinctive physical features on individuals (e.g., reproductive structures, deformities such as notched ear or shortened tail, eartags or radiocollars, sarcoptic mange symptoms). We then tallied the number of grid cells with kit foxes and the total estimated number of kit foxes detected. Finally, we estimated the total number of SJKF in Bakersfield using a simple proportional equation: