Density and population viability of coastal marten: a rare and geographically isolated small carnivore

Study area

We surveyed the northernmost population of coastal martens along the central Oregon coast within the 125 km² Oregon Dunes National Recreation Area (hereafter, “Oregon Dunes”). Coastal forests within the Oregon Dunes consisted of a narrow north-south strip along the margin of the Pacific Ocean bounded by two large rivers to the north and south (i.e., Siuslaw and Coos Rivers), Highway 101 to the east, and bisected by the Umpqua River, which is 600-m wide at the confluence with the Pacific Ocean (Fig. 1). Much of the forested area was the result of recent expansion over the last 70 years coincident with stabilization of near-coast beaches by European beach grass (Amophila arenaria) into mounded fore dunes, which limited sand deposition and facilitated vegetation expansion into previously shifting open sand (Christy, Kagan & Wiedemann, 1998).

Coastal dune forests grew on nutrient poor sandy soils (Christy, Kagan & Wiedemann, 1998), and they were dominated by young (<70 years-old) shore-pine (Pinus contorta contorta) and Sitka spruce trees (Picea sitchensis). The sub-tree canopy was dense, extended to >2.5 m in height, and it was dominated by willow (Salix hookeri), Pacific waxmyrtle (Myrica californica), salal (Gaultheria shallon), and slough sedge (Carex obnupta) on seasonally flooded sites, and berry-producing ericaceous shrubs (e.g., evergreen huckleberry Vaccinium ovatum, salal) on seasonally dry sites (Christy, Kagan & Wiedemann, 1998). Coastal forests differed substantially from inland forests east of Highway 101 in vegetation age, structure, composition, and their vertebrate communities (Eriksson, 2016). Inland forests were a mix of young (i.e., 0–80 years) and mature (i.e., >80 years old) Douglas-fir (Pseudotsuga menziesii) and Sitka spruce forests. Mature forests on federal lands were primarily managed as Late Successional Reserves to protect habitat for northern spotted owls (Strix occidentalis) and marbled murrelets (Brachyramphus marmoratus, Davis et al., 2015). Forests in the Oregon Dunes supported a high diversity of vertebrates, including several predators and competitors of martens (e.g., gray foxes Urocyon cinereoargenteus, coyotes Canis latrans, cougars Puma concolor; Eriksson, 2016).

To distinguish vegetation cover from open sand, we used airborne light detection and ranging data collected at 1-m resolution. We defined vegetation cover as >40% cover of pixels >1 m in height within a 100-m circular radius moving window of each pixel. This process produced a raster layer that smoothed small gaps in vegetation cover (i.e., sand gaps <30 m) that martens could presumably move through, but that excluded broad expanses of open sand that we assumed represented non-habitat for martens, particularly because of the presence of predators (Moriarty et al., 2015). The northern (i.e., north of Umpqua River, Fig. 1B) and southern (i.e., south of Umpqua River, Fig. 1C) study areas were comprised of 36.9 km² and 25.6 km² of vegetation cover, respectively.

Minimum and maximum temperatures in July and January were 10.1 °C and 20.3 °C and 3.2 °C and 10.2 °C, respectively. Annual precipitation averaged 176 cm, and occurred primarily between November and March (Western Regional Climate Center 1971–2016). Elevation within the study area ranged from eight to 80 m.

Live-capture and home range size estimates

We live-trapped and radio-marked martens from October to December 2015 using traps spaced approximately 1 km apart with some additional clustering of traps at <1 km (Fig. 1) using methods described in Moriarty et al. (2017) and Mortenson & Moriarty (2015). We fit adult martens (i.e., >two years old) with a VHF (Advanced Telemetry Systems, Minnesota, USA; 29 g) or GPS/VHF collar (Quantum 4000 Micro-Mini GPS collars, Telemetry Solutions, CA, USA; 41–44 g; or G10 snap technology GPS, Advanced Telemetry Systems, Isanti, MN, USA; 27 g; Table S1). Each individual marten was marked with a unique pattern of reflective tape attached to the antenna of the radio collar, which we used to resight marked individuals using black-LED remote cameras (Bushnell Aggressor, model: 119776; Bushnell Corporation, Overland Park, KS, USA; Fig. 2). Most martens (80%) were captured only once prior to the spatial mark-resight survey and then re-captured after this survey was complete when we removed collars in January and February 2016. All capture and handling procedures were approved by the USDA Forest Service’s Institute for Animal Care and Use Committee (USFS 2015–002) under an Oregon Department of Fish and Wildlife Scientific Take Permit (ODFW 119–15).

