Anthropogenic habitat fragmentation and modification can impact wildlife in many ways, from changes in nutrient resource quality and distribution, to increased exposure to pathogens and chemical pollution (Bradley and Altizer 2007, Becker et al. 2015, Murray et al. 2016). As such, anthropogenic stressors broadly have been associated with changes in the health, behavior and reproductive success of numerous vertebrate species (Ditchkoff et al. 2006, Lowry et al. 2013, Ryan and Partan 2014). While some species have been negatively affected by anthropogenic land-use changes, others have exhibited population increases and range expansions in the face of these same perturbations (Ditchkoff et al. 2006, Shochat et al. 2006). For example, large herbivores, such as mule and white-tailed deer, Odocoileus hemionus and O. virginianus, respectively, have successfully colonizedanthropogenic landscapes across North America (Conover 1995). In contrast, large carnivores occupy the highest trophic level, and therefore, exhibit exponentially lower population densities and have greater space requirements than their ungulate prey, making them disproportionately vulnerable to habitat disturbances (Woodroffe and Ginsberg 1998). Hence, anthropogenic habitat modifications may alter predator–prey relationships in non-linear ways, making successful conservation in mixed-use environments substantially more challenging.

The prey abundance hypothesis (Litvaitis et al. 1987, Palomares et al. 2001, Spong 2002, Broomhall et al. 2003) predicts that predator distribution should reflect that of their primary prey. This has been demonstrated at the largest spatial scales for a number of predator–prey relationships (Stoner et al. 2018). For example, ungulate density is a primary driver of tiger Panthera tigris abundance across several reserves in India (Karanth et al. 2004). Additionally both Canadian lynx Lynx canadensis and coyotes Canis latrans are more likely to use habitats where their prey, snowshoe hares Lepus americanus, are most dense (Murray et al. 1994).

However, an increasing number of studies suggest that habitat selection by predators and their prey show greater heterogeneity at progressively finer spatial scales (Hopcraft et al. 2005, Fuller et al. 2007, Kohl et al. 2018). These differences are attributable to stronger selection for micro-site features associated with foraging; that is, herbivores select for habitats where they can find adequate forage within the relative safety of their evolved anti-predation strategies. Conversely, for predators, this means selection of habitats or landscape features where prey can be encountered while facilitating their hunting style (Kauffman et al. 2007). For stalk and ambush predators, this typically means some combination of vegetative cover or rugged terrain.

Pumas Puma concolor are one of the most widely-distributed, large, stalk and ambush predators in the Americas (Pierce and Bleich 2003). At continental scales, puma distribution reflects that of their ungulate prey and historic eradication campaigns. However, they are also adaptable and resilient as evidenced by their reappearance in parts of the American Midwest (Thompson and Jenks 2010, LaRue et al. 2012) and their presence on the periphery of numerous urban areas (Vickers et al. 2015). Pumas occupy diverse ecosystems including tropical forests (Aranda and Sánchez-Cordero 1996, Garla et al. 2001), deserts (Choate et al. 2018), dry forests (Núñez et al. 2000, Lewis et al. 2015) and wetlands (Rodgers and Pienaar 2018). Despite differences in climate, prey composition and plant communities in these ecosystems, the literature on puma habitat selection has consistently identified vegetative or topographic cover as the prevailing factors predicting suitable habitat (Laing and Lindzey 1991, Dickson and Beier 2002, Burdett et al. 2010, Ordenana et al. 2010), which is thought to facilitate their stalk and pounce predation style.

