Ecological effects of fear: How spatiotemporal heterogeneity in predation risk influences mule deer access to forage in a sky‐island system

Abstract Forage availability and predation risk interact to affect habitat use of ungulates across many biomes. Within sky‐island habitats of the Mojave Desert, increased availability of diverse forage and cover may provide ungulates with unique opportunities to extend nutrient uptake and/or to mitigate predation risk. We addressed whether habitat use and foraging patterns of female mule deer (Odocoileus hemionus) responded to normalized difference vegetation index (NDVI), NDVI rate of change (green‐up), or the occurrence of cougars (Puma concolor). Female mule deer used available green‐up primarily in spring, although growing vegetation was available during other seasons. Mule deer and cougar shared similar habitat all year, and our models indicated cougars had a consistent, negative effect on mule deer access to growing vegetation, particularly in summer when cougar occurrence became concentrated at higher elevations. A seemingly late parturition date coincided with diminishing NDVI during the lactation period. Sky‐island populations, rarely studied, provide the opportunity to determine how mule deer respond to growing foliage along steep elevation and vegetation gradients when trapped with their predators and seasonally limited by aridity. Our findings indicate that fear of predation may restrict access to the forage resources found in sky islands.

In the Basin and Range Province of southwestern North America, mountain ranges dominated by conifers are widely separated by desert and shrub steppe regions to form "sky-island" landscapes, where a suite of species is postulated to have isolated for millennia (Brown, 1971). The relatively greater ecological diversity of sky islands is driven primarily by the elevation gradient that provides increasing precipitation, plant and animal diversity, complex terrain, and refuge from extreme environments (McCormack, Huang, & Knowles, 2009). Although sky islands may serve as refuges of high-quality habitat during periods of harsh conditions in the surrounding desert, local ungulate populations are contained with their predators and limited by high temperatures during the summer months that coincide with fawning.
Using the normalized difference vegetation index (NDVI) of primary production across different time periods and vegetation types, it is possible to quantify the connections between short-term increases in NDVI or "green-up," animal movements and habitat use, and the timing of reproduction (Hamel, Garel, Festa-Bianchet, Gaillard, & Coté, 2009;Pettorelli et al., 2011). Species such as mule deer (Odocoileus hemionus) are highly dependent on succulent forage to meet energy requirements (Hurley et al., 2014). They are also responsive to NDVI and green-up in terms of large-scale movements , habitat use (Marshal, Bleich, Krausman, Reed, & Andrew, 2006), and the timing of reproduction (Stoner, Sexton, Nagol, Bernales, & Edwards, 2016). We examined trade-offs between the intensity of cougar use and use of growing forage within sky-island habitats of the Mojave Desert. We further examined seasonal changes in mule deer habitat selection with respect to green-up, changes in the intensity of cougar use, and the timing of reproduction.
Mule deer populations typically pursue higher quality forage during green-up period(s) (Hebblewhite, Merrill, & McDermid, 2008) and then greatly reduce their foraging and energy expenditure during times of low vegetation growth (Alldredge, Lipscomb, & Whicker, 1974). In the temperate regions of the species' range, winter is the time of low forage quality, while green-up occurs during spring and summer (Smith, Krausman, & Painter, 2015). Parturition occurs in early to mid-summer in many populations (Butler et al., 2009;Smith et al., 2015) and appears to be related to the health of the female before gestation (Bowyer, 1991;Haskell et al., 2008).
This gives pregnant mule deer many months to improve nutritional condition both before parturition and during lactation, potentially increasing fawn survival (Lendrum, Anderson, Monteith, Jenks, & Bowyer, 2014;Monteith et al., 2014). In the Mojave Desert, however, green-up occurs primarily in the winter and spring, and summer is the time of plant desiccation and therefore low forage quality and availability (McKee et al., 2015). This suggests that deer which forage in greening areas earlier in the year would have a selective advantage not seen in more temperate regions. Unless desert mule deer are able to benefit from winter green-up, the time available to put on weight before parturition and during lactation may be significantly restricted.
Cougars (Puma concolor) are a major predator of mule deer within desert ecosystems, and fawns may represent a significant proportion of prey (Logan & Sweanor, 2001). In the isolating circumstances of a sky island, mule deer can neither migrate nor emigrate to safety, essentially making them a resident population under consistent predation risk (Hebblewhite & Merrill, 2009). Mule deer typically reduce risk of predation by early detection and outrunning predators, a strategy that should favor use of open shrub rather than forested habitats available on sky islands (Schmidt & Kuijper, 2015). Females, however, may move to the greater cover of forested areas during parturition to decrease the visibility of fawns and provide thermal protection (Marshal et al., 2006), potentially trading greater forage quality for reduced predation risk and/or heat stress.
In this study, we propose that sky-island mule deer modify their foraging and fawn-hiding behaviors to correspond with the timing and availability of plant resources, and that these adjustments can be measured by an analysis of seasonal resource selection functions (RSF).
Specifically, we predict female mule deer will demonstrate a pattern of foraging that includes the use of green-up areas during winter, and that seasonal changes in habitat variable coefficients will reflect a early summer parturition. Additionally, we predict that birth dates of sky-island mule deer will occur significantly earlier than in populations in more temperate climates, where plant resources are available later in the year. We further propose changes in intensity of cougar use will measurably modify mule deer foraging behavior. Specifically, we predict mule deer will reduce use of growing areas as intensity of cougar use increases especially during summer parturition. In a mule deer RSF model, the trade-off between green-up and intensity of cougar use will generate a negative interaction term and the positive relationship between increasing forage quality and female mule deer occurrence will diminish as cougar occurrence increases.

