An integrated spatial capture–recapture approach reveals the distribution of a cryptic carnivore in a protected area

Abstract

Quantifying animal abundance, density, and distributions affords the opportunity to understand the effects of landscape structure and change on species of conservation interest, but estimating these parameters can be difficult for rare and cryptic species. Noninvasive sampling methods, such as remote cameras or scat DNA, can mitigate the challenges of studying rare and cryptic species while also minimizing effects on species of conservation interest or concern. Data derived from these methods can be integrated into robust, contemporary quantitative methods, including spatial capture–recapture (SCR) models, which provide a hierarchical framework for jointly estimating animal abundance, distribution, and, consequently, density. Herein, we developed an integrated SCR model to estimate the abundance, density, and distribution of a large carnivore of conservation interest, the cougar (Puma concolor), in a rugged and remote protected area, Yosemite National Park, California, USA. We combined spatial encounter data from DNA-based individual identification of scats with detection count data derived from remote cameras to estimate cougar density and detection probability in Yosemite. We further estimated how cougar density and detection probability varied as a result of vegetation, topography, anthropogenic and natural linear features, and survey effort. Using data collected in 2019 and 2020, we estimated a median of 31 (SD = 3.96, 95% credible intervals = 24–39) cougars in Yosemite across the two years, with higher densities associated with productive, vegetated areas. We found detection probability by scent detection teams was higher for females than males and positively correlated with survey effort, proximity to trails, and distance farther from roads and streams. Our study illustrates the utility of noninvasive survey methods that yield individual identities in rugged and remote environments, where capture and handling of cryptic, low-density animals is logistically challenging and cost prohibitive. Integrated modeling approaches, as used here, allow ecologists to leverage empirical data using a robust quantitative framework that can effectively address conservation objectives. Through this work, we demonstrate the importance of large, contiguous, and heterogeneous ecosystems to the ecology of wide-ranging species that occur in dynamic landscapes.

INTRODUCTION

Though humans have sustainably managed land for millennia, rapid increases in human-induced industrialization and resource extraction since the turn of the 20th century have profoundly altered the structure and dynamics of many contemporary landscapes (Díaz et al., 2019). Habitat loss and fragmentation resultant from changes in land use and anthropogenic development can disproportionately affect certain taxa, such as carnivores (Crooks et al., 2011). Many carnivores have experienced range contractions, extirpation, and extinction due to human persecution, habitat loss and fragmentation, and nonconsumptive human activities (Estes et al., 2011; Ripple et al., 2014). Although carnivores vary in their resilience to change, anthropogenic activities can negatively affect their behavior, space-use, and demography, often with negative implications for population persistence (Green et al., 2018; Heinemeyer et al., 2019; Suraci et al., 2021). Declines in large carnivore populations may induce trophic cascades and disrupt the regulation of prey, subordinate predators, and disease, pathogens, and parasites with broader implications for ecosystem function and resilience to change (Crooks & Soulé, 1999; Lundgren et al., 2022; Monk et al., 2022; Prugh et al., 2009). Consequently, understanding changes in abundance and distribution of carnivores can inform their conservation and inform the restoration of ecosystems.

Quantifying the abundance and distribution of rare species, including carnivores, is often integral to their conservation. Yet, it is difficult to make inferences about the status of rare species, given they can exhibit cryptic behaviors (e.g., nocturnality), physical attributes (e.g., camouflage), or demographic attributes (e.g., low densities) that limit their detection or observation (McDonald & Thompson, 2004). Such difficulties are compounded in rugged and remote environments that are difficult to access, and where traditional, hands-on methods are challenging to execute (Jones, 2011). Noninvasive survey methods, particularly those that yield species-specific detections and individual identities via genetic materials, are increasingly being implemented to study cryptic populations and have increased capacity to evaluate population demography at broad spatial scales (Long et al., 2008; Pauli et al., 2010; Quinn et al., 2019; Smith et al., 2021). Recent developments in quantitative methods have improved the utility of noninvasive data by integrating data collected using multiple methods to maximize inferences and increase the precision of population parameter estimates (Clare et al., 2017; Jiménez et al., 2022; Linden et al., 2018; Sun et al., 2019). Though considerations should be made about the spacing and extent of sampling efforts (Fleming et al., 2021; Sun et al., 2014), noninvasive methods have also proven useful for studying cryptic and rare species, as they can produce parameter estimates similar to those derived from live-capture data, maximize efficiency of data collection at broad spatial scales, and decrease impacts on study animals (Ruprecht et al., 2021).

One cryptic species, the cougar (Puma concolor), is the most widely distributed wild mammal in the western hemisphere (Hornocker & Negri, 2010) and occurs in a variety of biomes (Blake & Gese, 2016; Dellinger et al., 2018, 2020; Zeller et al., 2017). Despite their widespread distribution, cougars can be difficult to study at broad geographic scales because they exhibit large home ranges, wide-ranging movement patterns, and low population densities (Lendrum et al., 2014; Murphy et al., 2019; Ruprecht et al., 2021). Logistical challenges in rugged or remote portions of their distribution may further limit inferences about cougar populations. Though conflicts between individual cougars and humans intermittently occur, cougars typically avoid human activities and development (Kertson et al., 2011; Knopff et al., 2014; Suraci et al., 2020; Wilmers et al., 2021). Anthropogenic change has isolated some contemporary cougar populations and resulted in limited gene flow, reduced genetic diversity, and decreased connectivity among populations (Gustafson et al., 2019). Conversely, areas with lower human density, development, and fragmentation are often associated with higher genetic diversity, survival, and genetic connectivity of cougars, and are thought to provide demographic and genetic sources for isolated and fragmented populations (Burdett et al., 2010; Gustafson et al., 2019; Miotto et al., 2014). For example, spatially and genetically isolated cougar populations in southern California were granted special species status and were proposed for candidacy under the California Endangered Species Act in 2020, but populations are thought to be relatively stable elsewhere in California (Center for Biological Diversity and Mountain Lion Foundation, 2020). Despite this, little is empirically known about the abundance or density of cougars in these remote and less-developed areas, and whether they actually support higher cougar densities than developed or fragmented landscapes (e.g., Miotto et al., 2014).

