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

Trends in population abundance can be challenging to quantify during range expansion and contraction, when there is spatial variation in trend, or the conservation area is large. We used genetic detection data from natural bear rubbing sites and spatial capture-recapture (SCR) modeling to estimate local density and population growth rates in a grizzly bear population in northwestern Montana, USA. We visited bear rubs to collect hair in 2004, 2009—2012 (3,579—4,802 rubs) and detected 249—355 individual bears each year. We estimated the finite annual population rate of change 2004—2012 was 1.043 (95% CI = 1.017—1.069). Population density shifted from being concentrated in the north in 2004 to a more even distribution across the ecosystem by 2012. Our genetic detection sampling approach coupled with SCR modeling allowed us to estimate spatially variable growth rates of an expanding grizzly bear population and provided insight into how those patterns developed. The ability of SCR to utilize unstructured data and produce spatially explicit maps that indicate where population change is occurring promises to facilitate the monitoring of difficult-to-study species across large spatial areas.


Introduction
Although bear populations are challenging to monitor, bears possess an attribute that facilitates sampling them. Bears routinely rub on trees and other objects vigorously, often leaving hair behind. The most likely explanation for this ubiquitous activity is that it is a form of marking behaviour and functions as chemical signaling to other bears 1 . This tendency to rub results in the deposition of hair from individual bears at multiple sites. Because bears often use trails to travel, rub objects along trail networks can be located visually and therefore represent entirely passive detection; the sampling process has no effect on bear behaviour. While bears with a previous history of capture and handling have lower detection rates at baited hair corrals than bears that have not been trapped 2 , no such effect has been found at bear rubs 3 .
Grizzly bear and black bear (U. americanus) marking behaviour is well documented, if imperfectly understood 4,5,6,7 . Bears of all sex and age classes rub throughout the active season with variable frequency 8 . Evidence of rubbing is found throughout occupied habitat off-trail as well as along maintained and wildlife trails and by roads. Bears do not appear to discriminate much in selecting objects to mark. Bear hair has been found on a wide variety of natural (e.g. trees of all sizes and taxa, stumps, rocks) and man-made objects (e.g. sign and fence posts, utility poles, bridges, cabins) 8,9,10 . Female and subadult grizzlies rarely rub during the breeding season (May-June), however, their detection probability increases after low rates in spring and approaches male rates by late summer 3,8,11,12 . Many studies have demonstrated that genetic detection through hair collection is a highly effective way to sample bear populations; hair collection from a grid of baited hair corrals has been used extensively to study bear population abundance throughout the world (e.g. 13,14,15,16, ). Kendall et al. 3,8 added hair collection at natural bear rub sites as a secondary sampling approach to a primary method using an array of baited hair corrals. Because there was only partial overlap in individuals detected at hair corrals and rubs, the addition of this second sampling approach increased sample coverage and population abundance estimate precision while abundance estimates themselves were not significantly changed 17 . Results of subsequent studies that used concurrent sampling at hair corrals and bear rubs also found that the addition of rub site sampling was a cost-effective way to increase sample size and estimate precision 11,15,18 . While detection rates per sample site are higher at baited hair corrals than at unbaited rub sites 3,8,17,19,20 it is possible to sample a similar or larger portion of the population at rub sites than at baited corrals 11,18 . Sampling at rubs can also be more cost effective than at baited hair corrals. In most studies using both methods, more field crew time was put into sampling at baited corrals than at rubs even when rubs greatly outnumbered corrals. For example, the ratio of field crews assigned to hair corral vs. rub sampling was over 4:1, while the ratio of individuals detected at corrals and rubs ranged from 1:1 to 2.5:1 for the studies conducted by Kendall et al. 3,8,18 . However, the overall efficiency of rub tree sampling depends on monitoring higher numbers of rub trees to offset the lower per site sampling efficiency of rub trees 19,20 . Therefore, rub tree sampling is most viable in areas with greater access, such as maintained trails.

