The fear induced by predators on their prey is well known to cause behavioural adjustments by prey that can ripple through food webs. Little is known, however, about the analogous impacts of humans as perceived top predators on the foraging behaviour of carnivores. Here, we investigate the influence of human-induced fear on puma foraging behaviour using location and prey consumption data from 30 tagged individuals living along a gradient of human development. We observed strong behavioural responses by female pumas to human development, whereby their fidelity to kill sites and overall consumption time of prey declined with increasing housing density by 36 and 42%, respectively. Females responded to this decline in prey consumption time by increasing the number of deer they killed in high housing density areas by 36% over what they killed in areas with little residential development. The loss of food from declines in prey consumption time paired with increases in energetic costs associated with killing more prey may have consequences for puma populations, particularly with regard to reproductive success. In addition, greater carcass availability is likely to alter community dynamics by augmenting food resources for scavengers. In light of the extensive and growing impact of habitat modification, our study emphasizes that knowledge of the indirect effects of human activity on animal behaviour is a necessary component in understanding anthropogenic impacts on community dynamics and food web function.

1. Introduction

Anthropogenic disturbance can cause shifts in biotic community dynamics, generally through the loss of native species, introduction of novel species or artificially enhanced populations of native generalists [1]. These changes are characterized most often by quantifying the population size or presence of particular species. However, behaviour-mediated interactions are predicted to have equal or greater impacts on animal populations than purely numerical mechanisms [2–4]. Animal behavioural responses to anthropogenic disturbances have the potential to be cryptic but powerful drivers of ecological change in modified habitats [5].

Large carnivores are widely recognized to be sensitive to human disturbances owing to slow life cycles, large space requirements, direct persecution by humans and avoidance of human-dominated areas, often resulting in their decline or extirpation. Reduced large carnivore density and occupancy in some developed areas has resulted in both mesopredator release [6,7] and overpopulation of primary consumers [8,9]. Yet, some large carnivores do persist in modified landscapes, but alter their behaviour to reduce interactions with humans [10,11]. Owing to the strong regulatory influence large carnivores can exert on their competitors and prey, changes in their behaviour are likely to contribute to whole community responses to anthropogenic disturbances [12].

Hunting and foraging are costly behaviours for carnivores, but are often assumed to be optimized so that individuals gain the most energy (or other limiting nutrients) for the least effort [13]. Animals that choose to expend more energy than what is perceived to be optimal may be responding to risks that increase the long-term pay-off of certain energy-expensive behaviours through decreased chance of non-starvation mortality [14]. Hunting in high-risk habitats can therefore create a risk-foraging trade-off in which animals sacrifice efficient foraging to compensate for increased vigilance and risk avoidance. Giving-up density studies support that an animal's perceived risk is inversely correlated with food intake, suggesting that time spent consuming a food item reflects the fear experienced by forager [15].

Top carnivores generally kill large-bodied prey species, which require a high initial energetic cost during the hunting stage [16], but provide high energy gain during consumption. However, because carnivores are constrained by gut capacity, solitary carnivores can only maximize their caloric yield by repeatedly returning to feed on a prey carcass. In developed habitats, carnivores can be particularly vulnerable to risk-foraging trade-offs because disturbance-induced carcass abandonment can result in food loss owing to scavenging [17] or decomposition [18]. Prey consumption time can therefore be limited by external forces that reduce carnivore access to a carcass. Anthropogenic disturbances can ultimately reduce the net caloric gain carnivores receive from consuming large prey [19] by displacing carnivores from kill sites and decreasing their prey consumption time at kills. If perceived risk increases with human disturbance, the magnitude of human impact should be a predictor of foraging efficiency and consumption time.

We examined puma (Puma concolor) behavioural changes associated with perceived risk at kill sites with increasing housing density levels and investigated the relationship between risk-avoidance behaviours and kill rates in disturbed areas. We hypothesized that disturbance would displace pumas more often in highly developed areas, reducing overall prey consumption time and increasing kill rates. In the long term, we anticipate that more frequent risk-avoidance behaviours will increase puma kill rate and subsequently alter interactions with prey, competitors and scavengers.

2. Material and methods

Our research was conducted in the Santa Cruz Mountains, which lie in the Central Coast region of California. We captured 30 pumas from 2008 to 2013 and fitted them with GPS/radio telemetry collars (IACUC no. WILMC1011 Model GPS Plus 1 or 2 D, Vectronics Aerospace, Berlin, Germany). Collars were programmed to record locations every 4 h, and location data were downloaded remotely via UHF once a month. We used a custom cluster generation algorithm integrated in the Geographical Information Systems program ArcGIS (v. 10; ESRI, 2010) using the programming languages R (v. 2.1.3.1; R Development Core Team, 2010) and Python (v. 2.6; Python Software Foundation, 2010) to identify groups of locations in which each location was within 100 m of the cluster centre and 6 days of another location in the cluster (for full details on the algorithm, see [10]). We field-investigated clusters in reverse chronological order from their time of formation using the most recently downloaded GPS data. Clusters were investigated within 30 days of their first recorded locations. At investigated clusters, we recorded whether a kill was present, and if so, the species, age and sex of the kill (if identifiable).

