Roe deer face competing risks between predators along a gradient in abundance

. Mortality rates and patterns are fundamental demographic traits for understanding the dynamics of populations of large herbivores in different environments. Despite the ongoing recovery of large carnivores in Europe and North America, few European studies on ungulate mortality are available from areas where both large carnivores and human hunters are present. We applied known fate models to estimate cause-specific mortality rates and Cox proportional hazard models to estimate the effects of environmental covariates on mortality risks of 330 radio-collared roe deer ( Capreolus capreolus ) (1995–2005) along a gradient in roe deer abundance in south-eastern Norway. The study area is characterized by the presence of human hunters, Eurasian lynx ( Lynx lynx ), red foxes ( Vulpes vulpes ) and occasionally wolves ( Canis lupus ). The main mortality causes were: hunter harvest, predation by lynx, predation by foxes (on fawns) and others (including wolves, dogs, diseases, vehicle collisions and accidents). The individual risk of roe deer being killed by lynx or by foxes was differently affected by covariates. In keeping with the specialist foraging behavior of lynx, predation risk by lynx decreased with increasing roe deer abundance. Conversely, consistent with the opportunistic habits of red fox, the risk of being preyed upon by foxes, tended to increase with increasing roe deer abundance, although the pattern was not so marked. Human hunters did not adjust their killing rate to changing roe deer abundance and annually harvested between 11 % and 28 % of the population according to different sexes and age classes.