Study area
Approximately 85% of the Northern Continental Divide Ecosystem (NCDE) is in public ownership, where grizzly bear conservation is a management objective. It The climate and vegetation of the NCDE are shaped by the uplifted fault-blocks which formed the Rocky Mountains and by Pleistocene era ice ages 21 . Higher elevations receive more precipitation than lower slopes and valleys with annual precipitation ranging from over 170 cm in the higher elevations on the west side of Glacier NP to less than 38 cm on the eastern border. The Pacific maritime-influenced climate west of the Continental Divide is moister and more temperate than the colder and drier continental climate on the eastern side. The rugged terrain creates an array of habitats from alpine tundra and subalpine forests to moist mountain valleys, foothill groves, and semiarid grasslands in the rain shadow of large mountain masses 22 . Vegetation is characterized by coniferous forest, shrub fields and alpine tundra in the mountains, mixed deciduous-coniferous woodlands and herbaceous meadows in the valleys, and aspen (Populus tremuloides) stands, prairie grasslands, and agricultural fields along the eastern boundary 3 .

Field methods
To identify bear rubs, we closely examined objects that included one or more of the following characteristics; smoothed and/or discoloured bark and branch stobs, bare ground at the base of trees/poles, entrenched bear trails leading to rubs from maintained trails, and bear bite or claw marks. We attempted to get broad and even geographic distribution of areas searched for rubs and did not prioritize areas thought to have higher densities of bears. In lower elevation areas where there were fewer trees for bears to rub on, utility poles and fence lines provided sampling opportunities. To prevent damage to pack stock and their cargo, we used doublestranded smooth wire on trees that appeared to have been bumped by horse packs. We separated the twisted ends of smooth wire slightly to make them more effective in snagging hair. We visited rubs at least twice annually. All hair deposited on the wire was collected during each visit. The hair found on rubs during the initial visit each year was not analyzed and did not contribute to detection histories because it could have been deposited the previous year. A flame was passed under each barb and wire after sample collection, and hair between wires was removed to prevent contamination when the next bear rubbed. All hairs from a 4point barb or wire end were considered one sample.

Genetic analysis
In 2004, we attempted genetic analysis on all hair samples collected. During 2009-2012, we subsampled hair for analysis to limit genotyping costs. When we collected more than one hair sample during a visit to a rub, we genotyped the two samples with the most and/or largest follicles to maximize the amount of DNA in the sample. If sample quality was equal, we selected non-adjacent samples to increase the probability of detecting more than one individual per site. If there were distinct differences in hair colour between samples, we selected the best 2 samples from each, if available.

Spatial capture-recapture modeling
We included several trap covariates to describe detection. We calculated security level based on the location of a trap, where National Park Service =10, Forest Service =7, State and private forest lands =3, and all other =1. We used the standard deviation of curvature within a 250-m buffer of the trap covariate based on Ironside et al. 23 , which found it to be the most robust measure of terrain ruggedness. We also evaluated the effects of percent canopy cover based on the cover layer from Landfire, using the most recent data for that year. We updated the 2004 layer to reflect reductions in canopy cover due to some logging and extensive fires between the 2000/2001 imagery and the 2004 dataset.

Rub sampling effort
Sampling primarily occurred in forested habitats, however, bears' predilection to rub on fence and utility poles ensured sampling opportunities in prairie and other treeless habitats and agricultural areas. In 2004, search effort to identify rubs was lower along the eastern edge of our study area because there were few trails or little bear sign to guide our efforts. Sampling in this area was primarily in riparian corridors. Anecdotal information also indicated that few bears ranged east of the Rocky Mountain foothills at that time. After 2004, we expanded the number of monitored rub objects along the eastern foothills and plains as bears began to reoccupy former range. Starting in 2009, additional rub objects on the eastern edge were, therefore, added to the monitoring network. Rubs continued to be identified and added as sampling sites during early 2010 and then they generally remained in our monitoring network for the duration through 2012. A small number (<1%) of trees blew down or could not be relocated during our study; however, other monitored rubs were typically < 1 km away and in some cases, nearby new rubs were added to the monitoring network to replace them. As a result, sampling distribution and effort was nearly identical 2010-2012.

