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.
Methods
Study area
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).
Juvenile captures
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).
Vegetation surveys
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).
Results
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.
Table 1.
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).
Figure 3.
Ridge specific Kaplan–Meier annual survival rates for black-tailed deer fawns in the Mendocino National Forest from 2009 to 2013. Confidence intervals (95%) are shown for survival rates for each ridge.
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).
Table 2.
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.
Figure 4.
The relative hazard of summer bear predation as bear abundance (camera detections per night × 30) increases at low and high oak forage levels. Low oak forage is the 25% quartile and high oak forage is the 75% quartile of oak forage biomass and 95% confidence intervals are shown for each line.
Discussion
Nutrition and predation both played important roles in the survival of juvenile black-tailed deer during their first year of life. The nutrition hypothesis received the most support in models of all mortality causes, while we found little support for the predator abundance hypothesis in summer. We also found support for an interaction between forage availability and predation risk from black bears (the most common predator). Together with limited evidence of a density limited population, we interpret these results to show that juvenile survival was affected by simultaneous effects of both predation and nutrition, especially during the first month of their life.
The importance of nutrition was shown by the results that juveniles with higher birth weights and in areas with more oak forage were more likely to survive until the end of summer. Nutrition may have even been even more important since we captured many fawns when they were already older than 2 days and possibly missed early fawn mortalities that may have been linked to nutrition. Oak forage was likely linked to better summer nutrition because it made up the majority of the summer diet and oak leaves had lower tannins and higher protein content than Ceanothus spp., the next most common shrub observed in the diets of deer in our study area. Protein is critical for early growth of juveniles and is just as important for summer nutrition as digestible energy (Parker et al. 2009). However, an alternate explanation could be that oak habitat provided more hiding cover for juveniles from predators, but it is unlikely that this benefit would have extended later in the summer after juveniles began following their mothers.
Overall, however, diet quality seemed to play a limited role in explaining the risk of mortality. Despite large differences in diet composition (Supplementary material Appendix 1) black-tailed deer in our study area seemed able to maintain similar diet quality (as measured by fecal N and DAPA) among fawning areas and between ridges. Black-tailed deer are selective foragers that prefer forage with maximum energy and minimal plant toxins (McArthur et al. 1993), and had diverse diets throughout our study area (Supplementary material Appendix 1). Diverse diets have also been linked to better nutritional condition in other ungulates (Parikh et al. 2017), and other herbivores are adept at exploiting the best forage available, often in the face of predation risk (Camp et al. 2017). There also seemed to be no carryover effects of nutrition from the preceding winter since there was no detectable difference in body condition of adult females across the study area at the beginning of the summer.
Despite the support for the nutrition hypothesis, there was little population level evidence that juvenile mortality was strictly a result of a food limited deer population. Determining the ultimate consequences of nutrition and predation on populations requires assessing the relationship of a population to its carrying capacity K (Errington 1946, Bowyer et al. 2005). We used established indices of population level processes (Bowyer et al. 2005, Pierce et al. 2012) to determine the degree of resource limitation present in our strongly declining black-tailed deer population (lgr; = 0.82 ± 0.13; Marescot et al. 2015). Despite a high variation in relative density between fawning areas, deer density did not explain variation in juvenile mortality. Estimates of pregnancy and fetal rates were higher (Marescot et al. 2015) than averages reported for mule and black-tailed deer across their distribution (Forrester and Wittmer 2013), results not typically observed in food-limited populations (Gaillard et al. 1998). Since estimates for average fetal rates reported in Forrester and Wittmer (2013) were mostly from high-density populations with an average lgr; = 0.99, they should provide a benchmark for equilibrium populations. The mean body condition of adult females in early summer was close to ‘good’ (mean rBCS = 2.8 on a 1–5 scale) at a time when fat reserves should be near the lowest point for the year (White 1992, Parker et al. 2009).
Predator abundance did not explain much variation in mortality risk. Predator abundance alone was not in the top summer mortality risk models, and predator abundance was not significant when it was present in top models with other variables. However, predation was still an important mortality source. Predation during the first month of life was the most significant cause of annual mortality. Bear predation was the largest overall source of juvenile mortality, although it was almost all concentrated on juveniles <30 days old, a pattern also found in other ungulates (Vreeland et al. 2004, Gustine et al. 2006, Griffin et al. 2011). The risk of being killed by a bear was unrelated to the birth weight of juveniles, a measure typically used as an index of nutritional status at birth (Monteith et al. 2014). Birth weight was likely not related to risk of bear predation because juveniles were still relying on hiding from predators during the first weeks of life (Geist 1981), and were unlikely to evade predators once detected despite differences in birth weight or overall nutritional condition.
The importance of predation was also shown by spatial variation in bear predation that resulted in significant differences in yearling recruitment. Higher bear predation on the FH7 ridge led to lower yearling recruitment despite higher birth weights on this ridge, which, together with the lack of evidence of food limitation, is evidence of additive summer mortality from bear predation. Bear predation was also linked to an interaction between forage availability and predation risk. The interaction between combined abundances of bears and coyotes and oak forage availability was in the top summer survival models even though the interaction was not significant. We did find a significant interaction between bear abundance and the biomass of oak forage. The risk of bear predation was generally lower in areas with more oak forage, but at low bear abundance the risk of bear predation was almost six times higher in areas with low oak forage (25% quartile) than high oak forage (75% quartile). However, as bear abundance increased, the predation risk dropped in areas with less oak forage and was comparable to risk in areas with more oak forage.
Our study was observational and so we were not able to determine the mechanism driving this interaction. One likely explanation for high bear predation risk in areas with low bear abundance and low oak forage is a spatial mismatch between oak forage and preferred bear foods while another is differences in the hunting methods of individual bears. Bears heavily forage on vegetation in spring (Greenleaf et al. 2009), and herbaceous forage was more abundant in the fawning areas on the FH7 ridge compared to the dry oak savannah that characterized fawning areas on the M1 ridge. Most bears have been shown to encounter juvenile ungulates opportunistically while consuming other forage resources (Zager and Beecham 2006, Bastille-Rousseau et al. 2011), and as a result bears on FH7 could have been more likely to encounter fawns. More oak resources on the M1 ridge in the fall likely sustained higher overall bear populations, and as a result the M1 ridge had a higher bear abundance but a lower risk of bear predation while the FH7 ridge had less bears, less oak forage and higher bear predation risk. Alternatively, a small subset of bears may have focused on hunting juvenile deer in certain areas (Zager and Beecham 2006, Bastille-Rousseau et al. 2011), and if more bears on the FH7 ridge were specializing on hunting ungulates, this could have caused the observed interaction. It is unlikely that the interaction between bear abundance and oak forage was simply because oak shrubs provided hiding cover (Shallow et al. 2015), due to the fact that bear predation was high in fawning areas with large amounts of cover from other shrubs, such as mountain whitethorn ceanothus C. cordulatus, that juveniles often used as hiding cover.
Unlike bear predation, the risk of coyote predation was related to the birth weight of juvenile deer. After accounting for the amount of herbaceous forage in a fawning area, juveniles with higher birth weights were less likely to be predated by coyotes. Since coyotes predated both neonate and older juveniles, older juveniles in better condition may have been able to escape coyotes (Lingle et al. 2008). It is also likely that coyotes were switching to preferred small mammal prey such as black-tailed jackrabbits Lepus californicus and California ground squirrels Otospermophilus beecheyi in habitats with more herbaceous vegetation (Hamlin et al. 1984, Hurley et al. 2011). Our findings that coyote predation was strongly related to birth weight and available forage matches previous research showing that coyote predation on juvenile deer was mostly compensatory (reviewed by Forrester and Wittmer 2013).
Combined, our findings provide evidence that juvenile black-tailed deer, like other mid-sized ungulates, are affected by both top–down and bottom–up forces (Hopcraft et al. 2010). Factors related to nutrition explained much of the variation in mortality risk, and coyote predation and mortality from all causes appeared to be mostly influenced by nutritional status through the birth weight of juveniles. These results support findings from other studies showing the importance of forage in explaining variation in juvenile survival (Shallow et al. 2015). Bear predation showed how top–down forces may also be simultaneously affecting ungulate populations and how predation may be a source of additive mortality. Bear predation has been shown to be a source of additive mortality for other ungulates (Griffin et al. 2011) and to cause the decline of deer populations even in areas with abundant forage (Monteith et al. 2014). Our results also complement findings that forage availability influenced vulnerability to predation for adult female deer in the study area (Forrester et al. 2015) and that low female survival due to predation was responsible for the declining population (Marescot et al. 2015). Understanding predator specific interactions and age specific mortality are key to understanding how predation and nutrition simultaneously affect mule deer populations. Meta-analysis similar to work in other ungulate communities (Hopcraft et al. 2010, Griffin et al. 2011) that accounts for predator identity will be vital to gaining deeper insight into the role of predation in juvenile survival and ungulate population growth rates across a range of predator diversity and abundances.
Acknowledgements
TDF thanks the Robert and Patricia Switzer Foundation, the UC Davis Graduate Group in Ecology, and the Stockton Sportsmen's Club. We thank D. Casady of the California Dept of Fish and Wildlife and our field crew for their dedicated efforts. Finally, we thank D. Kelt, A. Latimer, A. Sih, M. Hurley and J. Forbey for their comments that greatly improved this manuscript.
Funding – We received support from the California Department of Fish and Wildlife (contract no. P0880013), the California Deer Association and the Mendocino County Blacktail Association.
References
Appendices
Supplementary material (available online as Appendix wlb-00510 at < www.wildlifebiology.org/appendix/wlb-00510>). Appendix 1, 2.