We constructed a mixed-effects binomial logistic regression model to predict deer kill sites from all generated clusters. We chose to only use deer kills because deer are the preferred prey in our study area. In addition, small prey cannot be predicted accurately from location data in our study area owing to high variability in puma behaviours at small kill sites. The variables used to fit the model were behavioural characteristics associated with clusters, including number of night locations (NIGHT), a binary variable indicating greater than 1 day duration (BINARY), the harmonic mean distance of locations from the cluster centre within the cluster (HMDIST), the proportion of locations occurring at night (P.NIGHT), site fidelity measured by the proportion of points occurring within the cluster over the cluster duration (P.ACTIVE) and the farthest distance travelled during a cluster period (DIST). Total duration of a cluster (DUR) was excluded owing to high correlation (r > 0.7) with variables already used in the model. We allowed for random slopes and intercepts for individual pumas when fitting the best model. We also constructed a truncated model without behavioural variables that we expected to correlate with housing density in order to allow inference on behavioural influences on kill rates. This model excluded P.NIGHT, P.ACTIVE and DIST. We compared the truncated model and best-fit model to assess if we could confidently use the truncated model. We constructed receiver operating characteristic curves for both full and truncated models and calculated the area under the curve (AUC) to ensure that both models were similar in their predictive abilities. The AUC for the full model (AUC = 0.818) and the truncated model (AUC = 0.820) were nearly identical, and both support good discriminative ability [20]. We assigned each generated cluster as a deer kill or not a deer kill by applying the truncated model to all generated clusters.

We first investigated puma behavioural responses at the population level at four levels of housing density within 150 m of predicted kill sites: no housing, rural (greater than 0.0 and up to 0.062 houses per hectare), exurban (greater than 0.062 and up to 1.236 houses per hectare) and suburban (greater than 1.236 and up to 9.884 houses ha−1; [21]). Owing to the discrete nature of housing count data, each housing class was rounded up to the nearest number of houses, resulting in housing classes within 150 m of a cluster to be defined as 0 houses for no housing, one house for rural, two to nine houses for exurban and greater than nine houses for suburban. We used the 150 m buffer because this is the scale of development found to most impact puma hunting in our study area [10]. We constrained the response variables to only include behaviours we expected to be associated with risk aversion, which were narrowed down to P.NIGHT, P.ACTIVE, DUR, DIST and a final measure of prey consumption time (P.C.TIME) which was calculated as P.ACTIVE × DUR. We added the prey consumption time measure because it best reflects the energetic gain an animal receives from a kill. We tested the differences in behaviours at all predicted kill sites in different housing density categories using an ANOVA test and examined pairwise comparisons with Tukey's HSD test for all behaviours found to vary significantly with housing density.

To determine kill rates, we calculated the total number of deer predicted to be killed by each puma from the output of the deer kill prediction model using all generated GPS clusters. We divided the predicted number of kills by the total time each puma was actively collecting GPS data to obtain annual deer kill rates. We investigated sex-specific relationships between kill rates and average values for P.NIGHT, P.ACTIVE, DIST, DUR and P.C.TIME for individual pumas using univariate linear regressions owing to high correlation values (r > 0.7) among variables.

In order to assess the impact of housing on individual pumas, we calculated housing density within each puma home range. Puma home ranges were obtained using a local convex hull (LOCOH) home range estimator, where the 95% isopleth represented the home range boundary [22,23]. All housing points in the Santa Cruz Mountains were manually digitized from high-resolution satellite imagery. We calculated puma home range housing density as the number of houses per km2. We tested for the relationship between individual puma kill rates and home range housing densities using univariate linear regression.

3. Results

(a) Behavioural shifts