| Study area
The Desert National Wildlife Refuge (DNWR) of southern Nevada encompasses 6,540 km 2 (Figure 1). While the forested, montane environments form the "islands" in a matrix of lowland desert, the complex of desert and ranges creates the sky-island landscape. Figure 1 depicts the extent of forested areas on the refuge; however, it was beyond the scope of this paper to strictly delineate the boundaries of sky-island effects. Since both cougar and deer movements extended beyond forested environments, resulting in landscape-level effects up to the high elevation range, we defined the entire study area as being within the sky island (McCormack, et al., 2009). Precipitation is highly variable (30-year mean = 11.8 cm, SD = 7.8, range: 1.7-37.5) (source: Western Regional Climate Center). Desert shrub is dominant from 800 to 1,800 m and characterized by creosote (Larrea tridentata), blackbrush (Coleogyne ramosissima), and saltbush (Atriplex spp.) associations which often include Mohave yucca (Yucca schidigera) and Joshua tree (Y. brevifolia

| Predictor variables
NDVI: We used NDVI to estimate vegetation biomass, as an index of forage abundance. NDVI is a satellite-derived measure of the difference between visible (red) light absorbed by vegetation and the reflectance emitted in the near infra-red spectrum and is proportional to standing biomass (Pettorelli et al., 2011). Reflectance estimates were derived from the Moderate-resolution Imaging Spectroradiometer sensor data available from the National Aeronautics and Space Administration. NDVI was averaged within a 500 × 500 m 2 area (Stoner et al., 2016). NDVI rate of change (NDVIR): We used the rate of change of NDVI between two time periods, two weeks apart as a measure of changes in forage availability (Marshall et al., 2006). Vegetation type: A two-level categorical variable consisting of tree-or shrub-covered F I G U R E 1 Terrain, vegetation, and burn map and location of Desert National Wildlife Refuge areas was derived from the Southwestern Regional Gap Analysis Project. Shrub environments potentially provided more stalking cover for cougars as well as more forage than areas dominated by sparse tree cover. Therefore, we modeled both NDVIR and vegetation type to address the potential relationship between these variables. Season: We divided our data set according to three distinct seasonal time periods, based on periods biologically relevant for deer behavior, before analyses. Although there are benefits to running a single model with a season covariate, we believed the clarity of this approach reduced potential confusion in interpreting multiple three-way interactions necessary in a single-model approach. February-May (spring) was the period of greatest vegetation growth; June-September (summer) included high temperatures, plant desiccation, and greatest reliance on water sources; and October-January (winter) was the period of greatest precipitation and winter green-up. We used the interaction terms (cougar rsf [crsf defined below] × NDVI) and (crsf × NDVIR) to measure the effects of our indexes of vegetation biomass and changes in forage availability, respectively, on mule deer's response to cougar intensity of use. Distance to closest water source: Water sources were derived from USGS and Nevada Department of Wildlife data and ground-checked for occurrence of a year-round water supply.
Distance to closest burned area: Areas burned, regardless of vegetation type, were delineated using data measured by LandSat sensors  (Table 1). PC1 accounted for 75.8% of the variation in a PCA of these four variables with positive PC1 values representing the common covariance of increasing slope, ruggedness (VRM), and elevation and decreasing viewshed. Curvature variables were uncorrelated with the other terrain variables and therefore were not included in the PCA. PC1 represents the expected correlations between topographic variables which arise due to geomorphological processes that generate relatively scale-independent relationships between slope, valleys, and drainage areas (Montgomery & Dietrich, 1992). Increases in elevation from the base of a mountain are typically associated with increased slope and increased erosion of hills, gullies, and valleys which culminate near the tops of drainages at all scales (Istanbulluoglu, Yetemen, Vivoni, Gutiérrez-Jurado, & Bras, 2008) and cause increased vector measured ruggedness. Ruggedness is necessarily negatively related to viewshed due to topographic obstruction of view. To decouple the intercorrelated terrain variables but retain the ability to interpret different components of topography, we borrow a method from morphometrics and allometry (Bookstein, 1989) and remove the effects of the covarying variables by generating measures of relative terrain shape in their respective units with the residuals from regressions of each terrain variable against PC1. Slope, VRM, and Viewshed Residuals: Having accounted for 75.8% of the variance with PC1, we quantified residual slope, VRM, and viewshed relative to the common covariance of PC1 with separate linear regressions of each variable against PC1. These residual values represent the relative slope, VRM, or viewshed, which exceeds, either positively or negatively, that which is accounted for by PC1. Including PC1 and all four residual terrain variables in the linear model would create a perfect representation error. Therefore, residual elevation was excluded based on its relatively low PC loading and the difficulty of interpreting deviations in elevation.

| Mule deer captures
Capture techniques for all animals were performed under guidelines for the use of live animals (Sikes & Gannan, 2011) and followed pro-

| Cougar captures
Cougar captures occurred from October 2010 through May 2012. We spread trapping efforts across the study area and supplemented trapping with hound/tracking surveys for the more remote regions. We employed any of the following three methods depending on animal safety (e.g., ambient temperatures, terrain) and logistical considerations: Hounds pursued the cougar until it sought refuge in a tree or cliff (Hemker, Lindzey, & Ackerman, 1984); foot-hold snares were placed at kill sites or along cougar travel routes (Logan, Sweanor, Smith, & Hornocker, 1999); or, when access enabled transport, we set cage traps baited with deer carcasses (Bauer, Logan, Sweanor, & Boyce, 2005). Cougars were immobilized with a combination of 2 mg/kg ketamine and 0.2 mg/kg medetomidine (Kreeger, Raath, & Arnemo, 2002) and held for less than one-half hour while being fitted with a GPS-satellite collar (Telonics Gen4). Collars were placed on adult (>3 years) or subadults (1.5-3 years old) only. Collars acquired a location once every four hours starting at 12 a.m. for 24 months and were equipped with an automatic release mechanism and mortality sensor.

| Cougar intensity of habitat use or cougar RSF
We randomly selected 6,050 locations from a total of 7,672 collected over the two-year period and used a 1:1 ratio of used versus available points in a binary logistic model to determine the cougar resource selection function (RSF) or intensity of use within each season. A priori candidate models (for both cougar and deer, below) were derived based on ecological assumptions of terrain structure (hiding/es-  Hanley & McNeil, 1982).
Home ranges were derived for each cougar each season using a 95% Gaussian kernel density estimator and smooth cross-validation to estimate bandwidth. We combined individual cougar home ranges by season to determine the perimeter outlines of the seasonal home ranges and used these outlines to constrain random points. Although bighorn sheep (Ovis canadensis nelsoni) are a commonly exploited alternative prey species in this system, our preliminary analysis indicated bighorn occurrence did not contribute strongly enough for inclusion as a variable in the candidate models.

| Estimation of fawn birth dates
We placed cameras (Bushnell) at 35 water sources throughout the study area. Of these, 10 consistently photographed mule deer. From these 10, we selected two cameras within each of the four following habitat types along an elevation gradient: desert shrub, Joshua tree, pinyon-juniper, and ponderosa pine associations. Fawn birth dates were estimated from appearance of spotting in coats and nursing behavior (Anderson & Wallmo, 1984).

| Biomass and green-up
Mean biomass (NDVI) had the lowest value in January, increased to a peak in April and May, then declined continuously to a low in the winter months (Figure 2). The mean two-week change in NDVI or green-up (NDVIR) was greatest from late winter to early spring with a peak in March. Subsequently, mean NDVIR declined and became negative with plants desiccating in May and June, increased again with modest green-up in late summer, then declined in October. The green-up during most summer months, though positive, was relatively low (Figure 3).

| Cougar RSF or intensity of use model selection
We captured 4 female cougars (2 adults, 2 subadults) and 1 sub-  Table 2). Within the winter season, we found two models contained reasonable evidence (delta AIC < 2.0) for explaining habitat selection. The first contained all variables except curvature planform, and the second contained all variables except NDVI. We chose the best model based on performance (AUC = 0.881, 95% CI = 0.869-0.893, likelihood χ 2 = 2,254.5, p < 0.001. Deviance = 0.818, Pearson χ 2 = 1.12). Although averaging across reasonable models when using AIC is an option, it would preclude our ability to critically evaluate the relative contribution of individual coefficients from logistic regression (Cade, 2015). We used the highest-rated models to estimate cougar RSF values (cougar intensity of use) within each respective season.

