A hierarchical modeling approach to predict the distribution and density of Sierra Nevada Red Fox (Vulpes vulpes necator)

Abstract

Carnivores play critical roles in ecosystems, yet many species are declining worldwide. The Sierra Nevada Red Fox (Vulpes vulpes necator; SNRF) is a rare and endangered subspecies of red fox limited to upper montane forests, subalpine, and alpine environments of California and Oregon, United States. Having experienced significant distribution contractions and population declines in the last century, the subspecies is listed as at-risk by relevant federal and state agencies. Updated information on its contemporary distribution and density is needed to guide and evaluate conservation and management actions. We combined 12 years (2009–2020) of detection and nondetection data collected throughout California and Oregon to model the potential distribution and density of SNRFs throughout their historical and contemporary ranges. We used an integrated species distribution and density modeling approach, which predicted SNRF density in sampled locations based on observed relationships between environmental covariates and detection frequencies, and then projected those predictions to unsampled locations based on the estimated correlations with environmental covariates. This approach provided predictions that serve as density estimates in sampled regions and projections in unsampled areas. Our model predicted a density of 1.06 (95% credible interval = 0.8–1.36) foxes per 100 km² distributed throughout 22,926 km² in three distinct regions of California and Oregon–Sierra Nevada, Lassen Peak, and Oregon Cascades. SNRFs were most likely to be found in areas with low minimum temperatures and high snow water equivalent. Our results provide a contemporary baseline to inform the development and evaluation of conservation and management actions, and guide future survey efforts.

hierarchical modeling, occupancy modeling, remote cameras, Royle–Nichols model, Sierra Nevada Red Fox, species distribution model, Vulpes vulpes necator

Humans have contributed to global declines of many carnivore species. The distributions and abundances of many carnivores have been negatively affected by overexploitation, habitat loss, depletion of prey populations, direct persecution, and climate change (Ripple et al. 2014). Carnivores can have profound effects on the ecosystems they inhabit, and contracted distributions and decreased populations, including extirpations, can have cascading effects within local assemblages and ecosystems (Beschta and Ripple 2009; Prugh et al. 2009; Estes et al. 2011). For example, decreased abundance of large carnivores can trigger trophic cascades, where the numbers of smaller or subordinate carnivores and herbivores ordinarily limited by larger carnivores increase (Crooks and Soulé 1999; Prugh et al. 2009). Additional examples of trophic cascades include declines in bird populations and altered vegetation structure resulting from increased density of subordinate carnivores and prey species (e.g., Estes and Palmisano 1974; Beschta and Ripple 2009). The conservation status of carnivores can be difficult to determine because they tend to be rare, elusive, and often naturally occur at low densities (Gese 2001; Carbone and Gittleman 2002; Long et al. 2012). This poses a challenge to monitoring efforts, making it difficult to detect changes in population demographics or evaluate conservation status, and can obscure the efficacy of management actions.

It is important to delineate contemporary geographic distributions and estimate spatially explicit densities of species in order to identify habitat relationships, determine important areas to prioritize for conservation, and establish baselines for informing future management decisions. Quantifying distributions and densities of carnivores can provide a better understanding of how their populations may be affected by environmental change (Guisan and Zimmermann 2000; Parmesan and Matthews 2005; Elith and Leathwick 2009). In particular, global climate change threatens many carnivore species. Climate change can decrease the access of carnivores to prey or increase the cost of acquiring prey (Gormezano and Rockwell 2013; Pagano et al. 2018). Increased temperatures can also elevate physiological stress, increase the opportunity costs of foraging in extreme environments, or reduce reproductive success (Creel et al. 2016; Woodroffe et al. 2017; Rabaiotti and Woodroffe 2019). Climate change can alter interspecific interactions, potentially benefiting species better adapted to novel conditions (e.g., generalists) and decreasing the persistence of specialized species (Laliberte and Ripple 2004; Zielinski et al. 2017; Jensen and Humphries 2019; Peers et al. 2020). By understanding the current distribution and density of species located in areas sensitive to a changing climate, we can identify important habitat refugia and target management efforts to prevent further species losses.

Species inhabiting montane habitats may be particularly susceptible to adverse effects resulting from a changing climate as a consequence of upslope movement of novel plants and animals, inevitably resulting in range loss (e.g., Dullinger et al. 2012). It is valuable to identify and better understand the current distributions and habitat preferences of species, as such knowledge may inform their conservation. This is especially true for the Sierra Nevada Red Fox (Vulpes vulpes necator; hereafter SNRF), a rare subspecies that inhabits the high elevation, contiguous Cascade and Sierra Nevada mountain ranges (Grinnell et al. 1937; Aubry 1983; Aubry et al. 2009; Sacks et al. 2010). This subspecies is one of three montane red fox subspecies restricted to the western United States, including the Cascade Red Fox (V. v. cascadensis) and the Rocky Mountain Red Fox (V. v. macroura). These montane red foxes have several adaptations for surviving in snow (Grinnell et al. 1937; Aubry 1983; Cross et al. 2018), including thick winter pelage, smaller body size, and small toe pads covered by dense fur (Grinnell et al. 1937; Aubry 1983; Fuhrmann 1998).