Of 703 field-investigated clusters, 208 were classified as deer kill sites. The other remaining clusters included 66 non-deer kills (e.g. raccoons and house cats) and 429 non-kills (often bed sites). Our best-fit binomial logistic regression model to predict deer kills included NIGHT, BINARY, HMDIST, P.NIGHT and P.ACTIVE. The truncated model included NIGHT, BINARY and HMDIST. Neither the random intercept nor a random slope was included in the best-fit full model or the best-fit truncated model. We used the truncated model to predict 1537 deer kills from 8523 generated clusters (figure 1). Figure 1. (a) Study area and predicted kill sites in relation to housing density. Box shows area (b) in which kill site examples (c,d) are shown. Kill site 1 (c) belongs to puma 13F, whose home range is in the top quartile of housing densities among female pumas. Kill site 2 (d) belongs to 28F, whose home range in the bottom quartile of housing densities among female pumas. Grey labels depict location times during the first day at the kill site and black labels depict location times during the second day. Note that at kill site 1, the puma has made a kill close to development but spends the majority of time away from the kill, whereas at kill site 2, which has no nearby development, the puma stays within the vicinity of the kill. (Online version in colour.)

At predicted kill sites, females had lower P.ACTIVE (F = 67.7, ), higher DIST (F = 16.0, p < 0.001) and shorter P.C.TIME (F = 44.2, ) as housing density increased (figure 2; example shown in figure 1). In suburban habitat, female P.ACTIVE was 36% lower, DIST was 31% higher and P.C.TIME was 42% lower than in no housing areas. Both males (F = 19.3, p < 0.001) and females (F = 144.4, ) were more nocturnal (higher P.NIGHT) with increasing housing density at kill sites. Males did not show any responses to housing density concerning time spent at kill sites. Identical analyses using only confirmed kills supported each of the reported trends for predicted kills. Figure 2. Behaviours that vary with housing density at predicted kill sites. Pairwise comparisons from Tukey's HSD tests reported in superscripts, where different letters represent a statistically significant difference. Error bars represent two standard errors from the mean. Sample sizes of housing classes at female kills are: no housing, 719; rural, 83; exurban, 186; suburban, 71. Sample sizes of housing classes at male kill sites are: no housing, 389; rural, 32; exurban, 42; suburban, 15. (Online version in colour.)

(b) Deer kill rates

Male average home range size was 163.0 ± 7.7 s.e. km2 with 15.6 ± 0.8 s.e. houses km−2. Female average home range size was 53.8 ± 2.1 s.e. km2 with 25.5 ± 1.3 s.e. houses km−2. Males had an average deer kill rate of 43.7 deer yr−1, whereas females killed on average 67.3 deer yr−1. Male deer kill rates were not correlated with any of our variables of interest (P.ACTIVE, P.NIGHT, DIST, DUR or P.C.TIME), nor with home range housing density (p = 0.9, r2 = 0.005; figure 3). Conversely, female deer kill rates showed a strong positive and linear correlation with home range housing density within its observed range (p = 0.0003, r2 = 0.745; figure 3). Females with home ranges in the top quartile of housing density killed 36% more deer per year (81.2) than females in the bottom quartile of housing density (59.7). Female kill rates were also negatively correlated with average fidelity to kill sites (P.ACTIVE; p = 0.05, r2 = 0.322). Figure 3. Kill rates of individual male and female pumas in response to housing density per hectare within their 95% LOCOH home range. Regression line and 95% CIs for female kill rates in relation to housing density are shown. Housing density coefficient β 1 = 91.4. (Online version in colour.)

4. Discussion

Our estimate of average male kill rates (43.7 kills yr−1) stayed constant across housing densities and was comparable to previously reported values described by Knopff et al. (35 ungulates yr−1, [24]) and Anderson & Lindzey (47 kills yr−1, [25]). However, female kill rates increased positively, strongly and linearly with housing density. Although female kill rates in lower housing density areas (59.7 kills yr−1) were comparable to previously estimated mule deer kill rates for solitary adult females (52.5 kills yr−1, [25]) and females with kittens (62.4 kills yr−1, [26]; 68.1 kills yr−1, [25]; 57.2 kills yr−1, [24]), female kill rates in the highest quartile of home range housing densities were substantially higher (81.2 kills yr−1). This 36% increase in kill rates between the top and bottom quartiles indicates that development may exert a significant energetic cost associated with hunting behaviour.

Hunting deer requires large energetic investments in the stalking, subduing and killing stages for pumas [16]. We have documented a sizable increase in female kill rates that we expect represents higher energetic costs for females in developed landscapes. Although these costs do not appear to influence adult survival (C. C. Wilmers 2014, unpublished data), impacts on reproductive success possibly make development-interface zones sinks for the puma population. Anecdotally, we have observed that the tagged female living in the most developed habitat in our study area has lost at least three litters in the last 3 years, one of which was confirmed as abandonment (C. C. Wilmers 2014, unpublished data). The three other females living in less developed portions of our study area for which we have also documented at least three dens have had the majority of their litters survive. Although there are many stressors in a developed landscape that might influence kitten survival, we expect that higher energetic costs from increased hunting may contribute to this pattern.