| Mule deer habitat selection
Female mule deer used a wide variety of terrain within the study area ( Figure 5). During spring, mule deer responded positively though weakly to biomass (NDVI), and positively to green-up (NDVIR). Mule deer were associated with an increasing likelihood of cougar occurrence, suggesting similar habitat use (Table 3).
Interestingly, cougar intensity of use interacted strongly and neg- Distance to water was similar to springtime occurrence. Cougar habitat use again followed the shift toward shrub-covered areas, with risk greater in the shrub relative to tree-covered regions (Table 4).

| Estimation of fawn birth dates
We viewed over 28,700 photos over 36 months incorporating the study period (including 3 months before and after). From the estimated ages of fawns viewed, we determined birth dates within the relatively lower elevation desert shrub and Joshua tree associations to be from late May to early June, and within the higher elevation pinyon-juniper and ponderosa pine associations to be from early June to early July.

| D ISCUSS I ON
Based on mule deer association with NDVI and NDVI rate of change, we found little support for the hypothesis that female mule deer follow available green-up during the winter months, but instead found deer used this strategy primarily during the spring. In this sky-island population, a narrow four-month period stands in contrast to other western populations of mule deer, for which foraging on growing plants continues through the summer (Bowyer, 1991). The restricted availability of quality forage for lactating females may negatively impact nutritional requirements (Parker, Gillingham, Hanley, & Robbins, 1999), fawn survival Parker, Barboza, & Gillingham, 2009), and therefore population dynamics (Cook et al., 2004;Garrott, White, Bartmann, Carpenter, & Alldredge, 1987). A small but perhaps important increase in available NDVIR corresponded with an increase in use of burned areas, suggesting deer attempted to prolong access to growing forage into the summer months (Peek, Riggs, & Lauer, 1979). Greater use of burned areas continued into the winter months, further supporting the inference that these areas provide greater forage availability (McCullough, 1969).
Southwestern forest interiors generally provide lower forageplant density than that available in shrub-covered areas (Altendorf, Laundre, Gonzales, & Brown, 2001). As such, the summer shift from deer using primarily shrubs to both shrub and forested areas coincided with a declining use of available growing vegetation. The shift toward use of habitat with higher tree cover potentially provided shade and thermal buffering as well as cover to avoid detection of fawns during parturition . This habitat shift, along with the diminished availability of growing vegetation, TA B L E 2 Cougar seasonal habitat use candidate models effectively began a period of low-quality forage use in summer months. Our data indicated the timing of parturition occurred primarily in early to mid-June. This was similar to the timing in other populations for which the availability of growing forage generally occurs during spring and summer (Bowyer, 1991;Monteith et al., 2014). This implies Mojave Desert mule deer have not synchronized parturition to the relatively earlier growing season found on sky islands, resulting in a shortened period of forage availability during lactation (Bowyer, 1991). Alternatively, lactating females may have been able to take advantage of the two-month green-up occurring after summer rains. However, our data, which we analyzed over 4month periods, detected no such relationship. Distance to water sources greatly declined in summer, and the need for water may also explain a shift in habitat use (Longshore, Lowrey, & Thompson, 2009;Ordway & Krausman, 1986). However, because water is widely distributed throughout the Sheep Range, with 28 of 38 permanent sources occurring within shrub habitat, water use does not likely explain the shift toward forested areas.
Although we recognize our data are correlative, we found consistent support for our hypothesis that cougar intensity of use influenced both mule deer foraging and habitat use. There was a broad overlap between cougar and mule deer habitat use during spring, as both species used a sympatric range of terrain, elevation, and vegetation types (Hebblewhite & Merrill, 2009). We found the negative effects of cougar occurrence on the ability of mule deer to use greening areas were greatest during spring. However, as indicated by a positive association between mule deer occurrence and NDVIR, deer nonetheless benefitted from spring vegetation growth.
Mule deer may be willing to trade the risk of cougar predation for the seasonally greater benefit of growing vegetation occurring in spring (Festa-Bianchet, 1988). This inference is intuitive, as the greater availability and wider distribution of growing areas in the spring months necessarily reduced the relative patchiness of resources, and therefore risk, to foraging animals.
Cougar intensity of use, although still slightly greater in shrubcovered areas, increased in forested areas during the summer as cougars presumably focused predation activities near fawning sites (Pierce, Bleich, & Bowyer, 2000) and/or increased elevation to reduce heat load ( Figure 6). The seasonal concentration of both mule deer and cougar locations within the forested habitats in summer  (Bowyer, 1991). These findings reinforce our interpretation that although growing vegetation was available during winter, perceived predation risk precluded consistent use of these resources by deer (Barnier et al., 2014;Creel, Winnie, Maxwell, Hamlin, & Creel, 2005 shown to be an important limiting factor in female mule deer nutrition (Parker et al., 1999). As a result, vital nutritional reserves for pregnant females may be suppressed by an increase in cougar occurrence, potentially resulting in reduced reproductive success in Mojave sky-island landscapes occupied by cougars (Parker et al., 2009).

ACK N OWLED G EM ENTS
We thank USFWS refuge manager Amy Sprunger, biologist Dr.

CO N FLI C T O F I NTE R E S T
None declared.

DATA AVA I L A B I L I T Y
All data used for this manuscript will be available to the public within the