Predicted distribution of a rare and understudied forest carnivore: Humboldt marten (Martes caurina humboldtensis)

Study area

We collected data throughout coastal northern California and Oregon. The Humboldt marten is considered to occur in four Extant Population Areas (EPAs), which were created using minimum convex polygons around clusters of marten detections, but excluded clusters with smaller numbers of detections (<5) or detections >5 km from other detections (USFWS, 2019). As such, our surveys included both the recognized EPAs (Central Coastal Oregon, Southern Coastal Oregon, California-Oregon Border, and Northern Coastal California; Fig. 1) but also extended between these designated boundaries to include the historic range (USFWS, 2019).

Surveys in California occurred in both near-coastal and montane areas (Klamath Mountains, California Coast Range) that received substantial precipitation (100–300 cm annual precipitation) with cooler (7–10 °C) temperatures and drier summers dominated with fog and low cloud moisture (Rastogi et al., 2016). Forest types included a mix of coniferous and hardwood with a spatially-extensive shrub understory and dominant tree species included redwood (Sequoia sempervirens) along the coast and Douglas-fir (Pseudotsuga menziesii) in the mountains (Whittaker, 1960).

Surveys in Oregon similarly occurred in both near-coastal and montane areas (Oregon Coast Range) where dominant forest types included Sitka spruce (Picea sitchensis) and shore pine (Pinus contorta) along the coast and western hemlock (Tsuga heterophylla) slightly inland (Franklin & Dyrness, 1973). The Sitka spruce zone was characterized by a wet and moderately warm maritime climate with average annual temperatures of 10–11 °C, average annual precipitation of 200–300 cm, and frequent fog and cloud cover. The western hemlock zone, which was often co-dominated by Douglas-fir, was somewhat cooler (7–10 °C average annual temperature) and drier (150–300 cm annual precipitation) with fairly extensive summer fog and low cloud cover (Dye et al., 2020).

Common conifer species intermixed and included western hemlock, Port Orford cedar (Chamaecyparis lawsoniana), and western redcedar (Thuja plicata). Hardwood trees included tanoak (Notholithocarpus densiflora), giant chinquapin (Castanopsis chrysophylla), coastal live oak (Quercus agrifolia), canyon live oak (Q. chrysolepis), California bay (Umbellularia californica), red alder (Alnus rubra), bigleaf maple (Acer macrophyllum), and Pacific madrone (Arbutus menziesii). Dominant shrubs throughout the study area included salal (Gautheria shallon), evergreen huckleberry (Vaccinium ovatum), Pacific rhododendron (Rhododendron macrophyllum), and red huckleberry (V. parvifolium).

Marten locations

We used spatially-referenced Humboldt marten locations collected between 1996 and 2020 in California and Oregon. We excluded locations occurring in areas that were modified by fire or timber harvest after the location date and prior to 2016, the date represented by our vegetation data. If multiple locations occurred within a 500-m × 500-m grid cell, we spatially-thinned locations to randomly include one in each cell, attempting to achieve spatial independence for modeling (Kramer-Schadt et al., 2013). Priority for location retention from highest to lowest was: (1) rest and den locations from telemetry (Linnell et al., 2018; Delheimer et al., 2021); (2) locations from scat dog detection surveys (Moriarty et al., 2018; Moriarty et al., 2019); and (3) locations from baited camera and/or track plate surveys (Slauson, Baldwin & Zielinski, 2012; Barry, 2018; Gamblin, 2019; Moriarty et al., 2019). We used presence-only data because surveys that occurred prior to 2014 were often missing detection histories from non-detection (e.g., absence) locations.

For the data for which the authors were responsible, our protocols were reviewed and approved by the USDA Forest Service Research and Development Institutional Care and Use Committee (permits 2015-002, 2017-005) or Humboldt State University Institutional Care and Use Committee (permit 16/17.W.05-A). We obtained Scientific Take Permits for hair snares and samples collected through the Oregon Department of Fish and Wildlife (ODFW 119-15, 128-16, 033-16, 109-19, 107-20). Older verified survey data were provided by the US Fish and Wildlife Service with no additional information.

Modeling approach