Males did not alter their kill rates or prey consumption time at kills with increasing housing density. Because male life histories are constrained by requirements to defend much larger territories [27], this is perhaps not surprising. We found that male pumas have home ranges that are approximately three times as large as female home ranges on average. Male pumas are also known to spend significantly more time performing scent-marking behaviours than females [28]. We found that males have lower DUR at kills than females by 7.2 h on average (males = 2.86 days ± 0.06 s.e., females = 2.56 days ± 0.05 s.e.; t = 3.83, d.f. = 1000, ), probably owing to their need to patrol and defend their home range boundaries from encroaching males. Therefore, because males already tend to leave their kills early, they may be less influenced by chronic disturbance. In addition, male home ranges are characterized by much lower overall housing densities, indicating that males may exhibit risk-avoidance behaviours at the landscape scale rather than at the kill site scale.

Higher deer kill rates by females in response to increased housing density appear to be driven by a behavioural shift to a lower proportion of time spent at kill sites over the consumption period. Although females did not alter their total duration spent at clusters, their overall prey consumption times declined owing to a lower proportion of time spent at kills, indicating reduced utilization of carcasses at higher housing densities. Other possible explanations for female increase in deer kill rate are not supported by our understanding of puma energetics and reproduction. Deer in our study have lower detection rates in more developed habitats (C. C. Wilmers 2014, unpublished data), therefore variation in deer activity or abundance is unlikely to explain these patterns. In addition, it is unlikely that increased kill rate could be a result of greater reproductive activity in high housing density areas because females in our study area avoid anthropogenic development when denning [10]. We conclude that behavioural risk avoidance is a substantial contributor to female prey consumption time and hence hunting patterns, due to our observation that housing density is associated with decreased prey consumption time. Both food loss and increased movement as a result of these behavioural shifts may contribute to observed increased kill rate in human-modified habitats.

An increase in ungulate carcasses left by female pumas may impact the biotic community by providing additional carrion subsidies to scavengers. By leaving their kills for longer periods of time in more developed areas, female pumas might create greater opportunity for scavenging by mesopredators and birds. Subordinate predators often scavenge kills of apex predators when kills are abandoned or not guarded [29], and carrion can form a large proportion of their diets [30]. Pumas are known to be important sources of food subsidies to mesopredators through carcass abandonment [31]. Our results suggest that mesopredator release may occur not only through the well-documented pathway of apex predator extirpations, but also via behaviour changes in extant apex predators leading to increased food provisioning. The presence of scavengers can exacerbate this pattern by reducing apex predator prey consumption time via food loss [17].

5. Conclusion

The results presented here have bearing on human-modified systems globally. Behavioural responses are often overlooked as ecosystem drivers in modified systems, overshadowed by population declines and extirpations. However, many species are able to persist in developing landscapes, but in an altered behavioural state. Our findings suggest a strong, perceivable impact of observed human-induced behavioural change on species interactions instigated by the presence of development. Risk aversion behaviours that result from anthropogenic disturbances are likely to restructure predator–prey interactions in a variety of contexts, given the large effects risk has been shown to have on foraging across taxa. Behaviour-mediated interactions are powerful forces in biotic systems, often playing an even more impactful role than consumptive interactions. A greater focus on behaviour-mediated effects of habitat alteration can further expand our understanding of community-level processes in human-modified systems.

Ethics statement

The study was carried out under the IACUC approval no. WILMC1011.

Data accessibility

Because pumas are a specially protected species in California, the location data are sensitive and have not been made accessible. However, the data can be available on request by contacting Justine A. Smith at jsmith5@ucsc.edu.

Acknowledgements We thank the California Department of Fish and Wildlife, Cliff Wylie, Dan Tichenor and Troy Collinsworth for their assistance in helping us capture pumas with hounds. We thank landowners who have allowed us to capture pumas and investigate kills. We thank P. Houghtaling for running our field team and V. Yovovich, Y. Shakeri, C. Fust, S. McCain, J. Kermish-Wells, L. Hibbler and dozens of undergraduate field and laboratory assistants for their contributions to data collection, entry and management.

Author contributions

J.A.S. conceived the study, carried out the data analysis and drafted the manuscript. Y.W. assisted in project development and edited the manuscript. C.C.W. conceived and designed the overall study, supervised data collection, advised on data analysis and edited the manuscript.

Funding statement

Funding was provided by NSF grant nos. 0963022 and 1255913 , as well as by the Gordon and Betty Moore Foundation , The Nature Conservancy , Midpeninsula Regional Open Space District , UC Santa Cruz and the Felidae Conservation Fund .

Conflict of interests

We have no competing interests.

Footnotes