Integrating resource selection into spatial capture-recapture models for large carnivores

Wildlife managers need reliable methods to estimate large carnivore densities and population trends; yet large carnivores are elusive, difficult to detect, and occur at low densities making traditional approaches intractable. Recent advances in spatial capture-recapture (SCR) models have provided new approaches for monitoring trends in wildlife abundance and these methods are particularly applicable to large carnivores. We applied SCR models in a Bayesian framework to estimate mountain lion densities in the Bitterroot Mountains of west central Montana. We incorporate an existing resource selection function (RSF) as a density covariate to account for heterogeneity in habitat use across the study area and include data collected from harvested lions. We identify individuals through DNA samples collected by (1) biopsy darting mountain lions detected in systematic surveys of the study area, (2) opportunistically collecting hair and scat samples, and (3) sampling all harvested mountain lions. We included 80 DNA samples collected from 62 individuals in the analysis. Including information on predicted habitat use as a covariate on the distribution of activity centers reduced the median estimated density by 44%, the standard deviation by 7%, and the width of 95% credible intervals by 10% as compared to standard SCR models. Within the two management units of interest, we estimated a median mountain lion density of 4.5 mountain lions/100 km (95% CI1⁄4 2.9, 7.7) and 5.2 mountain lions/100 km (95% CI1⁄4 3.4, 9.1). Including harvested individuals (dead recovery) did not create a significant bias in the detection process by introducing individuals that could not be detected after removal. However, the dead recovery component of the model did have a substantial effect on results by increasing sample size. The ability to account for heterogeneity in habitat use provides a useful extension to SCR models, and will enhance the ability of wildlife managers to reliably and economically estimate density of wildlife populations, particularly large carnivores.