INTRODUCTION
European wild ungulates represent one of the best studied groups of mammals from the point of view of demography; however the best data come from just a few long-term studies. These are mainly from islands or fenced areas where bottom-up effects dominate (e.g., red deer [Cervus elaphus] on the island of Rum, Coulson et al. 1997, Soay sheep [Ovis aries] on St Kilda, Clutton-Brock et al. 1992, Milner et al. 1999 in French reserves, Gaillard et al. 1993Gaillard et al. , 1997; see Gaillard et al. 1998a for a review). However, in most regions a range of predators (both human and non-human) usually occur, especially as large carnivore recovery is proceeding in many parts of Europe and North America. As survival rates might vary greatly and might strongly contribute to variation in population growth rate in such environments (Nilsen et al. 2009a), comparative studies of ungulate mortality patterns in such areas are important for our ability to sustainably manage ungulates where both large predators and human hunters are present.
European roe deer are the smallest ungulate indigenous to Europe, are found throughout the continent, from Mediterranean to boreal areas, in human-dominated and wilderness areas, and are prey for many carnivores, ranging from red foxes (Vulpes vulpes) to wolves (Canis lupus). They have been widely studied in many areas (mainly without predators or human hunting), resulting in good knowledge of the range of demographic parameters that can be expected, making it possible to recognize shifts in parameters in areas where predation and human hunting coexist. Roe deer are characterized by low sexual dimorphism, which is expected to lead to similar life-history traits and particularly similar survival for the two sexes. Despite this, lower male adult survival has been documented in several roe deer populations (see Gaillard et al. 1998b for a review). While several studies did not find any sexual bias in fawn survival (Andersen and Linnell 1998, Jarnemo et al. 2004, Panzacchi et al. 2008, males were found to be at higher risk of neonatal predation from foxes in one study (Aanes and Andersen 1996).
The roe deer population we focus on occurs along a population abundance gradient (due to differences in climate and habitat suitability across the study area), and coexists with several predators; the most important are Eurasian lynx (Lynx lynx), red foxes and human hunters. Wolves and free-ranging dogs (Canis familiaris) occur sporadically in the study area and have been documented killing fawns only in a few instances. Previous studies from this study system have suggested a strong potential for top-down control on the roe deer population (Nilsen et al. 2009a, b, Gervasi et al. 2012. Wherever ungulates coexist with large predators, predation accounts for most of the mortality among cervids in multi-predator systems (see Jędrzejewski et al. 2011 for a review). However, mortality risk from different predator species varies in relation to age class, sex and season. For a given predator, the impact of predation also varies according to the abundance at which both predators and prey occur (e.g., Sinclair 1989, Messier 1991. For example, Eurasian lynx, when sympatric with roe deer, are roe deer specialist (Jędrzejewski et al. 1993, Jobin et al. 2000, Odden et al. 2006b and their kill-rates show only limited variation across a wide range of prey abundance (Breitenmoser and Haller 1993, Okarma et al. 1997, Molinari-Jobin et al. 2002, Nilsen et al. 2009b. Consequently, assuming that lynx density remains stable, lynx predation could be expected to be higher at low roe deer abundance and in areas with unfavorable environmental conditions (Melis et al. 2009, Melis et al. 2010). On the other hand, red foxes have been well documented to prey upon fawns (Cederlund and Lindströ m 1983, Linnell et al. 1995, Aanes and Andersen 1996, Jarnemo et al. 2004, but due to their generalist and opportunistic feeding behavior, they have been observed to specialize only when roe deer occur at high abundance and the higher encounter probability makes prey switching behavior profitable (Panzacchi et al. 2008). European hunters harvest roe deer for a wide range of reasons, focusing on meat, trophies, or population control under different circumstances. Therefore, their selection of age and sex classes and relative effects will vary with different ecological and cultural circumstances. According to Ginsberg and Milner-Gulland (1994) trophy hunting usually shows extreme selection for adult males in ungulates. In contrast lynx have not shown any selection for the two sexes when hunting roe deer (Andersen et al. 2007). In Norway, harvest management goals for roe deer mainly aim to maintain stable populations and hunting success for roe deer is low and not strongly regulated by quotas (Grøtan et al. 2005, Mysterud and Østbye 2006, Melis et al. 2010); therefore we expected that hunting mortality would be similar across the abundance gradient.
In this paper we estimated annual causespecific mortality rates of roe deer according to sex and age class in south-eastern Norway (1995)(1996)(1997)(1998)(1999)(2000)(2001)(2002)(2003)(2004)(2005) along a gradient in roe deer abundance. We tested the following predictions: (1) the risk of fox predation for roe deer fawns is positively related to roe deer abundance; (2) the risk of lynx predation is negatively related to roe deer abundance; (3) whereas the risk for a given roe deer of being preyed upon by lynx and fox does not vary for the two sexes, the risk of being harvested by hunters is higher for males and adults.

Index of roe deer abundance
The number of roe deer harvested annually in each municipality included in the study area was divided by the area of forested habitat and was used as an index of roe deer abundance. This type of index has been previously used for roe deer in Norway (Herfindal et al. 2005, Mysterud and Østbye 2006, Nilsen et al. 2009a, Nilsen et al. 2009b, Melis et al. 2010 and it correlates well with other indices of population abundance (Grøtan et al. 2005, Mysterud andØstbye 2006). Three independent lines of arguments suggest that roe deer hunting bag statistics can serve as a crude measure of roe deer abundance in our study system. First, roe deer hunting quotas are rarely filled (quota filling 27% 6 18.9% [mean 6 SD]; see Grøtan et al. 2005) and therefore the number of harvested roe deer is likely to reflect changes in roe deer density rather than being an artifact of quotas (Grøtan et al. 2005, Mysterud and Østbye 2006, Melis et al. 2010. Second, annual variation in the number of traffic killed roe deer generally correlates well with harvest bag statistics at the county scale (r . 0.55 for 12 of 13 counties; see Grøtan et al. 2005). Finally, the number of roe deer sighted at supplementary feeding sites in the study area and the number of roe deer harvested in the same areas showed strong correlation (r ¼ 0.94, n ¼ 5 years; see Grøtan et al. 2005). To further assess the suitability of this index in our study system, we fitted cox proportional hazard models (see Statistical methods) and used model selection (Burnham and Anderson 2002) to investigate which variables among sex, age and roe deer abundance better described the variation in the risk for a given roe deer of being harvested by a hunter (Appendix). The first ranked model contained the main effect of age (fawns and yearlings vs. adults), sex and their interaction (Table 1), but none of the top-ranked models retained the effect of roe deer abundance, showing that there is no bias in our index. The risk for a given roe deer of being harvested by hunters was higher for adult roe deer compared to fawns and yearlings, and was higher for males than for females. As might be expected, the highest risk of harvest mortality was found for adult males.

Study area
The study area encompasses an environmental gradient (roughly north-south) in the counties of Hedmark, Akershus and Østfold in south-eastern Norway. The northern portion of the study area is characterized by several river valleys at around 200-300 m, separated by hills reaching to 700-800 m. The forest is mainly composed of Norwegian spruce Picea abies and Scots pine Pinus sylvestris and most of it has been logged and regenerated throughout the last 100 years. The roe deer abundance in this portion is generally lower than in the southern one (0.01-0.65 individuals harvested annually/km 2 of forested area). The southern portion of the study area includes patches of deciduous forest, represented mainly by birch Betula spp. and the landscape is more human-modified, with the forest fragmented by cultivated land and water bodies. The altitude is not higher than 300 m. Here roe deer occur at higher abundance (0.10-2.50 individuals harvested annually/km 2 of forested area). Throughout the study area roe deer coexist with moose (Alces alces), mountain hare (Lepus timidus) and forest birds such as black grouse (Tetrao tetrix) and capercaillie (Tetrao urugallus). The abundance of lynx, including dependent offspring, inside the study area has been estimated to be 0.3 individuals/100 km 2 in the northern part of the study area and 0.4 individuals/100 km 2 in the southern part (Odden et al. 2000, Odden et al. 2006a). The lynx population is subject to an annual recreational harvest, such that numbers were more or less stable during the study period. Single wolves were only occasionally present. Red foxes were widespread and their abundance was approximately 3 times higher in the southern portion of the study area (based on snow tracking indices; Panzacchi et al. 2008). Roe deer hunting takes place every year and the hunting season starts for bucks on the 10th of August and for adult females, yearlings and fawns on the 25th of September and ends for all sexes and age classes on the 23rd of December. During the study period the climate was relatively constant, with a 10-year sequence of mild winters. There is, however, an underlying effect of climate as the environmental gradient along which our study was conducted is largely a product of a climate gradient (mild coastal to more continental inland climate zones).

Deer capture and radio-telemetry
Between 1995 and 2005 a total of 330 roe deer were captured and equipped with radio-transmitters (Televilt Int.) . Roe deer were captured in winter at artificial feeding sites using box traps, drop nets or canon nets. The animals were manually restrained (without use of anesthetic), aged, sexed, ear-tagged and radio-collared before being released. No complications were ever detected as a result of collaring, although two animals died during capture related activities and were therefore removed from the analyses. We distinguished between fawns (,12 months) yearlings (between 12 and 24 months), and adults (.24 months). Fawns were easily recognized based on size and tooth eruption patterns when captured during winter, and many were originally captured as neonates. Most yearlings were unambiguously classified on the basis of having been caught earlier as fawns. However, during winter capture we felt comfortable in separating yearlings (20-23 months of age at this time) from adults based on a range of morphological criteria including size, general shape, tooth wear, antler development, absence of following fawns, appearance of nipples. The roe deer capture was conducted by experienced people who had frequent opportunities to calibrate their classification from capture of known age animals and from watching the animals at the feeding sites. In summer roe deer fawns were captured soon after birth and provided with expanding radio-collars. The radio-transmitters were mortality-sensitive and the status of the animals was checked at least once a week, but usually more often. If a mortality signal was detected, the animal was immediately located and the cause of death was assessed directly in the field or, in dubious cases, autopsies were made at the Norwegian Veterinary Institute, Oslo.

Mortality causes
We knew the cause of mortality of 228 animals, after censoring the ones where the radio transmitter stopped working, to avoid bias towards human-caused mortality. Of these we knew the age class at death (adults, n ¼ 105; yearlings, n ¼ 28; fawns, n ¼ 95) and the sex (males, n ¼ 127; females, n ¼ 101). We defined the main causes of mortality as predation by lynx, fox, wolf, and domestic dog, hunter harvest, disease / starvation, car (vehicle collisions), accident (drowning, falling from a cliff or into a ditch), grass-cutter and stillbirth (fawns only). For the analysis of cause-specific mortality rates the mortality causes were grouped in 4 main categories: lynx, fox, hunting and other. All statistical tests were twotailed, with significance value set at 0.05.

Estimation of cause-specific mortality rates
Radio-telemetry records were converted into monthly encounter histories with staggered entry (White and Burnham 1999) for all the 330 radiomarked animals. Each month an animal could be classified as 'alive', 'dead' or 'censored' (i.e., 'lost' due to radio failure or other reasons; Cooch and White 2006 ). The encounter history of an individual started in May of the year of capture (coded as 'censored' until month of capture, according to Pollock's staggered entry design; Pollock et al. 1989) and continued throughout 12 months. The encounter history file also included information about the year, the sex and the age class of each individual. If an individual was captured as a fawn, it would be reclassified as a yearling in May of the subsequent year and as an adult in May when entering its third year of life, and it would then keep adult status throughout its life. Similarly, an individual captured as a yearling would be reclassified as an adult in May v www.esajournals.org of the subsequent year and then keep its adult status throughout. We initially explored the possibility of applying the non-parametric cumulative incidence function (NPCIFE; Heisey and Patterson 2006), which has recently been used to estimate cause-specific mortality rates for willow ptarmigans (Sandercock et al. 2011) and wolves (Liberg et al. 2012). However, if sample sizes are small to moderate at the left tail of the age distribution (i.e., left truncation), estimates of survival might be biased downwards (Woodroofe 1985, Tsai 1988. Initial exploration of our data indeed suggested this to be the case for fawns, because of the low sample size and the high frequency of fox predation during the first month after capture. We thus used the Heisey and Fuller estimator (Heisey and Fuller 1985) to calculate cause-specific mortality rates of roe deer of different sex-and age classes with the software Micromort 1.3 (Heisey and Fuller 1985). All other analyses were performed using the R 2.11.1 software (R Development Core Team 2010). To examine how age, sex and roe deer abundance affected the cause-specific mortality risk, we applied Cox proportional hazard models (Lunn and McNeil 1995, Heisey and Patterson 2006, Murray 2006, and stratified according to cause of mortality as described in Heisey and Patterson 2006. Model selection was based on AICc (Burnham and Anderson 2002). Prior to model selection, we assessed to which extent the assumption of proportional hazards were met with model diagnosis based on scaled Schoefield residuals. Cox proportional hazard models were implemented with functions in R add-on library Survival (Therneau 2012), and model diagnosis with the cox.zph function in the same software library (Fox 2002).

RESULTS
The annual cause-specific mortality rates according to age class and sex are presented in Table 2. The mortality rate due to lynx predation was very similar among age classes (0.13-0.15), and there were no statistically significant differences among age classes and sexes. Fox predation was exclusively directed towards fawns and accounted for a mortality rate of 0.26 (60.041 SE) among this age class. No sex difference was detected in fawn mortality due to fox (Table 2). There was a statistically significant difference in harvest mortality between adult and yearling females. When considering differences between sexes, adult males had three times higher rates caused by harvest than adult females (0.23 6 0.048 vs. 0.08 6 0.024) ( Table 2). Annual survival rates (sexes pooled) were 0.27 6 0.043 for fawns, Table 2. Cause-specific mortality rates with confidence intervals (CI) as estimated using the Heisey and Fuller method (1985) based on monthly intervals for radio-marked roe deer in south-eastern Norway (1995)(1996)(1997)(1998)(1999)(2000)(2001)(2002)(2003)(2004)(2005). The estimates whose CI do not overlap across age classes and sexes are shown in bold.  Table 2). We then examined how cause-specific mortality risks were affected by sex, age and roe deer abundance. In support of our first prediction, the risk of fox predation increased with increasing roe deer abundance (Table 3a, Fig. 1a) in a curvilinear manner. Although there was uncertainty in the results of the model selection, the best performing model containing the linear and quadratic effect of roe deer abundance was about 2.6 times more likely than the second ranked model. Also our second prediction was support-ed, as the risk of lynx predation declined as roe deer abundance increased (Table 3b, Fig. 1b). In addition fawns had a higher risk of lynx predation than older individuals (Table 3b, Fig.  1b). The best performing model described above was 2.2 times more likely than a model where sex was also included. Consequently, our third prediction was also supported as we did not find support for differences in mortality risks between the sexes for lynx and fox predation, whereas the hunters indeed selected adult males over adult females (Table 1). In Fig. 2 we report for comparison the survival rates of roe deer for the three age classes estimated by this study and other published studies on roe deer mortality Table 3. Estimates of the first ranked model (as selected by AIC c ) explaining the risk of (a) being preyed upon by fox and of (b) being preyed upon by lynx for roe deer in south-eastern Norway (1995)(1996)(1997)(1998)(1999)(2000)(2001)(2002)(2003)(2004)(2005) (1995)(1996)(1997)(1998)(1999)(2000)(2001)(2002)(2003)(2004)(2005). (b) Prediction lines with confidence intervals (shaded areas) for the first ranked model describing roe deer fawns (black line) and older (yearlings þ adults, dotted line) predation risk by lynx as function of roe deer abundance in south-eastern Norway (1995)(1996)(1997)(1998)(1999)(2000)(2001)(2002)(2003)(2004)(2005).
v www.esajournals.org conducted in areas with and without predation and / or human hunting. This figure suggests that when both human harvest and predation occur, the estimates of roe deer survival are considerably lower than when harvesting and/or predation are absent.

DISCUSSION
In our study area, where roe deer occur across a spatial gradient in abundance and are sympatric with several predators (mainly lynx, foxes and humans) we found that: (1) the risk of lynx predation was negatively related to roe deer abundance; (2) the risk of fox predation was positively related to roe deer abundance; (3) the risk of being preyed upon by lynx and foxes was similar for the two sexes. Moreover, the risk of being harvested by hunters did not vary with varying roe deer abundance, which supports the assumption that there is no relationship between these two variables and therefore that harvest statistics is a good index of roe deer abundance.
The risk of being harvested by hunters was higher for adult males than adult females and for yearling females than adult females.
The main causes of mortality for the three age classes and two sexes were predation and hunting. Gaillard et al. (1998a) showed that populations of large herbivores most often have a high and stable adult survival, while the survival of juveniles is generally lower and more variable in time, depending on several densitydependent and density-independent factors. Our results generally support these theories, but the estimated survival rates of yearlings and adults were very much lower than those commonly reported for large herbivores in Europe (e.g., Clutton-Brock et al. 1992, Gaillard et al. 1993, Coulson et al. 1997, Milner et al. 1999, Focardi et al. 2002, Cobben et al. 2009). Another published study on roe deer survival conducted in Tredozio (Italy) found different results (Focardi et al. 2002). In this population the annual estimate of fawn survival (0.38 6 0.07) was higher than ours (0.28 6 0.043), and also the estimate of adult Fig. 2. Annual survival rates of roe deer in presence of both large predators and human hunters, only large predators or hunters, in absence of predation and human hunting in Southern Norway (present study), Bavarian Forest National Park in presence and absence of lynx (Heurich et al. 2012), Tredozio (Focardi et al. 2002), Storfosna (Cobben et al. 2009), Chizé (Festa-Bianchet et al. 2003 and Trois Fontaine (Gaillard et al. 1993). The results from Trois Fontaine include only prime aged individuals.
v www.esajournals.org survival was much higher (0.90 6 0.07) compared to ours (0.51 6 0.036) (Fig. 2). Adult survival probabilities similar to the one recorded in Tredozio have been reported for two enclosed populations in France which were not subjected to predation (0.85 for males and 0.95 for females; Gaillard et al. 1993;Fig. 2). However, in the Italian study radio-collared animals were not harvested and the most important roe deer predator, the Eurasian lynx, was absent. Therefore predation was mainly due to foxes and directed towards fawns, which might explain why the estimates of adult survival were similar to the ones from sites without large predators. In our study, both human hunting and lynx predation are amongst the most important agents of mortality and therefore our estimates of annual adult survival were lower. The annual survival of adult roe deer was comparable, for example, to the values found for a population of black-tailed deer Odocoileus hemionus in Washington State (USA) (Bender et al. 2004), subjected to both human harvest and natural predation (by the mountain lion Puma concolor L.). The only study conducted in Europe on roe deer mortality in presence of both human hunting and lynx predation (Heurich et al. 2012) found that, after recolonisation of the study area by lynx, prime age roe deer survival was almost 20% lower than in absence of large predators (0.61 against 0.79).
As expected, lynx did not appear to select among sexes and age classes whereas predation by foxes did not differ between sexes but was exclusively directed towards fawns. Our result, in accordance to those found by Jarnemo et al. (2004) and Panzacchi et al. (2008), did not support the results of a higher predation risk from foxes on male fawns with respect to female fawns (Aanes and Andersen 1996). Moreover, for female roe deer human harvesting focused more on yearling than adult. Hunters killed more males than females; therefore trophy hunting seems to be a driver of prey selection by humans in our study area.
The inverse density-dependence in the risk of lynx predation on roe deer has been previously described by Jędrzejewska and Jędrzejewski (1998) in the multipredator system of Białowieża Primeval Forest (Poland). The annual mortality rate due to foxes was over 26%. This value is rather comparable to the estimate of 42% fawn mortality by foxes from an island in central Norway (Aanes and Andersen 1996). Our results reflect the different feeding tactics of lynx and foxes: lynx are efficient roe deer specialist (Breitenmoser and Haller 1993, Linnell et al. 1996, Okarma et al. 1997, Molinari-Jobin et al. 2002 and are not expected to substantially change their kill rate with varying roe deer abundance, and therefore will only have a limiting, or anti-regulatory effect. This has been clearly shown in the same study area by Nilsen et al. (2009b), although the functional response curve in that analysis was also strongly influenced by season and differences in lynx social status. Conversely, foxes are typical generalist predators with opportunistic diets and are likely to specialize in predating fawns only when their occurrence in the environment makes it worthwhile to actively spend time in their search. This supports previous findings in this study area, showing that foxes adjusted their functional response according to habitat characteristics and profitability of both main and alternative preys (Panzacchi et al. 2008(Panzacchi et al. , 2009). However, the fact that the number of foxes was three times higher in the southern part of the study area might have contributed to this pattern.
The annual survival rates that we report in this study are much lower than those previously documented in roe deer studies where predation and / or hunting are absent (Gaillard et al. 1993, Cobben et al. 2009, Heurich et al. 2012Fig. 2). Although juvenile survival rates were generally lower than those of older individuals also in our study populations, the estimated survival rates were much lower than those reported elsewhere from predator free areas. This clearly will have both a range of practical management and life-history implications (Nilsen et al. 2009a). Our study design does not allow us to formally explore the additive vs. compensatory nature of predation and hunting. However, by comparing the rates observed in our study to those previously reported for the species across a range of conditions (our mortality rates for all other causes were similar to those for total mortality in other study sites) it is highly likely that both predation and hunting were largely additive in our study. This implies that the recovery of large carnivores, especially lynx ), across Europe is going to induce a change in roe v www.esajournals.org deer mortality schedules, increasing total mortality, especially for yearling and adult females. If not accounted for this could cause declines in harvested roe deer populations, especially in marginal areas (Melis et al. 2009(Melis et al. , 2010. This study underlines how important the impact of predation from human and mammalian predators can be for roe deer survival. Therefore, while we acknowledge how crucial the studies conducted in simplified model environments have been for advancing our understanding of roe deer demographics (Gaillard et al. 1998a, Cobben et al. 2009), we feel that it is important to increasingly supplement these studies with others conducted under more complex conditions if we are to provide data relevant for management and conservation.

SUPPLEMENTAL MATERIAL APPENDIX
(2) preyed upon by fox; (3) preyed upon by lynx. DAIC c ¼ difference in AIC c between the best and the actual model; x i ¼ Akaike's weights, i.e., normalized likelihoods of the models. In bold is reported the first ranked model.