Using bear rub data and spatial capture-recapture models to estimate trend in a brown bear population

Field methods

We sampled the NCDE grizzly bear population in a 33,300-km2 area dominated by the rugged and remote terrain of the Rocky Mountains (Fig. 1a). We collected hair at bear rubs sites identified by looking for the presence of hair on trees and other objects while searching along trails, roads, and utility and fence lines. Hair deposition was a result of natural behaviour; no attractant was used to draw bears to survey routes or to encourage rubbing. Rubs were fitted with 3 to 4 40-cm lengths of 4-point barbed wire to facilitate hair deposition8. We monitored the rubs we located in 2004 and 2009–2012 (see Supplementary Information online for additional details).

Figure 1

(a) Northern Continental Divide Ecosystem study area in northwestern Montana, USA where growth rates of an expanding grizzly bear population were estimated during 2004–2012. State-space specified for the spatial capture-recapture models extended 45 km beyond the area where bear rubs were sampled. (b) Locations of naturally-occurring bear rub sites where the grizzly population was sampled by collecting hair for genetic determination of individuals 2004, 2009, 2010–2012. Maps were created from the USGS National Elevation Dataset58 using ArcMap 10.2.

We quantified rub sampling effort as the cumulative number of days between the first collection visit (after the initial visit to clear hair that may have been deposited the previous year) and last collection of the year for each rub sampled per year (Table 1). Sampling began in early June in all years, ended in early September in 2004, and most sampling ended in early October 2009–2012.

Table 1 Sampling effort and number of individual grizzly bears detected at natural bear rub sites monitored in the Northern Continental Divide Ecosystem in northwestern Montana, USA during 2004, 2009–2012.Genetic analysis

Hair samples were analyzed by a laboratory that specialized in genotyping low quantity and quality DNA (Wildlife Genetics International; http://www.wildlifegenetics.ca). We used the G10J marker to distinguish between grizzly bear and black bear (U. americanus) samples21. We added 6 additional nuclear microsatellite loci to determine individual identity of the grizzly bear hair samples (G1A, G10B, G1D, G10H, G10M, and G10P)22,23,24 and the amelogenin marker to identify sex25. We confirmed individuals identified using this 7-locus genotype by extending the genotypes to 9 additional loci (G10C, G10L, CXX110, CXX20, Mu50, Mu59, G10U, Mu23, and G10X22,23,24,26,27, and confirmed sex with a second proprietary test developed by Wildlife Genetics International. Sub-selection of samples for analysis is described in Supplementary Information online.

Spatial capture-recapture modeling

We used maximum-likelihood spatial capture-recapture (SCR) models28,29 to estimate local and overall density and annual population growth rates from spatial detection histories. SCR models specify a spatially explicit link between a summary of each individual’s location, henceforth their activity center, and locations where they may be detected, here, the sampled bear rubs. We treated the 5 years (2004, 2009, 2010, 2011, 2012) of detection data as independent “occasions” and ran separate models for males and females. For each year we denoted the coordinates of each rub by the vector xjt for years t = 2004 and 2009–2012 and for j = 1, 2, …, nt rub locations and the activity center for individual i as si, a 2-dimensional coordinate. SCR models parameterize the probability of detection, p(xjt,si), of an individual at each rub by a function of distance between the rub and the putative activity center. For each year and sex, we used the half-normal model30, where probability of detection at a location is assumed to decline with distance from the activity center following a half-normal distribution such that:

$$p({x}_{jt},{s}_{i})={p}_{0}exp(-dist{({x}_{jt},{s}_{i})}^{2}/(2{\sigma }^{2}))$$

(1)

The parameters we estimated were baseline detection probability at the activity center, p0, and the scale of movement parameter, σ, related to the extent of space used by individuals during the period of sampling. In addition, SCR models contain a parameter representing density, D, estimated from the observed encounter history data. Model parameters were estimated by maximum likelihood (MLE31) using the R package oSCR31,32. For all parameters of the model the standard error (SE) was from the Hessian of the likelihood function evaluated at the MLEs.

We developed our model set based on our a priori knowledge of grizzly bear ecology and theories about factors influencing bear behaviour. The hypotheses we considered modeled various effects on each of the 3 primary parameters; density (D), baseline detection probability (p0), and the scale of movement parameter (σ, Table S1 online). Effects were modeled on each parameter on a link-transformed scale, with an intercept and effects modeled additively. We modeled temporal variation in D by allowing full year-specificity, log(Dt) = βt, or a linear trend log(Dt) = β0+β1Year where Year was the centered year (2008 = 0). Density and σ were modeled on the log-scale and baseline detection probability was modeled on the logit scale. We included duration of sampling interval, linear effects of date, a trap-specific behavioural response, and landscape covariate effects in all models for p0 and evaluated hypotheses that included year-specific and quadratic effects of Julian date. We fit a behavioral model because rubbing functions as communication among individuals and thus prior marking at a rub may change the likelihood that site will be rubbed again by that or other individuals in the future. We calculated rub-specific landscape detection covariates to describe ruggedness, security level (protected land ownership), and the linear and quadratic effects of percent canopy cover within 250-m of a rub33 (see Supplementary Information online for additional details). We added these covariates to test the hypotheses that detection would be higher; 1) in rugged areas because bears might be more likely to walk on trails, 2) where there was higher security, and 3) where trees were of moderate to high density. We included hypotheses that σ was constant or year-specific. We used the sample size-adjusted Akaike’s Information Criterion (AICc) and AICc weights34 to evaluate relative support for each candidate model. Models with cumulative weights up to 1.0 were considered supported.