We programmed GPS collars to collect locations separated by 5 min, and only included locations in our analyses with predicted errors <30 m and time periods where data were collected for >72 consecutive hrs (details in Moriarty et al., 2017). We located individuals with VHF-only collars at least twice per week. We only used VHF locations where the variance of x and y was <400 m determined with Location of a Signal (Ecological Software Solutions LLC). We estimated 99% local convex hull home ranges, discarding 1% of the furthest dispersed locations (Lyons, Turner & Getz, 2013), using the t-LoCoH package in R. Local convex hulls were constructed using 35 neighboring locations (k = 35), regardless of time between locations (s = 0; Lyons, Turner & Getz, 2013; R Core Team, 2017). These parameterizations best reflected marten space use in our study (Moriarty et al., 2017) by limiting the formation of multiple activity centers, and they provided a smoothed outer contour boundary.

Spatial mark-resight

We monitored 31 sites for 39 consecutive nights in December 2015 and early January 2016 along a linear transect (henceforth, “SMR transect”) that overlapped the area occupied by radio-tracked martens in the northern study area. Each site was distributed (mean ± 1 standard deviation) 311 ± 91 m apart within vegetation cover with a minimum goal of four sites accessible to each female (Sun, Fuller & Royle, 2014). Each site consisted of one remote camera placed 0.5 m high, 2–4 m from attractants (i.e., bait, olfactory lure) and a measuring strip, such that the camera field of view was centered on the attractants (Fig. 2). We placed bait, in the form of ∼250 g of chicken and ∼100 g of strawberry jam, at each site during setup and replaced it at each site on 3 visits. Visits occurred every (mean ± 1 standard deviation) 8.2 ± 1.5 days after setup. We reviewed photos for each site and occasion, which we defined as the 24-hour day, to identify marked or unmarked individuals. We had (mean ± 1 standard deviation [range]) 32.8 ± 67.0 [1–637] photos available per site and occasion to identify marked or unmarked individuals. We censored photos that we were unable to determine if the marten was marked or unmarked (<1% of all photos).

We estimated the density of martens from our photographic data using a generalized spatial mark-resight model (Whittington, Hebblewhite & Chandler, 2018; Code in Article S1) based on the methods of spatial capture-recapture (Efford, 2004; Royle & Young, 2008). Generalized spatial mark-resight models combine the latent processes that generate both the capture (i.e., marking) and resight data (i.e., cameras) to estimate the number of latent activity centers (s_i) within the study area (Whittington, Hebblewhite & Chandler, 2018). We defined our study area as a discrete state-space S of a 100-m grid within a 5-km buffer around camera stations, excluding cells in the Pacific Ocean.

We defined live-capture data as the binomially distributed random variable ycap_ij representing the number of times that marten i was captured in trap j as a function of the probability of capture (pcap_ij) and the number of nights that trap j was open (Ktrap_j): ${\displaystyle {ycap}_{ij}\sim Binom\left({Ktrap}_j,{pcap}_{ij}\right).}$ We hypothesized that the probability of capture would vary by sex and the distance between the location of trap j and the activity center of marten i: ${\displaystyle {pcap}_{ij}={p\mathit{0}cap}_i\times e^{\left(-d_{ij}^2∕2{\sigma }_k^2\right)}}$ where the average probability of capture p0cap_i was modeled as a function of the sex of each marten (logit(p0cap_i) =β₀ + β₁ × sex_i), a half-normal decay function where d_ij is the distance between the trap and the latent activity center of individual i (s_i), and the standard deviation of a bivariate normal distribution reflecting space-use varying by sex (σ_k). We parameterized σ separately for each sex (See Article S1).