Of the three montane red fox subspecies, the SNRF is thought to have experienced the largest population decline during the 20th century (Sacks et al. 2010). SNRFs were historically distributed in many upper montane, subalpine, and alpine areas in California and Oregon (Bailey 1936; Grinnell et al. 1937; Sacks et al. 2010), but their current extant distribution is thought to be substantially reduced and more fragmented than their historical distribution (Quinn et al. 2022). A suite of potential stressors prompted a federal status of Endangered for the Sierra Nevada Distinct Population Segment of the SNRF in 2021 (U.S. Fish and Wildlife Service 2021). Stressors included low genetic diversity and subsequent inbreeding depression (Sacks et al. 2010; Quinn et al. 2019, 2022), hybridization with nonnative red fox subspecies, and the effects of climate change (e.g., decreased prey availability and increased interspecific competition). Describing the current geographic distribution, habitat associations, and environmental relationships of SNRF could help inform future conservation actions (Sierra Nevada Red Fox Conservation Advisory Team 2021).

Previous research has indicated that minimum winter temperatures, precipitation, and land cover are important factors influencing the occurrence of SNRF (Cleve et al. 2011; Quinn et al. 2018). However, these previous efforts considered only portions of the historical range and implemented presence-only modeling approaches that did not include nondetection data or the effects of detection probability in predictions of occurrence. Estimates using presence-only analyses can overestimate the effects of covariates on the occurrence of a species (Yackulic et al. 2013), thereby limiting inferences on species distribution. Modeling distribution and density using contemporary occurrence and detection/nondetection survey data from across the distribution of species or subspecies can provide a more comprehensive prediction of their distribution and habitat relationships. Incorporating nondetection data alongside detection data in an occupancy modeling framework can also inform site selection for future monitoring efforts.

We modeled the spatial distribution and potential density of SNRF in California and Oregon. We compiled detection/nondetection, and detection-only SNRF survey data throughout its historical range. Data were collected through systematic and opportunistic carnivore surveys, as well as telemetry locations from radio- and GPS-collared animals. We incorporated these data into a species distribution model using a hierarchical implementation of the Royle–Nichols model (Royle and Nichols 2003) that allowed for inclusion of opportunistic presence-only data. Using this analytical framework, we: (1) investigated relationships among SNRF density and habitat covariates; (2) identified detection covariates that can increase efficacy of surveying for and monitoring this subspecies in the future; and (3) predicted the contemporary distribution and density of this subspecies throughout its historical range.

Materials and Methods

To investigate the distribution and density of SNRF in California and Oregon, we used SNRF data collected from 2009 to 2020 within the Level II Ecoregion 6.2: Northwestern Forested Mountains/Western Cordillera region (Commission for Environmental Cooperation 2009), which roughly equates to the historical range of the SNRF. We included three types of data: (1) infrared trail camera survey data; (2) telemetry locations of SNRFs; and (3) observations of red foxes not associated with a camera deployment or telemetry (e.g., confirmed locations collected opportunistically or during other survey efforts, such as roadkills, genetically identified scats, and live captures; Fig. 1). We necessarily assumed that any red fox detected in the mountainous regions of California and Oregon was of the Sierra Nevada subspecies and hereafter are referred to as SNRF. Genetic verification to subspecies requires genetic samples (e.g., scats, hair, urine, tissue) and analyses using nuclear markers that require high-quality DNA, infrequently available from red fox samples (Quinn et al. 2019). Previous analysis of the genetic detections encompassed by our data set showed red fox populations in the Pacific coastal states are highly structured, with extant SNRF populations in the historical range retaining a distinct genetic signature despite some recent introgression from low-elevation red foxes (Quinn et al. 2017, 2019, 2022; Quinn and Sacks 2022). It is possible that our assumption may have resulted in some low-elevation dispersing foxes belonging to other subspecies being falsely identified as an SNRF. However, to the extent that genetically verified detections and photographic detections of SNRFs were spatially correlated (Fig. 1, Supplementary Data SD1), such misclassifications were assumed to be rare and inconsequential at the scale of our distribution modeling (Quinn et al. 2022).

Fig. 1.

Map of Sierra Nevada Red Fox (Vulpes vulpes necator; SNRF) detections and survey effort in California and Oregon, United States, in relation to the distinct population segment (DPS) designated by the U.S. Fish and Wildlife Service (2021). In (A) the gray-filled squares indicate a detection of SNRF on camera and the open squares show the locations of camera deployments that did not detect SNRFs. In (B) the gray-filled circles and triangles indicate a detection of a SNRF from telemetry locations or a detection of SNRF through other methods, respectively. Telemetry studies were located around Lassen Peak and Sonora Pass. In (C) the gray-filled diamonds indicate a genetically verified detection of a SNRF from high-quality genetic samples (e.g., scats, hair, urine, tissue) and analyzed using nuclear markers (Quinn et al. 2017, 2019, 2022; Quinn and Sacks 2022). In (A), (B), and (C), the gray dashed line is the 30-km buffer encompassing all survey locations used to define the state space for our modeling efforts, and the black dots show locations of major cities for reference. The DPS designated by the U.S. Fish and Wildlife Service is indicated in dashed black line within the southernmost portion of the historical range in both figures (U.S. Fish and Wildlife Service 2021).

SNRF detection and nondetection data collection with infrared trail cameras