In the Sierra Nevada of California, cougars are thought to be declining and occur within and interact with a diverse ecological community that includes a suite of other species of conservation interest (Thorne et al., 2006). Cougars forage primarily on ungulates, such as mule deer (Odocoileus hemionus), which are thought to be in decline throughout portions of their range due to rapid climate and land use changes in western North America (Mule Deer Working Group, 2021). Though cougars can affect prey populations by directly predating on them, they also affect prey, including ungulates or intraguild competitors, by changing their behavior and space-use patterns (Gaynor et al., 2022; Ruprecht et al., 2022). Cougars can also affect co-occurring species of conservation concern. For example, federally endangered species, such as Sierra Nevada red fox (Vulpes vulpes necator; U.S. Fish and Wildlife Service, 2021), and protected species, such as fisher (Pekania pennanti; U.S. Fish and Wildlife Service, 2020), and Sierra Nevada bighorn sheep (Ovis canadensis sierrae; U.S. Fish and Wildlife Service, 2007), have significantly declined from historical distributions and currently occur at low densities throughout the central Sierra Nevada. Consequently, these species are either proposed to be reintroduced or are actively being reintroduced to restore historical ecological communities in this region (Sierra Nevada Red Fox Conservation Advisory Team, 2021; U.S. Fish and Wildlife Service, 2007). Cougars can predate each of these species, but could also influence their behavior, foraging, and space-use, which could inhibit the success of proposed restoration efforts (Ernest et al., 2000; Gabriel et al., 2015; Gammons et al., 2021; Quinn et al., 2018; Rominger et al., 2004). Clarifying the density and distribution of cougars, including how landscape structure influences their occurrence, may both inform their conservation and mitigate potential negative effects on other species of conservation concern.

Herein, we developed an integrated spatial capture–recapture (SCR) model to estimate the abundance and distribution of cougars within a study area that encompassed Yosemite National Park, California, USA, and the surrounding area (hereafter, “Yosemite”). We used noninvasive methods to estimate cougar abundance and distribution, which included DNA from scats collected by detection dog teams and detection–nondetection data collected by remote cameras. Additionally, we wanted to understand how human development, topography, hydrology, and vegetation affected the detectability of cougars and the distribution of their activity centers in Yosemite. Through this work, we demonstrate the utility of noninvasive methods, particularly in logistically challenging areas, to understand the demography and distributions of a cryptic species of conservation interest. Further, we offer a broader context of the role that protected areas serve in the regional distribution of a wide-ranging carnivore in the heterogeneous landscape of western North America.

METHODS Study area

We surveyed for cougars within and immediately adjacent to the 3074-km² boundary of Yosemite National Park, approximately 2850 km² of which is designated wilderness, in the central Sierra Nevada of California, USA (latitude: 37.884, longitude: −119.529; Figure 1). A small portion of survey effort (i.e., <5%) occurred outside of Yosemite in U.S. Forest Service Stanislaus, Humboldt-Toiyabe, and Inyo National Forests. Climate in Yosemite is characterized by cool, wet winters and warm, dry summers, with the majority of average annual precipitation (93.3 cm) falling as snow (88.9 cm) (Cleland et al., 2007). Topography in Yosemite is complex, with elevation ranging from 600 to 4000 m above sea level. Vegetation varies by elevation, with oak woodlands and chaparral predominant at lower elevations, transitioning to mixed-conifer forests at moderate elevations and subalpine forests at higher elevations (Mayer & Laudenslayer, 1988). Above approximately 3300 m, the landscape is characterized as alpine and large areas are comprised of talus and scree slopes. Research in Yosemite was permitted by the U.S. National Park Service under permit number YOSE-2019-SCI-0080.

FIGURE 1. Study area in and around Yosemite National Park, California, USA, where we incorporated remote camera data and scat detection data in an integrated, hierarchical spatial capture–recapture model to evaluate the distribution and density of cougars (Puma concolor) in a 9234-km2 state-space. Sampling occurred in 2019 and 2020 and resulted in 121 genotyped cougar scats and 1340 remote camera images of cougars. We also depict the distribution of spatial covariates, including normalized difference vegetation index (NDVI), slope, roads, trails, and rivers, which were incorporated as effects on detection probability and the distribution and density of cougar activity centers.

Field methods Detection dog team study design and surveys

From 2 August to 14 October 2019 and 8 July to 30 October 2020, three detection dog teams from Rogue Detection Teams surveyed for cougars and co-occurring carnivores within Yosemite (Appendix S1: Figure S1). Teams were comprised of a human handler and a dog who completed laboratory and simulated field trials. Each handler and dog also had 2–10 years of experience searching for cougars during previous research efforts in other study systems. Detection dog teams focused on a suite of species, prioritizing surveys to detect cougars and two regionally rare species, Sierra Nevada red foxes, and fishers. Teams also collected samples from sympatric carnivores whose scats were similar to these target species, including gray foxes (Urocyon cinereoargenteus), coyotes (Canis latrans), and Pacific martens (Martes caurina). Maintaining dog focus and enthusiasm for rare species can be increased by having multiple targets, allowing a reward for detection of target species (Richards et al., 2021).

To determine where detection dog teams would focus their efforts, the survey area was split into 10.4-km² hexagonal cells following the multispecies monitoring approach developed by the California Department of Fish and Wildlife (Sierra Nevada Red Fox Conservation Advisory Team, 2021). Due to rugged terrain, low road density, and wilderness status within Yosemite, we implemented a modified sampling design to increase efficiency and mitigate access issues, while also surveying a sufficient portion of Yosemite to inform estimates of cougar abundance and distribution. We used a “clover” design where each clover was comprised of seven 10.4-km² hexagonal cells (Appendix S1: Figure S1). Each detection team camped in the center cell of a clover, or as near to the center as access or terrain allowed, for a week and surveyed all six cells surrounding the clover center. Twenty clovers were sampled in 2019 and 23 clovers were sampled in 2020, and the average distance between nearest clover centers was 10.63 km (±1.5 km, 9.14–13.86 km) apart in 2019 and 9.15 km (±0.01 km, 9.14–9.16 km) apart in 2020 (Appendix S1: Figure S1). Detection teams surveyed between 6–8 h or over a linear distance of 5–15 km within each surveyed cell. Handlers did not follow specific survey routes and were encouraged to follow natural travel corridors, such as ridgelines, saddles, drainage bottoms, and wildlife trails to maximize the likelihood of detection and survey effort over a gradient of habitat types. Detection teams recorded survey routes using a GPS track to determine the distance traveled, which was subsequently used to quantify survey effort. When a scat was encountered, the handler recorded the GPS location of the scat, photographed the scat with a Sharpie marker for scale, and placed the entire scat in a paper bag. Paper bags were subsequently sealed in 3.8 L resealable plastic bags containing approximately 60 mL of desiccant beads to absorb moisture. At the end of each survey day, approximately 2 mL of each collected scat was placed in a 15 mL centrifuge tube with at least 8 mL of 95% ethanol, where it was stored until laboratory analysis.