Our modeling approach included Humboldt marten locations, biotic and abiotic predictor variables, and randomly generated pseudo-absence points. We used a minimum convex polygon (MCP) around Humboldt marten locations buffered by 10 km to define the modeling region (Fig. 1B). We chose a 10 km buffer because it approximated the upper quartile of daily marten movement (Moriarty et al., 2017). We projected our model to available vegetation data from Gradient Nearest Neighbor (GNN) data supplied by the Landscape Ecology, Modeling, Mapping and Analysis lab (Bell, Gregory & Davis, 2020; Bell et al., 2021), which included the coastal and Klamath level-3 eco-provinces (U.S. Environmental Protection Agency, 2013). We removed urban areas and water from the background data (Davis et al., 2016). We summarized the range, average, and standard deviation for each variable within the modeling region and study area (Table 1, Fig. 1).

Biotic variables

Biotic variables in our models included forest structure and composition, forest age, canopy cover, OGSI, percent pine, percent mast, and predicted shrub cover, as described below.

We used the 2016 version of GNN (Ohmann & Gregory, 2002) to incorporate forest structure variables including forest age, canopy percent cover, OGSI, and percent pine. Forest age was the basal area-weighted age based on field-recorded or modeled ages of dominant and codominant trees. Canopy percent cover was calculated using the Forest Vegetation Simulator (Crookston & Stage, 1999). Our OGSI index ranged from 0–100 was based from 4 elements: density of large diameter live trees per hectare, density of large diameter snags per hectare, percentage of downed wood greater than 25 cm in diameter, and an index of tree diameter diversity computed from tree densities in different diameter classes (Davis et al., 2015). For live trees and snags, “large diameter” was dependent on forest type and was defined for twelve vegetative zones, each zone with a unique minimum diameter threshold (i.e., ranging 50–100 cm for live trees, 50–75 cm for snags; Davis et al., 2015); see Item S1 for more information on integration of the OGSI variable into our model.

We created a variable called “percent pine”, which was the combined percentage of total basal area of shore pine, Jeffreyi pine (P. jeffreyi), and knobcone pine (P. attenuata) from GNN. This variable was included because martens have been detected in shore pine communities in the Oregon Central Coast population (Linnell et al., 2018; Eriksson et al., 2019), and in areas with serpentine soils characterized by sparse cover of Jeffreyi and knobcone pine, stunted tree growth, and dense shrub understories (Kruckeberg, 1986; Safford, Viers & Harrison, 2005; Harrison et al., 2006; Slauson et al., 2019). We visually inspected the congruence of the serpentine soil layer created by the US Fish and Wildlife Service (Schrott & Shinn, 2020) with our percent pine layer, confirming overlap between the two variables.

Humboldt martens have been associated with dense shrub cover throughout their range (Slauson, Zielinski & Hayes, 2007; Moriarty et al., 2019). Salal and evergreen huckleberry appear particularly important, as the berries of each occur in Humboldt marten diets and provide food for marten prey species (Eriksson et al., 2019; Manlick et al., 2019; Moriarty et al., 2019). We modeled probabilities of species occurrence of salal and evergreen huckleberry, creating the model for evergreen huckleberry following methods published for salal and other shrub species (Prevéy, Parker & Harrington, 2020; Prevéy et al., 2020). We related locations to contemporary (1981–2010) bioclimatic variables from the AdaptWest project (Wang et al., 2016) to depict the probability of species occurrence (1–100%). Humboldt marten diet is dominated by animals (e.g., passerines, ground squirrels) that feed on berries and mast and Humboldt martens also directly consume berries (Slauson & Zielinski, 2017; Eriksson et al., 2019; Manlick et al., 2019). The “mast” variable represented hardwood tree and shrub species that produce nuts, seeds, buds, or fruits eaten by wildlife and was estimated using the 2016 GNN layer as the percent of total basal area comprised of tanoak, giant chinquapin, coastal live oak, canyon live oak, and California bay.

Abiotic variables