We defined camera resighting data as the Bernoulli distributed random variable ycam_ijk representing whether or not the previously live-captured marten i was resighted at camera station j on occasion k as: ${\displaystyle {ycam}_{ijk}\sim Bern\left({pcam}_{ijk}\right)}$ where pcam_ijk is a function describing the average daily rate of detecting martens on camera. Similar to the capture data, we hypothesized that the average daily rate of detection would vary by sex and the distance between the station and their latent activity center (${pcam}_{ijk}={\mathit{p0cam}}_{ijk}\times e^{\left(-d_{ij}^2∕2{\sigma }_k^2\right)}$). We also hypothesized that the time since bait addition would influence the rate of detecting martens on camera, so we added a variable to test the effect of days since baiting at station j on occasion k(days_jk) = (logit(p0cam_ijk) = δ₀ + δ₁ × sex_i + δ₂ × days_jk).

We modeled activity center locations using a non-homogeneous Poisson point process in S to examine whether locations of marten activity centers in the Oregon Dunes were associated with percent vegetation cover. We calculated the percent vegetation cover in each grid cell g in S, and used an intensity function following Royle et al. (2014) to model the probability of an individual being in grid cell g(p_g) as: ${\displaystyle p^g={\mu }_g∕EN}$ where μ_g is a function of an intercept (α₀), the linear effect of vegetation cover (α₁), and the size of the grid cell (area_g; μ_g = area_g × e^{a₀+a₁×vegetation cover_g}), and is divided by the expected number of martens in the study area (EN). We incorporated telemetry data to increase the precision of our estimates for the movement parameters and the location of marten activity centers (Royle et al., 2014; Sollmann et al., 2013). Telemetry locations of martens were modeled as being generated from a bivariate normal movement model with a mean of s_i and a precision of $\frac{1}{{\sigma }_k^2}$ (See Article S1). Martens can travel to any point in their home range within one hr (Moriarty, Epps & Zielinski, 2016; Moriarty et al., 2017). Accordingly, we only used locations that were at least one hr apart to ensure independence of telemetry locations (Sollmann et al., 2013).

We followed Whittington, Hebblewhite & Chandler (2018) in their treatment of the sightings of unmarked individuals; the detections of unmarked individuals at camera station j on occasion k (nU_jk) was modeled as: ${\displaystyle n{U}_{jk}\sim Sum\left(y{u}_{1jk},y{u}_{2jk},y{u}_{3jk\dots }y{u}_{ijk}\right)}$ ${\displaystyle y{u}_{ijk}\sim Bern\left({pcam}_{ijk}\right)}$ where the number of sightings of unmarked individuals was modeled to be generated from a latent Bernoulli process of resight probability (yu_ijk ∼ Bern(pcam_ijk)), based on the same probabilities of resighting as defined previously (pcam_ijk). Unmarked individuals were seen infrequently on our cameras (n = 14 sightings throughout the duration of the study), and identified as being present at the camera for a single bout typically <19 min. Thus, it was highly unlikely that more than 1 unmarked individual was captured on our cameras per day. The code for our SMR model can be found in Article S1.

We fit our models using data augmentation (Royle & Dorazio, 2008; Royle & Young, 2008) and the Markov-Chain Monte Carlo (MCMC) methods of JAGS (Plummer, 2003) with the jagsUI package (Kellner, 2014) in R v. 3.4.3 (R Core Team, 2017). We used uninformative prior distributions for all parameters (See Article S1). We calculated estimates from 3,000 MCMC samples, taken from three chains run for 10,000 iterations, thinned by five, following a burn-in of 5,000. We assessed model convergence by examining trace plots and $\hat{R}$ values for parameter estimates (Gelman & Hill, 2007; Gelman et al., 2014). All $\hat{R}$ values were <1.1, indicating chain convergence. We estimated the density of martens in our northern study area by determining the number of martens with estimated activity centers located in vegetation cover in the state-space, excluding open sand.

To evaluate our assumption that home range sizes were similar, we compared home range sizes in the northern and southern study areas using a general linear model with two parameters: sex and study area. We interpreted test statistics from this model and lacking any significant differences, we assumed density could be extrapolated to estimate population sizes (Moriarty et al., 2017).

Population viability