Remote camera surveys

In 2019 and 2020, we deployed remote cameras to augment scat collection efforts. Obtaining detection–nondetection data that both overlapped with detection team survey effort and occurred outside of detection team survey areas provided additional data to increase the precision of estimates of activity center distributions (e.g., Ferreras et al., 2021). We deployed remote cameras (Reconyx Hyperfire, Holmen, WI, USA) at 85 locations throughout Yosemite, focusing on natural corridors (e.g., recreational and wildlife trails), areas where cougar sign (e.g., scrapes) was observed, and vegetation structure that captured the range of vegetation conditions, ranging from productive, continuously vegetated areas in lower elevation portions of Yosemite to largely nonvegetated areas in alpine portions of Yosemite (Figure 1). Each remote camera was deployed on a wildlife trail, or similar linear feature, within 100 m of a randomly generated location. Cameras were mounted to a north-facing tree immediately adjacent to the trail and angled at approximately 45° to parallel (0°) with the trail to maximize the likelihood of capturing photographs of animals moving through the field-of-vision. Remote cameras were not baited, but a small amount of olfactory lure (Gusto, Minnesota Trapline Products, Pennock, MN, USA) was placed at the camera station upon deployment. Remote cameras were programmed to take a burst of five photographs each time the camera was triggered. The average distance between the nearest remote cameras was 2.24 km (±1.72 km, 0.17–11.82 km). We managed and tagged photographs to species using Colorado Parks Wildlife (CPW) photo software (Ivan & Newkirk, 2016) in conjunction with the photo viewer software IrfanView (https://www.irfanview.com/). The database was structured such that two independent observers viewed and tagged each photograph. If observers disagreed on the species' assignment, a third independent observer reconciled the discrepancy, resulting in more accurate species identification (Ivan & Newkirk, 2016). We grouped photographs to a daily observation interval, where a species was considered “detected” at a remote camera on a given date if at least one photograph of the species was obtained on that date.

Laboratory analyses Overview

We extracted DNA, identified species, and genotyped cougar samples at the Mammalian Ecology and Conservation Unit at the University of California, Davis. We extracted DNA from scats using the QiaAmp Stool Kit (Qiagen, Inc.) according to the manufacturer's instructions except that we eluted in 50 μL of buffer. Identification of individual cougars required a two-step process, including species identification and individual and sex genotyping.

Species identification

For species identification, we attempted to sequence each DNA extract at a fragment of the cytochrome b gene using primers RF14724: 5′-CAACTATAAGAACAT-TAATGACC-3′ and RF15149: 5′-CTCAGAATGA-TATTTGTCCTC-3′ (Perrine et al., 2007). We conducted PCR reactions in 11 μL total volume and consisted of 2 μL of DNA extract, 1.1 μL each of 10× PCR buffer, 2-mM dNTPs, and 25-mM MgCl₂, along with 0.28 μL each of 20 μM forward and reverse primers, 0.2 μL of 5 U/μL Taq (AmpliTaq, ThermoFisher Scientific) and 4.84 μL water. The thermal cycling conditions were 94°C for 3 min, then 45 cycles of 94°C for 30 s, 50°C for 30 s, and 72°C for 45 s, followed by 10 min at 72°C. Dye terminator sequencing reactions were performed for each PCR product for each primer using Applied Biosystems reagents with products sequenced on an ABI 3730 capillary sequencer (Applied Biosystems, Foster City, CA, USA). Sequences were aligned manually in Sequencher version 5.4.6 and compared against a database of orthologous sequences from all carnivore species occurring in the study area. When sequences did not match a carnivore, we used the basic local alignment search tool (BLAST) on the nucleotide database in GenBank (https://blast.ncbi.nlm.nih.gov/Blast.cgi) to obtain matches to other species, which was usually a prey species, such as deer. Given it is highly unlikely canine and human observers confused herbivore scat for carnivore scat, we presumed that prey species amplifying in the sample represented diet material. Therefore, we assumed that extracts sequenced as “deer” likely represented cougars, or other carnivores, particularly in cases where the scat was suspected to come from a cougar during sample collection in the field based on scat size and shape. In such cases, cougar DNA may have been present in the sample but was obscured by more abundant prey DNA.

Individual and sex genotyping

We attempted to genotype DNA extracts that were from cougars using a panel of 95 single-nucleotide polymorphisms (SNPs) designed specifically for individual and sex identification in cougars (Buchalski et al., 2022). Because the panel only reliably amplified felid DNA (and included markers that differentiate among felid species), we also included samples species-typed in the previous step as “deer” (i.e., samples presumed to be cougar scats containing deer as prey). We used Fluidigm's 96.96 Dynamic Arrays integrated fluidic circuits (IFCs) in conjunction with the Biomark HD system to conduct a pre-amplification reaction followed by a genotyping reaction according to the manufacturer's instructions except that we altered the pre-amplification and genotyping amplification steps to accommodate low-quantity, low-quality DNA samples. Specifically, for the pre-amplification step, we used 3 μL DNA (as opposed to 1.25 μL) and 19 PCR cycles (as opposed to 14 cycles) and diluted the PCR product 1:10 (as opposed to 1:100). For the genotyping step, we used seven additional cycles (41 vs. 34 post-stepdown cycles). We called genotypes using Fluidigm SNP Genotyping Analysis Software version 4.5.1. We excluded genotypes with >15% missing data assigned individuals based on pairwise comparisons of genotypes. The assay was designed for high resolution, such that pairwise estimates of allele sharing between sample genotypes results in a bimodal distribution, with a lower mode at approximately 75% allele sharing corresponding to pairs of genotypes from different individuals (including siblings and other close relatives) and a higher mode at 100% allele sharing, corresponding to pairs of genotypes from the same individual. Based on inspection of this distribution, we considered genotypes matching at >95% of alleles to be the same individual, and genotypes matching at <85% to be distinct individuals. Only in rare cases did two sample genotypes match at >85% but <95% of alleles, which typically occurred when one or both sample genotypes contained close to 15% missing data (the maximum we allowed). Therefore, we assumed these cases reflected high genotyping error rather than improbably similar genotypes.

Hierarchical modeling Overview

We developed a hierarchical, integrated SCR model to estimate the abundance, distribution, and density of cougars in Yosemite (Ferreras et al., 2021; Linden et al., 2018; Royle & Young, 2008; Figure 2). In this integrated SCR model, we incorporated scat detections, which provided locations of known individuals, and remote camera detections of unknown individuals (sensu Ferreras et al., 2021) to increase the precision of activity center estimates and improve density estimates of cougars. We used a nested approach, where we modeled the effects of survey year, sex of individual cougars, and spatial covariates on detection probability, analogous to finer scale space-use and probability of detecting cougars during detection team surveys (i.e., within the 3 × 3 km or 9 km², grid cell x), and modeled the effects of vegetation on the occurrence of activity centers, the latent variable s_i that is analogous to broader scale space-use of cougars and their distribution and density throughout the model state-space (i.e., within the 9 × 9 km or 81 km², grid cell c). To account for cougars whose activity centers may occur on the boundary of the survey area, we buffered our survey area by 22 km (i.e., approximately twice σ of male cougars; Hornocker & Negri, 2010; Royle et al., 2014). This resulted in a discrete state-space of 9234 km² throughout which we predicted cougar distribution and density. All analyses were completed in R version 4.0.2 (R Core Team, 2020).

FIGURE 2. Conceptual diagram of data collection and the joint analysis of spatial capture–recapture (SCR) and count data collected by remote cameras through a hierarchical, integrated SCR model. On the left, we depict a hypothetical state-space in which we collected cougar (Puma concolor) scats (hypothetical individuals represented by different colors) via detection dog team surveys and photographs via remote cameras. These data are subsequently incorporated into an integrated SCR model that estimates the effects of the landscape on detection probability, the distribution of cougar activity centers (circles), and, thus, the total number of estimated cougars throughout the state-space. On the right, we depict the structure of our integrated SCR model. We show data sources (green rectangles) that consequently informed estimates of population size N (orange square). Blue nodes represent observed quantities of two observation processes: (1) detections of unique cougars from detection dog team surveys (yxit) and (2) count data (yjt) of nonidentified individuals collected at remote cameras. White nodes represent latent variables and model parameters, including: λxit, SCR encounter rate; Λjt$$ {\Lambda}_{jt} $$, remote camera encounter rate; si, activity centers; zi, latent binary variable describing the membership of individual i in the population; gk, sex-specific detection function based on movement parameter σ. Conceptual diagram derived from Ferreras et al. (2021) and Barry et al. (2021).

Spatial covariates

We incorporated covariates in our analysis of cougar distribution and density to estimate the effects of landscape covariates on both detection probability and distribution of activity centers of cougars throughout the defined state-space. We included survey effort, topography, and distance to natural and anthropogenic linear features as covariates on detection probability. We determined year-specific survey effort within each 9-km² grid cell by calculating the total distance traveled by detection dog teams during surveys in each cell. Given cougar space-use can be affected by topography (Dellinger et al., 2018; Zeller et al., 2017), we also included the effect of slope at this scale. We derived slope data at a 30-m resolution from the US Geological Survey's (USGS) topographic dataset (U.S. Geological Survey, 2019a). We rescaled slope values using the exactextractr package's “exact_extract” function (Baston, 2022), which allowed us to calculate the mean of all slope pixels that occurred within each of the finer scale grid cells (i.e., 9 km²). We also incorporated the effect of distance to three types of linear features, recreational trails, roads, and perennial rivers and streams, on detection probability. We included recreational trails and roads to represent human activity and development (Kertson et al., 2011; Smereka et al., 2020) and included perennial rivers and streams to represent riparian corridors that prey may use for travel or foraging (Smereka et al., 2020; Zeller et al., 2017). Recreational trail, road, and major river spatial layers were obtained from the US Department of Interior's Integrated Resource Management Application Data Store (https://irma.nps.gov/DataStore). We also incorporated perennial rivers and streams described by the USGS National Hydrography Dataset (U.S. Geological Survey, 2019b), which included rivers and streams that occurred within our study area, but which were outside the Yosemite boundary and consequently not present in the Yosemite dataset. We calculated the minimum distance between the grid cell centroid and the closest of each of these linear features using the “gDistance” function in the “rgeos” R package (Bivand et al., 2017).

We calculated normalized difference vegetation index (NDVI) values to predict the distribution of cougar activity centers throughout Yosemite within each 81-km² cell. NDVI represented variation in vegetation greenness and is a proxy for primary productivity (Pettorelli et al., 2005). We considered primary productivity a proxy for the potential distribution of cougar prey, such as mule deer (O. hemionus), other ungulates, and small mammals, which occur in productive, vegetated areas that provide sufficient forage (Karandikar et al., 2022; Torstenson et al., 2006). Further, productive vegetated areas may provide ambush cover for cougars, which often use tall grasses, shrubs, and other vegetation cover to stalk and successfully capture prey (Williams et al., 2014). NDVI values were calculated by NASA at a resolution of 250 m every 16 days (Masek et al., 2006). In Google Earth Engine, we used these values to create a composite NDVI raster that reflected the average vegetation greenness during the survey periods in 2019 and 2020. We rescaled NDVI values using the exactexractr package's “exact_extract” function, which allowed us to calculate the mean of all NDVI pixels that occurred within each of the activity center grid cells (i.e., 81 km²).

Integrated SCR model

In our hierarchical formulation, spatially explicit detection counts y in grid cell x for individual i in year t (y_xit) were modeled as random variables that were a function of the location of an individuals' latent activity center (s_i), such that Pr(y_xit = 1|s_i):[Image Omitted. See PDF] [Image Omitted. See PDF]where y_xit is a latent random variable with a distribution defined by average encounter rate λ_xit and $$$ {r}_0 $$$ is a gamma-distributed random variable used to define the negative binomial distribution.

We modeled the average encounter rate λ_xit for individual i in grid cell x in year t as:[Image Omitted. See PDF]where average encounter rate is a function of detection probability p0_xit, a sex-specific detection function g_k, the Euclidean distance between grid cell x to the activity center of individual i (dist_xi²), and the Bernoulli distributed random variable z_i, which indicates the probability the individual is part of the population. We modeled sex of all individuals as the Bernoulli distributed random variable, sex_i ~ Bernoulli(Ψ_sex), where sex was estimated as coming from the population-level sex ratio (Ψ_sex).

The detection function g_k described how the sex-specific encounter rate of cougars decreased as a function of the distance between their activity center and the location of a known location (i.e., scat sample). We modeled g_k with a Gaussian encounter probability such that:[Image Omitted. See PDF]where capture probability of an individual of sex k is a function of the sex-specific movement parameter σ_k, which represents the distance to the latent activity center and was modeled as a random variable with a uniform distribution.

We modeled the probability of detection p0_xit as:[Image Omitted. See PDF]where the probability of detecting individual i in grid cell x in year t is a function of: an intercept (β₀), grid cell and year-specific survey effort by detection dog teams (β₁), sex of the cougar (β₂), year of the survey (β₃), slope within the grid cell (β₄), distance between the grid cell centroid to the nearest trail (β₅), distance between the grid cell centroid to the nearest road (β₆), and distance between the grid cell centroid to the nearest perennial river (β₇).

To model variation in cougar density, we used an inhomogeneous point process to estimate the distribution of activity centers s_i within our study area (Borchers & Efford, 2008). We modeled the expected density for each grid cell c in our discrete state-space as the intensity of a point process conditional on a linear model of spatially varying covariates:[Image Omitted. See PDF]where α are regression coefficients for the log-linear model estimating the expected density of cougars (Royle et al., 2014). The expected density of cougars in each given grid cell c was a log-linear function of NDVI within the grid cell and the area of the cell:[Image Omitted. See PDF]where expected density is a product of an intercept (α₀), vegetation productivity (α₁), and the area of the pixel.

The final component of the SCR model involved defining the distribution of cougar activity centers. For a basic SCR model having constant density, such that activity centers are distributed uniformly throughout the state-space, the probability of an activity center located in any given grid cell would be 1/c. As we were modeling variation in density conditional on spatial covariates, the probability was a ratio of the intensity function at a given grid cell, conditional on the coefficients of the linear model and the spatial covariate values within that grid cell, and the summed intensity function across all grid cells in the discrete state-space, such that:[Image Omitted. See PDF]Using data collected in 2019 and 2020, we calculated a single abundance estimate, which represents the total number of individuals N alive during the study period and was estimated as:[Image Omitted. See PDF]

The second dataset we incorporated into this integrated model was comprised of all remote camera detections of cougars. We modeled these count data under the conditional formulation (Chandler & Royle, 2013) such that:[Image Omitted. See PDF] [Image Omitted. See PDF] [Image Omitted. See PDF]In this formulation, ycam_jt are the latent spatial encounter histories, which were modeled as a function of remote camera detection probability p_jt on day j in year t with a negative binomial distribution.

Model fitting and assessment

The integrated SCR model was fit using the Markov chain Monte Carlo (MCMC) methods of JAGS version 4.2.0 (Plummer, 2003). We standardized each covariate to have a mean of zero and a SD of one. We tested for correlations among spatial covariates using Pearson's correlation coefficient, and all covariates had correlations <|0.6|. We used uninformative prior distributions for all estimated parameters. We calculated parameter estimates from 9000 MCMC samples, taken from 3 chains run for 50,000 iterations, thinned by 10, following a burn-in of 100,000 iterations. We assessed model convergence by examining trace plots and $$$ \hat{R} $$$ values for convergence (Gelman et al., 2013; Gelman & Hill, 2006). We present significance of covariates as percent probabilities from posterior distributions, which we calculated as the percent of posterior draws greater or less than zero, depending on the sign of the median value. We present the median estimate (±SD [±95% credible intervals]) for our results.

RESULTS Data collection

Detection teams surveyed 1287 linear km from 2 August to 14 October 2019 and 1853 linear km from 8 July to 30 October 2020. The teams collected 757 scats in 2019 and 1148 scats in 2020 for a total of 1905 scats. From these 1905 samples, we successfully typed 1583 (83.1%) to species, including 98 cougar scats (Appendix S1: Table S1). We then attempted to genotype the 98 cougar samples along with 37 suspected cougar samples (species-typed as deer) to individual. From these 135 samples, we successfully genotyped 121 samples, including 96 (98%) of those known to be cougar and 25 (68%) suspected cougar samples, and identified 44 unique cougars (25F, 19M). We located an average of 2.80 (±2.65) scats from each individual during the study period (Appendix S1: Figure S2). We deployed 42 remote cameras from 24 July to 4 November in 2019 and 43 cameras from 1 July to 21 October in 2020. These cameras collected 102,131 images, including 1340 images of cougars at 32 remote camera stations representing 91 unique detection events.

Integrated SCR model Cougar abundance and density

Using data collected in 2019 and 2020, we estimated a median of 31 (±3.96 [24–39]) cougars during the two-year sampling period within the 3027-km² boundary of Yosemite. This estimate included all age classes (i.e., juveniles, sub-adults, adults), as we were unable to distinguish cougar age from scat data. We estimated a median of 11 (±2.63 [5–15]) males and 21 (±2.97 [15–26]) females (Table 1, Figure 3) in Yosemite. In the broader 9234-km² state-space, we predicted there was a median of 84 (±13.89 [60–112]) cougars, with a median of 27 (±5.61 [18–38]) males and 56 (±12.37 [35–81]) females. Predicted density estimates were similar within the broader state-space and within the Yosemite boundary, with a mean of 0.75 (±0.78 [0.01–2.42]) and 0.88 (±0.89 [0.02–2.58]) cougars per 81-km² grid cell, respectively. Similar to abundance estimates, mean predicted density of female cougars in each 81-km² grid cell was higher than male cougars in both the broader state-space (females: 0.51 ± 0.55 [0.00–1.79]; males: 0.24 ± 0.26 [0.00–0.79]) and within the Yosemite boundary (females: 0.58 ± 0.66 [0.00–1.99]; males: 0.29 ± 0.29 [0.01–0.81]). The movement parameter $$$ \upsigma $$$ varied by sex, with male $$$ \upsigma $$$ estimated as 11.29 km (±1.65 [8.81–15.17]) and female $$$ \upsigma $$$ estimated as 6.30 km (±0.64 [5.25–7.77]). We found a slight bias in the estimated sex ratio of the population, with more females estimated to be in the population than males (0.32 ± 0.08 [0.19–0.49]; Table 1).

TABLE 1 Posterior parameter estimates, including mean, SD, 95% credible intervals, and median estimates, from hierarchical, integrated spatial capture–recapture model.

Parameter	Mean	SD	2.50%	Median	97.50%
Nyose	32	3.96	24	31	39
N	85	13.89	60	84	112
α0	−5.19	0.24	−5.71	−5.18	−4.78
α1	1.05	0.27	0.56	1.04	1.61
β0	−5.73	0.47	−6.67	−5.72	−4.83
β1	0.61	0.11	0.41	0.61	0.82
β2	−0.73	0.42	−1.57	−0.73	0.07
β3	0.38	0.19	0.02	0.38	0.75
β4	0.02	0.13	−0.24	0.02	0.28
β5	−2.39	0.45	−3.31	−2.37	−1.56
β6	0.17	0.13	−0.08	0.17	0.41
β7	0.23	0.30	−0.37	0.23	0.79
σfemales	6354.22	644.04	5252.67	6297.02	7768.26
σmales	11,471.21	1650.61	8813.99	11,289.85	15,174.85
Ψsex	0.33	0.08	0.19	0.33	0.49

Note: Included are the estimated abundance of cougars within Yosemite National Park, CA, USA (N_yose), estimated abundance of cougars in the discrete state-space (N), intercept (α₀), and effect of normalized difference vegetation index (α₁) on activity center densities, detection probability intercept (β₀), the effects of survey effort (β₁), sex (β₂: male), year (β₃), slope (β₄), and distance to trails (β₅), roads (β₆), and perennial rivers (β₇) on detection probability, sex-specific movement parameters in m (σ), and the sex ratio of the study population (Ψ_sex) where males are assigned a value of 1 and females are assigned a value of 0.

FIGURE 3. Predicted density and distribution of cougar (Puma concolor) activity centers in 81-km2 grid cells estimated using a hierarchical, integrated spatial capture–recapture model. The black line indicates the boundary of Yosemite National Park, whereas the study area state-space extended beyond this to account for cougars that may occur on or near the boundary of Yosemite. Darker grid cells indicate areas of higher predicted cougar density (i.e., more predicted activity centers) and lighter colors indicate relatively lower predicted cougar density. We estimated a median of 31 (±3.96 [24–39]) cougars and a median of 11 (±2.63 [5–15]) males and 21 (±2.97 [15–26]) females within the 3027-km2 boundary of Yosemite. In the broader 9234-km2 state-space, we estimated a median of 84 (±13.89 [60–112]) cougars, with a median of 27 (±5.61 [18–38]) males and 56 (±12.37 [35–81]) females. Density estimates were similar within the broader state-space and within the Yosemite boundary, with a mean of 0.75 (±0.78 [0.01–2.42]) and 0.88 (±0.89 [0.02–2.58]) cougars per 81-km2 grid cell, respectively. Predicted density of female cougars in each 81-km2 grid cell was higher than male cougars in both the broader state-space (females: 0.51 ± 0.55 [0.00–1.79]; males: 0.24 ± 0.26 [0.00–0.79]) and within the Yosemite boundary (females: 0.58 ± 0.66 [0.00–1.99]; males: 0.29 ± 0.29 [0.01–0.81]).

Covariate effects on detection probability and activity centers

Several covariates affected the probability of detecting cougars (Table 1, Figure 4). The probability of detection was positively associated with survey effort, which had a 100% probability of increasing detection probability (β₁ = 0.61 ± 0.10 [0.41–0.82]), which corresponded to an approximate 0.5% increase in detection probability for each 1 km increase in effort within a detection grid cell. We also found that detection probability was greater in the second year of surveys, with a 98% probability of the second year having higher detection probability (β₃ = 0.38 ± 0.19 [0.02–0.75]), which corresponded to an approximate 1.5% increase in the probability of detecting cougars. Though the credible intervals of both covariates overlapped zero, distance to roads and streams was positively related to detection probability, with a 91% probability of encountering cougars farther from paved roads (β₆ = 0.17 ± 0.13 [−0.08 to 0.41]) and a 77% probability of encountering cougars farther from perennial rivers (β₇ = 0.23 ± 0.30 [−0.37 to 0.79]). Conversely, distance to trails had a 100% probability of decreasing detection probability, suggesting cougars were more likely to be detected in areas closer to trails (β₅ = −2.37 ± 0.45 [−3.31 to −1.56]). Each additional 500 m from trails corresponded to a 2% decrease in predicted detection probability, and predicted detection probability declined to approximately zero when the average distance from trails was 10 km or greater. Sex also affected detection probability, with a 97% probability of male cougars having lower detection probability than females (β₂ = −0.73 ± 0.42 [−1.57 to 0.07]). Though the probability of detecting female cougars was approximately two times the probability of detecting male cougars, this only corresponded to a 2.1% increase in predicted detection probability. Slope had no discernable effect on the probability of detecting cougars (β₄ = 0.02 ± 0.13 [−0.24 to 0.28]). We found positive effects of vegetation productivity on the distribution of cougars; there was a 100% probability of cougar activity centers being more likely to occur in areas with greater values of NDVI (α₁ = 1.04 ± 0.26 [0.56–1.61]). Predicted cougar densities in state-space grid cells with high NDVI values (NDVI = 0.60) were approximately 30 times greater (2.8 [1.1–4.8] cougars/cell) than grid cells with low NDVI values (NDVI = −0.05; 0.09 [0.01–0.20] cougars/cell).

FIGURE 4. Effect of covariates on the detection probability (i.e., fine-scale space-use) of cougars (Puma concolor) in and around Yosemite National Park in California, USA. We show the median effect (square symbol), 80% credible intervals (thick line), and 95% credible intervals (thin line) of the covariate effects on the detection probability of cougars. Gray covariates indicate the 95% credible intervals overlap zero, whereas black covariates indicate the 95% credible intervals do not overlap zero. Values greater than zero suggest a positive relationship between the covariate and space-use, whereas values less than zero suggest a negative relationship between the covariate and detection probability.

DISCUSSION

Effective conservation often relies on the ability to estimate the distribution and abundance of species of interest. Here, we leveraged two noninvasive sampling approaches and an integrated SCR model to estimate the density and distribution of cougars in a rugged, remote, and logistically challenging landscape. During the two-year sampling period, we estimated a median of 31 cougars occur in Yosemite and 84 cougars in the broader landscape encompassing Yosemite, and predicted a positive effect of productive vegetation on the distribution and density of cougar activity centers. Further, we identified survey-specific and spatial covariates that influenced the likelihood of detecting cougars from detection dog team surveys. Despite working in rugged and undeveloped terrain, we were successful in estimating cougar abundance and density, and our findings corroborate previous research highlighting the efficacy of noninvasive detection dog teams in elucidating the demography and ecology of cryptic species (Davidson et al., 2014; Long et al., 2007; Petroelje et al., 2021). Though we did not conduct a formal comparison between traditional and integrated SCR models, our quantitative approach joins a rich body of literature demonstrating the utility of integrated approaches in estimating animal abundance and density (e.g., Morin et al., 2022; Ruprecht et al., 2021, 2022; Tourani et al., 2020; Twining et al., 2022) while also collecting ancillary data on locally rare species of interest (Appendix S1: Table S1, Figure S3). Through this work, we demonstrated the utility of noninvasive sampling approaches to elucidate important population parameters and provide foundational information about the population of a species of conservation interest.

The probability of detecting cougars during detection team surveys was influenced by survey methodology, cougar demography, and landscape structure. Greater survey effort increased the likelihood of detecting cougars and we found that detection probability was higher in the second year of surveys. This could be the result of the increased efficiency and efficacy of detection teams that occurred after the initial survey in 2019 (e.g., Richards et al., 2021). Detection probability was lower for male cougars than for female cougars, likely due to males exhibiting higher movement rates and larger territories than females (Benson et al., 2019). Due to these intersexual differences in space-use, any given survey route that intersects a single male cougar's putative activity center may intersect several female activity centers, resulting in a higher likelihood of detecting female cougars during surveys. Cougars were more likely to be detected near trails and were less likely to be detected near roads. Cougars may use linear features for travel, but are also wary of human activity (Riley et al., 2021). Trails may facilitate efficient movement across the landscape, while also being lower traffic and proximal to contiguous habitat, which allows cougars to avoid humans when they co-occur on trails (Benson et al., 2021). Alternatively, roads are relatively uncommon in Yosemite (Figure 1) and are associated with the highest rates of human use and density in the study area. Thus, cougars may have avoided these areas given their perceived association with human occurrence (e.g., Kertson et al., 2011). Despite the bulk of survey effort occurring in backcountry areas and off-trail, higher detection probability near trails may also reflect the challenges of implementing detection team surveys in rugged or remote environments. Cougars were more likely to be detected farther from perennial streams and rivers, though this effect was statistically weak. Cougars are thought to be associated with waterways given their association with prey activity (e.g., Zeller et al., 2017), but recent studies demonstrate ungulates may avoid certain sites (e.g., feeding or drinking areas) during temporally risky periods to minimize overlap with predators (Candino et al., 2022; Gaynor et al., 2022). It is also possible that there was no effect of hydrological features on the probability of detecting cougars due to their ubiquity, changes in persistence due to drought, or accumulation of dead and downed trees due to fire suppression and beetle kill. Though previous work shows cougars can effectively use steep, rugged terrain (Dunford et al., 2020), we found slope had no discernable effect on the probability of detecting cougars (Table 1, Figure 4). This could result from a difference in the scale of effect for microsite topography, where we were not able to detect any effect of slope on cougar space-use at this scale. Alternatively, changes in vegetation structure at differing slopes, or combined effects on prey distribution, may be more important to shaping cougar space-use (Peterson et al., 2021).

Cougar activity centers were more likely to occur within productive, vegetated areas, and cougar density was positively associated with increased vegetation productivity (Table 1). These results corroborate previous research that shows cougars are often associated with contiguous vegetation cover where they can pursue prey (Williams et al., 2014) and avoid humans or other large predators (LaBarge et al., 2022; Nickel et al., 2021; Smith et al., 2022). We posited productive vegetation may correlate with the occurrence of prey, such as mule deer, which often occur in vegetated areas that provide sufficient forage (e.g., Torstenson et al., 2006). However, mule deer were ubiquitous in Yosemite and were detected at 76% of camera stations during the two-year sampling period. Alternatively, it is possible that productive vegetation may provide cover and crypsis for these ambush predators and increase their likelihood of successfully acquiring prey (Allen et al., 2015; Andruskiw et al., 2008; Lehman et al., 2017; Pierce et al., 2004). Composite covariates, including NDVI, are helpful for interpreting the broad effects of vegetation on population-level parameters such as animal density and activity center distribution. However, we acknowledge that composites can limit inferences about habitat–wildlife relationships, particularly for practitioners interested in finer scale space-use. In the future, pairing broadscale, population-level patterns of distribution and density with finer scale space-use data (e.g., GPS or VHF-collar data) could clarify drivers of cougar demography and space-use at multiple scales.

Though others have suggested cougars in protected areas occur at higher densities than those in fragmented or developed landscapes, our cougar density estimates in Yosemite (i.e., 0.88 [0.02–2.58] cougars/81 km² or approximately 1.1 cougars/100 km²) were similar estimates throughout the species' distribution, including in California (Allen et al., 2015 [0.67 cougars/100 km²]; Dellinger et al., 2018 [0.87–1.42 cougars/100 km²]), more broadly throughout western North America (Russell et al., 2012 [3.6–5.6 cougars/100 km²]; Davidson et al., 2014 [2.3–5.2 cougars/100 km²]; Loonam et al., 2021 [3.2–6.5 cougars/100 km²]; Ruprecht et al., 2022 [1.7 cougars/100 km²]), and in South America (Paviolo et al., 2009 [0.3–2.9 cougars/100 km²]; Zanón-Martínez et al., 2016 [0.5–9.3 cougars/100 km²]). It can be difficult to robustly compare parameter estimates among studies given density estimates are often derived using varied field (i.e., live-capture vs. noninvasive) and analytical (e.g., SCR vs. space-to-event) methods, result from variable survey efforts, and occur in different study systems and during different time periods (Murphy et al., 2022). However, cougar densities in Yosemite were similar to estimates derived from live-capture data collected from 1983 to 1992 in the nearby Sierra National Forest (0.87 cougars/100 km² in summer, 1.42 cougars/100 km² in winter; Dellinger et al., 2018). Though cougar densities in Yosemite were not significantly greater than estimates elsewhere in their range, recent work has shown that these montane populations in California are important sources of dispersing individuals, and consequent genetic variation, for cougar populations in fragmented or developed areas (Ernest et al., 2000, 2003; Gustafson et al., 2019; Robinson et al., 2008). Thus, the conservation of a broader landscape mosaic that encompasses protected areas and private, state, and federal lands may support the persistence of cougars over time. Further, collaborative efforts to increase habitat restoration and connectivity, regardless of protection status, may elicit positive outcomes for species and ecological communities of conservation interest (Puri et al., 2022).

Cougars occur within a broader community of co-occurring species, including competitors and prey. Cougars occupy a unique trophic role as apex predators that can consume large-bodied prey, but are simultaneously subordinate to other large-bodied carnivores (Elbroch et al., 2015; Elbroch & Kusler, 2018). Consequently, cougars are sensitive to both direct competition, including intraguild predation or kleptoparasitism, and indirect competition, such as exploitation of shared prey species (Brunet et al., 2022; Lendrum et al., 2014; Ruprecht et al., 2022). For example, large carnivores can suppress cougar success by directly predating kittens (Elbroch et al., 2020), usurping and monopolizing carrion (Prugh & Sivy, 2020), and forcing cougars to increase kill rates of prey (Elbroch et al., 2015). Further, few studies have investigated the effects of co-occurring smaller bodied carnivores (e.g., Appendix S1: Table S1) on cougars despite the potential effects of resource competition and scavenging on cougar foraging patterns and density (Prugh & Sivy, 2020; Ruprecht et al., 2022; Smith et al., 2023). Future research efforts to understand the demography and ecology of cougars, or species of conservation concern, will benefit from incorporating information on species interactions to better understand the effects of competition and resources on populations.

There was little co-occurrence between cougars and other species of conservation concern during our study. While cougars are a primary predator of Sierra Nevada bighorn sheep elsewhere (Gammons et al., 2021; Rominger et al., 2004), we did not detect bighorn sheep in Yosemite during our two-year survey period. Further, Sierra Nevada red foxes were only detected at four remote camera stations during the study period, and only four of the 1905 scats collected were identified as Sierra Nevada red fox (Appendix S1: Table S1). All Sierra Nevada red fox detections occurred near the northern perimeter of Yosemite, where cougars were not detected on remote cameras and cougar scat detections were sparse (Appendix S1: Figures S1 and S3). Similarly, only 15 of 1905 scats collected were identified as fishers, and fishers only co-occurred with cougars at one remote camera station. However, fisher occurrence was highest in areas of high predicted cougar density (Figure 3; Appendix S1: Figure S3), which could affect fisher recovery and recolonization in this landscape. Nonetheless, it is difficult to speculate about the role of cougars in influencing the populations of these locally rare species given data limitations. Realistically, the interspecific interactions of cougars with co-occurring carnivores and prey species are likely complex and would benefit from further investigation to clarify community dynamics. Regardless, this work provides foundational information on the abundance and distribution of cougars in Yosemite, which will be useful in continued studies of the wildlife community in the region.

Our work joins a growing body of research that demonstrates the utility of combining noninvasive methods into integrated analyses to estimate the demography and distributions of rare and cryptic species of conservation interest. Contemporary quantitative methods allow applied ecologists to integrate different types of data, including those collected using noninvasive methods, into models that consequently provide more precise or robust estimates of population parameters of interest (e.g., Ferreras et al., 2021; Furnas et al., 2018; Linden et al., 2017; Ruprecht et al., 2021). Such estimates are imperative to inform assessments of population status and trends of species of conservation interest, but can also inform the conservation of co-occurring species negatively affected by the density and distribution of predators or competitors. Rapid environmental change requires the application of novel and effective methodologies, and incorporating effective survey methods with robust analyses will maximize information gleaned from limited resources and support the conservation of species and ecological communities.

AUTHOR CONTRIBUTIONS

Sean M. Matthews, David S. Green, Sarah L. Stock, Heather Mackey, Breeanne K. Jackson, Benjamin N. Sacks, B. Heath Smith, and Jennifer Hartman were involved in study conception and design. Sean M. Matthews, David S. Green, Sarah L. Stock, Heather Mackey, and Breeanne K. Jackson coordinated funding. Benjamin N. Sacks and Stevi L. Vanderzwan completed and supervised laboratory analyses. Sean M. Matthews, David S. Green, Mike A. McDonald, Heather Mackey, Breeanne K. Jackson, Tessa R. Smith, B. Heath Smith, and Jennifer Hartman collected and/or processed data. David S. Green completed data analyses, and David S. Green and Marie E. Martin created figures. All authors contributed to previous versions of the manuscript, and reviewed and approved the final version of the manuscript. Other than Marie E. Martin, David S. Green, Sarah L. Stock, Benjamin N. Sacks, and Sean M. Matthews, authors are listed in alphabetical order of their surname. Sarah L. Stock, Benjamin N. Sacks, and Sean M. Matthews are the senior authors of the manuscript.

ACKNOWLEDGMENTS

This project was funded by U.S. Fish and Wildlife Service, Yosemite Conservancy, and Yosemite National Park. We are grateful to the project partners and many field scientists who were instrumental in collecting data for this effort, including Crystal Barnes, Nicky Bunn, Sylvia Doyle, Dylan Dunn, Ryan Evans, Michelle Gilmore, Dan Gusset, Sarah Hecocks, Corrina Kamoroff, Jake Lammi, Elliot Lozano, Suzie Marlow, Vanessa Roy, Sean Smith, Alex Studd-Soika, and Sarah Sugarman. Phil Johnston trained field staff to identify cougar sign, which considerably improved camera placement in the second year of the study. The authors thank the associate editor, Matthew Mumma, and two anonymous reviewers for their constructive feedback, which improved this manuscript.

CONFLICT OF INTEREST STATEMENT

The authors declare no conflicts of interest.

DATA AVAILABILITY STATEMENT

All data and code are available from Zenodo: https://doi.org/10.5281/zenodo.8197227.