<figure><img src="image_1.png" /></figure>

ELSEVIER

Contents lists available at [ScienceDirect](https://www.sciencedirect.com)

Global Ecology and Conservation

journal homepage: [www.elsevier.com/locate/gecco](https://www.elsevier.com/locate/gecco)

<figure><img src="image_2.png" /></figure>

<figure><img src="image_3.png" /></figure>

Check for
updates

# Cumulative mortality effects on roe deer population dynamics in the boreal forest: Searching for pathways of population persistence

Cécile A.E. Carpentier<sup>a, *, 1</sup>, Marco Heurich<sup>a, b, c, 2</sup>, Olivier Gimenez<sup>d, 3</sup>, Olivier Devineau<sup>a, 4</sup>, John D.C. Linnell<sup>a, e, 5</sup>

## ARTICLE INFO

### Keywords:

- Additive mortality

- Matrix model

- Population dynamics

- Predator-prey

- Roe deer (*Capreolus capreolus*)

- Wildlife management

## ABSTRACT

Ungulate populations can exhibit various growth patterns, which are influenced by factors such as predation and resource availability. Favourable environments can lead to irruptive growth, resulting in resource depletion. However, additional pressures from predation, and hunting can potentially impact population development leading to declines or even local extinctions. This study uses simulation models to explore the potential impact of multiple mortality sources on roe deer populations. We develop an age-structured, two-sex demographic matrix model for roe deer, which we parameterise with empirical demographic estimates obtained from published studies in Norway. We develop scenarios to assess the influence of mortality sources such as hunting, predation by lynx and red foxes, and environmental stochasticity on roe deer population dynamics. When simulating favourable environments without predation, roe deer populations tended to erupt due to the species' rapid reproductive capacity. However, additional sources of mortality, such as predation or harvest, lead to severe population declines, and even to quasi-extinction, especially when occurring in combination. Environmental stochasticity such as periodic severe winters with heavy snowfall reduces the growth rate and population densities even further. On the other hand, accounting for some form of spatial heterogeneity through immigration and refuges stabilised populations, with a reduced risk of quasi-extinction. Our results provide meaningful insights into the properties of this system allow implications for the management and identify areas where further exploration is needed.

<footnote><sup>a</sup> Department of Forestry and Wildlife Management, Faculty of Applied Ecology, Agricultural Sciences and Biotechnology, University of Inland Norway, Koppang, Norway</footnote>

<footnote><sup>b</sup> Bavarian Forest National Park, Grafenau, Germany</footnote>

<footnote><sup>c</sup> Wildlife Ecology and Conservation Biology, Faculty of Environment and Natural Resources, University of Freiburg, Freiburg, Germany</footnote>

<footnote><sup>d</sup> CEFE, Univ Montpellier, CNRS, EPHE, IRD, Montpellier, France</footnote>

<footnote><sup>e</sup> Norwegian Institute for Nature Research, Lillehammer, Norway</footnote>

<footnote>* Corresponding author.
E-mail addresses: cecile.ae.carpentier@gmail.com, cecile.carpentier@inn.no (C.A.E. Carpentier).</footnote>

<footnote n="1">¹ ORCID ID: 0009-0007-6092-524X</footnote>

<footnote n="2">² ORCID ID: 0000-0003-0051-2930</footnote>

<footnote n="3">³ ORCID ID: 0000-0001-7001-5142</footnote>

<footnote n="4">⁴ ORCID ID: 0000-0002-7625-2816</footnote>

<footnote n="5">⁵ ORCID ID: 0000-0002-8370-5633</footnote>
# 1. Introduction

The ongoing recovery of large herbivores and large carnivores across the European continent which started during the 20th century (Apollonio et al., 2017; Chapron et al., 2014; Crétois et al., 2021; Linnell et al., 2020) has led to an urgent need for knowledge on their predator-prey interactions and population dynamics as a prerequisite for guiding their management in human dominated landscapes. The two main predator-prey systems that are emerging are the system with wolves (*Canis lupus*) feeding on moose (*Alces alces*), red deer (*Cervus elaphus*), roe deer (*Capreolus capreolus*), and wild boar (*Sus scrofa*), and the Eurasian lynx (*Lynx lynx*) system where roe deer are the major prey (in addition to localised systems where chamois *Rupicapra rupicapra* and reindeer *Rangifer tarandus* can be important). These systems have very different properties (Gervasi et al., 2012). In general, there has been far more research into group-living predators' (e.g. wolves) relationships with prey compared to solitary predators' (e.g. lynx) impact on prey. Furthermore, the role of roe deer as a prey for large carnivores is much less studied than that of moose or red deer despite roe deer being the most widespread wild ungulate in Europe.

Roe deer have long been a model species for studying large mammal population dynamics, with research spanning multiple European countries (e.g., studies from the UK, Denmark, France, Sweden and Norway; Gill et al., 1996; Hewison and Gaillard, 2001; Strandgaard, 1972; Gaillard et al., 1993; Hewison et al., 1996; Focardi et al., 2002; Kjellander et al., 2004; Andersen and Linnell, 1997, 1998; Cobben et al., 2009). These studies were central in the last part of the 20th century when researchers were trying to tease apart the roles of behavioural (social-regulation via dispersal), demographic (density-dependent) and environmental (weather, habitat) factors that influence population dynamics. This was especially relevant in the context of studying reproduction, and managing irrupting, high-density populations in areas considered "optimal" roe deer habitat (i.e., Central Europe). These studies typically took place in simplified ecosystems (islands, fenced areas, or highly human modified landscapes) where the roles of humans and mammalian predators as causes of mortality were limited. However, as the research agenda developed, priorities need to shift towards more applied questions, such as assessing the impact of hunting and recovering predator populations (Chapron et al., 2014; Crétois et al., 2021; Melis et al., 2009). Therefore, it is becoming necessary to also study roe deer population dynamics in more complex systems where they face human-induced mortality (hunting and traffic collisions) and predation by various predators (red foxes (*Vulpes vulpes*) preying on fawns, Eurasian lynx, and wolves targeting all age classes). In such settings, the research question change focus to studying mortality from different sources, preventing population declines, or managing recovering roe deer populations in the margins of their distribution (Torres et al., 2012). Furthermore, research has demonstrated that roe deer exhibit considerable mobility through natal dispersal and migration (Cagnacci et al., 2011; Hewison et al., 2021). Although highly variable across populations, this has a significant potential to influence population dynamics (Linnell et al., 1998) as does spatial heterogeneity in demographic rates (Pettorelli et al., 2003). In addition, periodic severe winters can greatly impact roe deer survival (Grøtan et al., 2005), particularly through increased predation and starvation (Cederlund and Lindström, 1983).

Studying complex predator-prey systems is challenging and expensive because the need to focus on cause-specific mortality has often required telemetry methods. The last few decades have seen several studies that have collected parameter estimates from multiple populations and quantified significant impacts of both human and mammalian predation (Andersen et al., 2007; Gervasi et al., 2012; Heurich et al., 2012; Jarnemo and Liberg, 2005; Lone et al., 2014; Melis et al., 2013). Analysis of long-term population indices in Sweden has also shown that predation can have dramatic and large-scale effects on roe deer population development (Andrén and Liberg, 2015, 2023; Gervasi et al., 2012). Moreover, the accumulation/addition of different sources of mortality on roe deer has shown that it can reach fairly high levels, substantially impacting potential population growth (Melis et al., 2010).

However, despite this accumulation of empirical demographic data, we are not aware of any attempt to explore the dynamics of these complex predator-prey systems through simulations and population projections. By exploring the impacts of various sources of mortality on roe deer populations in a simulation environment, it is possible to improve our understanding of the potential role of top-down mechanisms on population dynamics, providing a better understanding of the impact of different wildlife management strategies and predation pressures. Simulation studies are essentially explorative tools that represent a link between empirical studies that collect parameter estimates in the field and the production of testable hypotheses about the properties of the system that can be tested in future studies.

In this work we use published estimates of demographic parameters from multiple short-term studies of roe deer reproduction and survival across a gradient of Norwegian landscapes (Andersen and Linnell, 2000; Melis et al., 2013) as the basis for a stochastic population projection model (Caswell, 2000). This model allows us to simulate the trajectories of model roe deer populations over the medium term, projecting them under different predation and mortality scenarios to investigate their potential effects. Existing published data (Melis et al., 2010) indicate that mortality rates for roe deer from cumulative sources can be remarkably high. Therefore, in addition to exploring the relative and additive effects of major mortality factors, we also examine how previously unconsidered factors (like immigration (Linnell et al., 1998), the existence of individual spatial heterogeneity in risk exposure (Gorini et al., 2012), and periodic severe winters (Grøtan et al., 2005)) could potentially modify the population trajectories. Specifically, we model the development of roe deer populations over 25 years, considering different combinations of mortality sources (red fox predation, lynx predation and human hunting), under scenarios where extreme snow events could occur, increasing the overall mortality. Additionally, we simulated population outcomes after 25 years with varying immigration rates and exposure heterogeneity of predation risk. Because roe deer management is concerned with both over-abundance and decline, we explored what drives populations to irruption or quasi-extinction, and what influences changes in population abundance and population growth rate. Our specific objectives include: (1) Integrating published demographic rates into a modelling environment to explore their implications for population development over time and assess the effect of cumulative mortality factors, (2) Explore how additional spatio-temporal factors
such as immigration and spatial heterogeneity in predation risk influence the properties of the system, and (3) Identify important areas for further work using more complex modelling approaches.

## 2. Materials and methods

### 2.1. Study site and parameter estimates

We develop and parameterise an age-structured, two-sex demographic model for a roe deer population. We consider the influence of different mortality sources including hunting, predation (lynx and red fox) and environmental stochasticity (extreme snow events) as different scenarios of population development.

This model focuses on demographic parameters from roe deer populations, where lynx and roe deer have been studied since the 1990's in southeast Norway through a series of projects at different sites (Gervasi et al., 2014; Melis et al., 2009; Nilsen et al., 2009; Panzacchi et al., 2009). To parametrize the model, we mainly used demographic estimates from published data from the former Hedmark, Akershus and Østfold regions (now Innlandet and Viken; Table A.1 in Appendix). In the absence of estimates for reproduction for the study region, we used values from a study on the Norwegian island of Storfosna from 1991 to 1994. We used reproductive estimates including the average litter size per 2-year-old and adult female and the breeding proportion of 2-year-old and adult females (Andersen and Linnell, 2000).

The mortality rate estimates are based on a 10-year telemetry study in southeast Norway from 1995 to 2005 (Melis et al., 2013 - Table 2). We used different cause-specific mortality estimates for each age and sex group. The main cause-specific mortality estimates were used, as grouped in Melis et al. (2013), such as red fox predation on fawns, lynx predation and hunting mortality on adults, in addition to other mortality (including wolf predation, diseases, traffic collisions and natural accidents).

Extreme snow events have also been found to negatively affect roe deer population survival and growth rate (Grøtan et al., 2005), but because such events are episodic it was not possible to obtain empirical estimates of their occurrence. To include more realistic environmental stochasticity into our model, we considered plausible scenarios introducing about two extreme snow events happening within 25 years (i.e. frequency of event of 0.08), which caused a 30% additional mortality across all age classes using a binomial distribution. These plausible values are based on a crude estimate we obtained from a visual examination of the inter-annual fluctuations observed in regional hunting statistics [https://www.ssb.no - Statistisk sentralbyrå - Small game and roe deer hunting] correlated with weather data (including that used by Grøtan et al. (2005)).

There is currently no data on the extent to which habitat heterogeneity affects roe deer vulnerability to lynx predation (see Panzacchi et al. (2010) for a roe deer - red fox case) although there is clear evidence of the potential for this based on differential habitat selection by lynx and roe deer (see Bouyer et al., (2015), Bouyer et al. (2015), Bouyer et al. (2015)). Because these effects are also likely to vary with landscape structure, we have modelled their potential effects using an arbitrary, but plausible, range of values.

### 2.2. Model structure

The roe deer population was modelled using a matrix population model structured by age and sex (Caswell, 2000). We used a birth-pulse model with a post-breeding census, which means that the population census is done right after a breeding pulse (i.e. birth event) (Fig. 1).

The model is separated into 7 different stages: 3 stages for males and 4 stages for females of different age classes. For both males and females, we considered the *juvenile stage* for the first year from birth to 12 months ($J_m$ and $J_f$ respectively), and the *yearling stage* from 13 to 24 months ($Y_m$ and $Y_f$ respectively). For females, we separated the *primiparous stage* ($P_f$) for females between 25 and 36 months old, when their first reproduction can take place, and the *adult stage* from 37 months and older ($A_f$). We considered males to be in their *adult stage* from 25 months and older ($A_m$). The total population abundance at time $t$ is described by a vector $n_{(t)}$ whose elements represent the abundances of each stage. Each time step in the simulation represents a year of projection for the population. To project the population over time, this vector $n_{(t)}$ is multiplied by the population projection matrix called $A$, containing survival rates and

<figure><img src="image_4.png" /><figcaption>**Fig. 1.** Roe deer life cycle with 7 stages, 4 stages for females: Juvenile, Yearling, Primiparous and Adult stages ($J_f$, $Y_f$, $P_f$ and $A_f$ respectively) and 3 stages for males: Juvenile, Yearling, and Adult stages ($J_m$, $Y_m$ and $A_m$ respectively) and the different Survival rates ($S_i$, with $i$ respectively for each stage) and Fertility rates ($F_i$, with $i$ respectively for each stage).</figcaption></figure>
reproductive rates to get the population abundance $n_{(t+1)}$ at the next time step: $n_{(t+1)} = A_{(n_t)} * n_{(t)}$ With

$$A_{(n_t)} = \begin{pmatrix} 0 & \rho F_Y & \rho F_p & \rho F_A & 0 & 0 & 0 \\ S_{J_f} & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & S_{Y_f} & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & S_{P_f} & S_{A_f} & 0 & 0 & 0 \\ 0 & (1-\rho)F_Y & (1-\rho)F_p & (1-\rho)F_A & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & S_{J_m} & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & S_{Y_m} & S_{A_m} \end{pmatrix} \quad \text{and} \quad n_{(t)} = \begin{pmatrix} n J_{f_{(t)}} \\ n Y_{f_{(t)}} \\ n P_{f_{(t)}} \\ n A_{f_{(t)}} \\ n J_{m_{(t)}} \\ n Y_{m_{(t)}} \\ n A_{m_{(t)}} \end{pmatrix}$$

Every year, each stage survives to the next with a survival probability called $S_i$ in the matrix (with i the stage at year t). In the model, the survival probability for each stage is calculated from various sources of additive mortality, which is calculated as $S_i = 1 - (m_{i1} + m_{i2} + m_{i3} + ...)$ with $m_{xi}$ the cause-specific mortality rate. Due to the presence of both sexes in the model, the offspring produced need to be allocated between male and female newborns based on the binomial distribution for sex ratio. The proportion of female offspring is denoted as rho ($\rho$). In the matrix, the reproduction is shown from the female yearling stage due to the post-breeding census. The population is then censused after reproduction, which means that, reproductive individuals also include yearling females that have survived to the next year and will then reproduce for the first time. The fertility rate for each reproductive stage ($F_i$) is calculated considering the proportion of females of each reproductive stage that survive until the next census ($S_i$), the proportion of females of each stage giving birth ($p$) and using the fecundity rate ($f$). The fecundity rate is calculated using the average number of offspring per female ($k$), the number of reproductive males and females ($N_m$ and $N_f$ respectively) and the maximum number of females a male can reproduce with ($h$, here $h = 10$, see parameter table in Appendix, Table A.1).

$$F_i = S_i \times p \times f \quad \text{with}$$

$$f = k \times \min_{1} \left( N_m / \left( \frac{N_f}{h} \right) \right)$$

We assumed a closed (no immigration or emigration), homogeneous (same initial vital rates among all individuals in the same stage) and well-mixed (equal opportunities to interact with others) population. The model does not explicitly consider dispersal or immigration, nor does it include any additional sources of mortality besides predation and hunting (as per Melis et al., 2013), which are considered fully additive, as is likely to be the case in the study from which we extracted the estimates. We explore the impact of varying the relative importance of each mortality source on the population development as a scenario (explained under the section “scenarios”). We implement a linear density-dependent function on fawn mortality to prevent exponential population growth. This function adds extra mortality as density increases, thereby reducing the recruitment of new individuals. We chose the parameter strength of the function such that the population growth slows until it reaches the carrying capacity. We set carrying capacity to 40 ind./km², which represents a reasonably high density, close to what has been observed on the Norwegian Island of Storfosna without predation and hunting (Andersen et al., 1995) and in the Swedish study site of Bogesund after hunting was stopped (Kjellander et al., 2004). In our model, the population is simulated on a 10,000 km² area, which is a meaningful scale to embrace predator-prey dynamics and movements of roe deer (Linnell et al., 1998). The initial density was set to 10 individuals per square kilometre.

## 2.3. Simulations

We projected the roe deer populations to explore how populations are affected by different mortality rates and sources. The matrix population model enables the projection of populations and the prediction of fluctuations in population size over time. This allows us to determine population density, growth rate, and the proportions of populations that reach quasi-extinction or the carrying capacity, called irruption. We project the population over 25 years and ran 500 stochastic simulations, with slightly different parameter values (see Table A.1 in Appendix for stochasticity). Each stochastic simulation starts with a variation of 10% around the initial number of individuals. To vary the proportion of individuals in each stage, we used the Dirichlet distribution (which is continuous and whose probability vector is equal to 1). The population projection matrix is set with a variation around the vital rate estimates. The uniform distribution was used for the proportion of breeding females, so the proportion vary randomly within the confidential interval of the estimate. The binomial distribution was used to set the value of the sex ratio at birth considering the number of new individuals born. The normal distribution was used to draw the average litter size per female (with a minimum of 1) and a truncate normal distribution for the mortality estimates. Hence, the parameter value stays within 0 and 1 (see parameter table in appendix, Table A.1).

When the population is getting close to extinction, we have set a quasi-extinction threshold to stop the simulation after the population is considered quasi-extinct. For all simulations, we consider different thresholds to represent the state of the population. We consider **quasi-extinction** to be a state where a population goes below one-twentieth of its initial size (which corresponds to a density of 0.5 ind./km²). A population is considered in a state of **irruption** if it reaches the carrying capacity allowed by the density dependence function (i.e. 40 ind./km²). Everything between these two states is then in a quasi-steady state, which we call “stable”.

## 2.4. Scenarios

We called a scenario the simulation of a roe deer population exposed to one or more specific sources of mortality. We performed simulations across five main scenarios, each including a specific set of mortality sources. Our study includes two types of simulations: One explores the influence of different mortality factors on population trajectories based on empirical field-based parameter estimates,
and the other assesses the potential effects of immigration and refuges against predation on the population state after 25 years. For the latter approach, two gradients were used, one to take into account spatial heterogeneity creating refuges from predation and the other to represent an increasing amount of immigration into the population.

In the first type of simulations, the different additive mortality sources were switched on or off to expose populations to one or several sources, including red fox predation, lynx predation and hunting harvest in a first instance. We then added periodic extreme mortality events such as those that could be caused by deep snow events to the same set of scenarios as mentioned previously.

In the second type of simulations, we performed an additional set of simulations, including a theoretical rate of immigration and an additional adult survival to represent spatial refugees from predation risk. In one set of simulations, we introduced a fixed number of individuals into the population each year, evenly split between male and female yearlings, to represent potential immigration. This immigration was simulated across a gradient of number of added individuals to account for a range of immigration rates. Simultaneously, we reduced lynx predation on adult males and females equally to account for the effect of potential spatial refuges. We simulated this predation reduction along a gradient as well, representing different degrees of lowered predation pressure. For each combination of yearling immigration (shown as an additional percentage of the initial population) and additional adult survival (shown as percentage reduction of lynx predation mortality), we ran 500 stochastic simulations as in the previous set of scenarios. After 25 years, we established the proportion of stochastic simulations in which the population was in each state (irruption, stable or quasi-extinction as defined above).

All simulations were run in program R (v4.2.2; R Core Team 2022). We used the package *gtools* (Bolker et al., 2022) to draw the Dirichlet distribution in R, and the *truncnorm* package (Mersmann et al., 2023) for the truncated normal distribution. The R files to run the model and the data produced by the model are available via the GitHub repository (https://github.com/CecileAEC/MPM_roedeer).

<figure><img src="image_5.png" /><figcaption>**Fig. 2.** a, b, c, d, e: Roe deer population density [N/km²] (left column), population growth rate [‰] (middle column), and proportion of simulations [%] (right column) in different states (i.e. quasi extinction, stable and irruption) at each time step of the 25-year simulation under different scenarios: for each row, the pictograms indicate the additive mortality sources included in the model. For both left and middle column plots, the dark green line shows the median of the stochastic simulations, the dark green area represents the 25th and 75th quantiles, the light green area represents the 5th and 95th percentiles, and the dotted line represents the minimum and maximum of the simulations at each time step. The black dotted line on the growth rate plots shows the stable growth rate (lambda = 1).</figcaption></figure>
# 3. Results

## 3.1. Additive effects of mortality factors

When no predation or hunting was included in the model the average population density increased and reached a plateau at the pre-set carrying capacity (approx. 40 ind./km²) within approximately 10 years (Fig. 2a). The average population growth rate then fluctuated around 1, meaning that populations remained stable. In this scenario, none of the stochastic simulations resulted in populations reaching quasi extinction, while half of them reached irruption (i.e. carrying capacity). When only red fox predation was included in the model, the populations grew before stabilising at approx. 27 ind./km² and it took around 10 years to plateau and reach a stable growth rate of 1. In this scenario, none of the simulations resulted in populations reaching quasi extinction or irruption (Fig. 2b). When red fox predation was replaced by lynx predation, the average population density slowly decreased, but did not reach quasi-extinction by the end of the 25 years simulation period (Fig. 2c). In the scenario including both red fox and lynx predation, the populations started decreasing slowly by about 10% per year. In this scenario, 0.4% of the stochastic simulations resulted in populations reaching the quasi-extinction threshold within 15 years, although 62% reached this threshold within 25 years (Fig. 2d). In the final scenario, adding hunting mortality and predation from both red foxes and lynx led to faster population declines. On average, the population density decreased fast during the first 5 years and the population growth rates stabilised around 0.75 (i.e. 25% decrease each year). Quasi extinction was first reached by populations within 6 years, and all of them reached it within 10 years (Fig. 2e).

<figure><img src="image_6.png" /><figcaption>**Fig. 3.** a, b, c, d, e: Roe deer population density [N/km²] (left column), population growth rate [‰] (middle column), and proportion of simulations [%] (right column) in different states (i.e. quasi extinction, stable and irruption) at each time step of the 25-year simulation with 0,08 probability of extreme snow event under different scenarios: for each row, the pictograms indicate the additive mortality sources included in the model. For both left and middle column plots, the dark green line shows the median of the stochastic simulations, the dark green area represents the 25th and 75th quantiles, the light green area represents the 5th and 95th percentiles, and the dotted line represents the minimum and maximum of the simulations at each time step. The black dotted line on the growth rate plots shows the stable growth rate (lambda = 1).</figcaption></figure>
## 3.2. Effect of environmental stochasticity

For all scenarios, the model showed a higher proportion of the simulations experiencing lower population growth rates and densities. For the first scenario with neither predation nor hunting, the average growth rate stabilised around 1, but only 25 % of the populations reached irruption compared to 50 % in the previous scenarios when no extreme snow event was added to the model. This indicates that most populations stabilise at a lower density (Fig. 3a). In the scenario where only red fox predation occurred, no population reached quasi extinction, although more showed lower growth rates and densities due to extreme snow events (Fig. 3b). When we removed fox predation, but included lynx predation in the model, some populations reached quasi extinction within 13 years, but only 8 % reached it within 25 years (Fig. 3c). For the scenario with predation from both red fox and lynx, the first populations to reach quasi extinction did so within 5 years, and 93 % reached quasi extinction within 25 years (Fig. 3d), recalling that in the equivalent scenario without extreme snow events, only 50 % of the populations reached quasi extinction within 25 years (Fig. 2d). In the last scenario with hunting in addition to predation from red fox and lynx, the first populations reaching quasi extinction did so within only 2 years, and all of them reached quasi extinction within 10 years (Fig. 3e).

## 3.3. Stabilising effects of immigration and refuges

The gradient colour in the plots shows the proportion of the stochastic simulated populations that reached quasi extinction or remained stable. None of the simulations reached irruption in the range of our gradients, so we have truncated the y-axis up to 3 % additional yearlings each year to reveal more details at the lower end, comparing the proportion of the population states between quasi-extinction and stable states only.

In the scenario with lynx predation only, all populations were already stable before adding immigration or including refuges from lynx predation (Fig. 4c, Fig. 2c). In the scenario with both red fox and lynx predation, half of the populations reached extinction before adding extra immigration and including refuges from lynx predation (Fig. 2d). When reducing lynx predation on adult roe deer by 50 % or by adding only 0.5 % of the initial population as yearling immigrants (representing a net immigration of 5000 additional individuals each year), all populations remained at a stable state (Fig. 4d). In the last scenario where red fox predation and hunting are also included in addition to lynx predation, all populations went to quasi extinction (Fig. 2e). This remained the case even with completely relaxing the lynx predation (i.e. no lynx predation), meaning that red fox predation and hunting are sufficient to drive all populations to quasi extinction. Still, the proportion of populations at stable state increased when adding up to 25 % of the initial population as yearling immigration (which corresponds to a net immigration of 25000 individuals each year). If we completely relaxed the lynx predation, the required yearling immigration to prevent simulations from reaching quasi extinction could be reduced to approximately 20 % (Fig. 4e).

# 4. Discussion

This study demonstrates the potential negative cumulative impacts of predation from human and mammalian predators on the dynamics of roe deer populations simulated conditions reflecting the boreal forest ecosystems of southeast Norway. While high-density roe deer populations and irruptive growth rates dominated the early decades of their study across Europe, our models demonstrate the necessity of gaining a better understanding of predation and especially exploring the mechanisms that might buffer local populations against decline and quasi-extinction. Our simulated exploration of immigration and refuges from predation indicate the potential of these mechanisms.

<figure><img src="image_7.png" /><figcaption>**Fig. 4.** c, d, e: Proportion of roe deer population at quasi-extinction after 25 years over a gradient of yearling immigration and a gradient of lynx predation magnitude under the last 3 scenarios: c) lynx predation only, d) lynx and red fox predation, e) lynx and red fox predation, and hunting mortality. The colours represent the state proportion for each combination in the gradient. From dark red for 100 % of the simulated populations reaches quasi extinction to white where all simulations stay at a stable state.</figcaption></figure>
When the environment is favourable, with no predation or density-dependent food limitation, roe deer populations tend to irrupt due to the species ability to reproduce rapidly (Andersen and Linnell, 2000). However, this tendency is suppressed in less favourable environments with higher mortality due to predation and / or hunting. In our simulations, the hunting mortality shows a dramatic effect on roe deer populations. This is due to the combination of all the mortality factors together rather than only due to hunting mortality. Mortality from different sources greatly impacted the simulated roe deer populations. Our study indicates that even fawn predation by red fox alone can significantly stabilise or limit the population at a lower density. Indeed, roe deer fawns are an important part of the diet of red foxes during summer (Panzacchi et al., 2008), and fawn survival is strongly negatively dependent on red fox predation, which is directly linked to fluctuations in fox population density (Jarnemo and Liberg, 2005). The impact of lynx predation is even more important because lynx target different age classes (Andrén and Liberg, 2015; Heurich et al., 2016). By removing adults, and especially adult females with the greatest reproductive value, lynx predation exerts a greater pressure on roe deer mortality than fox predation on fawns (Gervasi et al., 2012). Some field studies have already shown the importance of lynx for roe deer mortality (Andrén and Liberg, 2015, 2023; Heurich et al., 2012). Moreover, lynx are stalking predators and are very efficient at preying on roe deer (Belotti et al., 2013). However, our model does not account for seasonal variations in predation rates, which may be linked to roe deer behaviour or habitat use, such as inter-seasonal migration (Cagnacci et al., 2011; Mejlgaard et al., 2013; Melis et al., 2009) which could induce more variability/stochasticity in the result but shouldn't change the overall population density.

It is well known that severe winters with heavy snowfall impact roe deer survival in different ways (Cederlund and Lindström, 1983; Grøtan et al., 2005). In our model, we implemented these extreme events to a medium/high additional mortality rate on the population to get an informative effect on the roe deer population. Although, it mainly added more variability in the output depending on when the snow event occurred during the 25-year runs. Most noticeably it drove a greater proportion of the simulated populations towards quasi-extinction under high mortality rates from multiple sources (i.e., predation and hunting) and particularly when populations were already at a low density. However, even in the absence of predation, mortality induced by extreme snow events seemed to greatly reduce the probability of irruptions. These results indicate that even limiting effects of mortality that are not necessarily density-dependent, such as the effect of red fox predation, can prevent populations from reaching levels where density dependence kicks in via food limitation.

Early research on roe deer commented that populations in many areas refrained from irruption, attributing this phenomenon to the regulatory effect of dispersal (Bobek, 1977; Strandgaard, 1972). Later research provided solid evidence for the role of density-dependence effects on demography in populations living under simple ecological conditions, where predation and hunting pressures are absent (Gaillard et al., 1998). However, in our simulations, population stability did not come from a density-dependent factor such as intraspecific competition, food availability, or disease. Instead, it arises from the combination of density-independent mortality sources including environmental and external factors such as weather events, hunting, and predation (set as a constant rate regardless of the density). The apparent stability was then reached via limitation, which imposes constraints on population growth, more than through regulation.

Population modelling and empirical data have shown the potential of spatial heterogeneity to stabilise predator-prey dynamics (Gorini et al., 2012). Although we included spatial heterogeneity in a straightforward manner in our model, adding net immigration and predation refuges stabilised the populations. The realism of these results will depend greatly on issues such as scale (i.e. dispersal distances) and how habitat heterogeneity differentially affects multiple mortality sources (Kauffman et al., 2007; Melis et al., 2013; Panzacchi et al., 2010). Resolving this requires new empirical analysis of how competing risks relate to fine-scaled habitat / landscape characteristics. Specifically, there is a need to explore how landscape characteristics and habitat contribute to lowering the risk of predation and how it affects other sources of mortality that depend on encounters (i.e. vehicle collisions, hunting, etc...). Additionally, the analysis of new empirical data could give a better understanding on how these factors influence the distance and direction of movements, such as natal dispersal.

The high frequency of simulations leading to decline and quasi-extinction raises the questions of whether the empirical vital rates (especially cause-specific mortality rates) are actually unbiased and realistic, and whether causes of mortality are purely additive. Our models considered the different sources of mortality as fully additive, which led to remarkably high mortality rates when all sources were included. We used data from field projects where the overall observed mortality rate was separated into different causes-specific mortality rates. The addition of all cause-specific mortality rates in our model was equal to the overall mortality rate from Melis et al., (2013), which is therefore not unrealistic as it reflects the same values as the original estimates. After careful consideration we have not been able to identify plausible mechanisms of compensation between hunting, lynx predation and fox predation. Even though there is little evidence for compensatory mortality in our study system (Andrén and Liberg, 2015), including a small amount of compensation would probably greatly impact the overall outcome of our model.

With respect to bias, the sample size was also high, >300 animals were studied, for periods from months up to multiple years. There is no evidence that collaring increases mortality (any cases where post-capture myopathy was suspected were excluded). Most animals older than 6 months were captured at winter feeding sites, which could have led to a bias. This potential bias could have positive implications, as feeding stations provide additional food and potentially increase survival rates (Milner et al., 2014). However, it could also potentially have negative implications, as prey are more easily predictable for lynx in these locations (Basille et al., 2009). However, in our study region almost all roe deer were associated with winter feeding sites. Therefore, while the risk of sampling bias cannot be excluded, we cannot see any obvious mechanism that would only lead to an upward bias in mortality. In our model, only the roe deer population was dynamic, and there was no feedback from roe deer population on the predator and human populations. This implies that predation (lynx and red fox) and hunting rates were considered constant throughout the 25-years simulation period disregarding any potential change in roe deer density. In other words, we assumed no functional responses and no lags in how predation or hunting rates respond to roe deer density.
Despite our model being a highly simplified representation of reality, some clear conclusions emerge. Firstly, the models suggest that severe population limitation, decline, or even quasi-extinction are highly probable outcomes when subject to multiple additive sources of mortality from predation and hunting based on field-based parameter estimates from our study ecosystem. This is supported by other empirical data (time series and cross population studies) on roe deer populations (e.g., Andrén and Liberg, 2015, 2023; Melis et al., 2010) exposed to lynx predation. The high impacts of lynx predation are mediated by the large proportion of adult females killed by lynx (Gervasi et al., 2012). Interestingly, the models show how density-independent effects can actually create stability (avoid irruption or quasi-extinction) via limitation. Combined, these insights show the contrast in how roe deer populations are likely to behave in real-world complex ecosystems as compared to the simple ecosystems where most real studies had previously been conducted. Secondly, the models show how episodic mortality events, like those that might be caused by extreme weather events or disease outbreaks, increase the risks of decline and quasi-extinction. Thirdly, the models show the potential of net immigration and spatial refuges from predation risk to buffer populations against decline and quasi-extinction. Although lynx have recolonised much of Scandinavia for decades, no roe deer population has been driven to local extinction although some have declined to lower densities (Andrén and Liberg, 2015, 2023; Crétois et al., 2021). This implies that ecological mechanisms like immigration or refuges, or others not included in our model such as source-sink dynamics, must allow populations to persist, even at low densities under real world conditions. Accordingly, there is a large scope for future studies (field-studies and simulations) to further explore these buffering persistence mechanisms. Among the likely mechanisms are the effects of roe deer dispersal and habitat refuges, as well as more complex functional and numerical responses from predators (also pointed out by Andrén and Liberg, 2023) and also factors inducing mortality in predators. Overall, we believe that there is considerable scope for exploring how spatio-temporal patterns in predator pressure and prey exposure as well as immigration / emigration of prey influence predator-prey relationships. While our models only hint at how these complex issues may influence prey dynamics, they do strongly indicate the need for their existence.

## Ethics statement

Not applicable: This manuscript does not include human or animal research.

## Declaration of Competing Interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

## Funding and acknowledgements

Funding for the first author's PhD scholarship was provided by the University of Inland Norway, and the analysis was based on data collected by the Scandlynx project using funding provided by the Research Council of Norway (grant numbers: 251112,281092 and 156810), the Norwegian Environment Agency, the Norwegian Institute for Nature Research, and multiple grants from the county wildlife management authorities in Innlandet and Viken counties over a 30-year period.

We would also like to thank Erlend B. Nilsen for his valuable feedback on the model and for the interesting discussion that followed.

## Appendix

<table><caption>**Table A.1**
List of parameter and value used in the matrix population model</caption><thead><tr><th>Parameter name</th><th>Description</th><th>Value Stochasticity</th><th>Unit</th><th>Study site (Reference)</th></tr></thead><tbody><tr><td colspan="5"><strong>Reproduction and recruitment</strong></td></tr><tr><td>rho.r</td><td>Proportion of female roe deer at birth</td><td>0.5</td><td>ratio (-)</td><td>(Assumed)</td></tr><tr><td>k1.r</td><td>Average litter size per primiparous female</td><td>2.14<br>N (μ = 2.14, σ = 0.38)</td><td>fawns (N)</td><td>Storfosna island, Norway (1993)<br>(Andersen and Linnell, 2000)</td></tr><tr><td>k2.r</td><td>Average litter size per adult female</td><td>2.18<br>N (μ = 2.18, σ = 0.73)</td><td>fawns (N)</td><td>Storfosna island, Norway (1993)<br>(Andersen and Linnell, 2000)</td></tr><tr><td>pr1.r</td><td>Breeding proportion of female roe deer 2-year-old</td><td>77.8 %<br>(64 % - 92 %)</td><td>proportion (-)</td><td>Storfosna island, Norway (1993)<br>(Andersen and Linnell, 2000)</td></tr><tr><td>pr2.r</td><td>Breeding proportion of female roe deer &gt;2-year-old</td><td>84.6 %<br>(77 % - 92 %)</td><td>proportion (-)</td><td>Storfosna island, Norway (1993)<br>(Andersen and Linnell, 2000)</td></tr><tr><td>h.r</td><td>Harem size, nb of female for one male</td><td>10</td><td>ind. (N)</td><td>(Assumed)</td></tr><tr><td colspan="5"><strong>Survival parameters</strong></td></tr><tr><td>m.jroef</td><td>Mortality rate for juvenile roe deer female (0–1)</td><td>0.175<br>N (μ = 0.175, σ = 0.056) *</td><td>proportion (-)</td><td>South-eastern Norway<br>(Melis et al., 2013) - Table 2</td></tr></tbody></table>

(continued on next page)
<table><caption>Table A.1 (continued)</caption><thead><tr><th>Parameter name</th><th>Description</th><th>Value Stochasticity</th><th>Unit</th><th>Study site (Reference)</th></tr></thead><tbody><tr><td>m.yroef</td><td>Mortality rate for yearling roe deer female (1–2)</td><td>0.241<br>N (μ = 0.241, σ = 0.093) *</td><td>proportion (-)</td><td>South-eastern Norway (Melis et al., 2013) - Table 2</td></tr><tr><td>m.proef</td><td>Mortality rate for primiparous roe deer female (2–3)</td><td>0.243<br>N (μ = 0.243, σ = 0.044) *</td><td>proportion (-)</td><td>South-eastern Norway (Melis et al., 2013) - Table 2</td></tr><tr><td>m.aroef</td><td>Mortality rate for adult roe deer female (&gt;3)</td><td>0.243<br>N (μ = 0.243, σ = 0.044) *</td><td>proportion (-)</td><td>South-eastern Norway (Melis et al., 2013) - Table 2</td></tr><tr><td>m.jroem</td><td>Mortality rate for juvenile roe deer male (0–1)</td><td>0.268<br>N (μ = 0.268, σ = 0.066) *</td><td>proportion (-)</td><td>South-eastern Norway (Melis et al., 2013) - Table 2</td></tr><tr><td>m.yroem</td><td>Mortality rate for yearling roe deer male (1–2)</td><td>0.116<br>N (μ = 0.116, σ = 0.089) *</td><td>proportion (-)</td><td>South-eastern Norway (Melis et al., 2013) - Table 2</td></tr><tr><td>m.aroem</td><td>Mortality rate for adult roe deer male (&gt;2)</td><td>0.162<br>N (μ = 0.162, σ = 0.049) *</td><td>proportion (-)</td><td>South-eastern Norway (Melis et al., 2013) - Table 2</td></tr><tr><td>redfox.r</td><td>Proportion of mortality due to red fox predation</td><td>0.263<br>N (μ = 0.263, σ = 0.047) *</td><td>proportion (-)</td><td>South-eastern Norway (Melis et al., 2013) - Table 2**</td></tr><tr><td>pred.jroef</td><td>Proportion of juvenile female predated by lynx</td><td>0.141<br>N (μ = 0.141, σ = 0.049) *</td><td>proportion (-)</td><td>South-eastern Norway (Melis et al., 2013) - Table 2</td></tr><tr><td>pred.yroef</td><td>Proportion of yearling female predated by lynx</td><td>0.131<br>N (μ = 0.131, σ = 0.083) *</td><td>proportion (-)</td><td>South-eastern Norway (Melis et al., 2013) - Table 2</td></tr><tr><td>pred.proef</td><td>Proportion of primiparous female predated by lynx</td><td>0.140<br>N (μ = 0.140, σ = 0.035) *</td><td>proportion (-)</td><td>South-eastern Norway (Melis et al., 2013) - Table 2</td></tr><tr><td>pred.aroef</td><td>Proportion of adult female predated by lynx</td><td>0.140<br>N (μ = 0.140, σ = 0.035) *</td><td>proportion (-)</td><td>South-eastern Norway (Melis et al., 2013) - Table 2</td></tr><tr><td>pred.jroem</td><td>Proportion of juvenile male predated by lynx</td><td>0.136<br>N (μ = 0.136, σ = 0.051) *</td><td>proportion (-)</td><td>South-eastern Norway (Melis et al., 2013) - Table 2</td></tr><tr><td>pred.yroem</td><td>Proportion of yearling male predated by lynx</td><td>0.213<br>N (μ = 0.213, σ = 0.113) *</td><td>proportion (-)</td><td>South-eastern Norway (Melis et al., 2013) - Table 2</td></tr><tr><td>pred.aroem</td><td>Proportion of adult male predated by lynx</td><td>0.152<br>N (μ = 0.152, σ = 0.048) *</td><td>proportion (-)</td><td>South-eastern Norway (Melis et al., 2013) - Table 2</td></tr><tr><td>hunt.jroef</td><td>Proportion of juvenile female harvested</td><td>0.111<br>N (μ = 0.111, σ = 0.051) *</td><td>proportion (-)</td><td>South-eastern Norway (Melis et al., 2013) - Table 2</td></tr><tr><td>hunt.yroef</td><td>Proportion of yearling female harvested</td><td>0.323<br>N (μ = 0.323, σ = 0.121) *</td><td>proportion (-)</td><td>South-eastern Norway (Melis et al., 2013) - Table 2</td></tr><tr><td>hunt.proef</td><td>Proportion of primiparous female harvested</td><td>0.076<br>N (μ = 0.076, σ = 0.030) *</td><td>proportion (-)</td><td>South-eastern Norway (Melis et al., 2013) - Table 2</td></tr><tr><td>hunt.aroef</td><td>Proportion of adult female harvested</td><td>0.076<br>N (μ = 0.076, σ = 0.030) *</td><td>proportion (-)</td><td>South-eastern Norway (Melis et al., 2013) - Table 2</td></tr><tr><td>hunt.jroem</td><td>Proportion of juvenile male harvested</td><td>0.115<br>N (μ = 0.115, σ = 0.047) *</td><td>proportion (-)</td><td>South-eastern Norway (Melis et al., 2013) - Table 2</td></tr><tr><td>hunt.yroem</td><td>Proportion of yearling male harvested</td><td>0.223<br>N (μ = 0.223, σ = 0.118) *</td><td>proportion (-)</td><td>South-eastern Norway (Melis et al., 2013) - Table 2</td></tr><tr><td>hunt.aroem</td><td>Proportion of adult male harvested</td><td>0.231<br>N (μ = 0.231, σ = 0.057) *</td><td>proportion (-)</td><td>South-eastern Norway (Melis et al., 2013) - Table 2</td></tr><tr><td>snow_frequency</td><td>Probability of this event to happen</td><td>0.08</td><td>frequency (-)</td><td>(Assumed)</td></tr><tr><td>snow_magnitude</td><td>How much mortality is induced from snow event</td><td>0.30</td><td>proportion (-)</td><td>(Assumed)</td></tr><tr><td>m.dens</td><td>Density dependence strength parameter</td><td>0.015</td><td>-</td><td>(Assumed)</td></tr></tbody></table>

<footnote>*Truncated Normal distribution between 0 and 1</footnote>

<footnote>**Pooled sexes estimates</footnote>
Data availability

All the data and R files are available on the GitHub repository: https://github.com/CecileAEC/MPM_roedeer

References

Andersen, R., Karlsen, J., Austmo, L.B., Odden, J., Linnell, J.D.C., Gaillard, J.-M., 2007. Selectivity of Eurasian Lynx lynx lynx and recreational hunters for age, sex and body condition in roe deer Capreolus capreolus. Article 4. Wildl. Biol. 13 (4). https://doi.org/10.2981/0909-6396(2007)13[467:SOELLL]2.0.CO;2

Andersen, R., Linnell, J.D.C., 1997. Variation in maternal investment in a small cervid; the effects of cohort, sex, litter size and time of birth in roe deer (Capreolus capreolus) fawns. Oecologia 109 (1), 74-79. https://doi.org/10.1007/s004420050060.

Andersen, R., Linnell, J.D.C., 1998. Ecological correlates of mortality of roe deer fawns in a predator-free environment. Can. J. Zool. 76, 9. https://doi.org/10.1139/z98-067.

Andersen, R., Linnell, J.D.C., 2000. Irruptive potential in roe deer: density-dependent effects on body mass and fertility. J. Wildl. Manag. 64 (3), 698. https://doi.org/10.2307/3802739.

Andersen, R., Linnell, J., Aanes, R., 1995. Rådyr. i Kult.: Slutrapport.

Andrén, H., Liberg, O., 2015. Large impact of Eurasian lynx predation on roe deer population dynamics. Article 3 PLOS ONE 10 (3). https://doi.org/10.1371/journal.pone.0120570.

Andrén, H., Liberg, O., 2023. Numerical response of predator to prey: dynamic interactions and population cycles in Eurasian lynx and roe deer. Ecol. Monogr. 94 (1), e1594. https://doi.org/10.1002/ecm.1594.

Apollonio, M., Belkin, V.V., Borkowski, J., Borodin, O.I., Borowik, T., Cagnacci, F., Danilkin, A.A., Danilov, P.I., Faybich, A., Ferretti, F., Gaillard, J.M., Hayward, M., Heshtaut, P., Heurich, M., Hurynovich, A., Kashtalyan, A., Kerley, G.I.H., Kjellander, P., Kowalczyk, R., Yanuta, G., 2017. Challenges and science-based implications for modern management and conservation of European ungulate populations. Mammal. Res. 62 (3), 209-217. https://doi.org/10.1007/s13364-017-0321-5.

Basille, M., Herfindal, I., Santin-Janin, H., Linnell, J.D.C., Odden, J., Andersen, R., Arild Høgda, K., Gaillard, J.-M., 2009. What shapes Eurasian lynx distribution in human dominated landscapes: selecting prey or avoiding people? Article 4. Ecography 32 (4). https://doi.org/10.1111/j.1600-0587.2009.05712.x.

Belotti, E., Cervený, J., Sustr, P., Kreisinger, J., Gaibani, G., Bufka, L., 2013. Foraging sites of Eurasian lynx (Lynx lynx): relative importance of microhabitat and prey occurrence. Article 2. Wildl. Biol. 19 (2). https://doi.org/10.2981/12-077.

Bobek, B., 1977. Summer food as the factor limiting roe deer population size. Nature 268 (5615), 47-49. https://doi.org/10.1038/268047a0.

Bolker B., Warnes G., Lumley T. (2022). *tools: Various R Programming Tools*. R package version 3.9.4, (https://CRAN.R-project.org/package=gtools).

Bouyer, Y., Gervasi, V., Poncin, P., Beudels-Jamar, R.C., Odden, J., Linnell, J.D.C., 2015. Tolerance to anthropogenic disturbance by a large carnivore: the case of Eurasian lynx in south-eastern eastern Europe. Article 3. Anim. Conserv. 18 (3). https://doi.org/10.1111/acv.12168.

Bouyer, Y., Rigot, T., Panzacchi, M., Moorter, B.V., Poncin, P., Beudels-Jamar, R., Odden, J., Linnell, J.D.C., 2015. Using zero-inflated models to predict the relative distribution and abundance of roe deer over very large spatial scales. Article 1-2 Ann. Zool. Fenn. 52 (1-2). https://doi.org/10.5735/086.052.0206.

Bouyer, Y., San Martin, G., Poncin, P., Beudels-Jamar, R.C., Odden, J., Linnell, J.D.C., 2015. Eurasian lynx habitat selection in human-modified landscape in Norway: effects of different human habitat modifications and behavioral states. Biol. Conserv. 191, 291-299. https://doi.org/10.1016/j.biocon.2015.07.007.

Cagnacci, F., Focardi, S., Heurich, M., Stache, A., Hewison, A.J.M., Morellet, N., Kjellander, P., Linnell, J.D.C., Mysterud, A., Neteler, M., Delucchi, L., Ossi, F., Urbano, F., 2011. Partial migration in roe deer: migratory and resident tactics are end points of a behavioural gradient determined by ecological factors. Article 12. Oikos 120 (12). https://doi.org/10.1111/j.1600-0706.2011.19441.x.

Caswell, H. (2000). *Matrix population models: Construction, analysis, and interpretation* (2. ed., 2. printing). Sinauer.

Cederlund, H., Lindström, E., 1983. Effects of severe winters and fox predation on roe deer mortality. Acta Theriol. 28, 129-145. https://doi.org/10.4098/AT.arch.83-9.

Chapron, G., Kaczensky, P., Linnell, J.D.C., von Arx, M., Huber, D., Andrén, H., López-Bao, J.V., Adamec, M., Álvares, F., Anders, O., Balčiauskas, L., Balys, V., Bedò, P., Bego, F., Blanco, J.C., Breitenmoser, U., Brøseth, H., Bufka, L., Bunikyte, R., Boitani, L., 2014. Recovery of large carnivores in Europe's modern human-dominated landscapes. *Sci.* **346**(6216), Artic. 6216. https://doi.org/10.1126/science.1257553.

Cobben, M.M.P., Linnell, J.D.C., Solberg, E.J., Andersen, R., 2009. Who wants to live forever? Roe deer survival in a favourable environment. Article 6. Ecol. Res. 24 (6). https://doi.org/10.1007/s11284-009-0600-y.

Crétois, B., Linnell, J.D.C., Van Moorter, B., Kaczensky, P., Nilsen, E.B., Parada, J., Rød, J.K., 2021. Coexistence of large mammals and humans is possible in Europe's anthropogenic landscapes. Article 9 *iScience* 24 (9). https://doi.org/10.1016/j.isci.2021.103083.

Focardi, S., Pelliccioni, E., Petrucco, R., Toso, S., 2002. Spatial patterns and density dependence in the dynamics of a roe deer (Capreolus capreolus) population in central Italy. Oecologia 130 (3), 411-419. https://doi.org/10.1007/s00442-001-0825-0.

Gaillard, J.-M., Andersen, R., Delorme, D., Linnell, J.D.C., 1998. Family effects on growth and survival of juvenile roe deer. Ecology 79 (8), 2878-2889. https://doi.org/10.1890/0012-9658(1998)079[2878:FEOGAS]2.0.CO;2.

Gaillard, J.-M., Delorme, D., Jean-Marie, B., Laere, G.V., Boisaubert, B., Pradel, R., 1993. Roe deer survival patterns: a comparative analysis of contrasting populations. J. Anim. Ecol. 62 (4), 778. https://doi.org/10.2307/5396.

Gervasi, V., Nilsen, E.B., Odden, J., Bouyer, Y., Linnell, J.D.C., 2014. The spatio-temporal distribution of wild and domestic ungulates modulates lynx kill rates in a multi-use landscape: Lynx prey selection in a multi-use landscape. Article 3 J. Zool. 292 (3). https://doi.org/10.1111/jzo.12088.

Gervasi, V., Nilsen, E.B., Sand, H., Panzacchi, M., Rauset, G.R., Pedersen, H.C., Kindberg, J., Wabakken, P., Zimmermann, B., Odden, J., Liberg, O., Swenson, J.E., Linnell, J.D.C., 2012. Predicting the potential demographic impact of predators on their prey: A comparative analysis of two carnivore-ungulate systems in Scandinavia. Article 2 J. Anim. Ecol. 81 (2). https://doi.org/10.1111/j.1365-2656.2011.01928.x.

Gill, R.M.A., Johnson, A.L., Francis, A., Hiscocks, K., Peace, A.J., 1996. Changes in roe deer (Capreolus capreolus L.) population density in response to forest habitat succession. For. Ecol. Manag. 88 (1-2), 31-41. https://doi.org/10.1016/S0378-1127(96)03807-8.

Gorini, L., Linnell, J.D.C., May, R., Panzacchi, M., Boitani, L., Odden, M., Nilsen, Erlend B., 2012. Habitat heterogeneity and mammalian predator-prey interactions: predator-prey interactions in a spatial world. Article 1. Mammal. Rev. 42 (1). https://doi.org/10.1111/j.1365-2907.2011.00189.x.

Grøtan, V., Sæther, B.-E., Engen, S., Solberg, E.J., Linnell, J.D.C., Andersen, R., Brøseth, H., Lund, E., 2005. Climate causes large-scale spatial synchrony in population fluctuations of a temperate herbivore. Article 6. Ecology 86 (6). https://doi.org/10.1890/04-1502.

Heurich, M., Möst, L., Schauburger, G., Reulen, H., Sustr, P., Hothorn, T., 2012. Survival and causes of death of European Roe Deer before and after Eurasian Lynx reintroduction in the Bavarian Forest National Park. Eur. J. Wildl. Res. 58 (3), 567-578. https://doi.org/10.1007/s10344-011-0606-y.

Heurich, M., Zeis, K., Küchenhoff, H., Müller, J., Belotti, E., Bufka, L., Woelfing, B., 2016. Selective predation of a stalking predator on ungulate prey. PLOS ONE 11 (8), e0158449. https://doi.org/10.1371/journal.pone.0158449.

Hewison, A.J.M., Angibault, J.M., Bideau, E., Vincent, J.P., Boutin, J., Sempéré, A., 1996. Annual variation in body composition of roe deer (Capreolus capreolus) in moderate environmental conditions. Can. J. Zool. 74 (2), 245-253. https://doi.org/10.1139/z96-031.

Hewison, A.J.M., Gaillard, J.-M., Morellet, N., Cagnacci, F., Debeffe, L., Cargnelutti, B., Gehr, B., Kröschel, M., Heurich, M., Coulon, A., Kjellander, P., Börger, L., Focardi, S., 2021. Sex differences in condition dependence of natal dispersal in a large herbivore: dispersal propensity and distance are decoupled. Proc. R. Soc. B: Biol. Sci. 288 (1946), 20202947. https://doi.org/10.1098/rspb.2020.2947.

Hewison, A.J.M., Gaillard, J.M., 2001. Phenotypic quality and senescence affect different components of reproductive output in roe deer. J. Anim. Ecol. 70 (4), 600-608. https://doi.org/10.1046/j.1365-2656.2001.00528.x.

Jarnemo, A., Liberg, O., 2005. Red fox removal and roe deer fawn survival-A 14 years study. J. Wildl. Manag. 69 (3), 1090-1098. https://doi.org/10.2193/0022-541X (2005)069[1090:RFRARD]2.0.CO;2.
Kauffman, M.J., Varley, N., Smith, D.W., Stahler, D.R., MacNulty, D.R., Boyce, M.S., 2007. Landscape heterogeneity shapes predation in a newly restored predator–prey system. Ecol. Lett. 10 (8), 690–700. https://doi.org/10.1111/j.1461-0248.2007.01059.x

Kjellander, P., Hewison, A.J.M., Liberg, O., Angibault, J.-M., Bideau, E., Cargnelutti, B., 2004. Experimental evidence for density-dependence of home-range size in roe deer (Capreolus capreolus L.): a comparison of two long-term studies. Oecologia 139 (3), 478–485. https://doi.org/10.1007/s00442-004-1529-z.

Linnell, J.D.C., Cretois, B., Nilsen, E.B., Rolandsen, C.M., Solberg, E.J., Veiberg, V., Kaczensky, P., Van Moorter, B., Panzacchi, M., Rauset, G.R., Kaltenborn, B., 2020. The challenges and opportunities of coexisting with wild ungulates in the human-dominated landscapes of Europe’s Anthropocene. Biol. Conserv. 244, 108500. https://doi.org/10.1016/j.biocon.2020.108500.

Linnell, J.D.C., Wahlström, L.K., Gaillard, J.M., 1998. From birth to independence: birth, growth, neonatal mortality, hiding behaviour and dispersal. In: Andersen, R., Duncan, P., Linnell, J.D.C. (Eds.), In The European Roe Deer: the Biology of Success. Scandinavian University Press, pp. 257–284.

Lone, K., Loe, L.E., Gobakken, T., Linnell, J.D.C., Odden, J., Remmen, J., Mysterud, A., 2014. Living and dying in a multi-predator landscape of fear: roe deer are squeezed by contrasting pattern of predation risk imposed by lynx and humans. Oikos 123 (6), 641–651. https://doi.org/10.1111/j.1600-0706.2013.00938.x.

Mejlgaard, T., Loe, L.E., Odden, J., Linnell, J.D.C., Nilsen, E.B., 2013. Lynx prey selection for age and sex classes of roe deer varies with season: age- and sex-dependent predation. Article 3 J. Zool. 289 (3). https://doi.org/10.1111/j.1469-7998.2012.00980.x.

Melis, C., Basille, M., Herfindal, I., Linnell, J.D.C., Odden, J., Gaillard, J.-M., Høgda, K.-A., Andersen, R., 2010. Roe deer population growth and lynx predation along a gradient of environmental productivity and climate in Norway. Article 2 Écoscience 17 (2). https://doi.org/10.2980/17-2-3314.

Melis, C., Jędrzejewska, B., Apollonio, M., Barton, K.A., Jędrzejewski, W., Linnell, J.D.C., Kojola, I., Kusak, J., Adamic, M., Ciuti, S., Delehan, I., Dykyy, I., Krapinec, K., Mattioli, L., Sågaydak, A., Samchuk, N., Schmidt, K., Shkvyrya, M., Sidorovich, V.E., Zhyla, S., 2009. Predation has a greater impact in less productive environments: variation in roe deer, Capreolus capreolus, population density across Europe. Article 6 Glob. Ecol. Biogeogr. 18 (6). https://doi.org/10.1111/j.1466-8238.2009.00480.x.

Melis, C., Nilsen, E.B., Panzacchi, M., Linnell, J.D.C., Odden, J., 2013. Roe deer face competing risks between predators along a gradient in abundance. Article 9 Ecosphere 4 (9). https://doi.org/10.1890/ES13-00099.1.

Mersmann O., Trautmann H., Steuer D., Bornkamp B. (2023). truncnorm: Truncated Normal Distribution. R package version 1.0-9. (https://CRAN.R-project.org/package=truncnorm).

Milner, J.M., Van Beest, F.M., Schmidt, K.T., Brook, R.K., Storaas, T., 2014. To feed or not to feed? Evidence of the intended and unintended effects of feeding wild ungulates. J. Wildl. Manag. 78 (8), 1322–1334. https://doi.org/10.1002/jwmg.798.

Nilsen, E.B., Gaillard, J.-M., Andersen, R., Odden, J., Delorme, D., van Laere, G., Linnell, J.D.C., 2009. A slow life in hell or a fast life in heaven: Demographic analyses of contrasting roe deer populations. Article 3 J. Anim. Ecol. 78 (3). https://doi.org/10.1111/j.1365-2656.2009.01523.x.

Panzacchi, M., Herfindal, I., Linnell, J.D.C., Odden, M., Odden, J., Andersen, R., 2010. Trade-offs between maternal foraging and fawn predation risk in an income breeder. Article 8 Behav. Ecol. Sociobiol. 64 (8). https://doi.org/10.1007/s00265-010-0941-0.

Panzacchi, M., Linnell, J.D.C., Odden, M., Odden, J., Andersen, R., 2009. Habitat and roe deer fawn vulnerability to red fox predation. Article 6 J. Anim. Ecol. 78 (6). https://doi.org/10.1111/j.1365-2656.2009.01584.x.

Panzacchi, M., Linnell, J.D.C., Serrao, G., Eie, S., Odden, M., Odden, J., Andersen, R., 2008. Evaluation of the importance of roe deer fawns in the spring–summer diet of red foxes in southeastern Norway. Article 5 Ecol. Res. 23 (5). https://doi.org/10.1007/s11284-007-0452-2.

Pettorelli, N., Gaillard, J.-M., Duncan, P., Maillard, D., Van Laere, G., Delorme, D., 2003. Age and Density Modify the Effects of Habitat Quality on Survival and Movements of Roe deer. Ecology 84 (12), 3307–3316. https://doi.org/10.1890/02-0602.

Strandgaard, H., 1972. The roe deer (Capreolus capreolus) population at Kalø and the factors regulating its size. In: Joensen, A.H. (Ed.), Danish review of game biology, 7.

Torres, R.T., Virgós, E., Panzacchi, M., Linnell, J.D.C., Fonseca, C., 2012. Life at the edge: Roe deer occurrence at the opposite ends of their geographical distribution, Norway and Portugal. Article 2 Mamm. Biol. 77 (2). https://doi.org/10.1016/j.mambio.2011.11.002.