To address human-bear conflicts, wildlife managers sometimes move bears large distances resulting in unnaturally large distances between detections35. To assess the impact of this on our density estimates, we used the best model for each sex from the initial set of candidate models, but allowed σ and p0 to vary by translocated individuals. These models would thus estimate the proportion of translocated individuals, defined as Ψtrans = Pr(translocated = 1)36.

To calculate trend, we used the most supported models. For models varying linearly with year, the trend parameter was estimated directly within the model. For models with year-specific density parameters, we characterize trend by the annual geometric mean rate of population change37, λ, from 2004–2012 (8 annual intervals) according to:

$$\lambda =({D}_{{\rm{2012}}}/{{\rm{D}}}_{2004})\wedge (1/8)$$

(2)

The MLE of λ was obtained by plugging in the MLEs of year-specific densities. The standard error (SE) of estimated rates of population change was computed using a delta approximation38. The SEs were used to obtain 95% Wald type confidence interval. These growth rate estimates apply to the entire state space. We also compared our estimates to the trend in effective population size, Ne, a metric referring to the number of bears that are contributing reproductively to the population, under a set of genetic assumptions. For this we used analysis of covariance as described in Pierson et al.39.

In SCR models, MLEs are based on the marginal likelihood of the encounter histories computed by averaging over all possible locations of the activity centers associated with detected individuals. This spatial region is called the “state space”. The state-space should be chosen to be sufficiently large such that an individual having an activity center on the periphery of the state-space has negligible probability of capture (40 p. 132). For the half-normal detection model we used, we concluded that 3 σ was sufficiently close to negligible, yielding a probability of encounter < 0.01 based on estimates from our model. In our initial model runs, male σ was 15 km so we defined our state space by a regular grid of points (4 km grid spacing) extending 45 km beyond any rub location (Fig. 1). The number of state-space points varied with locations of rubs monitored each year and ranged from 3,992 (2004) to 4,422 (2012). SCR models regard the activity center, si, as a latent variable that is estimated along with other parameters. We estimated the latent activity centers using the estimated best unbiased predictor41 which is the conditional mean of Si given the data, evaluated at the MLEs. Predicted density surfaces were generated by aggregating the estimated posterior distribution of each Si. To examine changes in local trend over time, we estimated the realized local population size (aka realized density) for each state-space point (40 sec. 5.11), using the MLEs of the model parameters. The estimated realized population size for any state-space point u is defined as the estimated posterior mean of the number of activity centers associated with point u. This is computed by adding (a) the sum of the posterior probability of each observed individuals’ activity center for u to (b) the expected number of individuals for that state-space point that went undetected. The former quantity, (a), is denoted by \(Pr({{\bf{s}}}_{i}={\bf{u}}|{{\bf{y}}}_{i},\theta )\) (40 ch. 6) where \(\theta \) represents the collection of all model parameters for a given model, and we plug-in the MLEs for \(\theta \), \(\hat{\theta }\) for this evaluation. Component (b) is denoted by \({n}_{0}({\bf{u}}|\theta )\) which is also a function of all model parameters (40 ch. 6). Thus, the expression for the estimated realized population size for any point u of the discrete state-space is given by:

$$\hat{D}({\bf{u}})=\{\mathop{\sum }\limits_{i=1}^{n}Pr({{\bf{s}}}_{i}={\bf{u}}|{{\bf{y}}}_{i},\hat{\theta })\}+E({n}_{0}({\bf{u}}|\hat{\theta }))$$

(3)

This calculation can be done using the model parameters for each year to produce \({\hat{D}}_{t}({\bf{u}})\), from which estimated local growth rates can be computed for each unit of the state space,

$$\lambda ({\bf{u}})={(\frac{{\hat{D}}_{2012}({\bf{u}})}{{\hat{D}}_{2004}({\bf{u}})})}^{1/8}.$$

(4)

Local growth-rate estimates were limited to the area that was sampled in 2004 as this area was sampled in all years. We also calculated growth rates above and below Highway 2, a known barrier to gene flow8, by summing the estimated density of the state space points within each region in each year and calculating λ from those regional estimates. Variance estimates are not available for these sub-zones, however, because variance was calculated for the full state space and cannot be partitioned.

In one year, we detected one female bear 190 km to the south of other detections of this individual. Although grizzly bears are capable of traveling long distances, having a female range this far between detections, i.e. almost 3 times the distance of the second largest female movement detected (64 km), was highly unusual and constituted the lone extreme movement in our data set. While we did not identify any errors that could explain this large movement and it was not associated with a known translocation, in spatial modeling it is common to remove outlier locations when they are not representative of the population and leverage results42. We, therefore, fit all hypothesized models with and without the southern location for this individual, but consider primary results to be those associated with models with the location removed. Results from complete models that include this location are presented in the online Supplementary Information.