Understanding how top–down and bottom–up effects influence population dynamics of ungulates is essential for effective management and conservation, and there is an emerging consensus that forage productivity and quality interact with predation to influence survival. From 2009 to 2013, we captured and monitored 135 neonatal black-tailed deer Odocoileus hemionus columbianus in coastal California to study possible interactions between forage and predation on survival. We estimated seasonal and annual survival rates, assessed the cause of all mortalities (n = 93), measured available forage and diet quality, estimated relative abundances of predators on summer range each year, and used remote sensing to quantify habitat types on winter range. We evaluated the relationship between fawn survival and environmental covariates with cumulative incidence and proportional hazards functions. Summer survival rates averaged 0.40 (SE = 0.05) across all years and the mean annual survival rate was 0.26 (SE = 0.04). We found that most juvenile mortality resulted from predation during summer, and spatial differences in summer survival persisted until recruitment. Variation in mortality risk from all causes was related to the birth weight of juveniles and available oak forage but not predator abundance. The risk of black bear predation, the single largest cause of mortality, was unrelated to birth weight and showed an interaction between bear abundance and the amount of oak forage. Additionally, the good body condition of adult females in spring and a lack of relationship between mortality risk and variation in deer density did not provide evidence for purely bottom–up limitation. Rather our results provide evidence that both bottom–up and top–down effects were influencing fawn survival in this declining population, and that predator identity and the timing of mortality affected the impact of predation.
Unexpected declines of mule Odocoileus hemionus and black-tailed deer O. h. columbianus across the western USA (Ballard et al. 2001, Forrester and Wittmer 2013) have highlighted widespread challenges in determining the role of top–down and bottom–up effects on game species in complex communities (Sinclair and Krebs 2002, Sinclair 2003, Owen-Smith and Mills 2008). Understanding how predation and forage interact to influence survival and population dynamics is of particular interest for ungulates because they are dominant herbivores (Augustine and McNaughton 1998) and important components of trophic cascades (Ripple et al. 2014). More recently a consensus has begun to emerge suggesting that populations of mid-sized ungulates are regulated by the combined effects of predation and forage availability (McNaughton et al. 1989, Sinclair et al. 2003, Hopcraft et al. 2010).
Understanding relative effects of predation and nutrition on juveniles is particularly important since highly variable juvenile survival and recruitment typically drive population dynamics in ungulates (Gaillard et al. 1998, Forrester and Wittmer 2013). It is also evident that both predation and forage availability have a larger effect on juvenile ungulates than adults (Gaillard et al. 2000). Untangling the effects of predation and forage availability on the survival of juvenile ungulates is complicated by the fact that predation is the primary source of mortality for juveniles in almost every community and species (Gaillard et al. 2000), including mule and black-tailed deer (Forrester and Wittmer 2013).
Predation is high for juvenile ungulates because juveniles may be preyed upon by up to four times as many predators as adults (Linnell et al. 1995), and even in species-poor mammal communities there is still at least one predator of juveniles (Moorter et al. 2009). Predation mortality predominately occurs in the summer months immediately following birth when juveniles are the most vulnerable (Gaillard et al. 2000). However, the true effect of predation on populations can be hard to discern because predation may have compensatory or additive effects on mortality (Ballard et al. 2001, Bowyer et al. 2005). Predation on juveniles has been found to cause declines in some mule deer populations even in areas with abundant forage (Monteith et al. 2014), while predator reductions in other areas have increased juvenile survival but did not change population growth rates (Bartmann et al. 1992, Hurley et al. 2011).
In contrast to predation, bottom–up effects on juveniles are typically moderated by maternal condition (Parker et al. 2009). The highest nutritional demands of the year for female ungulates occur during late pregnancy and lactation (Clutton-Brock et al. 1989). During this time females require forage with high amounts of energy and protein to support juvenile growth and to replenish body fat (Parker et al. 2009). Because nutrition effects are mediated through the mother and starvation may not be the ultimate cause of death, nutritional effects on juvenile survival are almost always more cryptic than predation (Pierce et al. 2012, Monteith et al. 2014). Lower forage availability has been shown to result in poor maternal body condition and lower juvenile survival (Shallow et al. 2015), while increased winter nutrition has been found to reduce mortality from both starvation and predation (Bishop et al. 2009). While bottom–up effects in ungulates can be cryptic, they can still be assessed in observational studies because juvenile survival and reproduction in ungulates show predictable patterns with increasing population density (Gaillard et al. 1998, Eberhardt 2002). In nutritionally stressed populations both fawn birth weight and maternal body condition should decline, even though adult survival will likely be unaffected (Parker et al. 2009). These and other population level indicators can be used to determine how close a population is to nutritional carrying capacity (K) and evaluate the role of predation and nutrition as ultimate causes of fawn mortality (Pierce et al. 2012).
Despite our understanding that most mid-sized ungulates are regulated by a combination of top–down and bottom–up effects (Hopcraft et al. 2010), predation and nutrition are still mostly treated as dichotomous. Here we present a study of the effects of forage availability and predation on juvenile survival from birth to recruitment in black-tailed deer in northern California. We considered three alternate hypotheses regarding the effects of nutrition, predator abundance and possible interactions between nutrition and predator abundance on juvenile survival in the summer. The nutrition hypothesis tested whether variation in forage availability, forage and diet quality, the density of deer populations in known birth areas, and birth weight contributed to juvenile mortality. Contrarily, the predator abundance hypothesis tested whether variation in predator occurrences across birth areas best explained differences in juvenile mortality. The predation and forage interaction hypothesis tested whether an interaction occurred between predator abundances and covariates that measured nutritional status to determine if top down and bottom up forces were both influencing juvenile mortality. We evaluated these hypotheses in relation to all juvenile mortality and for predator specific mortality in summer. Finally, we also evaluated whether habitat composition, elevation, summer nutrition, or weather were affecting winter survival of juveniles.
Our study area covered ∼1000 km2 in the Mendocino National Forest in the northwestern California Coast Range, and was composed of two major ridges (named M1 and FH7 for the forest roads traversing them; Fig. 1). The terrain in the area is rugged with sharp gradients in elevation, ranging from 500 m in valley bottoms to >2000 m on ridgetops. The climate is considered Mediterranean and >85% of all precipitation occurred from October through April. Snow cover was generally limited to elevations >1000 m and was irregular, particularly during mild winters.
Vegetation communities in the area were diverse and consisted of a mix of oak woodlands (Quercus spp.), dense chaparral and grasslands at low elevations, while mid elevations were mainly mixed-coniferous forests dominated by pine (Pinus spp.) and Douglas-fir Pseudotsuga menziesii. High elevation areas supported a mix of true fir forests (Abies spp.), patches of shrubs (mainly Ceanothus spp.) and scattered dry and wet meadows. Past logging and grazing left a mosaic dominated by even-aged conifers and grasslands dominated by exotic species, with occasional mature timber stands.
Black-tailed deer wintered in valleys and moved to high quality habitats at high elevations during summer to give birth and raise young (Allen et al. 2014). Predators included American black bears Ursus americanus, coyotes Canis latrans, bobcats Lynx rufus and pumas Puma concolor. Black-tailed deer were the only resident ungulate and pumas were the principal predator of adult deer (Marescot et al. 2015). Hunters were only allowed to harvest male black-tailed deer and poaching rates were likely low (Wittmer et al. 2014).
All handling procedures were approved by an Institutional Animal Care and Use Committee at the University of California, Davis (protocols 15 341 and 16 886) and adhered to guidelines established by the American Society of Mammalogists (Sikes and Gannon 2011). We captured juvenile black-tailed deer from mid-June to mid-July during 2009–2012 by driving along unpaved forest roads during daylight, by using spotlights to locate juveniles at night, and by scanning meadows and forest habitat with binoculars for post-parturition does to find hidden juveniles. We captured juveniles by hand or with handheld nets wearing new latex gloves for each capture to avoid scent contamination. Upon capture, we weighed, sexed and then fitted juveniles with a small colored and numbered plastic ID tag in one ear and a very high frequency (VHF) motion-sensitive transmitter (Sirtrack, Havelock North, New Zealand) in the other ear. Battery life of transmitters was one year. We estimated juvenile age in the field using a combination of hoof growth line measurements (Sams et al. 1996) for juveniles over approximately 4–5 days old and the status of the umbilical cord and standing/walking coordination for younger juveniles. Since age estimations were not precise we assigned each juvenile an age in weeks since birth (one or two weeks) and used this value for survival analysis. Juveniles were released near the capture site immediately after processing, which averaged approximately 10 min. Juvenile birth weight was estimated by regressing juvenile age against weight for all juveniles in a given fawning area and using the linear regression equation to predict weight at birth. Birth weight is regarded as an index of body condition of the mother and juvenile for ungulate species (Parker et al. 2009, Monteith et al. 2014).
Monitoring and mortality investigation
We monitored the status of juveniles daily from capture to mid-September and every 7–14 days from either the ground or air until June the year following capture, ending in June 2013. VHF transmitters switched to mortality signal after remaining stationary for 4 h. We investigated summer mortalities almost immediately following detection of a mortality signal (AVG = 1.1 days, SE = 0.25, n = 63), while inclement weather and limited access delayed investigations in winter (AVG = 24.2 days, SE = 8.1, n = 19) and prevented accurately assessing cause of mortality. Cause of mortality in summer was determined by site investigations using systematic criteria including disposition of the carcass, predator sign, evidence of caching, bite marks and blood (Atkinson and Janz 1994). If no obvious evidence of predation was found, we performed a field necropsy to assess if death was result from adenovirus (plasma found in body cavity), other diseases, or malnutrition (emaciated carcass; red jelly-like bone marrow). To confirm our field assessment, we also swabbed carcasses and tested samples for predator DNA following published protocols (Wengert et al. 2014). We considered deaths capture related if they occurred within a week of capture and our investigation could determine no other cause of mortality (including disease).
Population structure and body condition
As a part of a larger research project we also monitored 60 adult GPS collared female black-tailed deer captured on the same summer areas in the same years as juveniles, and estimated fetal and pregnancy rates (methods described in Marescot et al. 2015). From adult location data, we delineated birth areas (hereafter fawning areas) (n = 4) and winter ranges (n = 4) using 95% minimum convex polygons (Fig. 1). We used an average of 13 adult females to delineate each winter seasonal range and 11 adult females to delineate each summer seasonal range (winter ranges: Elk Cr. = 11, Upper Black Butte = 11, Lower Black Butte = 14, Grindstone = 16; summer ranges: Cold Spring = 9, Plaskett Meadows = 13, Coyote Rock = 10, Cherry Hill = 12) and all fawns were captured within summer ranges. We considered areas to be spatially separated if location data indicated gaps greater than twice the size of the average diameter of an adult home range and where geographic features created barriers to movement. DNA results confirmed significant female reproductive isolation between the two major ridges (Bose et al. 2017). During adult captures we also assessed body condition using a modified rump fat body condition scores (rBCS range from 1 to 5, Gerhart et al. 1996, Cook et al. 2010), and we tested for differences in body condition score among fawning areas using a one-way analysis of variance (ANOVA). We estimated deer densities associated with identified fawning areas using fecal DNA and capture–mark–recapture (Lounsberry et al. 2015).
We surveyed all fawning areas in the summer of 2010 and 2011 to quantify percent cover of deer forage types and to estimate biomass of shrubs, forbs and grasses. Forbs and grasses were also surveyed again in the summer of 2012. We conducted surveys along 100 m line transects with random starting points that we located on grids covering each fawning area with 1 × 1 km spacing. We estimated shrub cover and species composition using line-intercept surveys on each transect (Bonham 1989) and estimated shrub forage biomass using twig counts on three 1 × 3 m quadrats per transect (Shafer 1963). We estimated herbaceous biomass using the comparative yield (CY) and dry weight ranking (DWR) methods using 10 different 0.25-m2 quadrats per transect (Haydock and Shaw 1975, Jones and Hargreaves 1979). We identified shrubs to species, classified all small flowering plants as forbs, and categorized grasses as annual or perennial. We conducted 157 line transect surveys, conducted CY and DWR surveys on 1770 quadrats, and counted all twigs equal or smaller to typical deer browse diameter on 471 quadrats. We measured 100–200 browsed twigs to obtain the mean species-specific browse diameter for important deer browse (Ceanothus, Prunus, Arctostaphylos and Quercus species). We estimated habitat specific forage amounts for classification and assessment with Landsat of visible ecological groupings (CALVEG) cover types (Schwind and Gordon 2001); conifer, hardwood, mixed conifer and hardwood, shrub and herbaceous. We created habitat weighted estimates of forage by estimating the amount of forage for each habitat type per fawning area (forage g/m2 × habitat area), summing these values from all habitat types, and then dividing by the total area. Herbaceous biomass was variable among years so we calculated herbaceous forage values for each year, while shrub biomass was less variable and we calculated average shrub browse for all years combined.
Deer diet analysis
We collected deer pellets in fawning areas during the summers of 2010, 2011 and 2012 on transects that followed deer trails with randomly located starting points distributed across available habitat types. Diet composition was analyzed using microhistological analysis (Holechek et al. 1982, Leslie et al. 1983) and diet quality was indexed using fecal nitrogen and diaminopimelic acid (DAPA) (Hodgman et al. 1996). We used previous work (Wallmo 1981, Kie et al. 1984) and our dietary analysis to determine the most important shrubs for deer in our study area and estimated nutritional quality for these species from samples collected in mid to late summer, including crude protein, in vitro dry matter digestibility, and tannin analysis (Supplementary material Appendix 1 Table A1, A2). All diet analyses were performed by the Wildlife Habitat and Nutrition Laboratory at Washington State University. We assessed differences in diet among fawning areas using one-way ANOVA.
Winter range weather, elevation and habitat
We used the monthly Pacific Decadal Oscillation (PDO) and the El Nino Southern Oscillation (ENSO) values as an index for climate variation for the winter months of our study. Both indices are measures of surface sea temperatures that directly relate to temperature and precipitation trends in the Pacific Northwest and northern California. We also included the average monthly elevation on winter range of GPS collared deer to account for possible variation in temperature and precipitation related to elevation. Although we did not collar mother and offspring pairs we confirmed that winter juvenile locations fell within the designated winter ranges.
We created a forage availability index for winter ranges by estimating the total area of CALVEG vegetation types containing high-quality forage. High quality vegetation types included oak woodland, herbaceous meadows and shrub types with high quality forage (e.g. montane mixed chaparral) (Wallmo 1981, Livezey 1991). We included a variable for the available oak forage on summer range for each surviving juvenile to test for possible carryover affects from summer range.
Relative abundance of predators
We estimated the abundance of predators in summer fawning areas using temperature and motion-triggered cameras (Bushnell Trophy Cam, Bushnell, Overland Park, KS and Cuddeback Capture IR, Cuddeback, Green Bay, WI). We randomly sampled fawning areas by placing a 12–14 km2 grid of 1-km2 cells with a random starting point over the four fawning areas. We randomly selected 8–10 of the grid cell centers to deploy cameras each month, and randomly selected another 8–10 grid cell centers with replacement for each subsequent sampling period. We placed cameras at areas of animal activity (trails, closed roads, springs, etc.) within 100 m of randomly selected points (Rowcliffe et al. 2008). We used the average summer home range size (∼1 km2) of adult female black-tailed deer in our study area (Bose et al. 2017) as a grid cell to estimate predator use and abundance on the scale of a female deer home range (MacKenzie et al. 2005). We deployed cameras for three one-month periods beginning in mid-June and ending in September during 2010, 2011 and 2012, and deployed 275 cameras for 8980 trap nights over three summers.
We used program PRESENCE to model both probability of use and detection probability for predator species, but excluded pumas due to insufficient detections (puma densities of 0.68 pumas/100 km2 in the entire study area were comparatively low; Allen et al. 2015). Detection probabilities of predators did not differ among fawning areas, and probability of use often approached 1 so we used the monthly detection rate of predators ((no. predator detections/camera days) × 30) as an index of predator relative abundance. We estimated predator relative abundance separately for three one-month periods, mid-June to mid-July, mid-July to mid-August and mid-August to mid-September to account for possible variation over the summer season.
Modeling juvenile mortality risk
We defined summer separately for each juvenile deer as the period from capture until their last location on summer range and winter as the time from the first relocation on winter range until one year of age. The date of the initial mortality signal was used as the date of death or the date between the last live location and the first mortality signal if there was a gap of >3 days. We limited the Cox proportional hazards analysis to the three cohorts captured from 2010 to 2012 (n = 121, mortalitysummer = 63) since we did not collect covariate information for the 2009 cohort.
We used cumulative incidence functions (CIFs) that model cause-specific mortality while accounting for all other causes of mortality to estimate juvenile mortality and survival rates by month throughout the first year of life (Heisey and Patterson 2006). These functions are based on proportional hazards models (Cox 1972) and model the probability of a mortality from cause i occurring before time t.
We modeled juvenile mortality risk with Cox proportional hazards using the formula
where t is time as specified in the model (e.g. days since birth), h(t|Xj) is the hazard rate for the jth deer at time t, h0(t) is the baseline hazard, and the regression coefficients βx are estimated from the risk covariates Xj for the jth deer (Cox 1972, Therneau and Grambsch 2000). The βx are used to estimate hazard ratios (HR) that are similar to an odds ratio, where a HR of less than or greater than 1 represents a smaller or greater risk of death respectively. We considered a HR significant if the 95% confidence interval did not overlap 1.
We modeled survival as a function of age in days because we captured juveniles soon after birth (Fieberg and DelGiudice 2009). We used a delayed entry design and estimated risk beginning at birth but juveniles entered the analysis at the day of capture for summer survival and the day of arrival on winter range for winter survival. We censored animals from analysis after death, after the last day on summer range, or after survival to one year in the winter analysis (Hosmer et al. 2011). We adhered to the guideline of 8 mortalities per covariate recommended by Vittinghoff and McCulloch (2007). The maximum number of parameters for summer survival models (63 summer mortalities) was seven and we used one parameter for winter survival models. We tested the assumption of proportional hazards for covariates in the model using Shoenfeld residual plots (Grambsch and Therneau 1994). We assessed if outliers unduly affected the model by graphing DFBETA residuals (Cleves et al. 2010) and likelihood displacement values (Collett 2003) against analysis time. We used pairwise correlation coefficients to assess if covariates were highly correlated (correlation >0.5) and then chose the most biologically relevant covariate based on literature references.
We modeled whether fawning area or study year explained more variation in survival than the null model to determine if we needed to include these nuisance variables in our models. Nutrition covariates for summer survival models included the biomass of oak and herbaceous forage, fecal nitrogen and DAPA as measures of diet quality, forage biomass weighted by diet quality measures, relative density of deer in fawning areas, and estimated birth weight of juveniles. The predation covariate was the relative abundance of bears and coyotes, the predators that accounted for 77% of all predation mortalities. We also included covariates that may have been related to juvenile survival to control for confounding effects, including previous winter precipitation, spring precipitation, sex of juveniles and age in weeks of juveniles at capture.
Covariates for winter models included the monthly Pacific Decadal Oscillation (PDO) index, the total area of oak, herbaceous and shrub habitat types on winter range, the average monthly elevation of all collared female deer for each wintering area, and birthweight and summer range oak forage biomass. Due to lower numbers winter mortalities in our sample we only compared seven models with a single covariate. The PDO was chosen as a climate covariate because it indexes temperature and precipitation values together (Mantua et al. 1997), winter habitat types that contained important forage were chosen to model general forage availability, the average elevation of all collared deer per month in a given wintering area was used to control for the difference in elevation between the wintering areas (202–369 m average Jan. elevation), and birthweight and oak biomass on summer range were included to account for possible nutritional carryover effects.
We used our hypotheses to build an a priori model set for juvenile summer survival. Since the numbers of covariates that we could use were limited (Peduzzi et al. 1995, Vittinghoff and McCulloch 2007) we first selected a top nutrition model from a model set containing all nutrition variables, including interactions between forage and birthweight and deer density and birthweight ( material Appendix 2). We also compared models for sex, age in weeks at capture, and both covariates together to determine the top model for nuisance variables. We then combined the top nutritional and nuisance models with models for predator abundance and the best weather model. Finally we tested for an interaction between forage abundance and predator abundance in all models containing both variables (Supplementary material Appendix 2). We tested possible non-linear interactions using fractional polynomials (Royston and Sauerbrei 2004).
We evaluated the importance of nested models using criteria defined in Arnold (2010). We used Akaike information criterion adjusted for small sample sizes (AICc) to rank models and considered a model to be supported if the AICc score was <2 AICc from the next model (Burnham and Anderson 2002). We did not consider a model that only differed from the best model by one parameter with similar log-likelihood to be competitive (Burnham and Anderson 2002, Arnold 2010).
We modeled the cause specific risk of bear and coyote predation using cumulative incidence functions (CIFs) (Fine and Gray 1999) that models the CIF for cause i as the cumulative sub-hazard function for that cause alone. Covariate effects for cause i can be interpreted similarly to a Cox proportional hazards model. We tested assumptions of the CIF models with the same methods as the Cox proportional hazards models and ranked models using AICc (Burnham and Anderson 2002). We used a maximum of three parameters for bear CIF models and two parameters for coyote CIF models based on the number of cause-specific mortalities (Vittinghoff and McCulloch 2007). We modeled eight forage and diet quality covariates (oak forage, herbaceous forage, diet quality weighted oak and herbaceous forage, and the two fecal measures of diet quality) to determine the single best nutritional model. This model was compared to models of bear or coyote abundance (depending on the species being modeled), birthweight and age at capture in weeks, and a model with the interaction between forage and bear abundance (bear CIF only). We performed all statistical tests in STATA ver. 12.1 (StataCorp, College Station, TX).
Juvenile capture and monitoring
We captured 137 juveniles (72 females, 64 males, 1 unknown) during the summers of 2009–2012. Two juveniles were censored due to tag failure (n = 1) and capture related mortality (n = 1). The mean capture date was June 27 (SEamong years = 6.40 days, range 6 June–19 July), the mean age at capture was 4.8 days (range of 0–10), although the true mean age likely ranged from 2.9 to 6.5 days given sampling variation in age estimation (Grovenburg et al. 2014). The mean capture weight was 3.7 kg (SEall years = 0.08, range 2–7). Mean capture date differed significantly among years (ANOVA, F3,134 = 17.24, p < 0.001), but mean capture age (Kruskal–Wallis, Xdf 2 = 3 = 6.873, p = 0.076) and capture weight (ANOVA, F3,134 = 0.80, p = 0.493) were consistent.
Summer diet and fawning area vegetation
Deer diet composition and forage quality results are reported in detail in Supplementary material Appendix 1. Diet was averaged between years and was mostly composed of shrubs (Cherry Hill = 88%, Coyote Rock = 83.1%, Cold Spring = 85.8%, Plaskett Meadows = 53.6%), while forbs contributed only a small proportion (Cherry Hill = 2.1%, Coyote Rock = 3.9%, Cold Spring = 4.8%, Plaskett Meadows = 11.5%). Forbs may have been underestimated due to known issues with differential digestibility of forage types. Oak leaves composed most of the diet in summer in all areas except for Plaskett Meadows (Cherry Hill = 76.1%, Coyote Rock = 65.4%, Cold Spring = 73.6%, Plaskett Meadows = 21.8%). The amount of fecal nitrogen in deer pellets (Supplementary material Appendix 1) did not differ among fawning areas (ANOVA, F3,8 = 0.83, p = 0.51) or years (ANOVA, F2,9 = 2.52, p = 0.14). The amount of DAPA found in deer pellets (Supplementary material Appendix 1) did not differ among fawning areas (ANOVA, F3,8 = 0.14, p = 0.94), but did differ among years (ANOVA, F2,9 = 8.87, p = 0.007) showing an increase each year from 2010 to 2012.
Maternal condition, pregnancy and deer density
The average adult female rBCS score was 2.8 (SE = 0.37), and there were no differences in body condition among fawning areas (ANOVA, F3,54 = 0.91, p = 0.44). Pregnancy rates averaged 0.87 (SE = 0.05) and average fecundity was 1.9 juveniles per doe (Marescot et al. 2015). The relative density of black-tailed deer varied significantly among fawning areas (Lounsberry et al. 2015) (Cherry Hill = 42.1, Coyote Rock = 19.5, Cold Spring = 37.6, Plaskett Meadows = 41.1, deer per km2).
Temporal and spatial patterns in juvenile survival
A total of 93 juveniles died during our study, including 74 during summer. Summer survival rates for juveniles averaged 0.40 across all years (SEamong years = 0.05) and the mean annual survival rate was 0.26 (SE = 0.04). Summer (Min = 0.42; 2010, Max = 0.51; 2011) and annual (Min = 0.21; 2012, Max = 0.37; 2011) juvenile survival fluctuated among years, but the differences in summer survival were not significant (Cox hazards, probability of difference from 2010, p2011 = 0.51, p2012 = 0.57). Winter survival of juveniles averaged 0.63 across all years (SEamong years = 0.07).
Predation was the primary source of juvenile mortality (Fig. 2), and black bear predation was the largest single source of mortality (Table 1). The majority (61%) of total mortality and of predation mortality (69%) occurred within 30 days of birth. During summer, there were low numbers of mortalities assessed as unknown predators (5% of summer mortality) or unknown cause (7.7% of summer mortality). Only 22% of annual mortality occurred on winter range, and most known causes were attributed to predation. No winter mortalities were attributed to malnutrition but we could not assess the cause of mortality in most instances (unknown mortalities = 14 of 19 total). Summer survival did not differ among fawning areas (Cox hazards, probability survival different than Cherry Hill area, pCoyoteRock = 0.68, pColdSpring = 0.17, pPlaskettMeadows = 0.14) or winter areas (Cox hazards, probability survival different than Elk Creek area, pUpperBlackButte = 0.76, pLowerBlackButte = 0.34, pGrindstone = 0.15). However, differences in summer survival between the two ridges in the study area resulted in a difference in annual survival (Fig. 3, LR test, Xdf 2 =1 = 3.69, p = 0.05).
Summer mortality risk
We pooled data across years after confirming that there were no differences in summer survival among years. Correlation coefficients were >|0.5| between herbaceous forage and both overall shrub cover and Ceanothus species. We retained herbaceous forage for modeling because herbaceous forage is critical summer forage for mule deer (Wickstrom et al. 1984, White 1992), and dropped Ceanothus spp. since these species did not contribute much to summer diets (Supplementary material Appendix 1 Table A1). We dropped shrub cover because it was measured as vegetation composition and not specifically as hiding cover. All remaining environmental covariates met proportional hazards assumptions.
There were five models within 2 AICc of the top summer survival model. The estimated birthweight of juveniles and the amount of oak forage in the fawning area were included in all of the top models (Table 2) and were related to significantly lower risk of juvenile mortality. The effect size of birthweight was consistent among models, with a 0.5 kg increase in birthweight resulting in an average 28% (27–29% range) reduction in mortality risk (standard deviation of juvenile birthweight = 0.6 kg). The effect of available oak forage was consistent as well, with a 10% increase in the average available forage resulting in an average of 3% reduction in mortality risk (3–4.5% range). The notable exception was the predation × oak forage model where a 10% increase in oak forage reduced mortality risk by 7.5%. The age in weeks in capture was not a significant predictor of mortality in any models, but did reduce the model deviance compared to nested models when it was included (Table 2) and so we considered it to be an informative parameter (Arnold 2010). Age at capture showed a trend that juveniles caught in their second week of life were more likely to die during the summer than juveniles caught in their first week. Bear and coyote abundances were not a significant predictor of mortality in any of the top models, and the interaction between oak forage and bear and coyote abundance and birthweight and bear and coyote abundance were not significant in any model. However, models that included the oak forage and bear and coyote abundance interaction showed reduced deviance compared to nested models (Table 2), and also showed a stronger relationship between juvenile mortality and oak forage biomass compared to other models (Supplementary material Appendix 2). There were no non-linear interactions between bear and coyote abundance and oak forage.
Fawn mortality rates. Cause specific mortality rates, total mortality and survival rates calculated from cumulative incidence functions (CIF) for black-tailed deer fawns in Mendocino National Forest from 2009 to 2013.
The three models that best explained coyote predation risk were herbaceous forage (AICc = 146.85, wi = 0.33), birthweight (AICc = 147.22, wi = 0.27), and birthweight and herbaceous forage together (AICc = 148.01, wi = 0.18). The only covariate that showed a significant relationship with coyote predation risk was birthweight in the oak forage–birthweight model. Juveniles that were 0.5 kg higher in birthweight had a 34% lower risk of coyote predation. However, birthweight was not significant in the birthweight only model. The amount of herbaceous forage in a fawning area showed a trend for lower mortality risk, but was not a significant predictor of coyote predation (Supplementary material Appendix 2).
The top models explaining variation in predation risk from bears contained the ridge that the fawns occupied during the summer (AICc = 259.78, wi = 0.27), the amount of oak forage on a fawning area (AICc = 260.43, wi = 0.18), the forage–bear abundance interaction model (oak forage, relative abundance of bears, oak forage × bear abundance; AICc = 260.46, wi = 0.18), and the ridge–bear abundance interaction model (ridge, relative abundance of bears, ridge × bear abundance; AICc = 261.32, wi = 0.17). The other seven models were >2 AICc from the top model. Ridge was not significant in its own model, but in the interaction model juveniles on FH7 were up to 8 times more likely to die from bear predation on FH7 than M1 (sub-hazard ratio = 8.15, p = 0.023) after controlling for bear abundance and the interaction between bear abundance and ridge. The interaction between ride and bear abundance was marginally non-significant (p = 0.07). Oak forage was not significantly related to mortality risk from bears on its own, but in the interaction model juveniles in areas with more oak forage were 35% less likely to die from bear predation (sub-hazard ratio = 0.65, p = 0.016). The interaction between bear abundance and oak forage was also significant (sub-hazard ratio = 1.18, p = 0.043). In areas with lower bear abundance the risk of bear predation was higher in areas with less oak forage, while in areas with more bear abundance bear predation risk is similar for juveniles at all levels of oak forage (Fig. 4).
Summer mortality risk models. Top Cox proportional hazards models (<2 ΔAICc). Model deviance is -2 times the log likelihood (-2LL). Nutrition, predation and predation × nutrition correspond to the three hypotheses of fawn mortality, while Individ. refers to models with capture age included. Capture week is the age of juveniles at capture in weeks. Oak is the biomass of oak (Quercus spp.) forage and Predator is the relative abundance of bears and coyotes, both within a given fawning area. Variables that were significantly correlated with mortality risk are indicated by (+) and (-) showing increasing or decreasing mortality risk respectively, while (0) indicates no significance.
Winter mortality risk
All covariates in winter hazards models met proportional hazards assumptions. The single best model predicting winter mortality risk was the winter monthly PDO index (2.91 ΔAICc from the next nearest model). All other models (winter habitat, elevation and summer nutrition) did not perform better than the null model (Supplementary material Appendix 2). The PDO index ranged from 0.08 to -2.33 (x̄ = -0.81, SD = 0.56) in the winter months of our study, and the PDO is positively associated with warmer temperatures and lower precipitation and negatively associated with lower temperatures and higher precipitation (Mantua et al. 1997). Positive increases in the index (warmer winter temperatures and lower winter precipitation) were associated with a higher risk of mortality, and an increase of one standard deviation in the PDO index was associated with an 87% higher risk of winter mortality.