Infrared trail cameras (henceforth, ‘cameras’) were deployed widely in upper montane, subalpine, and alpine areas within the historical range of SNRF and in lower-elevation habitats for surveys targeting other carnivores throughout California and Oregon (Fig. 1A; e.g., Hatfield et al. 2020). These surveys were completed by many other researchers and, thus, survey protocols, goals, and durations varied. Although many of these projects were explicitly targeted toward SNRF, some survey data used in our analyses were focused on documenting the distribution of other carnivores occurring at lower elevations, e.g., Fisher (Pekania pennanti; Barry et al. 2021). By including survey data from regions not historically believed to be occupied by SNRFs, we were able to incorporate additional nondetection data to better understand the landscape covariates associated with their current distribution. In general, cameras were deployed in regions of interest, some with bait (e.g., chicken leg) and/or lure (e.g., Gusto; Minnesota Trapline Company, Pennock, Minnesota), and checked either repeatedly during the length of their deployment to refresh bait or lure, or visited only at the end of their deployment when cameras were being retrieved. Camera deployments often spanned multiple seasons. For each camera survey, we recorded the length of deployment and any days when cameras were not functioning properly (e.g., covered by snow, dead batteries, disturbed by wildlife, etc.), whether bait or scent lure were deployed as part of the survey, and if SNRFs were detected at each site.

SNRF capturing and collaring efforts

One female SNRF at Sonora Pass in 2015 and six females and one male in Lassen Volcanic National Park were captured between 2018 and 2020 using a combination of custom box traps, modified log cabin traps, and Tomahawk wire traps (107 × 38 × 51 cm; model 109.5; Tomahawk Live Trap Company, Tomahawk, Wisconsin). Upon capture, SNRFs were either anesthetized with a 5:1 mixture of ketamine hydrochloride (15 mg/kg) and xylazine (3 mg/kg) administered using a hand syringe or physically restrained by hand without sedation. We took tissue and hair samples for genetic analysis, marked individuals with unique colored ear tags or a Passive Integrated Transponder (Trovan Inc., United Kingdom) tag inserted into the scapular region, and fit them with a GPS radio collar, for example, KiwiSat303 (Sirtrack, Havelock North, New Zealand), Litetrack 130-150 (Lotek Wireless, Newmarket, Ontario, Canada), or TGW-4170-4 (Telonics, Mesa, Arizona). The collar compatible with Argos telemetry (KiwiSat303) was programmed to emit a signal every 45 s when activated; to prolong battery life, we scheduled collars to be active for one of every three days, cycling for 4 h on and 2 h off. Location acquisition rates for Argos telemetry depend on the orbital path of Argos satellites as well as local or environmental factors (e.g., topography, speed of movement), but on average we received 1–5 ‘high-quality’ fixes (location classes of two or three) for every 16 h of activity.

SNRF opportunistic detections

We collected opportunistic detections of SNRFs using hiking surveys, scent detection dog team surveys, visual observations, and hair snares within the elevation range of SNRF. Biological samples collected noninvasively at these locations, such as scat and hair, were molecularly confirmed to be from SNRFs through the course of other monitoring or population genetic studies (e.g., Hiller et al. 2015; Quinn et al. 2018, 2019, 2022). Samples were identified to species using a conserved segment of mammalian mitochondrial DNA with laboratory procedures for DNA extraction and sequencing described in detail in those previous studies. Briefly, we used previously published primers (Perrine et al. 2007) to PCR amplify and Sanger sequence a 354-bp portion of the mitochondrial cytochrome b gene and checked for >98% sequence match to previously accessioned SNRF haplotypes in GenBank NCBI (all detected sequences matched previously documented SNRF haplotypes 100%).

Spatial covariates

We created a state space to delineate the geographic range of our model. The results of species distribution modeling can be biased by the area of inference (Elith et al. 2010). Average home range size of SNRFs varies from 22 to 126 km² with a mean of 66 km² (Sierra Nevada Red Fox Conservation Advisory Team 2021). To limit the state space that we were modeling to areas where we had detection and nondetection data, we buffered all locations of the data we compiled by 30 km. We then created a sampling area within this state space by dividing it into a grid of 10 × 10 km squares to investigate the distribution and landscape associations of SNRFs in California and Oregon.