We assessed the risk of extirpation for a marten subpopulation over the next 30 years in the context of threats from human activities (e.g., trapping, roadkill; Gerber, Buenau & Vanblaricom, 2004). We estimated the maximum intrinsic population growth rate using a modified Euler-Lotka equation proposed by Skalski, Millspaugh & Ryding (2008) ${\displaystyle e^{ra}-e^{-M}{\left(e^r\right)}^{a-1}-{ml}_a=0,}$ where r is the maximum intrinsic growth rate, a is the age at first birth, m is the fecundity constant (number of female offspring/female/year), e^−M is the probability of survival, and l_a is the probability of survival to maturity (survivorship). We obtained a range of parameter estimates associated with the maximum reproductive output of two closely related species of North American martens (Martes americana, Martes caurina) in wild populations from the literature (Table 1). Estimates of r are sensitive to uncertainty in annual survival, e^−M; we estimated maximum intrinsic growth rate assuming average, high, and very high survival rates (0.7, 0.8, 0.9 respectively, McCann, Zollner & Gilbert, 2010) to obtain three values of r = 0.143, 0.205, 0.268 (Table 2). We used the intermediate value of r = 0.205 in our population projections, but we also implemented a stochastic element with σ = 0.06 such that the low and high estimates of r would bracket one standard deviation from the mean.

We simulated the dynamics of a population beginning at carrying capacity using initial values of the population size (K) equal to 20, 30, and 40 to illustrate how estimates of extirpation risk depend on our uncertainty about the current population size, assuming that immigration between the northern and southern study areas was infrequent due to a large barrier (i.e., Umpqua River; see population estimates in Results). The density-dependent population dynamics are given by the discrete theta-logistic model with an annual mortality component: ${\displaystyle N_{t+1}=e^{r\left(1-{\left(\frac{N_t}{K}\right)}^{\theta }\right)+\epsilon }{-H}_t,}$ where $\epsilon \sim N\left(0,\sigma \right)$, with σ = 0.06 based on the variation in our best estimate of r (Table 2), and mortalities resulting from trapping and road-kills as H_t ∼ Pois(λ). The rate parameter of the Poisson distribution, λ, defines both the mean and variance of the annual mortality through road-kills or trapping (H_t), which takes values of λ = 1, 2, or 3 martens in our models (36 martens harvested 1969–1995, 0–4/year; Verts & Carraway, 1998). We assumed a small density-independent harvest to illustrate how extirpation risk can be influenced by relatively low levels of human-caused mortality. We conservatively assumed a standard logistic population growth (θ = 1), but we also assumed that density-dependent declines in per-capita growth occurred at higher population densities (θ = 2), which is expected for long-lived mammals (Boyce, 1992). We simulated 1,000 population trajectories for each of three initial conditions (K = 20, 30, and 40), three stochastic human-caused mortality rates (λ = 1, 2, and 3), and two values of the strength of density dependence using theta (θ = 1, 2). Finally, we report observed mortalities during our study period. Where appropriate, we report results as mean ± 1 standard deviation.

Density and population size

We incorporated 79.3 ± 59.2 telemetry locations per individual into our SMR models. No marked individuals of the same sex were observed visiting the same camera station, and 1.1 ($\overline{x}$, range: 0.7–1.9) km and 4.7, 2.1 (n = 2) km was the furthest distance between camera station detections for females and males, respectively (Fig. 4).

We estimated marten density as (mean ± 1 standard deviation) 1.13 ± 0.15 individuals/km² (95% CRI [0.81–1.39]), or 9.75 ± 1.32 individuals within the SMR area. Assuming density was constant within vegetation in the 62.5 km² Oregon Dunes, we estimated a median population size of 42 (CRI = 30–51) north of the Umpqua River and 29 (CRI = 21–36) south of it. Sex did not have a significant effect on the probability of live-capture (Table 3), but female martens had a higher resight probability than males (Table 3). Days since baiting had a significant effect on resight probability; martens were more likely to visit baited cameras closer to a baiting event (Table 3). Percent vegetation cover had a significant effect on the distribution of activity centers (Table 3).

Population viability

We estimated that two or more annual human-caused mortalities on martens (e.g., trapping and road-kills) would lead to a substantial risk of extirpation, particularly at smaller population sizes (Figs. 5 and 6) and for θ = 1 (Fig. 5) relative to θ = 2 (Fig. 6). The likelihood of extirpation when θ = 1 for a population of 30 individuals, which approximated the average of our estimates for each study area, was 32% and 99% with two and three annual mortalities, respectively. The probabilities decreased to 1% and 60% when θ = 2 with two and three annual mortalities, respectively. The probability of extirpation increased to 89–100% and 65–100% for a population of 20 individuals with two or three annual mortalities.