The Humboldt marten (Martes caurina humboldtensis) has declined from over 95% of its historic range and currently occurs in just four extant population areas (EPAs). Prior to their listing under the Endangered Species Act, a conservation strategy was developed to identify key conservation needs for this species. This assessment identified an area near the California–Oregon (CA–OR) border as the second EPA in California, yet little was known about the overall distribution or habitat used by this population. This prompted our investigation to provide the first systematic survey of the CA–OR EPA and to assess habitat use using an occupancy modeling framework. Between 2017 and 2018 we surveyed 51 survey units in and adjacent to the EPA and detected martens at 20 units (39.2%). We found that occupancy by martens was most influenced by the amounts of low-elevation late-seral old-growth forest habitat, riparian habitat, and mid-seral forest habitat. The probability of occupancy by martens was greatest in low-elevation (< 800 m, Odds Ratio [OR] = 0.33, 95% Confidence Interval [CI] = 0.13–0.81) habitat and was positively associated with late-seral forest habitat at the 1,170 m home range scale (OR = 35.31, 95% CI = 1.30–958.07), riparian habitat at the 1,170 m home range scale (OR = 3.20, 95% CI = 1.01–10.1), and marginally by mid-seral forest habitat at the 50 m microhabitat scale (OR = 1.28, 95% CI = 0.95–1.73). Our findings identified habitat characteristics important for explaining the distribution of this understudied population, addressing two of the highest priority research needs identified in the Humboldt marten conservation strategy.
Key Points
Our surveys detected Humboldt martens in areas beyond the previously mapped California–Oregon extant population area, expanding the known distribution.
Martens were detected at 20 of 51 (39%) survey units in and adjacent to the California–Oregon extant population area, suggesting a patchy distribution.
Occupancy by martens was most influenced by low-elevation late-seral forest and riparian habitat (home range scale).
Introduction
Over the last few centuries, the global loss of biodiversity has occurred at an alarming rate (Estes et al. 2011, Segan et al. 2016), and this trend is particularly profound for rare species (Dirzo and Raven 2003). Rare species are inherently vulnerable to population declines due to their limited distributions and low abundances (Drever et al. 2012). Furthermore, the difficulties associated with studying elusive species can pose challenges in developing timely conservation initiatives (Martin et al. 2022). Understanding habitat use of at-risk species is an important first step in identifying key areas for management and recovery (Krausman 1999), yet lack of sufficient data is a common challenge in modeling habitat use for rare and elusive species (Hamilton et al. 2015, Todman et al. 2023).
The Humboldt marten (Martes caurina humboldtensis), also known as the coastal marten, is a subspecies of the Pacific marten (M. caurina) and is an example of a rare and elusive species for which knowledge of key population dynamics is lacking (Martin et al. 2022). The Humboldt marten is a medium-sized forest carnivore that historically occurred throughout the coastal forests of Oregon and northern California and has declined from > 95% of its historic range (Slauson et al. 2018, Moriarty et al. 2021). Signs of decline began to appear in the early 1900s due to unregulated and excessive trapping for their fur (Grinnell et al. 1937), while continued declines and lack of recovery following cessation of trapping has been attributed to extensive timber harvesting that followed throughout the late 1900s (USFWS 2015, Slauson et al. 2018). After 50 years without verifiable detections, the Humboldt marten was considered extirpated throughout its California range (Zielinski and Golightly 1996). However, in 1996 the subspecies was rediscovered in a remote portion of its historical range in northwestern California (Zielinski et al. 2001).
Contemporary surveys conducted throughout the historical range of the Humboldt marten in California and Oregon have identified four extant population areas (EPAs): two disjunct EPAs in Oregon along the central and southern Coast Range, and two disjunct EPAs in California, one in the northern Coast Range and the other farther inland near the California–Oregon border (CA–OR EPA; Slauson et al. 2018). Despite extensive survey efforts, there is still uncertainty about the exact distribution, population size, and habitat use of the few populations of Humboldt martens that remain (Moriarty et al. 2016, Slauson et al. 2019). In 2009, the northern coastal California EPA was estimated to contain fewer than 100 individuals (Slauson et al. 2009), and in 2018 the population size of the central coastal Oregon EPA was estimated at 71 individuals (95% Confidence Interval [CI] = 41–87; Linnell et al. 2018). Concerns over the persistence of this subspecies, known from only a few small and geographically isolated populations, prompted the listing of the Humboldt marten as Endangered under the California Endangered Species Act in 2018 (CDFW 2019) and Threatened under the federal Endangered Species Act in 2020 (83 FR 50574).
With Humboldt martens occupying < 5% of their historic range, it is critical to understand the habitat conditions important for supporting the few existing populations. Humboldt martens are considered habitat specialists and like other carnivores have large home ranges relative to their small body size (Lindstedt et al. 1986). Consistent with marten species across much of their North American range, Humboldt martens occur in structurally complex, late-seral and old-growth forests (Andruskiw et al. 2008, Kirk and Zielinski 2009, Thompson et al. 2012). This habitat type contains large trees and snags that are used for resting and denning, abundant prey resources, tree and shrub cover for protection from aerial predators, and downed woody debris near the forest floor that can improve hunting success (Andruskiw et al. 2008, Kirk and Zielinski 2009, Thompson et al. 2012).
Surveys of the northern coastal California EPA found that Humboldt martens were primarily associated with late-seral forest habitats, but they have also been detected in two low-productivity forest habitat types: shore pine (Pinus contorta) dominated coastal forest habitat found only on stabilized dunes (Linnell et al. 2018, Moriarty et al. 2019), and serpentine forest habitat found only on ultramafic soils (Slauson et al. 2019). The central coastal Oregon EPA persists entirely in young, coastal forest habitat (< 70 years old; Eriksson et al. 2019), and detections in serpentine forest habitat have occurred in both the southern coastal Oregon EPA and northern coastal California EPA (Moriarty et al. 2019, Slauson et al 2019). Collectively, these two low-productivity forest habitat types are endemic to their parent soil types and are limited to < 8% of the Humboldt marten's historic range (Slauson et al. 2019). Martens can persist in these two less productive forest habitat types, so long as key habitat types are available for supporting resting, denning, and prey resources (Slauson et al. 2007, Moriarty et al. 2016, Moriarty et al. 2021). However, martens do not occur in their structural analogs (i.e., forest habitat with small diameter or young trees) in the productive forest habitats that comprise the majority (> 90%) of their historical range where they have been largely extirpated (Slauson et al. 2018). Dense shrub cover, typically dominated by ericaceous species, is the most consistent habitat feature within these three distinct habitat types used by Humboldt martens (Slauson et al. 2018, Moriarty et al. 2019).
The occurrence of martens in these three distinct habitat types demonstrates the variation in habitat use among the EPAs (Slauson et al. 2007, Eriksson et al. 2019, Moriarty et al. 2021). This variation highlights the importance of using localized data to model habitat use that may be particular to each remnant population. With only a handful of verified detections in the CA–OR EPA, little is known about the habitat types that are most important for Humboldt martens in this area (Slauson et al. 2018). The first verified detection of a marten in the CA–OR EPA occurred in 2011, with subsequent surveys between 2012 and 2014 detecting martens at five additional locations (Slauson et al. 2018). No formal assessments of the distribution or habitat associations of martens in this EPA have been conducted to date, and these assessments have been identified as high priority information needs in the Humboldt marten conservation strategy (Slauson et al. 2018).
Our primary objective was to conduct the first systematic survey of the CA–OR EPA and provide a formal assessment of the habitat use and distribution of Humboldt martens in the least studied population. This population-level assessment provides new insight into the habitat types used by martens in the CA–OR EPA. Understanding habitat requirements for species of conservation concern is essential for developing effective management and conservation actions. Our study addresses one of the most important information needs identified in the Humboldt marten conservation strategy (Slauson et al. 2018).
Methods
Study Area
The CA–OR EPA is located primarily on federal lands managed by the Six Rivers and Siskiyou National Forests in northwestern California, just south of the Oregon border (-123°42′58″W, 41°53′41″N, Figure 1). The study area encompassed approximately 406 km2 and ranged from 27 to 48 km inland from the Pacific Ocean. The climate is characterized by warm, dry summers and cool, wet winters (3–30 °C, Jimerson 1989), with annual averages for precipitation of 237 cm and snowfall of 6 cm.
The study area was composed mainly of two habitat types used by Humboldt martens: serpentine forest habitats found on low-productivity ultramafic soils (17.0%), and productive forest habitats found on high-productivity soil types (83.0%; Soil Survey Staff 2022). The productive forest habitats were dominated by Douglas-fir (Pseudotsuga menziesii), incense-cedar (Calocedrus decurrens), Port Orford-cedar (Chamaecyparis lawsoniana), red fir (Abies magnifica), and white fir (A. concolor) plant associations (USFS 2018, CDFW 2021). Hardwoods, such as tanoak (Notholithocarpus densiflora), Pacific madrone (Arbutus menziesii), and canyon live oak (Quercus chrysolepsis) were also subdominant in the tree overstory. Ericaceous shrubs, such as evergreen huckleberry (Vaccinium ovatum) and salal (Gaultheria shallon), dominated the shrub layers in productive forest habitats. Serpentine forest habitats were dominated by Jeffrey pine (Pinus jeffreyi), knobcone pine (P. attenuata), and Douglas-fir plant associations. The dominant shrub species in serpentine habitats were huckleberry oak (Q. vacciniifolia), manzanita (Arctostaphylos spp.), bush tanoak (N. d. echinoides), and red huckleberry (V. parvifolium).
The study area was characterized by a mixture of forest seral stages. Using the Gradient Nearest Neighbor (GNN) vegetation structure map data (LEMMA 2017), the tree size class attribute data characterized seral stages based on quadratic mean diameter (QMD) and canopy cover, with early-seral stages represented by size class 0–3, mid-seral stages by size class 4, and late-seral stages by size class 5–6. Overall, early-seral stages (59.3%, size class 0–3) included 6.6% classified as unvegetated or in the shrub/seedling stage (size class 0–1, QMD 0–2.4 cm and canopy cover < 10.0%), 28.8% in the sapling/pole stage (size class 2, QMD 2.5–24.9 cm and canopy cover 10.0–24.9%), and 23.9% in the small tree stage (size class 3, QMD 25.0–37.4 cm and canopy cover 25.0–37.4%). Mid-seral forest habitat in the medium tree stage (size class 4, QMD 37.5–49.9 cm and canopy cover 37.5–49.9%) composed 17.3% of the study area, and late-seral forest habitat in the large and giant tree stages (size class 5–6, QMD ≥ 50.0 cm and canopy cover ≥ 50%) composed 23.5% of the study area (LEMMA 2017).
Detection Surveys
We used the Humboldt marten population monitoring protocol to survey for martens (Slauson and Moriarty 2014). This survey protocol is based on a 2 km systematic grid that covers the entire historical range. The 2 km distance between grid points is larger than the average radius of home ranges for male martens elsewhere in California (Moriarty et al. 2021), likely ensuring spatial independence from detecting the same individual at adjacent survey units. The survey period occurred during the latter half of the denning period (June–mid-August; Delheimer et al. 2021) to increase the likelihood of detecting resident adults rather than dispersing juveniles (Slauson and Moriarty 2014, Zielinski et al. 2015). At each point on the 2 km grid, we established a two-station survey unit: one placed on the grid point (station A) and the second placed 500 m away in a random direction (station B). In 2017, one remote camera station and one track-plate station were deployed within each survey unit. We randomly assigned either a track plate or remote camera to station A, and station B was assigned the alternative detection device. The Humboldt marten surveying protocol recommends the use of both remote cameras and track plates as both device types typically yield similar detection probabilities for martens (Gompper et al. 2006, Slauson and Moriarty 2014). However, we used remote cameras at all stations in 2018 due to the difficulties of deploying track plates in our study area and the similarities in detection events observed between device types within the survey units deployed in 2017.
Figure 1.
Study area and locations of survey units sampled in and around the California–Oregon Extant Population Area (CA–OR EPA) in northern California, USA, 2017–2018, depicting survey units with Humboldt marten detections (n = 20, closed circles) and non-detections (n = 31, open circles).
At stations with remote cameras, we used passive infrared-triggered cameras (Command Ops Pro; Browning Trail Cameras, Morgan, UT) programmed to take 8-shot photo bursts when triggered. Cameras were placed in metal security boxes to prevent damage from American black bears (Ursus americanus) and mounted to trees using lag bolts and straps. Bait was mounted < 0.6 m from the ground on a tree < 10 m away from the camera. Track plate stations consisting of a square Coroplast cubby closed on one end with hardware cloth were placed alongside a stabilizing structure (i.e., tree, stump, rocks). Aluminum plates were set inside the cubbies with toner and sticky contact paper to collect tracks when animals entered. Surrounding debris was placed along the sides and top to further stabilize the track plates, and bait was placed inside near the closed end. Each station was baited with two chicken drumsticks on the camera bait tree or in the back of the track plate and a sponge soaked in commercial trapping lure (Gusto; Minnesota Trapline Products, Pennock, MN) to attract martens (Baldwin and Bender 2008, Moriarty et al. 2018). The lure-soaked sponges were hung approximately 2 m above the ground in the tree or shrub nearest to the camera or track plate station. Once established, each station was deployed for a minimum of 21 days and revisited approximately every 3–5 days to replace bait, refresh the lure, and retrieve photographs from cameras or tracks from track plates. All survey methods were approved by the Cal Poly Humboldt (previously Humboldt State University) Institutional Animal Care and Use Committee (protocol 16/17.W.05-A).
Occupancy Modeling Approach
We used occupancy modeling to account for imperfect detection and to model the influences of habitat characteristics on the probability of occupancy by marten using our detection/non-detection data (MacKenzie et al. 2002). To create detection histories for each survey unit we first defined our survey occasion and then identified whether a marten was (1) or was not (0) detected during each occasion. Survey occasions were defined by each of the 3- to 5-day station check intervals, for a total of five survey occasions for each station. Detections were combined for both stations A and B to create a single detection history for each survey unit, including when track plates and cameras were used in 2017. To alleviate any concerns about different detection probabilities between cameras and track plates, there were no instances of a marten being detected only at the track plate and not at the camera station in 2017, therefore the resulting detection histories for survey units where both detection devices were used would remain unchanged if detection histories for track plates were constructed separately from stations with cameras. A survey unit was considered occupied if a marten was detected at either station using either method on at least one survey occasion.
We used a hierarchical modeling approach to develop and evaluate our candidate occupancy models by first modeling the detection process (p), and then using the top detection probability model in all occupancy models (Ψ). We used an information-theoretic approach to develop a candidate model set (Burnham and Anderson 2002) by first developing a set of a priori models representing alternative hypotheses of the most influential variables on the detection process and marten occurrence. Alternative a priori hypotheses were developed using variables known to influence habitat use in the three other Humboldt marten EPAs (Slauson et al. 2007, 2019; Moriarty et al. 2019, 2021), expert opinion, and hypotheses developed while conducting fieldwork in the study area ( Supplementary Table S1 (04_Gamblin_et_al._Supplementary_Material.pdf)).
Candidate Variable Selection
Twenty-three variables (3 detection, 20 occupancy) were considered for inclusion when developing candidate models ( Supplementary Table S1 (04_Gamblin_et_al._Supplementary_Material.pdf)). To evaluate the influence of survey-specific variables on detection probability, we included the variables survey month (June or July–early August) to account for temporal variation, and total survey duration (number of days) and station check interval length (number of days) to account for any effects of differences in overall survey duration. To account for potential heterogeneity in detection probability over the survey occasions, we considered both constant detection probability (p.) and occasion-specific (i.e., time-varying) detection probability (pt). For occasion-specific detection probability models, we incorporated the variable check interval length (check) to capture the realized differences in the number of days between when stations at each survey unit were checked.
We calculated a number of physical and biological variables to represent the habitat characteristics of the survey units ( Supplementary Table S1 (04_Gamblin_et_al._Supplementary_Material.pdf)). We used topographic and environmental variables from USGS, TIGER, and PRISM, including elevation, slope, road density, stream density, and precipitation ( Supplementary Table S1 (04_Gamblin_et_al._Supplementary_Material.pdf)). We used forest structure and composition variables from the GNN structure and composition datasets (LEMMA 2017): tree size classes (small, medium, and large), canopy cover, dominant/codominant conifer QMD, snag density, regionalized old-growth structure index (OGSI), late-seral old-growth forest (LSOG), mean forest ages, hard masting trees, coarse woody debris, and pine basal area ( Supplementary Table S1 (04_Gamblin_et_al._Supplementary_Material.pdf)). We generated shrub cover using modeled data published for available understory shrub species in the study area (Prevéy et al. 2022). We used the USDA Gridded Soil Survey Geographic Database (Soil Survey Staff 2022) and groups associated with gabbro and serpentinite soil types to identify serpentine habitat. All geographic information system (GIS) calculations were conducted in ArcMap 10.3 (ESRI 2015).
We evaluated each variable for inclusion in the candidate model set. Variables were excluded if there was incomplete GIS coverage in our study area, there was redundancy with other variables, or if they were inapplicable to our dataset. Excluded variables included slope, precipitation, small tree size classes, road density, snag density, coarse woody debris, forest age, serpentine, and pine basal area. Using this approach, we retained 14 (3 detection, 11 occupancy) variables ( Supplementary Table S2 (04_Gamblin_et_al._Supplementary_Material.pdf)). We evaluated correlations between the variables retained using the ‘corrplot’ package in RStudio (RStudio Team 2022). If a variable pair was highly correlated (correlation coefficient |r| ≥ 0.6), those variables were not included in the same model. We used the ‘car’ package in RStudio to test for collinearity among covariates within a single model by evaluating variance inflation factor (VIF) values (Zuur et al. 2010). Covariates with VIF ≥ 2 were removed from the model.
We evaluated the inclusion of sample units dominated by serpentine habitat prior to developing candidate models. We conducted an exploratory principal components analysis to compare detections in survey units located in low-productivity serpentine habitat (n = 4) to those located in high-productivity forest habitat (n = 16) ( Supplementary Table S3 (04_Gamblin_et_al._Supplementary_Material.pdf)). We found that 19 of the 20 candidate variables were significantly different between these unique habitat types ( Supplementary Table S4 (04_Gamblin_et_al._Supplementary_Material.pdf)). Due to the significant differences between these habitats and the small number of serpentine-dominated survey units, we excluded units dominated by serpentine habitat from the occupancy analysis (see Supplementary Material 3 (04_Gamblin_et_al._Supplementary_Material.pdf)). We reported the means and standard errors for the variables for survey units composed primarily of serpentine versus productive forest habitats separately and combined ( Supplementary Table S5 (04_Gamblin_et_al._Supplementary_Material.pdf)). We also reported the results of an exploratory occupancy analysis for all survey units combined ( Supplementary Table S6 (04_Gamblin_et_al._Supplementary_Material.pdf), Figure S1 (04_Gamblin_et_al._Supplementary_Material.pdf)).
Spatial Scale Optimization of Habitat Variables
Martens can exhibit habitat selection at multiple spatial scales (Slauson et al. 2007, Kirk and Zielinski 2009, Thompson et al. 2012). We used bivariate spatial scale optimization to identify the optimal spatial scale for each variable, which is a technique used to capture scale-dependent effects of habitat selection for martens (Shirk et al. 2014, Tweedy et al. 2019, Martin et al. 2021, Moriarty et al. 2021). We created six spatial scales represented by buffers around the 2 km grid point (station A) for each survey unit with radii of 50, 270, 500, 750, 1,170, and 3,000 m. The smallest spatial scale (50 m) represented fine-scale microhabitat types measured at the station level. The 270 m and 500 m scales represented within-home range (core area) scales (Tweedy et al. 2019, Slauson et al. 2019). The 750 m and 1,170 m scales represented the average female and male home range size, respectively (Moriarty et al. 2021). Our broadest spatial scale (3,000 m) incorporated landscape-level effects that may influence where martens position their home ranges within the surrounding area (Slauson et al. 2019). All occupancy models included only each variable's optimal spatial scale ( Supplementary Table S1 (04_Gamblin_et_al._Supplementary_Material.pdf)).
Candidate Models
We developed 11 candidate models for detection probability and 26 candidate models for occupancy to evaluate both additive and interactive effects of variables on the probability of occupancy ( Supplementary Table S2 (04_Gamblin_et_al._Supplementary_Material.pdf)). Due to the small sample size, we limited the total number of variables included in any occupancy model to ≤ 3 variables to reduce the risk of overfitting (Burnham and Anderson 2002) and maintain a ratio of ≥ 10 observations per estimated parameter. Models were fit using Program MARK (White 2001) and evaluated using Akaike's Information Criterion adjusted for small sample size (AICc). Models with ΔAICc < 2 units were considered to have substantial support (Burnham and Anderson 2002).
To interpret the relationship between each variable and marten occurrence or detection, we calculated odds ratios for variables present in models with substantial support. Odds ratios (hereafter reported as OR) were calculated by exponentiating the beta coefficients to estimate the influence of a one-unit shift on the odds of occurrence or detection. For variables where a one-unit shift was not biologically meaningful (i.e., 1 m elevation), we adjusted the Odds ratios (hereafter reported as OR) to reflect a scale appropriate to the range of the data by multiplying the beta coefficient by a more meaningful value (i.e., 100 m change in elevation) and exponentiating the adjusted beta coefficient. To evaluate the relative strength of each variable in the model set, we also calculated adjusted variable importance weights by taking the sum of AICc weights for models containing the variable and adjusting it relative to the number of models the variable appears in (Burnham and Anderson 2002). We created boxplots to visually examine the univariate relationship between the scale-optimized variables at detected and non-detected productive forest habitat survey units ( Supplementary Figure S2 (04_Gamblin_et_al._Supplementary_Material.pdf)).
Model Fit
Individual model fit was evaluated in program PRESENCE (MacKenzie and Hines 2006) using a parametric bootstrap goodness of fit test with 10,000 simulations. The goodness of fit test was used to generate an estimate of overdispersion, ĉ, to evaluate whether the top model adequately fit the data. The general approach for this method is to run the test on the global model. However, when the number of parameters in the global model is too large this results in reduced precision in the estimate of ĉ, which can make it difficult to detect lack-of-fit. We used the most parsimonious model to assess model fit, as that method is recommended when the global model has a large number of parameters (MacKenzie and Bailey 2004). The goodness of fit test generated an overdispersion estimate (ĉ) of 0.67 for the most parsimonious model, which is generally considered to reflect underdispersion (Cooch and White 2001). When ĉ < 1 it is recommended to set ĉ = 1 and proceed with model interpretation, and so we followed this guideline before interpreting parameter estimates (Cooch and White 2001).
Results
Occupancy Surveys
During June–August in 2017 and 2018 we surveyed 51 survey units (21 in 2017, 30 in 2018). Survey durations differed somewhat from the protocol, averaging 20 days (range = 14–28 days) and five survey occasions (range = 4–7 occasions). Stations with fewer than the recommended 21 days of survey effort occurred due to a nearby wildfire that required the removal of stations for safety concerns, or due to camera malfunctions. Survey durations were extended beyond 21 days at some stations to increase the chances of capturing hair samples for a complementary study.
Overall, martens were detected at 20 of 51 survey units (39.2% naïve occupancy; Figure 1). Martens were detected at a total of 24 of 102 stations (23.5%) across all two-station survey units, with only four survey units (20%) detecting martens at both (n = 8) stations and 16 survey units (80%) detecting martens at only one station (n = 16). At stations where martens were detected, detections occurred on an average of two survey occasions (range = 1–6 survey occasions). Mean latency to the first detection was six days (range = 1–13 days). Martens were detected at four of the nine survey units that were dominated by serpentine habitat (44.4% naïve occupancy). Martens were detected at 16 of 42 survey units dominated by productive forest habitat (38.1% naïve occupancy). Limited road access and remote terrain limited our ability to survey substantial portions of the eastern part of the CA–OR EPA; therefore, approximately half of the survey units we sampled occurred within the CA–OR EPA boundary and the rest were distributed along the western edge of the EPA boundary (Figure 1).
Table 1.
Beta estimates and odds ratios (OR) for the one top detection probability (p) and three top occupancy (Ψ) models for Humboldt martens monitored in northern California, USA, 2017–2018, along with associated standard error (SE) and 95% lower (LCI) and upper (UCI) confidence intervals. The optimal spatial scale (m) for each occupancy variable is included in the parameter name.
Occupancy Analysis
Of the 11 models for estimating detection probability, only one model had strong support (ΔAICc < 2; SupplementaryTable S2 (04_Gamblin_et_al._Supplementary_Material.pdf)). The top model for detection probability included survey month and total survey duration (Table 1), indicating these two variables accounted for sources of heterogeneity realized in the detection process. This model was used as the base detection probability model for all occupancy models.
The odds of detecting a marten during surveys conducted in July–early August were 281% greater than in surveys conducted in June (OR = 3.81, 95% CI = 1.31–11.10), after accounting for the effects of survey duration. The estimated detection probability for each survey occasion was 0.23 in June (95% CI = 0.12–0.38) and 0.53 in July–early August (95% CI = 0.34–0.71). For each additional survey day added to the mean survey duration of 20 days, the odds of detection increased by 14% (OR = 1.14, 95% CI = 1.02–1.28), after accounting for the effects of the month when the surveys were conducted (Table 1). The overall probability of detection using these month-specific detection probabilities for the typical 5-occasion (21 day) survey duration we used was 72.9% for June surveys and 97.7% for July–mid August surveys.
Of the 26 models evaluated for estimating the probability of occupancy by martens, three models showed substantial support (ΔAICc < 2; Supplementary Table S2 (04_Gamblin_et_al._Supplementary_Material.pdf)). The top-ranked model included the variables elevation (Elev) and mid-seral forest habitat (SC_Med). The second most competitive model included the variables riparian habitat (Stream) and late-seral forest habitat (LSOG), and the third most competitive model included an interaction between late-seral forest habitat and elevation (Table 1, Figure 2, Figure 3).
The amount of late-seral forest habitat and elevation had the greatest individual variable importance weights influencing the probability of occupancy by martens, followed by the amounts of mid-seral forest and riparian habitat (Table 2). The mean amount of late-seral forest habitat measured at the 1,170 m spatial scale was greater at survey units where martens were detected (mean = 46.0% [197.6 ha], SE = 1.8%, range = 35.0–58.8% [150.7–252.8 ha]) compared to units where they were not detected (mean = 35.8% [154.1 ha], SE = 2.5%, range = 16.3–66.0% [70.2–283.9 ha]; Table 2, Figure 2b). Using the beta estimates from the second-ranked model (Table 1), for every 5% (21.5 ha) increase in the amount of late-seral forest habitat at the 1,170 m scale, the odds of marten occurrence was 35.3 times greater (OR = 35.3, 95% CI = 1.3–958.0; Figure 3d). Martens were not detected in high-productivity survey units with < 150 ha (35% composition) of late-seral forest habitat at the 1,170 m spatial scale.
Martens were detected more often at survey units located at lower elevations (mean = 582 m, SE = 36.9 m, range = 362–858 m; survey units with no detection: mean = 964 m, SE = 67.3 m, range 458–1,655 m; Table 2). Using the beta coefficients from the best-supported model (model 1 in Table 1), a 100 m increase in elevation resulted in a 67.1% decrease in the odds of occurrence by martens (OR = 0.33, 95% CI = 0.13–0.81; Figure 3a). The influence of elevation and the amount of late-seral forest habitat on occupancy by martens also appeared to be interactive as one of the top models included their interaction (model 3 in Table 1). Most marten detections occurred at survey units with greater amounts of late-seral forest habitat located at the lowest elevations (Figure 2b). There was a 69.4% decrease in the odds of occurrence by marten for every 100 m increase in elevation (OR = 0.30, 95% CI = 0.21–0.40, Figure 3e) when using the beta coefficients from the interactive model (model 3 in Table 1) and the interacting variable at its mean value. Using the beta coefficients from the interactive model, for every 5% (21.5 ha) increase in the amount of late-seral forest habitat at the 1,170 m scale, the odds of marten occurrence increased by 198% (OR = 2.98, 95% CI = 2.88–3.08; Figure 3f).
Figure 2.
The habitat values associated with Humboldt marten detections (n = 16, closed diamonds) and non-detections (n = 26, open circles) in northern California, USA, 2017–2018, for the scale-optimized habitat variables present in the top three occupancy models: (a) elevation at the 50 m scale (Elevation_50) and mid-seral forest habitat at the 50 m scale (Size Class Medium_50), (b) elevation at the 50 m scale (Elevation_50) and late-seral forest habitat at the 1,170 m scale (LSOG_1170), and (c) riparian habitat at the 1,170 m scale (Stream_1170) and late-seral forest habitat at the 1,170 m scale (LSOG_1170).
The mean amount of mid-seral forest habitat measured at the 50 m spatial scale was slightly greater at survey units where martens were detected (mean = 17.0% [0.14 ha], SE = 5.7 %, range = 0.0–87.5% [0.0–0.69 ha]) than where they were not detected (mean = 13.2% [0.10 ha], SE = 3.4%, range = 0.0–62.5% [0.0–0.49 ha]; Table 2, Figure 2a). Using the beta coefficients from the top model (model 1 in Table 1), a 5% (0.04 ha) increase in mid-seral forest habitat at the 50 m spatial scale resulted in a 28.4% increase in the odds of marten occurrence (OR = 1.28, 95% CI = 0.95–1.73; Figure 3b).
Riparian habitat at the 1,170 m spatial scale was more abundant at survey units where martens were detected (mean = 1.55 km∙km-2, SE = 0.09 km∙km-2, range = 0.75–1.96 km∙km-2) than where they were not detected (mean = 1.17 km∙km-2, SE = 0.09 km∙km-2, range = 0.16–2.06 km∙km-2; Table 2, Figure 2c). Using the beta coefficients from model 2 (Table 1), a 100 m∙km-2 increase in the amount of riparian habitat resulted in a 220% increase in the odds of marten occurrence (OR = 3.20, 95% CI = 1.01–10.1, Figure 3c). No martens were detected at high-productivity survey units with < 0.75 km∙km-2 of riparian habitat at the 1,170 m spatial scale.
Figure 3.
Probability of occupancy (Ψ) by Humboldt martens in northern California, USA, 2017–2018, along with associated 95% confidence intervals for habitat variables in the top three occupancy models (AICc < 2) while holding the other variables present within the model at their average values. The top model depicts Ψ as a function of (a) elevation (Elev) and (b) the amount of size class medium trees (SC_Med) present at the 50 m scale (additive model, denoted by +). The second best-supported model depicts Ψ as a function of (c) riparian habitat (Stream) and (d) the amount of late-seral old-growth (LSOG) habitat present at the 1,170 m scale (additive model, denoted by +). The third best-supported model depicts Ψ as a function of (e) elevation (Elev) at the 50 m scale and (f) the amount of late-seral old-growth (LSOG) present at the 1,170 m scale (interactive model, denoted by *).
Table 2.
Adjusted variable importance weights for variables in the occupancy model set for Humboldt martens monitored in northern California, USA, 2017–2018. Variable weights were calculated as the sum of Akaike's Information Criterion weights (AICc) for models containing the variable relative to the number of models the variable appeared in. Variables are listed in decreasing order of importance. The average (x̄) values for each scale-optimized variable at detection and non-detection productive forest habitat survey units are reported along with associated standard error (SE).
Discussion
This study provides the first systematic survey of the CA–OR EPA and addresses two of the key information needs identified in the Humboldt marten conservation strategy: 1) to determine the distribution of martens in the CA–OR EPA, and 2) to identify habitat types that most influence the distribution of marten in this area. Martens were detected both in and adjacent to the previously mapped EPA boundary, suggesting the population was distributed more broadly than initially predicted and reported in the Humboldt marten conservation strategy (Slauson et al. 2019). We suspect that the distribution of this population may exist most significantly to the south, east, and southwest of the area we surveyed, based on the presence of similar habitat conditions to where most martens were detected during our efforts. Overall, occupancy of habitat by martens was most influenced by productive forest habitats located at lower elevations, with greater amounts of late-seral forest and riparian habitat at the home range scale (1,170 m), and to a lesser extent by greater amounts of mid-seral forest habitat at the microscale (50 m).
The amount of late-seral forest habitat at the home range scale and elevation collectively had the greatest influence on the occupancy of productive forest by Humboldt martens. The importance of late-seral forest for this population was consistent with habitat selection by martens in the larger California population of Humboldt martens (Slauson et al. 2007) and elsewhere for Pacific martens (Buskirk and Ruggiero 1994, Kirk and Zielinski 2009, Delheimer et al. 2019). Humboldt martens can occur at all elevations present within their historical range, from sea level to approximately 1,500 m (Slauson et al. 2018), yet martens in the CA–OR EPA primarily occupied low-elevation areas. However, the CA–OR EPA is located further from the coast than most of the northern coastal California EPA and the two Oregon EPAs, and it occurs in a more xeric climate than the other EPAs. The CA–OR EPA is one of the most inland locations where Humboldt martens have been found within their historic range, and these low-elevation (< 800 m) forest habitats may provide mesic microclimatic conditions that support more productive habitat for this more inland EPA.
The amount of mid-seral forest habitat and riparian habitat were present in the top two occupancy models, suggesting that occupancy of lower elevation sites by martens may be further influenced by the inclusion of more riparian habitat. Similar to the two Oregon EPAs (Eriksson et al. 2019, Moriarty et al. 2021), we found that Humboldt martens in the CA–OR EPA used areas associated with greater amounts of mid-seral forest habitat. However, the influence of mid-seral forest was only significant at the microscale (50 m) which represented < 1% of a typical marten home range. With such a small amount of habitat represented by the 50 m scale, this association may reflect microhabitat use rather than the influence of mid-seral forest on home range occupancy in the CA–OR EPA.
We used stream density as an indicator of the amount of riparian habitat, as riparian zones often support increased vegetation productivity and truffle production, leading to higher densities of prey (Doyle 1990, Waters et al. 2001). Riparian areas are important foraging areas for martens (Zielinski 2014), and these areas provide mesic microenvironments for thermoregulation that can be especially important during the warmest periods of the year. Riparian habitat was also positively associated with Humboldt marten occurrence at the core area scale (500 m radius) in broader habitat modeling efforts (Slauson et al. 2019), although its influence was much less than the amount of late-seral forest habitat in widespread productive forest habitats and the amount of serpentine habitat in the limited distribution of low-productivity habitats. The importance of riparian habitat may increase with distance from the coast or other dominant orographic features, such as major river valleys, as key habitat elements for Humboldt martens (e.g., dense and spatially extensive ericaceous shrub cover) are influenced by factors such as moisture and summer fog, which are less prevalent further inland.
Martens select resources at multiple spatial scales, therefore habitat models accounting for this scale-dependency can identify stronger relationships between resources and animal occurrences than single-scale models (Shirk et al. 2012). We tested six spatial scales (50–3,000 m) that have been used in prior analyses of habitat use by Humboldt martens (Slauson et al. 2019, Moriarty et al. 2021). However, the use of the smaller scales (50–270 m) departed from those theorized or demonstrated to influence home range scale habitat selection. Thompson et al.'s (2012) review of scale-specific habitat use by martens in North America found that habitat selection was strongest at the landscape scale across many studies, suggesting a robust connection between home range composition and individual fitness. Two of the most influential habitat variables in our analyses, late-seral forest and riparian habitat, were consistent with this home-range scale pattern of importance for key resources, while elevation and mid-seral forest habitat showed scale-specific optimization at the smallest microhabitat scale (50 m).
While elevation was statistically optimized at the 50 m scale, it was only marginally more significant than larger spatial scales. Moreover, nearly all topographic variables had the strongest statistical differences at the smallest spatial scales, raising further questions about their biological relevance. Finally, the interaction between elevation and late-seral forest habitat suggested that lower elevation late-seral forest at the home range scale was most influencing site occupancy by marten rather than the elevation of a small portion (< 1%; 50 m scale) of the home range.
The significance of mid-seral forest habitat at the 50 m scale may represent patterns of within-home range use, but because the scale represents < 1% of a marten home range, its biological relevance for home range selection and composition is questionable. While martens, like most animals, select resources at multiple spatial scales, they do not exhibit selection at all spatial scales at the same time (Mayor et al. 2009). Selection of resources to incorporate into a home range to provide for an animal's year-round resource needs may happen once in an individual's life, while selection of specific habitat types at the microscale may happen on a daily or hourly basis while they are foraging (Rettie and Messier 2000, Mayor et al. 2009). Therefore, it is critical to identify and constrain the selection of spatial scales for evaluation in multi-scale habitat modeling to those that the dataset is capable of addressing. In our study, we compared the portions of the study area occupied by martens to those not occupied by martens, essentially comparing where marten home ranges occurred versus where they did not. The spatial scales most relevant for modeling resource influence on home range occupancy should therefore be constrained to those representing significant portions of the study area (e.g., core areas, the entire home range, or the larger landscape area encompassing the home range). Although recent examples of modeling with spatial scale optimization for Humboldt martens include all six spatial scales (Slauson et al. 2019, Moriarty et al. 2021) and we sought to follow these methods, it may have been more appropriate to exclude the use of the smaller spatial scales (50–270 m), as these did not match the scales of habitat selection we were explicitly modeling. We recommend that the spatial scales used in multi-scale habitat analyses carefully evaluate scales of habitat selection that the study design and dataset can address and select only spatial scales for consideration that are relevant to the specific research objectives.
Although the majority of Humboldt marten detections in the CA–OR EPA occurred in high-productivity low-elevation forest habitats, four marten detections also occurred in low productivity serpentine forest habitats. This confirms that the two distinct habitat types present in the CA–OR EPA that are used by Humboldt martens elsewhere are also used by martens in this population. However, despite the large amount of serpentine habitat present in the broader region around the CA–OR EPA, previous research suggests the use of serpentine forest habitat may depend on its spatial juxtaposition to areas with large patches of late-seral productive forest (Slauson et al. 2018). The significant structural and compositional differences in the tree characteristics, primarily age and size classes/seral stages, between high-productivity and low-productivity forest habitat used by Humboldt martens have prompted researchers to assess characteristics for these distinct habitat types separately (Slauson et al. 2007). Our exploratory analysis of the differences in characteristics of the locations where martens were detected in each of these habitat types confirmed the stark differences between these habitat types ( Supplementary Table S4 (04_Gamblin_et_al._Supplementary_Material.pdf)). Our limited sample size for survey units dominated by serpentine habitat (n = 9) precluded our inclusion of these unique areas in this analysis. However, these data will be valuable when combined with larger samples for areas dominated by low-productivity serpentine habitats.
This study represents the first stage of determining the spatial extent of martens in this population and provides a timely assessment of habitat use in this area. We provide evidence that martens in the CA–OR EPA primarily occupy productive forest habitats located at low elevations and composed of large amounts of late-seral forest, mid-seral forest, and riparian habitat. In addition, some martens in the CA–OR EPA also occupy low-productivity forest composed of serpentine habitat. The CA–OR EPA has been affected by multiple recent wildfires since the completion of our surveys (USFS 2020), providing an opportunity to assess the short-term influence of mixed-severity wildfires on this population. Nearly all of the EPA burned between 2018 and 2023. Our surveys provide a pre-fire baseline of occupancy of habitat by martens in the CA–OR EPA that can be used to compare their distribution and post-fire habitat use, and to evaluate the effects of fire severity on post-fire occupancy patterns. Managers can help maintain and promote the expansion of Humboldt martens in and around the CA–OR EPA by using our results to prioritize the maintenance and restoration of habitat management areas that are composed of: 1) large patches of low-elevation (< 858 m) late-seral forest habitat (> 197.6 ha within 1,170 m radius areas), 2) large amounts of riparian habitat (> 1.55 km∙km-2 within 1,170 m radius areas), and 3) adjacent areas of low-productivity serpentine habitat.
Acknowledgements
We thank all those who contributed and supported us through this project. We thank T. Bean, B. Devlin, D. Barton, S. Hart, A. Benn, B. Carniello, K. Wright, and the many volunteers who contributed their time to the project. We also thank the U.S. Fish and Wildlife Service, the U.S. Forest Service, and the Cal Poly Humboldt (previously Humboldt State University) Sponsored Programs Foundation for their financial support.
© 2024 The Authors. This open access article is licensed under a Creative Commons Attribution 4.0 International License [ https://creativecommons.org/licenses/by/4.0/].
Conflict of Interest
The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.
Data Availability Statement
The datasets generated during the study are available from the corresponding author upon reasonable request.
Supplementary Materials
Supplementary materials are hosted online by BioOne at https://doi.org/10.3955/046.097.0404.s1.
Author Contributions
HELG: Conceptualization, data collection, writing—original draft, visualization, validation, formal analysis. KMS: Conceptualization, data collection, writing—review and editing, visualization, validation, formal analysis. MSG: Conceptualization, writing—review and editing, validation, supervision, funding acquisition.