SNRFs are high-elevation specialists (Aubry et al. 2009) and, consequently, their distribution is thought to be strongly correlated with cold temperatures, snow–and more specifically–areas of increased primary productivity associated with upper montane, subalpine, and alpine environments (e.g., montane meadows). Thus, we incorporated minimum temperature, snow water equivalent (the amount of water contained within the snowpack), and normalized difference vegetation index (NDVI) as model covariates to estimate relationships between landscape characteristics and presence of SNRFs to predict their distribution and density. Daily minimum temperature data were acquired from the PRISM Climate Group at Oregon State University at a resolution of 800 m (PRISM Climate Group, Oregon State University, http://prism.oregonstate.edu). Daily snow water equivalent data were accessed from the Snow Data Assimilation System (SNODAS; National Snow and Ice Data Center) at a resolution of 1,000 m. NDVI values measure the variability in landscape greenness, providing a proxy for primary productivity (Pettorelli et al. 2005), and were calculated by NASA at a resolution of 250 m every 16 days. We aggregated NDVI data to 1,000 m using the raster package’s v. 3.4-10 ‘aggregate’ function by calculating the average of a 3 × 3 moving window and ‘resample’ function using bilinear interpolation (Hijmans 2021) in R v.4.0.2 (R Core Team 2020). We calculated average values for each covariate for all days for all years (2009–2019) to create a composite value that reflected the variation in each covariate. We standardized each covariate to have a mean of 0 and a standard deviation of 1 to indicate areas in California and Oregon that had comparatively high or low values of each covariate. We tested for correlations among minimum temperature, snow water equivalent, and NDVI across all grid cells within the sampling area; Pearson’s correlation coefficients were <|0.6|.

Detecting SNRFs is difficult due to their being rare, elusive, and persisting in areas that are challenging to access. An additional objective of this work was to determine how to improve detection probabilities when surveying for SNRF. Bait and lure were used at some of the camera deployments to increase detection probabilities of foxes. We investigated whether the presence of bait or scent lure, and the slope, aspect, and land cover at the camera site might influence the probability of detection. We extracted slope and aspect at the location of each camera from the USGS 1-Arc Second National Elevation Dataset. We cosine-transformed aspect (cosine(radians(aspect))) to use this variable as a measure of north–south angle of the camera location (a measure of sun exposure; North = 1 and South = −1). We standardized slope and aspect to have a mean of 0 and a standard deviation of 1. We also extracted the land cover classes at each camera deployment from the National Land Cover Dataset (Dewitz 2019) and determined whether or not each camera was deployed on open ground (where vegetation accounted for <15% of total cover) or in forest to estimate the effects of local habitat type on detection probability.

Species distribution and density modeling

We developed a species distribution model fit in a hierarchical framework to predict the distribution and density of SNRFs across their historical range. We coded a hierarchical model because detection data were collected from systematic remote camera deployments as well as from opportunistic and additional systematic survey efforts (e.g., roadkill collections, scat detection surveys, telemetry, etc.). We used a Royle–Nichols estimation (Royle and Nichols 2003) hierarchically linked to a spatially replicated occupancy modeling framework (MacKenzie et al. 2002). A main assumption of this approach is that the probability of detection at a specified location is a function of local density, and is described as:

$${\displaystyle \begin{array}{l}{\displaystyle P_i=1-{(1-r)}^{N_i},}\end{array}}$$

where P_i is the detection probability at location i, r is the binomial sampling probability that an individual of the focal species is detected, and N_i is the density (number per area) of species N at site i. Thus, as local density increases, so does the probability of detecting an individual at that location. We used this equation as a foundation for our modeling, and then modified it to predict the distribution and density of SNRFs throughout California and Oregon. We estimated a single potential distribution and density of foxes that reflected the 12 years of data that were collected. Due to data limitations, we were unable to investigate temporal changes in this distribution and density across years or effects of seasonality.

We estimated the distribution and density of SNRFs in 2,365 grid cells within the sampling area of our state space by identifying the grid cell j that each camera i was deployed in. We then modeled detection and nondetection data ycam of SNRF at camera i as:

$${\displaystyle \begin{array}{l}{\displaystyle yca{m}_i\sim \mathrm{B}\mathrm{e}\mathrm{r}\mathrm{n}\mathrm{o}\mathrm{u}\mathrm{l}\mathrm{l}\mathrm{i}(P_i)}\end{array}}$$

$${\displaystyle \begin{array}{l}{\displaystyle P_i=1-{(1-r_i)}^{N_{grid(i)}},}\end{array}}$$

where P_i is the average detection probability, r_i is the binomial sampling probability that a fox is detected during camera i’s deployment, and N_grid(i) is the predicted number (# per 100 km²) of SNRFs in the grid cell j where camera i was located. Data on the temporal deployments of cameras were limited, so we were unable to investigate the effects of seasonality on distribution, density, or detection probability. We predicted that lure, bait, slope, aspect, the length of the survey, and land cover type (open vs. forest) would influence the probability of detection. Thus, we modeled the logit-linear binomial sampling probability r at camera deployment i as:

$${\displaystyle \begin{array}{l}{\displaystyle \mathrm{L}\mathrm{o}\mathrm{g}\mathrm{i}\mathrm{t}(r_i)=\beta 0+\beta 1\ast lur{e}_i+\beta 2\ast bai{t}_i+\beta 3\ast slop{e}_i}\\ +\beta 4\ast aspec{t}_i+\beta 5\ast day{s}_i+\beta 6\ast ope{n}_i+\beta 7\ast fores{t}_i,\end{array}}$$

where the probability of detecting a SNRF is a function of an intercept (β0), effects of lure and bait (β1, β2), effects of the slope and aspect of where the camera was deployed (β3, β4), the number of days the camera was deployed (β5), and whether the camera was deployed on open ground or in forest (β6, β7). Bait, lure, and land cover type were coded as binary variables indicating whether or not bait or lure were applied during camera surveys or if the camera was deployed in an open or forested part of the landscape. Slope, aspect, and days deployed were coded as continuous variables and standardized to have a mean of 0 and a standard deviation of 1.

To determine habitat associations and extrapolate the predicted distribution and density of SNRFs throughout our modeled state space, we included habitat covariates associated with each grid cell j. We modeled the log-linear predicted number of SNRFs N per 100 km² grid cell j as:

$${\displaystyle \begin{array}{l}{\displaystyle N_j\sim \mathrm{P}\mathrm{o}\mathrm{i}\mathrm{s}\mathrm{s}\mathrm{o}\mathrm{n}({\lambda }_j)}\end{array}}$$

$${\displaystyle \begin{array}{l}{\displaystyle \mathrm{L}\mathrm{o}\mathrm{g}({\lambda }_j)=\mathrm{a}0+\mathrm{a}1\ast sno{w}_j+\mathrm{a}2\ast NDV{I}_j+\mathrm{a}3\ast mintem{p}_j}\\ +\mathrm{a}4\ast sno{w_j}^2+\mathrm{a}5\ast NDV{I_j}^2+\mathrm{a}6\ast mintem{p_j}^2,\end{array}}$$

where fox abundance N in each grid cell j is modeled as a function of an intercept (a0), snow water equivalent (a1), NDVI (a2), minimum temperature (a3), and the quadratic effects of these factors (a4, a5, a6). We included quadratic effects to test the hypothesis that for each of these main effects there might also be a threshold above or below which habitat is either suitable or no longer suitable for SNRFs.

Using the predicted density of SNRFs in each grid cell N_j, we used the step function (a test for whether or not N_j ≥ 0) to calculate the binary variable z indicating whether or not ≥1 SNRFs were estimated to be in each grid cell j as:

$${\displaystyle \begin{array}{l}{\displaystyle z_j=step(N_j-1),}\end{array}}$$

where z_j = 1 when a grid cell is estimated to have ≥1 SNRFs, and otherwise z_j = 0. We then identified the grid cells j that contained additional nonsystematically collected fox detections yopp in grid cell j as:

$${\displaystyle \begin{array}{l}{\displaystyle yop{p}_j\sim \mathrm{B}\mathrm{e}\mathrm{r}\mathrm{n}\mathrm{o}\mathrm{u}\mathrm{l}\mathrm{l}\mathrm{i}(z_j).}\end{array}}$$

Thus, grid cells where SNRFs were detected that were not collected from a camera were coded to be occupied by SNRFs and were modeled accordingly. This variable can be explicitly used to estimate the probability of occupancy by SNRFs in each grid cell. Full model code is presented in Supplementary Data SD2.

Predicting the range and density of SNRFs in California and Oregon

We used the posterior distributions of z_j to predict the distribution of SNRFs in California and Oregon. We calculated the mean occupancy probabilities of z_j and applied focal statistics on this predicted surface to calculate the mean occupancy of SNRFs within a window of 3 × 3 grid cells surrounding each grid cell j using the raster package v. 3.4-10 (Hijmans 2021) in R v. 4.0.2 (R Core Team 2020). We then determined the windows throughout California and Oregon that were predicted to have ≥0.3 probability of occupancy by SNRFs and used kernel smoothing to determine the predicted range. We used 3 × 3 grid window and a ≥0.3 threshold after visually inspecting predicted distributions using windows of 1 × 1, 3 × 3, and 5 × 5 and thresholds from 0.1 to 0.5 at 0.1 intervals. We determined that the 3 × 3 grid window and ≥0.3 threshold excluded small areas that were insufficient to support a population of SNRFs and best reflected the data from available surveys. To predict localized densities of SNRFs, we extracted the spatially explicit predicted densities in the Lassen Peak region, the Sierra Nevada Range, the Oregon Cascades, within the USFWS designated distinct population segment (DPS), and throughout California and Oregon.

Model fitting and assessment

We fit our distribution and density model using the Markov chain Monte Carlo (MCMC) methods of JAGS v.4.2.0 (Plummer 2003). We used uninformative prior distributions for all estimated parameters, including a uniform distribution with a minimum of −20 and a maximum of 20 for intercepts and a normal distribution with a mean of 0 and precision of 0.000001 for fixed effects. We calculated parameter estimates from 15,000 MCMC samples, taken from three chains run for 50,000 iterations, thinned by 10, following a burn-in of 100,000. We assessed model convergence by examining trace plots and R values for convergence (Gelman and Hill 2006; Gelman et al. 2013). We determined significance of covariates as percent probabilities from posterior distributions, which we calculated as the percent of posterior draws greater or less than 0, depending on the sign of the median value. We present the medians and 95% credible intervals (CIs) for our results.

Results

We compiled data from 10,616 remote cameras in California and Oregon that were deployed between February 2009 and December 2019 (Fig. 1A). Cameras were deployed for an average of 43.5 ± 79.8 (SD) days, and 222 cameras (2%) detected SNRFs. We also collected 1,326 opportunistic detections of SNRFs from their sign (e.g., scat, urine, hair), and 8,908 GPS telemetry locations from eight collared SNRFs (Fig. 1B). Photographic, opportunistic, and telemetry locations of SNRFs were located an average of 2,085 ± 8,541 (SD) m from a genetically verified detection (Fig. 1), which in turn were associated with one of the three biological populations previously determined to be composed of ‘pure’ or admixed SNRF ancestry–Sierra Nevada, Oregon Cascades, and Lassen Peak (Quinn et al. 2022). SNRFs were detected at a mean elevation of 2,113 m above sea level (asl; range = 1,200–3,681 m asl) at cameras, 2,335 m asl (range = 1,431–3,360 m asl) at telemetry locations, and 2,589 m asl (range = 694–3,731) at opportunistic observations (Fig. 1).

All three MCMC chains converged and $\widehat{R}$ values were <1.01 for all posterior parameter estimates. We found strong effects of minimum temperature and snow water equivalent on SNRF density. The linear and quadratic effects of average yearly minimum temperature had a 100% probability of negatively influencing the density of SNRFs (Table 1). These results indicate that the optimal average annual minimum temperature for SNRFs was −1.2°C (Fig. 2A). There was a 100% probability that the linear effect of average yearly snow water equivalent had a positive effect on SNRF density and a 79% probability that the quadratic effect of average yearly snow water equivalent had a negative effect on SNRF density (Table 1). The predicted density of SNRFs increased exponentially with an increase in average yearly snow water equivalent (Fig. 2B). The linear term for average yearly NDVI had an 82% probability of positively influencing SNRF density, and the quadratic effect had a 63% probability of negatively influencing SNRF density (Table 1). Increasing values of average yearly NDVI (i.e., primary productivity) were associated with greater SNRF density (Table 1; Fig. 2C).

Table 1.

Posterior parameter estimates from the Sierra Nevada Red Fox (Vulpes vulpes necator; SNRF) species distribution model. The predicted effects of covariates on SNRF density and the probability of detecting foxes at remote camera stations are presented on the log and logit scales, respectively. Posterior parameters for covariate effects with >90% probability of influence are indicated in bold. CI = credible interval; NDVI = normalized difference vegetation index.

	Parameter	Mean	SD	Low CI	Median	High CI
Density	Intercept	−3.808	0.533	−4.895	−3.778	−2.837
	Snow water equivalent	0.557	0.175	0.227	0.554	0.913
	Snow water equivalent2	−0.019	0.023	−0.067	−0.018	0.024
	Minimum temperature	−3.694	0.788	−5.254	−3.671	−2.231
	Minimum temperature2	−1.167	0.290	−1.765	−1.151	−0.642
	NDVI	0.276	0.304	−0.318	0.278	0.860
	NDVI2	−0.052	0.149	−0.348	−0.052	0.235
Detection	Intercept	−2.152	0.276	−2.700	−2.150	−1.618
	Lure	0.086	0.305	−0.507	0.083	0.687
	Bait	−0.643	0.281	−1.190	−0.645	−0.087
	Slope	−0.239	0.100	−0.439	−0.239	−0.044
	Aspect	−0.004	0.077	−0.153	−0.004	0.146
	Survey days	0.200	0.051	0.102	0.200	0.301
	Barren land	0.356	0.488	−0.646	0.375	1.281
	Forested land	−0.109	0.202	−0.504	−0.110	0.288

Fig. 2.

Predicted effects of minimum temperature (A), snow water equivalent (B), and normalized difference vegetation index (NDVI) (C) on Sierra Nevada Red Fox (Vulpes vulpes necator) density in California and Oregon, United States. The solid lines indicate the median effect of each of these covariates and the dashed lines indicate the 95% credible interval around this effect, while holding the other habitat covariates constant.

Our results also suggested that survey methods and site characteristics influenced the probability of detecting SNRFs on our cameras. The intercept for detection probability was small (Table 1), and many of the covariates we included in our model influenced detection probability. Using bait had a 99% probability of decreasing the detection probability of SNRFs (Table 1). Survey duration had a 100% probability of increasing detection probability (Table 1, Fig. 3); that is, the longer the survey, the more likely a camera was to detect SNRFs in locations where they occurred. The slope at camera sites had a 99% probability of negatively influencing detection probability; SNRFs were more likely to be detected at cameras deployed on flatter rather than steeper slopes (Table 1). Placing cameras in areas of open ground had a 77% probability of increasing the detection probability of SNRFs, and placing cameras in forests had a 70% probability of decreasing detection probability of SNRFs (Table 1). Aspect of camera site and including lure at survey sites did not influence detection probability (Table 1).

Fig. 3.

Predicted effect of the number of sampling days on the cumulative detection probability of Sierra Nevada Red Foxes (Vulpes vulpes necator) at cameras placed on open ground and on a slope of 0°. Here, the solid line indicates the predicted median effect and the dashed lines indicate the predicted 95% credible interval around this effect.

The predicted occupancy of SNRFs varied throughout California and Oregon, but was generally limited to upper montane, subalpine, and alpine regions with high amounts of snow water equivalent and low minimum temperatures (Fig. 4; Supplementary Data SD3 and SD4). The predicted geographic distribution of SNRFs in California and Oregon was 22,926 km² (Supplementary Data SD5), comprised of 12,238 km² in California, and 10,687 km² in Oregon (Fig. 5). The model predicted a median density of 1.11 (95% CI: 0.76–1.48) individuals per 100 km² in the 10,687 km² Oregon Cascades region, 0.96 (0.65–1.289) in the 9,077 km² Sierra Nevada region, 1.14 (0.89–1.42) in the 3,161 km² Lassen Peak region, and 1.16 (0.79–1.55) in the 4,908 km² U.S. Fish and Wildlife Service DPS.

Fig. 4.

The predicted mean (A) and standard deviation (B) of Sierra Nevada Red Fox (Vulpes vulpes necator; SNRF) occupancy in California and Oregon derived from the output of our hierarchical formulation of Royle–Nichols (2003). Darker colors indicate comparatively higher values of each. The dashed gray line indicates the extent of our modeling region and the solid black outline indicates the historical range of SNRFs (Grinnell et al. 1937; Sierra Nevada Red Fox Conservation Advisory Team 2021). Regions within the historical distribution with a value of 0 in (B) indicate grid cells where foxes were detected. These files are included in Supplementary Data SD2 and SD3.

Fig. 5.

The predicted distribution of the Sierra Nevada Red Fox (Vulpes vulpes necator) in California and Oregon, United States, based on the results of our species distribution modeling efforts in reference to the historical distribution displayed in light gray (Grinnell et al. 1937; Sierra Nevada Red Fox Conservation Advisory Team 2021) and distinct population segment (DPS) designated by the U.S. Fish and Wildlife Service displayed as a dashed black outline (U.S. Fish and Wildlife Service 2021). The current predicted distribution that we calculated based on Royle–Nichols (2003) modeling and kernel smoothing is shown in dark gray and included as Supplementary Data SD3.

Discussion

By combining data from systematic camera surveys, telemetry locations, and opportunistic observations of SNRFs into a novel hierarchical model, we were able to predict landscape conditions associated with SNRF occurrence and density and refine survey methods that may increase detection probabilities during future surveys. This study represents the most comprehensive analysis of SNRF survey data to date, and provides a range-wide model of a putative contemporary distribution and density of SNRFs that can be updated and improved as additional data become available. Our results have important implications for the conservation and recovery of the SNRF by documenting range contractions and local extinctions compared to putative historical distributions, quantifying the importance of average yearly minimum temperature and snow water equivalent to their occupancy, identifying priority areas for future surveys, and offering an empirical baseline to measure future conservation actions.

The predicted contemporary distribution of the SNRF we calculated is 37% the size of their presumed historical distribution–22,926 km² in the current study compared to the historically defined distribution of 63,640 km² (Grinnell et al. 1937; Sierra Nevada Red Fox Conservation Advisory Team 2021). The predicted distribution that we calculated (10,686 km²) is 41% the size of the putative historical range in Oregon (26,348 km²) and 33% (12,239 km²) the size of the putative historical range in California (37,292 km²). Within California, the predicted range in the Sierra Nevada region is 37% the size of the putative historical range (a reduction from 24,776 to 9,077 km²), the Lassen Peak region was 55% the size of the putative historical range (a reduction from 5,710 to 3,161 km²), and our results suggest a local extinction of SNRFs around Mt. Shasta (a reduction from 6,806 to 0 km²). It is important to note that the putative historical distributions were estimated using different techniques than our study (i.e., expert opinion in the historical distribution estimation), and explicit or direct comparisons of these distributions should be done with care. Specifically, the differences in distribution that we present may not represent range contractions for the SNRF distribution in California and Oregon over time, but rather the historical distributions might have been overestimated, for example, the polygons mapped in California (Grinnell et al. 1937) were labeled at the time as the ‘assumed general range.’ Further, by using a threshold of ≥0.3 for the current distribution, there are additional areas with lower estimated occupancy that may not be included in the distribution that we calculated that are within the true distribution of the SNRF.

The differences in the SNRF distribution that we estimated based on contemporary detections and the distribution presented by Grinnell et al. (1937) may correspond to their purported substantial population declines in the last 100 years (Grinnell et al. 1937; Aubry 1983; Perrine et al. 2010; Sacks et al. 2010). The causes of historical SNRF declines are unclear but thought to have included trapping and poisoning in the last two centuries (Perrine et al. 2010; Sierra Nevada Red Fox Conservation Advisory Team 2021). Contemporary climatic warming may be further exacerbating these declines. The geographic distributions of other high-elevation species, for example, the Pacific Marten (Martes caurina) and Wolverine (Gulo gulo) have contracted during this time, which is believed to be a result of a changing climate (Laliberte and Ripple 2004; Parmesan and Matthews 2005). Thus, areas where SNRFs are known to presently occur may represent important refugia and should be prioritized for protection and monitoring, especially because they are naturally rare even in environments with the most suitable conditions (Grinnell et al. 1937; Aubry et al. 2009; Sacks et al. 2010).

Our predictions of SNRF distribution and density should be interpreted as putative because they are based on extrapolations from well-surveyed locations to locations with less survey effort. In regions such as the southern Sierra Nevada (e.g., Sequoia and Kings Canyon National Parks; Fig. 4), where survey effort was low and evidence of SNRF presence is currently lacking, predicted density should be interpreted as the potential density for SNRFs given current landscape conditions. The value of such predictions is in highlighting locations with suitable conditions for SNRFs that can be prioritized for subsequent surveys. Moreover, if surveys ultimately yield no evidence of fox presence, such locations may serve as potential reintroduction sites. Conversely, if SNRF numbers increase in the future, they could potentially occur over a broader set of habitat conditions than observed in the present study. Thus, just as model predictions should not be used as estimates of current density where foxes have not been confirmed to occur, these predictions should also not be used to limit recovery goals to the landscape conditions described in our study.

A critical caveat of our models is that we assumed all detections of red foxes on cameras were representative of ‘pure’ SNRF ancestry. Recent genetic analyses indicate that all SNRF populations have experienced low-level introgression from other red fox subspecies, with the highest levels of admixture observed in the Sierra Nevada (Quinn et al. 2022). Genetic admixture from low-elevation subspecies of red foxes or the misidentification of a nonnative or low-elevation red fox for a SNRF could have resulted in an inflated predicted distribution or inaccurate habitat associations (e.g., Aubry et al. 2017). The potential of a predicted distribution modeled using possible detections of other subspecies of red foxes emphasizes the magnitude of our estimated range contraction and local extinctions and the need for more intensive genetic analyses (Quinn et al. 2018).

We found that average yearly minimum temperature and snow water equivalent were correlated with the density of SNRFs. These results are consistent with other distribution models for Sierra Nevada and other montane red foxes (Cleve et al. 2011; Akins 2016; Quinn et al. 2018). The quadratic effect for average yearly minimum temperature indicated that SNRF density was highest in a narrow threshold of temperatures with a peak at −1.18°C (Fig. 2A), which corroborates previous research (Quinn et al. 2018). This indicates that areas with average yearly minimum temperatures in this range were correlated positively with SNRF density. Further, SNRF density was predicted to increase as snow water equivalent increased. We do not know the exact mechanism for this correlation, but one hypothesis is that deep snow may limit interspecific competition with species including the Coyote (Canis latrans) and Gray Fox (Urocyon cinereoargenteus; Sierra Nevada Red Fox Conservation Advisory Team 2021; U.S. Fish and Wildlife Service 2021). Given the effects of climate change in montane environments, and the likely importance of cool temperatures and snow for SNRFs, increases in daily minimum temperatures and decreases in snow may negatively affect the distribution of SNRFs. The correlation of NDVI with SNRF density was limited, possibly because of correlations with other variables (e.g., minimum temperature). Holding all other variables equal, however, we found that SNRFs were more likely to be found in areas of comparatively high NDVI. This result supports the hypothesis that subalpine and alpine meadows may be important habitat characteristics for SNRFs, likely because they are important habitats for small mammals, which they use as prey.

We also empirically identified how camera survey design influenced the probability of detecting SNRFs. We estimated a mean intercept for the probability of detection (on the normal scale) to be approximately 0.1. This low detection probability indicates that a single camera deployment is likely to be insufficient to effectively survey for SNRFs and that additional cameras may be needed to appropriately determine their presence or absence. Contemporary surveys for SNRFs have often relied upon placing cameras in upper montane, subalpine, and alpine areas where SNRFs, and wildlife more generally, are thought to move (e.g., saddles, gaps between peaks) and set to run throughout the winter. These surveys are challenging and time-consuming to complete because survey locations are remote and difficult to access, cameras are often buried or damaged during periods of heavy snowfall, and it can be unsafe to revisit stations during winter. Thus, it would be valuable to maximize detection probabilities and optimize effort while minimizing risks to researchers conducting surveys in winter conditions.

We found that the probability of detecting SNRFs was highest at cameras deployed on patches of open ground, on gentle slopes, and with longer survey duration (Fig. 3). Although we were unable to explicitly test the effects of seasonality on detection, it is possible that the relationship between survey duration and detection probability is a function of differences in detection probabilities by season–cameras that were deployed for longer periods spanned multiple seasons. Determining the optimal season for surveys would be valuable. Interestingly, we found that SNRFs were less likely to be detected at baited cameras. We caution, however, that our data set was a compilation of many studies that used different protocols, so it is also possible that differences in baiting methodology confounded this covariate. Further investigation may clarify whether the presence of bait influences the probability of detecting foxes. Similar to Akins (2016), we also found that deploying cameras in patches of open ground and on gentle slopes increased detection probabilities of SNRFs. This supports current best practices for camera placement for SNRF surveys (Sierra Nevada Red Fox Conservation Advisory Team 2021). The above caveats notwithstanding, our findings suggest that cameras deployed on gentle slopes in open areas and set to run for as long as possible may maximize detection probabilities of SNRFs. Although scent lure did not have significant effect on detection probability here, it is possible that in some study areas lure may cause animals to linger longer or attract them into the field of view of a camera.

Studying rare and elusive wildlife can be challenging yet essential given the increasing rate of species–and geographically and genetically distinct populations within a species–extinctions and the fragility of montane habitats to climate change. Detections of such species are often sparse, and it can be difficult to discern true absence from a lack of detection (Lobo et al. 2010). For rare species, such as the SNRF, this uncertainty can limit inferences of density and distribution, and make it difficult to determine appropriate conservation and management actions. Thus, for rare species, it may be important to use all available data as possible to inform conservation and management decisions. We demonstrated that a broadscale study that employs contemporary quantitative approaches can maximize limited detection data for elusive species such as the Sierra Nevada Red Fox (V. v. necator).

Supplementary Data

Supplementary data are available at Journal of Mammalogy online.

Supplementary Data SD1.—Frequency histograms of minimum temperature, snow water equivalent, and normalized difference vegetation index (NDVI) at (A) locations of remote camera, opportunistic, and telemetry detections of red foxes and (B) locations of genetically verified detections of Sierra Nevada Red Fox (Vulpes vulpes necator; SNRFs), used as spatial covariates to model the mean predicted SNRF occupancy in California and Oregon. Solid black lines indicate mean values and dashed black lines indicate standard deviation values.

Supplementary Data SD2.—JAGS code to fit a Royle–Nichols integrated distribution and density model.

Supplementary Data SD3.—Mean predicted Sierra Nevada Red Fox (Vulpes vulpes necator) occupancy in California and Oregon. Lighter colors and no color indicate comparatively higher predicted occupancy values.

Supplementary Data SD4.—Standard deviation of predicted Sierra Nevada Red Fox (Vulpes vulpes necator) occupancy in California and Oregon. Lighter colors and no color indicate comparatively higher standard deviation values of predicted occupancy.

Supplementary Data SD5.—Predicted Sierra Nevada Red Fox (Vulpes vulpes necator) distribution in California and Oregon.

Acknowledgments

This project was funded by U.S. Fish and Wildlife Service, Yosemite Conservancy, and Yosemite National Park. We are grateful to the project partners and many field technicians who were instrumental in collecting or contributing data to this effort, including Chris Stermer, Brett Furnas (Furnas et al. 2017), Jay Power, Brent Barry, Katie Moriarty, Ellen Brandell, Todd Calfee, Lacey Greene, Ann Kleinfelter, Jackie Leary, Cody Massing, Elsbeth Otto, Julia Runcie, Liz Siemion, Amy Sturgill, Jon Weissman, Maria Immel, Michael McDonald, Stephanie Eyes, Breeanne Jackson, Laura Pilewski, Rob Pilewski, Toren Johnson, Brian Malloure, Ryan Kelly, Lacey Green, Catherine Gumtow-Farrier, Daniel Gumtow-Farrier, Andrew Irvin, Michael Magnuson, Jeff McFarland, Marie Tosa, Jamie Bowles, Cascadia Wild, Oregon Department of Fish and Wildlife, California Department of Fish and Wildlife, and the U.S. Forest Service Sierra Nevada Carnivore Monitoring Program. We also thank Stephanie Eyes, Julia Runcie, Kalysta Adkins, Thomas Jung, and two anonymous reviewers for helpful comments on an earlier draft of our manuscript.