INTRODUCTION
Understanding patterns and drivers of abundance and density is a central tenant of ecology and wildlife management (Andrewartha and Birch 1954), and the estimation of these parameters has long challenged ecologists (Seber 1973, Jolly 1982. The ability to estimate the abundance and density of wildlife is particularly important in applied settings, where conservation or management actions are often assessed by measuring trends in these parameters. Capture-recapture (CR) methods for closed populations are the standard methodology for the estimation of animal abundance from fixed arrays of traps or other sampling devices (Borchers and Buckland 2002). However, some species such as large carnivores introduce additional complexities for the estimation of density because these species are wide ranging, occur at low densities, and are difficult to detect. These species commonly violate assumptions about geographic closure in the CR framework and make estimation of the effective sampling area difficult (Royle et al. 2013a). Additionally, individual heterogeneity in recapture probability (detection probability) over fixed arrays due to behavioral differences or differences in space use challenges conventional CR methods. Recent methodological advances in spatial capture-recapture (SCR) methods address these shortcomings of conventional CR methods by incorporating the spatial organization of individuals through the estimation of trapspecific capture probabilities (Efford 2004, Efford et al. 2009, Gardner et al. 2010, Royle et al. 2013a. SCR methods have recently been extended to accommodate unstructured spatial sampling (Thompson et al. 2012) where effort is variable across the study area and applied to the estimation of mountain lion density (Russell et al. 2012). Mountain lions (Puma concolor) in North America have slowly increased in number and expanded their range over the last several decades (Hornocker and Negri 2009). Wildlife managers need reliable methods to monitor mountain lion population trends to manage harvest, minimize conflicts with humans and balance mountain lion densities with ungulate management objectives. Traditional approaches to estimate mountain lion abundance have focused on marking and counting individual lions, a method that is labor-intensive and expensive, and often minimum counts of individuals are treated as a true census for management purposes (Robinson et al. 2015). The resources required by these traditional methods have limited the spatial scope and utility of the resulting estimates for population management (Stoner et al. 2006, Quigley andHornocker 2010), and minimum counts of known individuals underestimate true population sizes. Therefore, the SCR approach provides managers with an effective and economical alternative method to estimate mountain lion density. Royle et al. (2013b) provided an approach to incorporating covariates on density including information about habitat selection into SCR in the Bayesian framework. Mountain lions have been shown to have strong selective preferences for areas that offer cover, forest edges, and moderate slopes (Newby 2011). The assumption that all habitat within the study area will be used equally is unlikely to be true for mountain lions and most large carnivores. Therefore, we extend the Russell et al. (2012) model for estimating mountain lion density to include the integration of a previously existing resource selection model as a covariate to account for this differential habitat selection. For many studies, pre-existing information about habitat use could inform the distribution of activity centers without requiring ancillary data (e.g., GPS or telemetry data), or use of the methodologically intensive joint estimation framework. Here, we evaluate the effects of introducing an existing RSF into Bayesian SCR models by comparing models results with and without the RSF. Additionally, we incorporate the ability to include information about harvested animals (dead recoveries) into Bayesian SCR models, increasing the utility of the modelling approach for species subject to harvest.

Study area
The 2,625 km 2 study area was located in the southern Bitterroot watershed in western Montana, USA, primarily within Ravalli County and spanning portions of two mountain lion management units (LMU 250 and LMU 270; Fig. 1

Data collection
We overlaid a 5 3 5 km grid across the study area and assigned each cell a grid identification number. We randomly generated a list of grid cells and started search effort each day in the randomly assigned grid cell. We stratified sampling in this manner to ensure sampling was allocated across both the high and low quality habitat. Mountain lion hair, scat, and muscle samples were collected by trackers and houndsmen for genetic analysis to identify individuals. When a fresh mountain lion track was located, the houndsmen would release trained hounds to locate and tree the mountain lion. Tracks were backtracked and inspected to determine if the mountain lion was independent or associated with a family group, and group size was recorded. Sex was determined based on track v www.esajournals.org size and physical features of treed mountain lions. Sex was assigned by genetic analysis when sex could not be determined in the field or field staff was uncertain of sex. Muscle samples were collected from treed animals using biopsy darts fired from a CO 2 -powered rifle (Palmer Cap-Chur brand model 1200c). When older mountain lion tracks were located, a tracker or houndsmen would backtrack the tracks and collect any hair or scat samples along the tracks. All field crews used a Global Positioning System to record the length (in km) and location of their search effort. Harvest and management removals occurred during the sampling period and we used samples collected from harvested animals within the study area. In Montana, the hide and skull of all mountain lions publically harvested must be presented to Montana Department of Fish, Wildlife and Parks. During the mandatory check, officers collected a muscle sample from each harvested animal. We also collected harvest samples from all adjacent hunting districts to determine if animals marked within the study area may have moved out of the study area. Adjacent districts had similar harvest management regulations as the study area.
To estimate the density of independent mountain lions in the study area, we censored the dataset to include only samples from independent animals or the adult female of a family group. This eliminated multiple samples from within family groups, and eliminated all groups where only a subadult animal within the group was sampled. The average age that dependent offspring disperse and become independent of their mother is approximately 15 months of age (Sweanor et al. 2000, Robinson andDeSimone 2011), therefore our density estimates include all animals .15 months of age. We estimated the number of independent mountain lions rather than density of all mountain lions because harvest management quotas are based on the number of independent, subadult or adult-aged animals.
We performed genetic analysis of hair, scat, and muscle samples to identify the sex and individual identity of sampled mountain lions following methods described in Russell et al. (2012). We genotyped tissue samples using 20 variable microsatellite loci used previously in mountain lions (see Appendix A for details).

Spatial capture-recapture modelling
We followed the hierarchical model formulation described by Royle et al. (2009) applied to genetic capture-recapture data from unstructured spatial sampling. For our purposes we decomposed the DNA observations of individuals at particular sites during a sampling period into two components: a spatial point process that describes the distribution of animals in space and an encounter process that describes the capture of individuals in grid cells given their activity center. The spatial point process model allows for the spatial information from the locations of individual captures to be incorporated into the estimate of density or abundance. Within the SCR framework we assume that a population consists of n individuals, and that each individual, i ¼ 1, 2, . . . , n, in the population has a fixed activity center representing the center of the area occupied by the individual during the study period. Each individual moves about these activity centers, which are allowed to overlap, according to some distribution defined in the observation model. In addition, we assume that the farther an observer is from the animal's activity center the less likely the animal is to be detected.
For the encounter model, we constructed individual, cell-specific encounter histories for each time period of the study, y i,j,k for individual i; cells j ¼ 1, 2, . . . , J; and sample periods k ¼ 1, 2, . . . , K. In this study individual encounters did not arise from discrete trap locations (e.g., camera traps), but instead encounters could occur at any spatial location searched by trackers and houndsmen. To accommodate this unstructured spatial sampling, we used the center of sampling grid cells as conceptual traps following Russell et al. (2012). Previous simulations demonstrated little effect of grid cell size on lion density estimates using this design (Russell et al. 2012), and, more generally in SCR models (Sollmann et al. 2012). Because animals were also harvested and removed from the population during the sampling period, we defined 4 sampling periods so that harvested individuals' detection probabilities could be adjusted in later sampling periods (see below for dead recoveries). We selected 4 sampling periods (December, January, February and March-April) that roughly corresponded to monthly sampling intervals in v www.esajournals.org efforts to adjust detection probabilities of harvested animals and still allow for adequate sampling effort and detections per sampling interval. We assumed that individuals could have been encountered in all traps for each sampling period through search effort or harvest. We developed an encounter history for each animal during each sampling period where y i,j k ; Binomial (K, p i,j ) (Eq. 5.2.1, Royle et al. 2013a) and the value of y i,j,k ¼ 0 if the animal i was not encountered in the grid cell j on sampling occasion k and y i,j,k ¼ 1 if the animal i was encountered 1 or more times in the grid cell j on sampling occasion k. If an animal i was harvested during a previous sampling period, we set the encounter history to y i,j,k ¼ 0 for all grid cells for all subsequent sampling periods, and removed the individual from likelihoods within the Markov chain Monte Carlo estimation of parameters (i.e., the detection probability for that individual after known harvest did not influence the estimate of detection probability for the remaining live individuals).
Following Gardner et al. (2010) and Russell et al. (2012), we assume that cell-specific encounter probabilities for each i individual and j cell ( p ij ) are related to the Euclidean distance between an individual's activity center (s i ) and cell j, r is a scale parameter depending on the measurement units, and h is a shape parameter of the activity distribution (Gardner et al. 2010) Royle et al. 2013a). When h ¼ 1, this model describes a bivariate-normal distribution of animal activity, and a value of h ¼ 0.5 corresponds to an exponential activity distribution. We used a prior distribution for h ; U(0.5, 1). This model of detection probability assumes a baseline encounter rate a 0 which is the probability of encounter at the animal's center of activity (Royle et al. 2013a). We considered the following covariates as potentially modifying detection probability: (1) sex, and (2) log-transformed search effort within a given sampling period (Thompson et al. 2012). We did not record hunter search effort, so we assumed that hunter search effort was equal to the mean of our search effort during each sample period because both hunters and researchers utilized the same road network to search for tracks. This allowed us to estimate the effects of search effort, given a small number of samples were collected via harvest in grid cells with no recorded search effort. We further considered the interactive effects of distance and sex on detection probability, wherein capture probabilities differ by sex because females are expected to have smaller activity distributions than males (Hornocker and Negri 2009). In practice, the statespace is chosen to encompass the movements of all animals within the study area but the statespace does not determine the extent of animal movements as part of the density estimate, as one-half mean maximum distance moved (MMDM; Royle et al. 2013a) or similar measures do in traditional studies. The statespace is chosen to be larger than the trapping grid such that the model accounts for individuals with activity centers outside of the trapping grid, but whose activity range extends into the trapping grid. The exact size of the trapping grid buffer used to construct the statespace does not strongly influence density estimates under the assumption of a uniform distribution of activity centers (Russell et al. 2012). When spatial covariates are included, however, the size of the buffer may influence density, as different values of spatial covariates may be included within the buffer. We choose a buffer size of 10 km around our study grid and identified potential activity centers within this area every 2 km. Because our study grid was not a square, we applied the 10 km from the most extreme edges of our grid in each cardinal direction, such that the buffer exceeded 10 km in some areas. In one area the buffer zone was reduced to ,10 km to exclude areas that did not include wintering lion habitat. This resulted in a 5912 km 2 statespace.
Currently, methods exist in ''secr'' (Efford 2014) to incorporate dead recoveries into spatial capture-recapture models in the likelihood framework. We developed a method to incorporate information about harvested lions in a Bayesian framework by including a matrix, indicating whether the animal was potentially alive in each sampling period and therefore available to be captured or if the animal had v www.esajournals.org been harvested during a previous session. This matrix was then used to censor when animals were not available to be captured from the likelihood, thus preventing a bias in detection probabilities from harvested individuals not available for detection post-harvest.
We included information about winter mountain lion habitat quality by incorporating predictions from an existing resource selection function (RSF) as a density covariate. We incorporated values of the RSF into our estimate of density using the following equation: . We used the logit of continuous predictions rescaled from 0 to 1 from this RSF as a covariate in the SCR model. We evaluated the effects of including harvested animals in the analyses in two ways. First, to evaluate potential effects of bias in the detection process by introducing individuals that could not be detected after removal, we compared density and precision of estimates from a model that masked harvested animals from sampling periods post-harvest and a model that treated harvested animals as live throughout the sampling period. Second, to evaluate the effects of the harvest samples on the density estimates, we compared density estimates from a model that included only the live animal samples with a model that included both the live animal and harvest samples. In both these evaluations, we used the best model identified through our model selection process for comparisons.

Bayesian analysis by MCMC
We used Bayesian MCMC methods in the SCRbayes package for the R programming language to estimate the posterior distribution of the model parameters (https://sites.google.com/site/ spatialcapturerecapture/scrbayes-r-package; see Supplement). This approach uses data augmentation to add a sufficiently large number of all-zero (unencountered) capture histories to create a dataset of M individuals. We determined the data augmentation to be large enough when final posterior estimates of population size were not limited by the number of augmented unencountered individuals (see Royle et al. 2013a). In essence, we assume a uniform prior distribution on N, population size, from 0 to M, where M includes the unencountered animals. Here, we augmented with 1000 all-zero encounter histories. Following Russell et al. (2012), models were run for 30,000 iterations with the first 10,000 iterations discarded as burn-in, leaving 20,000 samples from the posterior distribution. Starting values for parameters were: r ¼ 1, h ¼ 0.75, ln(a 0 ) ¼ 0, b ¼ 0, w ¼ 0, and w sex ¼ proportion of individuals sampled that were male (0.40 for our sample). We used improper priors (À', ') for a 0 and all b parameters, (0, ') for r, (0.5, 1) for h, and (0, 1) for w and w sex . We assessed model convergence by examining posterior parameter-wise trace plots and histograms.
To estimate lion abundance (N) within the statespace and within the two management units of interest, we counted the number of activity centers within the statespace and within the management units. Since this estimate of abundance is linked explicitly to an area (either the statespace or a management unit), we calculate density (D) by finding the quotient of the N activity centers and the area of the statespace or v www.esajournals.org management unit.

Bayesian model selection and model goodness of fit
We evaluated 16 potential models (Appendix C) that we fit to the data in a Bayesian framework. We conducted a form of model selection by examining the posterior significance of the parameters in each model weighted by our prior knowledge of mountain lion biology. For the sex and effort effects on the detection probability, as well as RSF covariate on activity center distribution, 95% credible intervals that excluded zero would provide evidence of significance. Similarly, for the sex effects on the scale parameter of the activity distribution, we examined whether male and female values had significantly different 95% credible intervals. We did not, however, immediately exclude models that included apparently non-significant parameters, as many of these parameters derive support from previous ecological studies (i.e., differences in male and female activity range sizes). We evaluated the effects of introducing the RSF covariate by comparing models with without the RSF as a covariate on the distributions of activity centers. We compared the density and precision of estimates of the best model with and without the RSF covariate.
We evaluated model goodness of fit for the top models following methods described by Russell et al. (2012) including a standard Bayesian Pvalue approach (Gelman andShalizi 2011, Royle et al. 2013a). We tested the GOF of the encounter process separately from the GOF of the underlying spatial point process. For the encounter process, we calculated a discrepancy measure for the cell-specific individual encounter frequencies to compare posterior samples and new realizations of the data generated from the posterior distribution. We used the Freeman-Tukey statistic to construct a Bayesian P-value where n i is the (observed or simulated) encounter frequency conditional on s i (activity center) and e i is the expected value under the model. The Pvalue is the proportion of time D(obs) . D(posterior).
For the point-process, we examined model GOF by testing whether estimated activity centers were independently and uniformly distributed over the statespace. We calculated a Bayesian P-value based on the statistic, I ¼ (G -1) 3 s 2 / n, where G is the total number of grid cells, and n and s 2 are the mean and variance of the number of activity centers per grid cell. We compare I from the estimated posterior distribution of the point-process to simulations under complete spatial randomness. We did not apply the point-process GOF test to models with the RSF covariate on the distribution of activity centers, because we would not expect activity centers to be independently and uniformly distributed across space for these models.
As a final metric of model goodness of fit, we compared the observed and expected number of individuals captured for each model to holistically examine both the point-process and detection process. We calculated the expected number of individuals captured with where E s i represents the exposure probability of an individual with an activity center at s i and n s i is the number of activity centers estimated at s i . By computing these values for each MCMC iteration, we constructed a 95% confidence interval for the number of individuals captured given the complete process described by the model. An observed number of captures that fell outside this range would indicate poor model fit.

RESULTS
We searched for mountain lion sign over a total of 8382 km during 98 person-days from December 7, 2012 to April 2, 2013. Search effort was distributed across 85 of 105 grid cells, and 80% of the effort occurred before February 15, 2014 (Appendix D). Hunter effort was unquantified, but there were only 8 grid cells in which a harvest but no live recapture occurred. Animals were sampled in 35 grid cells, and individual grid cells contained 0-6 samples (Fig. 1). A total of 80 samples from independent animals were included in the analysis, and 62 unique individuals were identified. Three individuals were identified from hair samples, 4 individuals were identified from scat samples, 43 individuals were v www.esajournals.org identified from biopsy muscle samples, and 30 individuals were identified from harvest samples. There were 25 male mountain lions identified and 37 female mountain lions identified. Fifteen individuals were recaptured 2-4 times during the sampling period (13 animals captured 2 times, 1 animal captured 3 times, and 1 animal captured 4 times). Thirty of the 62 individuals were harvested (17 female, 13 male), and of these, 10 (6 female, 4 male) were previously sampled and identified. The detection of new genotypes scaled linearly with search effort throughout the entire sampling period (Appendix D). No animals detected within the study area were later detected in the harvest sampling conducted outside the study area, suggesting limited movement occurred during the sampling period.

Model selection and goodness of fit
We evaluated 16 candidate models. Across all models, the parameter estimates for effort and RSF were consistently positive with 95% credible intervals that did not include zero (Appendix C). The 95% credible interval for the effect of sex on detectability and the sex-specific scale parameters overlapped zero in all models, but the effect of sex and sex-specific scale parameters were in the expected direction (Appendix C). The Bayesian P-value for the encounter process GOF produced reasonable results for all models and did not aid discrimination between models. Our ad hoc GOF measure of the predicted number of captures showed all models plausibly described the combination of the underlying point process and capture of individuals because the number of individuals captured (n ¼ 62) falls within the expected range. Therefore, we present results of all 4 models that included effort and RSF, and we selected the full model that included detection covariates for search effort and sex, RSF-driven densities and sex-specific activity distributions for further evaluation (Table 1).
Estimates from the full model indicated monthly detection probabilities were higher in grid cells with more search effort and mountain lion activity center densities were higher in areas with larger RSF values (Table 2). Females had higher baseline detection probability. Males were more likely to be detected farther from their activity center than females (e.g., males had larger observed movements; Appendix E). Using this best model, we estimated a median of 226 mountain lions over the entire statespace of 5,912 km 2 (Fig. 2), corresponding to a median realized density (D) of 3.8 (6 1.02 SD) mountain lions/100 km 2 . We estimated the proportion of males in the population was 0.41 (95% CI ¼ 0.26-0.61). Extracting estimates from the two management units of interest, in HD 250 we estimate a median of 82 (95% CI ¼ 54, 141) mountain lions, corresponding to a median density of 4.5 mountain lions/100 km 2 and a median of 79 (95% CI ¼ 51, 137) mountain lions in HD 270, corresponding to a median density of 5.2 mountain lions/100 km 2 .

Effects of including habitat quality on density estimates
Including the RSF as a covariate on the distribution of activity centers reduced estimated v www.esajournals.org mountain lion densities across all models (Fig. 3, Table 3), and had mixed effects on extrapolated densities in the two management units of interest. Comparing the best model with and without the RSF covariate, we found estimated abundance and density was reduced and the effective sampling area was similar (Fig. 4).
Median density was 44% less in models with the RSF covariate (3.6 lions/100 km 2 ; 95% CI ¼ 2.4, 7.4) compared to the average of models  Fig. 2. The spatial densities of mountain lions/4 km 2 across the statespace in the Bitterroot Watershed of western Montana estimated from the best model assuming uniform distribution of activity centers (Effort þ Sex þ r sex ; A) and best model that estimated the distribution of activity centers as a function of the resource selection function (RSF, Effort þ Sex þ r sex þ RSF; B) for sampling conducted during winter 2012-2013. Fig. 3. Effects of including a pre-existing mountain lion resource selection function (RSF) as a covariate on the distribution of activity centers on estimated median mountain lion population density (mountain lions/100 km 2 ) the Bitterroot Watershed of western Montana during winter 2012-2013 from 16 candidate spatial capture-recapture models. The uniform model estimates are based on a prior assumption of uniform distribution of activity centers and did not include RSF as a covariate on the distribution of activity centers. Error bars represent 95% credible intervals.
v www.esajournals.org without the RSF covariate (5.2 lions/100 km 2 ; 95% CI ¼ 3.4, 8.8). The inclusion of the RSF improved precision by decreasing the standard deviation by 7% and decreasing the width of 95% credible intervals by 10%. Within the two management units of interest, median density was 17% less in the best model with the RSF covariate (4.5; 95% CI ¼ 2.9, 7.7) as compared to the same model without the RSF covariate (5.4; 95% CI ¼ 3.4, 9.2) in lion management unit 250 and 4% greater in the best model with the RSF covariate (5.2; 95% CI ¼ 3.4, 9.1) as compared to the same model without the RSF covariate (5.0; 95% CI ¼ 3.1, 9.0) in lion management unit 270.

Effects of including harvest on density estimates
Treating harvested individuals as live captures did not create a significant bias in the detection process by introducing individuals that could not be detected after removal. When samples from harvested individuals were treated as live captures, we estimated a median density of 3.9 (95% CI ¼ 2.6, 6.7) mountain lions/100 km 2 , which did not represent a significant difference from the model that masked these individuals from sampling periods after they had been harvested (3.8, 95% CI ¼ 2.6, 6.5). Similarly, we observed no significant differences in any of the parameter estimates from models fit with these two permutations of the data.
The dead recovery component of the model did have a substantial effect on results by increasing sample size. When captures from harvested individuals were excluded from the sample, the data set had 50 total captures of which 5 were recaptures. With this reduced data set, we estimated a median density of 6.8 (95% CI ¼ 2.7, 16.6) mountain lions/100 km 2 . This median density represented a 78% increase over the model fit with the complete data set (3.8, 95% CI ¼ 2.6, 6.5). Moreover, with a reduction in sample size, the precision of the estimate decreased substantially. The standard deviation increased by 277% and the 95% credibility interval width increased by 251%. The effects on our estimate of density manifested through poor estimates of the sex-specific parameters on both detection probability and the scale of the activity distribution. Removing samples from harvested individuals eliminated all recaptures of male individuals such that the sex-specific parameters could not be estimated.

DISCUSSION
Our results indicate that incorporating prior knowledge of animal habitat selection into spatial capture-recapture models may improve model fit and the precision of the abundance estimates. In this case, the model where the probability of an activity center being located in a grid cell was a positive function of an existing RSF reduced the overall estimate of abundance by 44%, the SD by 7% and the CI width by 10%, an important improvement in both biological realism (i.e., high quality habitats had higher density than lower quality habitats) as well as precision of estimates. This approach to SCR modelling that increases the precision of estimates may increase the applicability of SCR modelling as an applied tool for monitoring trends in population abundance and effects of v www.esajournals.org various management actions on these trends. The overall strength of the SCR approach is the ability to explicitly estimate density over a defined area and account for animals whose activity ranges overlap the periphery of the area surveyed. This is particularly important for large carnivore species that tend to be wide-ranging and violate traditional capture-recapture geographic closure assumptions. For example, we estimated the location of activity centers across the statespace and calculated a statespace density, then estimate the density across two defined areas of interest (the hunting districts) as a function of the number of activity centers within the area of the hunting districts. In this manner, animals whose activity ranges overlap the periphery of the hunting districts are accounted for. These strengths add ecological realism, and address the challenge of comparing estimates across studies. The differences in estimated density between the statespace and management units of interest highlight the fact that density estimates are sensitive to the area over which the density estimates are generated. Thus, the ability to generate spatially explicit spatial abundances as a function of underlying habitat quality through the SCR approach may improve applicability of extrapolated density estimates beyond a given study area, making the approach more relevant to wildlife managers making decisions for larger landscapes rather than distinct study areas. Further improvements and refinements to SCR methods, including pooling data across study areas (Howe et al. 2013) or data sources (Gopalaswamy et al. 2012) will continue to improve rigor and reliability of SCR methods for monitoring trends in wildlife population abundance.
Our methods provide a method of integrating harvest into Bayesian SCR models. Similar methods exist within the likelihood framework Efford 2008, Efford 2014), however harvest has not been previously included within the Bayesian framework. SCR models have been applied to other species, for example, wolverines (Gulo gulo) and Lynx (Lynx lynx; Royle et al. 2011, Blanc et al. 2013 with potentially open harvest seasons without explicitly considering effects of ongoing harvest. Given that most large carnivore species are harvested or subject to management removals, our dead-recovery approach is likely to Our median estimates of 4.5 and 5.2 independent mountain lions/100 km 2 in the management units of interest represent the estimated density of a hunted population living sympatrically with wolves, and are higher than previously published estimates of mountain lion densities (Hornocker and Negri 2009). There could be important methodological as well as ecological reasons for our high estimates. Most published mountain lion density estimates are based on intensive, multi-year marking, radiocollaring and population monitoring methodologies that generate a mean minimum population density estimate. These studies assume that all individuals within the study population are detected, an individual's presence in the population may be backdated based on age to account for their presence prior to detection, and the study area includes the annual ranges of radiocollared mountain lions (Hornocker and Negri 2009). However, even within these types of intensive radiocollaring studies differences in sampling and estimation methodologies make comparisons of mountain lion density across study areas challenging. Specifically, the inclusion of different sex-age classes and, crucially, differences in area over which density is calculated (i.e., winter range vs. annual range) differs between studies and obviously challenges comparisons. Using SCR methods, we do not assume perfect detection and we include transient individuals in the estimate.
Additionally, two important ecological factors may be contributing to a high mountain lion density in the study area. First, mountain lion harvest in the study area has been conservative during the past decade and the population has likely increased throughout this period of conservative harvest management. Second, the abundance and diversity of prey in this area, resulting partially from ungulate management practices, likely sustains a higher than average mountain lion population (Carbone and Gittleman 2002, Karanth et al. 2004). The potential for individuals to move into the study area to occupy territories vacated by harvested individuals could result in overestimating mountain lion density, as our estimate is designed as an estimate of density pre-harvest (i.e., on the first day of the sampling period). However, the study area was located within a watershed under uniform mountain lion harvest management, and source-sink dynamics are unlikely in this scenario. Further, we designed our sampling plan to minimize the potential for movement into the study area by minimizing the duration of sampling and by concentrating our search effort during the period prior to the hunting season opening to the general public (see Appendix D). Although subtle shifts in territory use may have occurred as nearby animals were harvested, this effect would positively bias our estimates of activity range and negatively bias our estimates of population density. Therefore, shifting territories during the sampling period is not likely to explain our high density estimates. Finally, while SCR models do not necessarily allow flexibility in violating the closure assumption (but there are open SCR models that can do so), the explicit integration of space combined with choice of effective study area boundaries can minimize potential problems with closure compared to non-spatial capture-recapture models (Royle et al. 2013a).
The DNA sampling methodology used in this study reduces the time and effort involved in capturing and handling animals, but does not provide information about the age, body condition or other individual characteristics of an animal. Additionally, the ability to distinguish between transient and resident individuals is limited. For some areas and species transients may constitute a major portion of the population, and methodological development may be required to adjust density estimates for these areas. For many other large carnivores, transient or dispersing individuals comprise a significant portion of the harvestable population (e.g., wolves in Alaska where 50% of the harvest were such animals; Adams et al. 2008). We suggest that for mountain lions including transient animals in the population estimate is appropriate because these animals are present, likely affect the dynamics of local ungulate populations, and are legally harvestable during the hunting season. We also expect over the long-term the number of transients moving into and out of the study area are roughly equal, resulting in a consistent effect v www.esajournals.org on both mountain lion density and ungulate populations.
In our study, like many other studies of large carnivores, detection probability was relatively low, 0.09/month for males and 0.21/month for females with average effort at the center of the home range of the lion. Despite being low, our detection rates were similar to other large carnivore studies; for example, detection probability in non-spatial capture-recapture ranged from 0.12 to 0.26 for grizzly bears in British Columbia (Boulanger et al. 2002), from 0.11 to 0.26 for Bengal tigers in non-spatial models in India (Karanth and Nichols 1998), from 0.01 to 0.04/night for tiger density estimated in a spatial capture-recapture model (Royle et al. 2009), and from 0.05 to 0.09/period for jaguars in Brazil (Sollmann et al. 2013). In our case, with mountain lions, our detection probabilities resulted in reduced precision of our parameter estimates, increased estimates of population abundance well beyond the raw number of individuals identified, and emphasizes the cryptic nature of mountain lions. Future studies of large carnivores could probably improve precision of SCR estimates by incorporating existing information regarding habitat quality (or if unavailable, radiocollaring and estimating a RSF or some other index of habitat suitability) or developing sampling plans that target recaptures of sampled individuals to improve the estimation of the detection function (i.e., deliberately resampling areas that have previously been sampled).
We recommend that the decision to approach future large carnivore studies using the SCR or traditional radiocollaring approach be made based on the questions and applications of interest. For monitoring long-terms trends in animal abundance, non-invasive sampling techniques for capture-recapture studies provide a fast and economical method to estimate the number of individuals and monitor trends in segments of a population across time. The SCR method for the analysis of capture-recapture data provides estimates of density using a repeatable methodology which makes comparisons across time and space possible. Studies seeking to estimate vital rates, assess space use, distinguish between residents and transients or understand cause-specific mortality may be best approached using traditional tracking methodologies or capture-recapture approaches.