Abiotic variables included temperature (°C), precipitation (cm), cloud cover (%), coastal proximity, percent slope, and topographic position index. We used 30-year normal PRISM variables of Average Annual Precipitation converted to cm and Maximum Temperature in August at an 800-m scale (1981-2010, PRISM Climate Group, Oregon State University, http://prism.oregonstate.edu, created 10/17/2019) as a proxy for maximum annual temperature. We explored annual data for temperature (2010–2018), but the available 4 km resolution produced artifacts in the model.

We created models with the variable Coastal Proximity, which uses PRISM data and combines coastal proximity and temperature advection influenced by terrain (Daly, Helmer & Quiñones, 2003) modified for the western United States (Daly et al., 2008). We derived percent slope and topographic position index from US Geological Survey digital elevation models. Topographic position index is an indicator of slope position and landform category; it is the difference between the elevation at a single cell and the average elevation of the user-defined radius around that cell (Jenness, 2006).

Scale optimization

Given that martens select habitat at multiple scales (e.g., broad-scale landscape features) and fine-scale features within home ranges (4th order selection; e.g., Minta, Kareiva & Curlee, 1999), we optimized the spatial scale of each variable included in the model. We smoothed variables using the extract function in package raster in R (Hijmans, 2020; R Core Team, 2020) with a radius of 50 m, 270 m, 742 m, and 1,170 m. Our smallest scale (50 m, 0.81 ha) provided local and fine-scale conditions. We assumed 270 m (20 ha) approximated the size of a Humboldt marten core area, similar to optimized scales of vegetation characteristics used in predicting conditions for marten rest structures elsewhere in California (Tweedy et al., 2019). The scale of 742 m (174 ha) represented an approximate female Humboldt marten home range size, calculated as the average of female home range estimates (173 ha) from two previous studies (Linnell et al., 2018; Data S1; PSW, 2019). Our broadest scale was based on the largest size of a Humboldt marten male home range (1,170 m, 428 ha, Data S1), assuming a male would overlap multiple females and could be interpreted as the smallest unit of population level selection (Linnell et al., 2018; PSW, 2019). We used individual univariate linear models (glm) for each spatial scale using our training location data and a random background sample of 9,600 points (25 times the location data) within the MCP at different locations than the Maxent generated pseudo-absence data (Data S2). Similar to prior examples (Wasserman et al., 2010; McGarigal et al., 2016; Zeller et al., 2017), we selected the scale for each variable that had the most extreme, and thus the most predictive, coefficient as well as the lowest Akaike’s Information Criterion (AIC) value. We also visually inspected the fit of each spatial scale using boxplots (Figs. S1–S3).

We provided boxplots to visually estimate whether our final variables were similar between all marten locations, thinned marten locations, available surveyed locations without detections (non-detection), and random locations (Fig. 2).

Predicted distribution

We used Maxent modeling software v3.4.1 (Phillips, Anderson & Schapire, 2006) to estimate the relative probability of Humboldt marten presence (Merow, Smith & Silander, 2013). Maxent uses a machine learning process to develop algorithms that relate environmental conditions at documented species’ presence locations to that of the surrounding background environment in which they occurred (Phillips & Dudík, 2008; Elith et al., 2011). We excluded variables with highly correlated predictors (—Pearson coefficient—>0.6), selecting the variable that was most interpretable for managers (Table S2). During this process, we considered the variance inflation (Table S3), which allows for evaluation of correlation and multicollinearity. Variance inflation factors equal to 1 are not correlated and factors greater than 5 are highly correlated as determined by (1/(1-R_i²)), where R_i² is squared multiple correlation of the variable i (Velleman & Welsch, 1981).

Within each model iteration, we selected the bootstrap option with 10 replicates, random seed, and 500 iterations. We trained our models using a random subset of 75% of presence locations and tested these using the remaining 25% with logistic output. We used the default of 10,000 pseudo-absence background samples. We varied the response functions to include linear, product, and quadratic features. We selected the “auto features” option for all runs, which allows Maxent to further limit the subset of response features from those selected by retaining only those with some effect.

Species distribution maps were produced from all models using the maximum training sensitivity plus specificity threshold, which minimizes both false negatives and false positives. We evaluated the AUC statistic to determine model accuracy and fit to the testing data (Fielding & Bell, 1997). The AUC statistic is a measure of the model’s predictive accuracy, producing an index value from 0.5 to 1, with values close to 0.5 indicating poor discrimination and a value of 1 indicating perfect predictions (Elith et al., 2006). We assessed variables using response curves, variable contributions, and jackknife tests. We used percent contribution and permutation importance to determine importance of input variables in the final model (e.g., Halvorsen 2013). Percent contribution can be more informative with uncorrelated variables (Halvorsen, 2013), while permutation importance provides better variable assessment when models and variables are correlated (Searcy & Shaffer (2016).

Because over-parameterized models tend to underestimate habitat availability when transferred to a new geography or time period, we used selection methods suggested by Warren & Seifert (2011). Maxent provides the option of reducing overfitting with a regularization multiplier that can be altered by the user to apply a penalty for each term included in the model (β regularization parameter) to prevent overcomplexity or overfitting (Merow, Smith & Silander, 2013; Morales, Fernández & Baca-González, 2017). A higher regularization multiplier will reduce the number of covariates in the model, becoming more lenient with an increased sample size (Merow, Smith & Silander, 2013). We did not include model replicates, an option in the interface, to output the required data (lambda file) and set output to logistic. We altered the Regularization Multiplier from 0.5 to 4 for each 0.5 increment (e.g., Radosavljevic & Anderson 2014).

We ranked candidate models using AIC corrected for small sample sizes (AIC_c; Burnham & Anderson 2002). We considered the model with the lowest AIC_c value to be our top model with those with ΔAIC_c<2 to be competitive models. For our top model, we generated predicted-to-expected (P/E) ratio curves for our model using only the testing data to evaluate its predictive performance, which was based on the shape of the curves, a continuous Boyce index (Boyce et al., 2002), and Spearman rank statistics. We used the predicted-to-expected curve to inform our suitability thresholds following Hirzel et al. (2006). We defined unsuitable in areas where the model performed equal to or poorer than random chance (P/E ≤ 1) with the lower 95% confidence interval of the P/E curve overlapping 0. For predicted suitable and highly suitable locations, we divided P/E and their respective 95% confidence values greater than 1, categorizing the lower half of data as suitable and the upper portion as predicted highly suitable.

Locations

We compiled 10,229 Humboldt marten locations collected during 1996–2020 (542 baited station, 263 detection dog team, 831 VHF telemetry, 8,537 GPS telemetry, 15 roadkill, and 41 others). Our GPS data represented locations taken every 2.5–5 min on 7 individuals within the Central Coast (Linnell et al., 2018), and we did not display those clustered data. After we spatially-thinned locations, 384 locations remained and were spread among Extant Population Areas: Central Coastal Oregon (n = 77 locations), Southern Coastal Oregon (n = 77 locations), California-Oregon Border (n = 33 locations), and Northern Coastal California (n = 192 locations) (Fig. 1). There were 5 locations that did not occur within boundaries of any EPA (USFWS, 2019). Location types included den or rest structure locations (18%), genetically verified scats or telemetry locations (32%), and baited camera or track plate locations (50%).

Thinned locations had similar medians and data distributions to the full location dataset, except for mast and precipitation where the medians were slightly lower for the thinned locations (Fig. 2). Non-detection locations had similar medians and data distributions to random locations, with the most notable difference between medians for salal (Table 1, Fig. 2). Differences between non-detection and random locations were likely due to clustered sampling efforts (Fig. 1B).

Distribution modeling

Our final model included 8 variables after excluding correlated variables (Tables S2, S3). Variables in our model were optimized at the home range spatial scale (1,170 m) except OGSI (50 m), but differences between scales were modest (Figs. S1–S3). Our top model had a Regularization Multiplier of 1.5. Predictor variables, in order of percent contribution, included a positive relationship with salal (23.3%), percent pine (22.5%), average annual precipitation (21.6%), canopy cover (18.7%), and mast (5.4%) followed by a negative relationship with average maximum August temperature (4.7%), percent slope (2.7%), and OGSI (1.1%, Table 2). Permutation importance was similar with the same top four variables highly contributing, but with a slightly modified order of percent pine (30.3%), average annual precipitation (25.3%), canopy cover (20.2%), and salal (15.5%; Table 2). The OGSI variable contributed least for both metrics.

We interpreted Maxent’s univariate response curves and provide the marginal plots as a supplemental figure (Fig. S4). Marten locations were correlated with both low and high amounts of canopy cover and percent slope (quadratic response, Fig. 3). Moderate amounts of canopy cover (e.g., 5–50%) appeared to be negatively correlated with marten locations. Predicted marten distribution was positively correlated with salal with some likelihood of a threshold at high values (Fig. 3), percent pine (Fig. 3), average annual precipitation (Fig. 3), and mast (Fig. 3). There was a negative correlation between marten locations and August temperature (Fig. 3) and a slightly negative to neutral relationship between marten locations and OGSI (Fig. 3).

The predicted versus expected curve of our final model delineated unsuitable areas as <14%, suitable areas as 15–30%, and predicted highly suitable at >30% predicted probability (Fig. 4, Data S3) with an AUC value on the test data at 92%. The model depicted southern Oregon and northern California as having the largest spatial extent for predicted marten distribution, including areas south of the current known distribution (Fig. 5, Data S3).