Discussion
Our point estimate for the annual rate of change, 1.043, was fairly high for an interior grizzly/brown bear population. Stable growth rates (λ ~ 1.0) have been reported in the past for a number of interior populations in the Rocky Mountains of the United States and Canada (e.g. 24,25 ). More recently, however, several populations have grown at rates comparable to those we found in the NCDE, likely a reflection of demographic recovery in response to conservation measures 26,27,28 . A study examining demographic mechanisms underpinning genetic assimilation of remnant groups of grizzly bears in the NCDE found that the increases in genetic diversity they identified resulted primarily from immigration 29 . It is likely that a combination of growth from within the remnant population coupled with bears moving in from higher density areas and the mathematics of small numbers in a local area (a change from 1 to 2 bears is a 100% increase) was responsible for the highest rates of pixel-specific growth we found in the NCDE (Figure 3). Quantifying local change rates using alternate measures such as gain or loss in density per pixel may be useful in smoothing this effect.
Bear species throughout the world exhibit rubbing behaviour; e.g., grizzly 30 , European brown bear 31,32 , Japanese brown bear 33 , American black bears 34 , and Asiatic black bears 15 . A growing number of studies have taken advantage of this behaviour to study bear demography (e.g., 8,11,35 ), dispersal 29 , and behaviour 7,36 . Developing a better understanding of which bears rub, why they rub, and where and when they do it will make this sampling method more powerful. The data obtained from genetic detection can inform multiple aspects of population status and mechanisms. For example, studies of habitat effects on abundance 38 , connectivity and dispersal 38,39 , the role of social learning 36 , genetic assimilation in remnant subpopulations 29 , and the relationship between effective and demographic population size in continuously distributed populations 40 , were based on the genetic samples collected for this and earlier genetic detection studies in the NCDE 3,8 .
Efficient sampling methods increase sample coverage while decreasing monitoring costs. The number of individual bears we sampled through hair collection at rub sites was an order of magnitude larger and more evenly spatially distributed, e.g. in the expansion area, than would typically be possible through live capture methods. Thus, we expect that genetic monitoring will be more sensitive to spatial and temporal changes in population demography than approaches that involve fewer individuals. Our genetic detection network yielded large sample sizes and precise estimates of population growth rate with 5 years of rub detection data (2004,(2009)(2010)(2011)(2012). Other approaches to obtaining similar spatial data, such as a grid of bear hair corrals, cost millions of dollars to implement across a region of this size 8 ; repeated sampling to obtain trend estimates using these methods across areas as broad as the NCDE has to date proved unsustainable. It may be possible to monitor populations such as ours with adequate precision with less sampling effort and at lower cost than in this study. The use of multiple detection methods 17 , clustered and stratified sampling designs 19,35 , adding telemetry data, or integrated population models that combine SCR data with lower resolution information 41,42 could further reduce the sampling and, thus, expense of SCR-based studies in the future. ~ linear ~ year + dur + jul + jul2 + beh + cover + cover2 + curv + sec ~ 1 a Model notation: D density (bears/1,000 km 2 ) p0 baseline detection probability σ sigma; spatial scale parameter related to the amount of space used by each individual ~ "a function of" year parameter is year-specific linear density is a linear function of time dur duration of sampling interval; number of days since previous sampling visit jul Julian day; linear effect of season on detection probability jul2 Julian day squared; quadratic effect of season on detection probability beh behavioural response; a visit of some bear at some bear rub changes its subsequent probability of detection at that rub cover percent of mean canopy cover within a 250 m radius of the trap cover2 percent of mean canopy cover squared curv standard deviation of terrain curvature within a 250 m radius of the trap sec security level as determined by land ownership policies at the trap: 10=NPS, 7=FS, 3= state/other, public, 1=private) year-specific detection probability p.beh behavioral response; a visit of some bear at some bear rub changes its subsequent probability of detection at that rub p.dur duration of sampling interval; days since previous sampling visit p.jul Julian day; linear effect of season on detection probability p.jul2 Julian day squared; quadratic effect of season on detection probability p.cover proportion of canopy cover