Pumas do not typically reside in core urban areas. Indeed, evidence suggests that although pumas will use near-urban areas for travel and hunting (Ordenana et al. 2010, Wilmers et al. 2013, Lewis et al. 2015, Benson et al. 2016), they tend to avoid direct contact with humans (Crooks 2002, Ordenana et al. 2010, Kertson et al. 2011b, Wang et al. 2015). Yet, peri-urban ecosystems may be infrequent but valuable hunting grounds for pumas because ungulate prey are often drawn to anthropogenic landcover types that pumas tend to under-utilize (Van Dyke et al. 1986, Dickson and Beier 2002, Burdett et al. 2010, DeAngelo et al. 2011, Kertson et al. 2011a). For example, parks, golf courses and residential areas attract deer to anthropogenic food resources, i.e. irrigated, early successional plant communities (DeStefano and DeGraaf 2003). Additionally, deer may frequent these sites because they may offer refuge from some native predators (Šálek et al. 2015).

Taken together it is unclear whether puma occupancy is better explained by the prevalence of undisturbed wildlands and stalking cover, or by prey abundance, especially given that anthropogenically disturbed habitats may harbor relatively higher prey densities (Coon et al. 2019). Hence our goal was to identify variables predicting seasonal mule deer and puma habitat occupancy in a highly fragmented landscape in the green and dry seasons. Based on previous research, we hypothesized that deer would occupy a broader range of habitats than pumas, including anthropogenically developed areas, whereas puma occupancy would largely be associated with natural vegetative cover. The prey-abundance hypothesis (Litvaitis et al. 1987, Palomares et al. 2001, Spong 2002, Broomhall et al. 2003) predicts that, as obligate carnivores and habitat generalists, puma distribution and habitat selection should also display a tight correlation with mule deer, their primary prey (Iriarte et al. 1990, Pierce and Bleich 2003, Allen et al. 2015). Likewise, we expected that occupied deer habitat would be broader in the dry season when they would likely have to range farther in order to find sufficient forage and then pumas, following the deer, would also exhibit this pattern.

Material and methods

Camera trapping

Cameras were set in a variety of open spaces throughout the study area in locations that were potentially suitable for medium to large wildlife such as pumas and black-tailed mule deer Odocoileus hemionus columbiana (hereafter: ‘deer'; Penrod et al. 2013, Allen et al. 2015). Though pumas avoid major highways, they do use fire roads and trails for movement (Van Dyke et al. 1986, Dickson et al. 2005, Nichols 2017). Likewise, deer are known to travel along roads and trails for movement and foraging (Ager et al. 2003, Bird 2019). For this reason, cameras were set on dirt and fire roads, hiking routes and game trails in areas with identifiable wildlife sign (scat, tracks, biologist sighting), without consistent vehicular traffic and ideally near trail junctions or water sources. All camera monitoring was passive – no exogenous (scent, audio) lures were used.

Figure 1.

Map of study region which includes camera placements (red dots) in the San Francisco Bay Area, California, USA. Grid cells with at least one camera are outlined.

Cameras were placed between 90 and 120cm above ground and angled down trails and unpaved roads to maximize the possibility of detecting medium to large wildlife. Most cameras were Bushnell Trophy or Natureview models (∼95%), though we also used some Reconyx Hyperfire and Moultrie M-class cameras. All cameras were set to take three photos after each motion trigger. Trigger speeds and delays between photos were less than 1s for all camera models. Sensor sensitivities were set on auto or low depending on the density of nearby vegetation.

Camera traps were set at a total of 416 locations and monitored between June 2013 and May 2017. In most cases, there was a single camera within a 1 km² grid cell at a given time (Supplementary material Appendix 1 Data A1), however, over time, camera sites were sometimes less than 1 km apart or camera malfunctions may have necessitated taking the camera down for extended periods of time. To avoid pseudo-replication bias between close cameras we overlaid a 1km² grid over our study area and grouped all data collected within a single grid cell (Fig. 1, 2; Sun et al. 2014, Wang et al. 2015, Macdonald et al. 2019). After grouping, we had 181 grid cells with at least one active camera during some part of the study period.

Figure 2.

Map indicating the 180 study grid cells of the 18 463 total grid cells in the study. Grid cells are labeled as having detected puma, deer, both or neither. Non-study cells that are majority covered by water bodies, agriculture, barren lands or human development are masked (dark grey) due to low sampling effort in these areas.

Both deer and pumas move approximately 5km per day (Ager et al. 2003). As such, it is possible that the same individual of either species may have been detected within multiple grid cells during a survey period (one week). This means that sites were not closed which violates an assumption of strict occupancy modelling (MacKenzie et al. 2017). However, for these models we are considering ‘occupancy' as ‘habitat use’ as suggested by Burton et al. (2012) which permits non-closed sites.

Data preparation

Camera trap photos were categorized using Camerabase software (The Botanical Research Institute of Texas) by trained volunteers and spot checked by staff. Presence–absence (binary) data was organized into matrices by grid cell and week of the study and split into two seasons in an attempt to minimize missing camera data and facilitate model convergence. Of the 181 sites, 162 were used for models of the green season habitat use and 141 were used for models of the dry season habitat use for each species.

Using data from all 181 sites, covariates were centered and scaled by 1 SD (i.e. z-score transformed). A number of variables were tested in our initial models including elevation, slope, distance to development, topographic position index, topographic ruggedness index, vector ruggedness measure and normalized difference vegetation index (NDVI; Yang et al. 2013) as well as several second-order variables. These variables were not included in the final global models because they were either colinear with other variables (Pearson's correlation; Supplementary material Appendix 1 Data A2) or prevented the models from converging.

Global models for puma and deer included three common variables hypothesized to affect both species' habitat choice (Table 1): distance to the nearest paved and trafficked road as a measure of distance to human disturbance (Ordenana et al. 2010, Benson et al. 2016, Wang et al. 2017), paved road density as a measure of the density of human disturbance (Angelieri et al. 2016), and forest cover – includ evergreen, mixed and deciduous forest (Nicholson et al. 1997, Dickson and Beier 2002, LaRue and Nielsen 2011). The global deer models also included shrub cover and distance to water (Nicholson et al. 1997), whereas global puma models included a simple estimate of deer abundance which is described in Table 1.

Table 1.

Description of variables used in occupancy models for pumas and deer. ‘Type’ identifies whether the variable was used to model occupancy (O) or detection (D). All continuous, occupancy variables were z-score transformed and all spatial variables were calculated in ArcGIS for the 1 km² grid cells. NASS refers to the National Agricultural Statistics Service which is a division of the US Department of Agriculture. For all data from NASS, a majority filter was applied to the Cropscape/Cropland layer. CLN refers to the Conservation Lands Network which is led by the Bay Area Open Space Council. Both mapping data sets (NASS and CLN) are at a resolution of 30 m.

For paved road density and distance to fresh water, we tested for bias in camera placement by comparing the average distance within each cell with the average distance for each camera within each cell using zonal statistic in ArcGIS. Two-sample T-tests for unequal variances confirmed that there were no differences between camera placements and average locations within a grid cell (n = 204; road density: t₄₀₆ = –0.25, p = 0.81; distance from water: t₄₀₆ = 0.16, p = 0.87). To match camera trap data, all spatial variables were summarized at a 1-km² resolution.

Occupancy maps

Using the best models for pumas and deer and data from the 181 grid cells, we generated maps depicting the probability of puma and deer occupancy across our study area in each season in the remaining 18 283 grid cells in Fig. 2. To do so, the spatial variables were scaled to the 1 km² grid in ArcGIS and standardized using the z-score transformation described for the occupancy models. Using coefficients from the occupancy models, we applied the following logistic formula as a transformation of linear predictors into a probability:

The results were displayed on the 1-km² grid using the natural break (Jenks) classification. This method was chosen over the ‘equal interval’ classification method in order to better cluster cells with similar values, and differentiate cells with relatively large differences. For better legibility, a smoother surface was overlaid on top of the grid. This surface was obtained by turning the 1-km² grid into points on which the Kriging interpolation method was applied. This method generates a prediction surface based on a set of points.

Some landcover types were not adequately sampled under the assumption that they would not accommodate our focal species, which violates a requirement of occupancy modeling (Burton et al. 2015, MacKenzie et al. 2017). For this reason, cells containing water bodies, agriculture, barren lands, and high, medium and low development were masked in all maps. To do so, we applied the zonal statistic tool using the ‘majority’ statistic to identify the most frequently occurring value in each cell (Fig. 2).

Results

We used data from 68 299 trap nights from 181 sites which were organized into weekly survey periods by site. The average site had 54 weekly survey periods (range: 3–470; SD: 65). Pumas were detected at least once in 68 sites (Fig. 2; green season: 59; dry season: 38). These included 432 surveys with puma detections (green season: 247; dry season: 185), and 6483 surveys without (green season: 3595; dry season: 2888). Deer were detected at 166 sites (Fig. 2; green season: 163; dry season: 142). In total, 2592 surveys had deer detections (green season: 1293; dry season: 1299), and 4505 did not (green season: 2539; dry season: 1866).

Seasonal puma occupancy was predicted by three equally-well-performing models (AICc within two points). Estimates and significance for each of these models plus the global model is reported in Table 2. We used estimates from the top model (lowest AICc) for maps because the three models were, for the most part, similar in regards to which variables were included and their beta estimates. Variables predicting occupancy of deer in the green and dry seasons was less consistent, with seven and nine models within two AICc points of the top model for each season respectively. For this reason, we chose to average these models and report the estimates of the averaged model in conjunction with the top and global models in Table 2. Likewise, estimates from the averaged model are used for maps.

Extrapolated occupancy maps

We used the occupancy variable estimates from the top puma models and the averaged deer models for each season. We defined high-occupancy as the top 33% of calculated occupancy probability. In the green season, high-occupancy areas spanned 2077.96km² (11.4% of the total mapped area) for pumas and 4548km² (25.1% of the total) for deer (Fig. 3). In the dry season, high-occupancy areas spanned 2833.84km² (15.6% of the total mapped area) for pumas and 6029km² (33.2% of the total) for deer (Fig. 3).

Table 2.

Model estimates. There were only three best models for pumas (within two AIC points) so all three models are reported. Deer, on the other hand, had several models within two AIC points of the best model, so the best model plus a fully averaged model are reported for each season. Estimates from the global models are shown for both pumas and deer. Abbreviations are defined in Table 1 and as follows: Col = colonization variables, Ext = extinction variables and int = intercept. The significance notation for these variables corresponds to the least significant p-value as compared to the default categories: p = 0.01–0.05*, p = 0.001–0.01** and p ≤ 0.001***.

Figure 3.

Extrapolated occupancy maps for pumas (top) and deer (bottom) in the green (Dec–May; left) and dry (June–Nov; right) seasons. The best model of puma occupancy was predicted by forest cover, distance to the nearest road and road density, and negatively correlated with deer abundance in the green season and with forest cover and distance to the nearest road in the dry season. The average of the best deer occupancy models was (non-significantly) correlated with shrub cover, distance to the nearest road and road density in both the green and dry seasons, as well as with distance to the nearest fresh water source in the dry season.

Discussion

Our first objective was to identify landscape features most closely associated with puma and mule deer occupancy. Deer occupancy was best predicted by shrub cover in the green season, though the positive correlation was not significant and the model's predictive power was low. Likewise, none of our variables performed well when predicting deer occupancy in the dry season, perhaps because deer in this region are ubiquitous in the habitats we surveyed. Puma occupancy during the green season was positively correlated with forest cover, distance to the nearest road and road density. Ironically, puma occupancy was negatively correlated with deer abundance in the green season. In the dry season, only forest cover and distance to nearest road predicted puma occupancy. Our second objective was to map the results of our occupancy analyses and visually predict habitat use of pumas and deer. Our results suggest that puma habitat is restricted to between 11 and 16% of our study area whereas deer habitat covers at least a quarter and up to a third of the study area, depending on season.

We expected that high probability puma and deer habitats would demonstrate a near one-to-one correspondence or that puma occupancy would be predicted by deer abundance, but we found that less than half of predicted puma habitat overlaps with predicted mule deer habitat in the San Francisco Bay Area. Puma occupancy was not significantly predicted by deer abundance in either seasonal model and, in our green season models, puma occupancy was actually negatively correlated with deer abundance. Simply put, puma habitat selection was not strongly predicted by deer. Rather than choosing habitats with abundant deer, pumas appeared to prefer habitats with adequate forest cover which likely facilitates pumas' stalking hunting style (Dickson and Beier 2002, Burdett et al. 2010, Kertson et al. 2011a, Wilmers et al. 2013, Angelieri et al. 2016, Benson et al. 2016). Additionally, given the positive correlation between puma occupancy and distance to road in both seasons, one other aspect of puma habitat choice is likely direct avoidance of humans (which is likely also positively correlated with deer abundance). Puma avoidance of humans is consistent with many other puma studies conducted in urban-interface settings (Crooks 2002, Ordenana et al. 2010, Kertson et al. 2011b, Wang et al. 2015).

Taken altogether, our results are inconsistent with the prey-abundance hypothesis which suggests that predator habitats should reflect the density of their primary prey. On the contrary, our data supports what Hopcraft et al. (2005) call the ‘ambush-habitat hypothesis,’ which proposes that predators choose habitats not based on prey density or encounter rates per se, but on other factors such as prey vulnerability and minimization of bodily risk. Pumas are not cooperative hunters, and so the element of surprise is critical to neutralize prey quickly to avoid injuries. In fact, puma kill attempts made > 25m from ungulates are rarely successful (Wilson 1984, Holmes and Laundré 2006, Nichols 2017), and injuries sustained while hunting can be a source of puma mortality (Elbroch et al. 2016). One way that pumas can maximize the element of surprise is to choose habitats where prey occur in smaller group sizes or at lower densities since large aggregations of herding animals tend to exhibit higher vigilance rates (Hunter and Skinner 1998, Lashley et al. 2014). For this reason, habitats with higher deer densities may be associated with lower hunting success rates (Kohl et al. 2018).

Previous work on pumas also supports the ambush-habitat hypothesis without testing it directly. Pierce et al. (2004) reported that puma kills of mule deer in California were more prevalent in substandard deer foraging habitats, suggesting that deer were less vulnerable to puma predation in their optimal foraging habitats. Nevertheless, prey-abundance and ambush-habitat hypotheses are not mutually exclusive in that they both predict that predators should occupy habitats that optimize the tradeoff between energy expenditure and resource acquisition (Balme et al. 2007). Furthermore, different individuals or populations of the same species may tend toward either of the two tactics. For example, studies in bobcats Lynx rufus and African lions Panthera leo give support for both the prey-abundance hypothesis (Litvaitis et al. 1987, Spong 2002) and the ambush-habitat hypothesis (Hopcraft et al. 2005, Davidson et al. 2012, Lewis et al. 2015) depending on the population.

Although pumas eat small prey and may do so more frequently in anthropogenically disturbed habitats (Kertson et al. 2011b, Smith et al. 2016, Coon et al. 2019), we focused on mule deer habitat use because deer comprise the majority of prey biomass for pumas in our study area and throughout western USA (Iriarte et al. 1990). Despite individual and seasonal variation, mule deer typically account for > 80% of the average puma's diet in our study area (Iriarte et al. 1990, Allen et al. 2014). Thus, even if deer abundance is not a good predictor of puma habitat use at the scales measured here, deer presence still is required for puma survival in this region, and cannot be discounted in any conservation planning for this species.

During the green season, road density is significantly, positively correlated with puma occupancy, but not strongly correlated with deer occupancy. Consequently, there is something attractive to pumas about this type of anthropogenic development, other than deer, during this period. We suspect that these habitats have more potential non-live-ungulate food resources, such as roadkill or urban-tolerant prey species (e.g. rats, squirrels), which have been shown to have a positive influence on body condition in pumas in near-urban environments (Bauer et al. 2005, Coon et al. 2019). This is further supported by the fact that we only find this correlation during the green season which is when reproduction, dispersal and foraging is at its peak for many smaller, alternative prey species (Lowry et al. 2013).

The increased availability of food resources in the green season may also explain why preferred puma habitat contracts during this season as compared to the dry season – pumas do not have to range as far in order to find available prey (Stoner et al. 2018). When we mapped predicted occupancy, we see the same pattern in deer: an expanded range in the dry season and contracted range in the green season. This is consistent with previous studies of urban deer, suggesting that they do not have to range as far for forage in the green season when food resources are higher quality and more abundant (Grund et al. 2002).

Unlike pumas, deer occupancy was not significantly correlated with any single predictor in either season even though we tested all variables typically associated with mule deer habitat use. Importantly, deer are ubiquitous throughout our study area, being commonly detected in all available landcover types in the study area which may explain our inability to model occupancy at the 1 km² scale. Our predicted occupancy maps corroborate this hypothesis in that much of our habitat map is moderately to highly likely to be occupied by deer. This pattern came as some surprise since some studies suggest that deer avoid roads and human infrastructure (Rost and Bailey 1979). However, these investigations largely come from remote habitats where deer are less acclimated to human disturbance, or may easily avoid roads. In the highly developed landscape we studied here, deer may be equally attracted to traditional shrublands as they are to disturbed areas, which often contain early successional forage and irrigated landscaping (Ng et al. 2008).

In addition to predicting occupancy, our models also suggest important variables influencing detection probability. Understanding factors affecting detection is critical for long-term monitoring. Here we found that year was an important factor in detection for pumas in the dry season and for deer in both seasons, whereas camera effort was only important for deer. The utility of year as a detection variable likely stems from variation in camera placement over time. Over the course of the study we acquired more cameras and began sampling more diverse habitats as we gained access to more properties. Cameras shifted in location to become more spaced out (i.e. less effort within single grid cells). The importance of effort for deer but not for pumas is interesting in that we expected a positive correlation between effort and puma detection given that pumas occur at substantially lower densities. However, despite the notable differences in density between the two species, species-specific differences in movements likely affect detection rates on static occupancy sampling grids, potentially compensating for the need for increased effort.

Some additional caveats should be considered when drawing conclusions from our study. These include the inability to distinguish between resident and dispersing animals (i.e. how long animals spent in the habitat where they were detected), or between sexes, which can impact predictions of occupancy (Beyer et al. 2010). However, we were interested in habitat use, including use as movement corridors, more so than strict occupancy and, as such, lack of individual identification should not affect our interpretation. Additionally, though we did not find that pumas select habitats based on deer abundance, that does not mean that deer abundance does not affect puma feeding habits. For example, leopards Panthera pardus preferentially consume wild prey in central India but in a scenario where wild prey decline (due to habitat loss, hunting, etc.), models suggest a dramatic increase in leopard consumption of livestock (Puri et al. 2020). This, of course, can have important implications for wildlife management.

Puma management implications

Habitat loss and fragmentation represent the greatest threats to large mammals, worldwide. Conservation of pumas in California is a primary example of this management challenge (Ernest et al. 2003) as freeways and housing developments are overlaid on extensive natural fragmentation stemming from inland waterways and widespread agricultural lands with little cover. Our findings suggest that the simple presence of ungulate prey will not be adequate to insure puma conservation. Although prey density is a widely accepted correlate of habitat quality for many carnivores, our results suggest that structural elements of habitat interact with prey abundance to influence predation success, which may influence life-history traits that ultimately determine population viability. In short, given the social and financial constraints on land conservation in mixed-use landscapes, management efforts for stalk and ambush predators such as pumas should prioritize habitat patches and corridors that maximize prey vulnerability over absolute abundance.