Demographic projections of the Apennine brown bear population Ursus arctos marsicanus (Mammalia: Ursidae) under alternative management scenarios

Abstract

Apennine brown bears are a very small, isolated population of central Italy, consisting of about 50 individuals and under a severe risk of extinction. We performed a population viability analysis (PVA) for this population, contrasting a deterministic model and an individual-based stochastic model, using a set of demographic parameters estimated for the same population during the last decade. We also built a set of simulated management scenarios, in which we compared the effectiveness of alternative conservation measures and assessed the susceptibility of the population to catastrophic mortality events. The deterministic model produced an estimate of the asymptotic population growth rate r = 0.001, corresponding to an asymptotically stable population. The stochastic model produced an estimate of r = − 0.013 (standard deviation = 0.103), corresponding to an annual population decrease of 1.3%, a 17% extinction risk in 100 years, an average population of 27 bears for non-extinct populations, and an average time to extinction of 81 years for those gone extinct. Extinction probability increased to more alarming levels (> 0.4) when at least one catastrophic event occurred during a 100-year period. Current vital rates of the population are not compatible with a more than negligible numerical increase, and this bear population is likely to remain small and exposed to a relatively high risk of extinction, if the average survival or reproductive rates do not increase. Management efforts aimed to increase food availability generated minimal to moderate variations in population growth rate and in the associated risk of extinction, whereas interventions meant to reduce adult female mortality were highly effective in increasing persistence probability. We propose that the general objectives of the action plan for the conservation of the Apennine brown bear for the incoming decade should explicitly contemplate quantitative demographic goals, focusing in particular on adult female and cub mortality.

Keywords:

• Extinction risk
• small populations
• individual-based stochastic simulation models
• viability analysis
• VORTEX

Introduction

Apennine brown bears (Ursus arctos marsicanus Altobello, 1921) are a very small, isolated and endemic population (Loy et al. 2008; Colangelo et al. 2012), consisting of about 50 individuals (Ciucci et al. 2015a) facing a severe risk of extinction (Ciucci & Boitani 2008). This population is the last remnant, genetically isolated portion of a formerly larger one, with an historical distribution along a wider range in the central Apennines of Italy (Ciucci & Boitani 2008; Benazzo et al. 2017). During the last 15 years, considerable resources have been invested in ecological research, monitoring and conservation actions on this population. These interventions have focused both at the administrative level, by fostering an interregional platform through the adoption of a renewed Action Plan (Anonymous 2011), and at the bear–human interface, by facilitating conflict resolution and mitigation of human-caused mortality (e.g. PNALM 2015; SLO 2017). Although the current Action Plan provides for general, quantitative objectives (i.e. a 25% increase in the overall Apennine brown bear population within 2020 and a 50% reduction in known bear mortality due to illegal causes, compared to the decade 2000–2010; Anonymous 2011), such goals were established without any detailed knowledge of the underlying demographic processes, which had been only minimally explored at the time.

Despite these renewed efforts, brown bears in central Italy still appear substantially concentrated in a core population within the Abruzzo, Lazio and Molise National Park (PNALM) and adjacent areas. Although a few bears are increasingly detected in the peripheral portions of the range, their demographic relevance has been negligible until recent years, as the core population is the only one stably hosting reproductive females (Ciucci et al. 2017). Accordingly, no evidence of population growth has been detected in the core distribution during the past 13 years (Gervasi et al. 2017). Whereas this could partly be explained by a relatively high bear density within the core population (i.e. 39.7 bears/1000 km²; Ciucci et al. 2015a), suggesting it is approximating carrying capacity, no other stable reproductive nuclei have been established in this time in the peripheral portions of the range, nor has the population’s distribution noticeably expanded beyond the historical occurrence of erratic bears (Ciucci et al. 2017). However, habitat suitability and connectivity at the landscape scale do not seem to be limiting (Falcucci et al. 2009; Maiorano et al. 2017). Therefore, the apparent failure of the population to expand could be due to intrinsic factors (i.e. a small number of reproducing females and a relatively low reproductive rate; Gervasi et al. 2017; Tosoni et al. 2017a,b) coupled with high levels of human-caused mortality (Falcucci et al. 2009). A substantial fraction of the adult mortality in Apennine brown bears occurs from causes related to humans. Over 50% (n = 110) of the adult bears reported as dead between 1970 and 2010 (i.e. the minimum known mortality) died because of poaching or poisoning, or in car/train accidents (Leonardo Gentile, National Park of Abruzzo Lazio and Molise, pers. comm.). In addition, a decreased demographic vigour of this isolated bear population cannot be discounted, as the long isolation time, the highly reduced genetic variability, an extremely limited effective population size and a high level of inbreeding all suggest that inbreeding depression could likely be in place (Benazzo et al. 2017).

Population viability analysis (PVA; Beissinger & McCullough 2002) is a useful but much debated tool to support managers and policymakers in the decision process for the conservation of endangered populations. By projecting the populations in the future, PVAs allow the estimation of population size, trends and probability of persistence based on the current state. Rather than being useful for providing estimates of these projections per se, PVAs are of practical value by illustrating the expected, relative outcome of alternative management scenarios (Beissinger & Westphal 1998). In these terms, PVAs allow management and conservation to be placed within their appropriate biological and ecological context (Reed et al. 2002). PVAs have been criticised due to the high number of parameters required to build a model with an acceptable degree of realism, and hence the high degree of residual uncertainty that is usually associated with their results (Coulson et al. 2001). However, PVAs rarely have the goal of providing accurate estimates of population status and trend, and more often they are instrumental to inform decision-making in a context of uncertainty. As a matter of fact, PVAs have been performed for both of the other two small and highly imperilled bear populations in western Europe (i.e. brown bears in the Cantabrian mountains and in the Pyrenees), using mainly simulated or literature-based demographic parameters (Wiegand et al. 1998; Chapron et al. 2003; Martinez-Cano et al. 2016). In all these cases, however, viability analyses emphasised the practical importance of maintaining high survival levels for the breeding segment of the population, while rejecting the hypothesis that nutritional stress could be the main cause of demographic stagnation (Wiegand et al. 1998; Chapron et al. 2003). Also, the PVA performed on the Cantabrian bear population confirmed that the use of annual counts of family units contained valuable information on the state of the whole bear population, thus supporting the continued use of such a monitoring tool in future years (Wiegand et al. 1998).

No PVAs have been previously conducted on the relict and isolated population of Apennine bears, due to the paucity of data on the most relevant demographic parameters (Ciucci & Boitani 2008). However, the recent availability of demographic information (Ciucci et al. 2015a; Gervasi et al. 2017; Tosoni et al. 2017a) provides for the first time the opportunity to develop a demographic projection model of the dynamics of this bear population. Following Beissinger and Westphal (1998), who cautioned against using too-complex models when performing PVAs on endangered species, we aimed to assess relative rather than absolute extinction risk, and extended projections over short- rather than long-term periods. To this aim, we projected the demographic dynamics of Apennine bears contrasting a deterministic, population-based approach, and a stochastic, individual-based one. Specifically, we identified the following objectives: (i) to estimate the expected growth rate, and the associated extinction risk, under current demographic parameters and conservation measures; (ii) to assess how alternative management scenarios, ultimately affecting either reproductive performances or human-caused mortality, could enhance the demography of the population and its probability of persistence; and (iii) to assess the extent to which the occurrence of catastrophic events should be expected to decrease population persistence. We focussed our analysis on the core population because this is currently the sole demographic source (Ciucci et al. 2017). Accordingly, we were particularly interested in evaluating future conservation investments that would most effectively facilitate demographic and range expansion of the core population beyond its current limits. Because effective monitoring of this population is crucial for a timely assessessment of the probability of population persistence in response to management interventions, we also interpret our PVA results to define efficient and realistic monitoring objectives for the future of this bear population.

Materials and methods

Study area

The core of the Apennine brown bear distribution roughly corresponds to the PNALM and its outer buffer zone (Ciucci et al. 2017), which cover about 1200 km². The area is located in the central Apennines (Italy) along a NW–SE direction, with elevation ranging from 400 to 2285 m. Average temperatures range from 2°C in winter to 20°C in summer, while snow cover generally extends from mid-December to March. Beech (Fagus sylvatica) and oak (Quercus cerris and Q. pubescens) cover about 60% of the area. Average human density is 14.6 inhabitants/km². Forest cutting in the PNALM is strictly regulated by the park authority and hunting of any kind is prohibited within the PNALM, although hunting with dogs is allowed in the external buffer zone (Maiorano et al. 2015). Additional details of the study area can be found elsewhere (Gervasi et al. 2008; Ciucci et al. 2015a).

Deterministic, population-based demographic model

We first explored the dynamics of the Apennine brown bear population with a deterministic model of its asymptotic behaviour under constant vital rates (Caswell 2001). Using the software ULM (Legendre & Clobert 1995), we implemented a Lefkovitch stage-structured matrix (Lefkovitch 1965), including five age classes for females and two for males, configured as follows:

$\begin{array}{ccccccc}0& 0& 0& 0& L\ast \rho \ast R& 0& 0\\ 1-M_c& 0& 0& 0& 0& 0& 0\\ 0& 1-M_f& 0& 0& 0& 0& 0\\ 0& 0& 1-M_f& 0& 0& 0& 0\\ 0& 0& 0& 1-M_f& 1-M_f& 0& 0\\ 0& 0& 0& 0& L\ast \rho \ast \left(1-R\right)& 0& 0\\ 0& 0& 0& 0& 0& 1-M_c& 1-M_m\end{array}$

where M_c is the annual cub survival of both sexes, M_f is the annual survival for females older than 1 year, M_m is the annual survival for males older than 1 year, L is litter size, ρ is the proportion of females reproducing each year and R is the sex ratio at birth.

Most of the parameters had been estimated for the Apennine brown bear population over the preceding 10 years, using data from non-invasive genetic sampling, live-tapping and radio-telemetry, and from direct observations of family units during summer (for methodological details see Ciucci et al. 2015a; Gervasi et al. 2017; Tosoni et al. 2017a). The only exception was the age at first reproduction, for which we lacked reliable information and had to resort to a literature-based parameter estimate. In particular, we assumed the age of first reproduction of female bears was 4 years, corresponding to the average value estimated in other brown bear populations in western Europe (Cantabrian Mountains in Spain: Wiegand et al. 1998; Eastern Alps: Groff et al. 2015; Scandinavia: Zedrosser et al. 2013). After estimating the asymptotic growth rate, we performed a sensitivity analysis to rank vital rates in order of their relative influence on the long-term population performance, based on which we successively designed management scenarios (see below). All parameter values are summarised in Table I.

Table I. Input parameters for the population projection model of Apennine brown bears (Ursus arctos marsicanus) in central Italy. Initial values were used to parameterise the reference scenario (demographic parameters estimated on the population between 2003 and 2014), whereas simulation intervals refer the range used for each parameter during the sensitivity analysis.

Parameter	Initial value	Simulation interval	Source
No. iterations	10,000	-	-
No. simulated years	100	-	-
Mortality cubs of the year (Mc, both sexes)	0.49	0.3–0.6	Gervasi et al. 2017
Female mortality (Mf, age > 1 year old)	0.08	0.02–0.15	Gervasi et al. 2017
Male mortality (Mm, age > 1 year old)	0.15	0.05–0.2	Gervasi et al. 2017
Type of reproductive system	Polygynous	-
Age of first reproduction	4	-	From literature
Reproductive senescence	25	-	Schwartz et al. 2003
Litter size (L)	1.9	-	Tosoni et al. 2017b
Sex ratio at birth (R)	50:50	-	From literature
Proportion of reproducing females ($\rho )$	0.24	0.15–0.3	Tosoni et al. 2017b
Initial population size	51	-	Ciucci et al. 2015a
Carrying capacity (K)	64	51–102	-

Stochastic, individual-based demographic model

We also ran an individual-based, stochastic population projection model using VORTEX (v. 10, Lacy & Pollock 2015) and the same demographic parameters illustrated above (Table I). Differently from the deterministic model, however, VORTEX used each parameter value to define a normal probability distribution from which individual survival and reproduction events were extracted. For reproductive senescence we referred to Schwartz et al. (2003), fixing it at 25 years, but we expected that this parameter had a minimal effect on the model’s performance since the estimated proportion of individuals > 25 years old did not exceed 5%, based on the stable age distribution of the population (see Results). We ran the model over a time interval of 100 years, allowing 10,000 iterations for each scenario, and setting an initial population size of 51 bears, as estimated for the Apennine brown bear population in 2014 (Ciucci et al. 2015b). In VORTEX, the carrying capacity function acts on population dynamics through a k-truncation, so that populations exceeding K are pushed back by a proportional increase in annual mortality (Lacy & Pollock 2015). Although carrying capacity in the bear core distribution has not been formally estimated, there are indications that the population is likely approaching it (Ciucci et al. 2015a; Gervasi et al. 2017). We conservatively assumed a carrying capacity of 125% with respect to the bear density estimated in 2014. We initially simulated environmental stochasticity according to three levels of increasing variability (coefficient of variation, CV = 0.1, 0.2, 0.3) of survival and reproductive parameters. Given the similarity of population trajectories under different levels of environmental variance, we used a fixed value of CV = 0.2 in all subsequent analyses, corresponding to an intermediate level of environmental stochasticity. We also explored the effects of catastrophic events on the probability of persistence of Apennine bears. Based on an estimated frequency of catastrophic events in vertebrate populations, defined as a 50% or higher reduction in survival rate, of about 0.14 per generation (Reed et al. 2003), we fixed the expected frequency of extreme mortality events to 1.4 every 100 years, according to a generation time of 10 years (Harris & Allendorf 1989; Gaillard et al. 2005). We simulated low, moderate and high mortality effects of catastrophic events by associating increasing levels of mortality (10–30% reduction of survival) to increasing frequencies of catastrophic events (0–10 events every 100 years).

To conduct a sensitivity analysis of the stochastic model, we first used Monte Carlo simulations to obtain 1000 datasets corresponding to an equal number of combinations among all parameters, each varying randomly within a pre-defined range (Table I). Then, we ran VORTEX for each dataset and regressed (multiple linear regression) all the estimated growth rates of the population (response variable) against the standardised input parameter values. By doing so, model coefficients represented the elasticity of the stochastic population growth rate to changes in the corresponding demographic parameter (Cross & Beissinger 2001).

We ran simulations over a period of 100 years, and evaluated models’ output in terms of: (i) average population growth rate (r); (ii) probability of population persistence (i.e. the proportion of simulated populations that survived); (iii) average time of extinction; and (iv) average size of the population at the end of the simulated time interval. While a period of a century for the analysis of viability is compatible with the brown bear life cycle, such a time scale may be of limited practical value for management and monitoring purposes. Thus, we extracted the projection model results at shorter time intervals (i.e. 10, 20, 30 years) to define practical demographic goals for monitoring the efficacy of future management efforts. Additionally, by estimating the population trajectories separately for the populations gone extinct and those which persisted, we estimated a probability distribution of persistence probability as a function of projected population size.

To account for alternative management scenarios in VORTEX, we first defined a reference scenario (scenario 1) corresponding to stable demographic parameters over the simulated time interval (Table I). In a second group of simulations (scenario 2), we mimicked the effects of conservation efforts meant to increase the availability of food resources, as these are expected to positively affect reproductive performances of the population, cub survival and carrying capacity (Newton 1998; Sibly & Hone 2002). Accordingly, we simulated a 5, 10 and 15% increase in the proportion of females reproducing each year and in K, and a 10% increase in cub survival. In a third group of simulations (scenario 3) we contemplated a more effective prevention of human-caused mortality, with the immediate effect of enhancing the survival of females of reproducing age: (i) a 25% reduction in total mortality of adult females, that therefore decreases from 0.08 (Gervasi et al. 2017) to 0.06; and (ii) a 50% reduction in mortality of adult females (i.e. from 0.08 to 0.04). Finally, we ran a mixed scenario (scenario 4), in which we contemplated an intermediate achievement both in limiting human-caused mortality (i.e. 10% reduction in adult female mortality) and in enhancing habitat productivity (10% increase in fecundity and carrying capacity).

Results

The deterministic model produced an estimate of the asymptotic population growth rate r = 0.001, corresponding to an asymptotically stable population of 59 bears in 100 years. At the equilibrium, cubs represented 17% of the population, subadults 23%, adult females 37%, adult males 18%, and senescent individuals 5%. Female generation time was 11.4 years. Population growth rate was most influenced by mortality of female bears of reproductive age, whereas both cub mortality and the proportion of females reproducing each year corresponded to much lower elasticity values (Table II).

Table II. Sensitivity analysis under the deterministic and stochastic projection models for the Apennine brown bear (Ursus arctos marsicanus) population in central Italy. Model parameters used for the analysis are cub survival (M_c), female survival (M_f), male survival (M_m), the proportion of females reproducing each year ($\rho$) and carrying capacity (K).

	Deterministic model		Stochastic model
Parameter	Sensitivity	Elasticity	Regression coefficient
Mc	−0.128	0.065	−0.014
Mf (1–3 years old)	−0.212	0.193	−0.019
Mf (> 3 years old)	−0.741	0.677	−0.037
Mm	−0.000	0.000	−0.001
$\rho$	0.271	0.065	0.016
K	-	-	0.001

The stochastic model produced an estimate of population growth rate r = − 0.013 (standard deviation, SD = 0.103), corresponding to an annual population decrease of 1.3%, a 17% extinction risk in 100 years, an average population of 27 bears for non-extinct populations, and an average time to extinction of 81 years for those gone extinct (Table III). The introduction of progressively higher levels of environmental stochasticity did not drastically alter the average growth rate of the population, even though its variance increased up to 30% of the reference value (Table III). The sensitivity of stochastic population growth rate to changes in adult female survival was about 2.5 times higher than that referring to cub survival (Table II), compared to a 10-fold difference in the deterministic model. Both adult male mortality and carrying capacity exhibited minimum values of elasticity and a very weak link with population growth rate.

Table III. Results of the stochastic population projection model of Apennine brown bears (Ursus arctos marsicanus) in central Italy, using VORTEX and under the reference scenario (i.e. constant demographic parameters estimated on the population between 2003 and 2014). Increasing levels of environmental stochasticity were defined by the coefficient of variation (CV) of both survival and reproductive parameters in the model. The average population growth rate (r) and its standard deviation (SD) are provided, along with measures of population persistence in 100 years.

	Population growth
Environmental stochasticity (CV)	r	SD	Persistence probability in 100 years	Mean time to extinction (years)	Mean population size in 100 years
0	−0.013	0.103	0.887	81	27
0.1	−0.013	0.107	0.854	81	27
0.2	−0.014	0.116	0.831	80	26
0.3	−0.016	0.130	0.787	78	25

The trajectories of extinct and persistent populations diverged markedly through time (Figure 1(a)). Only the probability distributions of bear population size within the first 10 years of simulation largely overlapped (Figure 1(a)). Also, the probability of extinction was expected to remain at low levels for the next 50 years, while it cumulated up to 17% during the second 50 years of the simulations (Figure 1(b)). Using the current extinction risk as a baseline, if population size increased to 60 bears in the next 10 years, extinction risk in 100 years would drop to 8%, whereas it would rise to 23% if the population should decline to 40 individuals (Figure 1(c)).

Figure 1. Extinction probability for the Apennine brown bear (Ursus arctos marsicanus) population in central Italy, under the reference scenario of the stochastic population projection model: (a) the average trajectories of persistent and extinct populations over a period of 100 years; (b) the probability of extinction as a function of the number of years; (c) the probability of extinction as a function of population size; vertical dashed line represents population estimated in 2014.

Under moderate but relatively frequent catastrophic events, extinction probability steadily increased to about 0.4 (Figure 2). Less frequent but more intense catastrophic events produced a steep increase in the extinction probability, so that a single mass mortality event (−50% in bear survival) every 50 years would be sufficient to raise extinction probability to 0.6 (Figure 2).

Figure 2. Probability of extinction for the Apennine brown bear (Ursus arctos marsicanus) population, as a function of the frequency and intensity of simulated catastrophic events. Results refer to the stochastic version of the population projection model, under moderate levels of environmental stochasticity.

Compared to the reference scenario, all simulated management scenarios benefitted the population to some extent (Table IV). However, management efforts aimed to increase food availability (scenario 2) generated minimal to moderate variations in population growth rate, and in the associated risk of extinction, compared to management interventions meant to reduce mortality (scenario 3; Table IV, Figure 3(a and b)). A 25% reduction in adult females’ mortality (i.e. a 2% increase in survival) produced an estimated positive growth rate of r = 0.018, which corresponded to a 99.8% probability of persistence and to an average population size of 57 bears in 100 years (Table IV). Further reductions in adult female mortality (i.e. an additional increase of 2% in annual survival) corresponded to a higher population growth rate (r = 0.040) and to no population gone extinct over a period of 100 years (Table IV).

Table IV. Results of the stochastic population projection model of Apennine brown bears (Ursus arctos marsicanus) in central Italy, under different future management scenarios. Parameters used to construct the scenarios included cub survival (M_c), carrying capacity (K), the proportion of females reproducing each year ($\rho$) and adult female survival (M_f). For each scenario the average population growth rate (r) and its standard deviation (SD) are provided.

			Population growth
N.	Scenario ^a		r	SD	Persistence probability in 100 years	Mean time to extinction (years)	Mean population size in 100 years
1	Reference (scenario 1)		−0.013	0.116	0.831	80	26
2	Increased food availability (scenario 2)	+5% $\rho$, K	−0.008	0.107	0.89.6	80	34
3		+10% $\rho$, K	−0.003	0.099	0.942	81	42
4		+15% $\rho$, K	0.001	0.095	0.958	80	48
5		−10% Mc	−0.008	0.106	0.899	82	33
6	Reduced mortality (scenario 3)	−25% Mf	0.018	0.078	0.998	83	57
7		−50% Mf	0.040	0.070	1.000	-	62
8	Combined food increase and mortality reduction (scenario 4)	+10% $\rho$, K; −10% Mf	0.012	0.097	0.981	81	53

^a In all scenarios, environmental stochasticity was set to CV = 0.2 (see Table III).

Figure 3. Projected extinction probabilities for the Apennine brown bear (Ursus arctos marsicanus) population in a period of 100 years, under different future management scenarios: (a) progressive increase in the availability of bear food, that will mainly affect productivity and cub survival; (b) increasing efforts at reducing human-related mortality. The parameters used to construct the above scenarios included cub survival (M_c), carrying capacity (K), the proportion of females reproducing each year ($\rho$) and adult female survival (M_f).

Discussion

From a retrospective point of view, our PVA indicates that the core Apennine brown bear population is in all likelihood in a condition of demographic stagnation or slow decline. This supports all recent estimates of population size (for 2003: Gervasi et al. 2008; for 2008: Gervasi et al. 2012; for 2011: Ciucci et al. 2015a; for 2014: Ciucci et al. 2015b), and it is also consistent with a capture–recapture-based estimate of the realised population growth rate during 2003–2014 (Gervasi et al. 2017). Accordingly, the concurrent outcomes of different methodological approaches and datasets strongly support the claim that the conservation efforts implemented in the past few decades have not yet promoted the desired numerical recovery of the population in the core of its distribution (Ciucci & Boitani 2008; Anonymous 2011). From a prospective point of view, simulations over the next 100 years indicate that the demographic stagnation of Apennine bears has a high probability of persisting unless the demographic rates of the population improve in the near future. Under current conditions, the extinction risk in the next 100 years is not trivial (between 11 and 21%, depending on environmental stochasticity), and there is a substantial risk that the core population will be markedly smaller than the current size in a few decades. Compared to similar applications on small bear populations (e.g. Wiegand et al. 1998; Chapron et al. 2003; Martínez Cano et al. 2016), the strength of our findings is that all the main model’s parameters have been formally estimated in recent years.

The use of both a deterministic and a stochastic population projection model produced complementary information regarding the expected trend and the persistence probability of Apennine bears. Deterministic matrix models do not account for temporal variation in vital rates (i.e. environmental stochasticity), nor for the random fluctuations in the realisation of the same parameters due to the small size of a population (i.e. demographic stochasticity; Caswell 2001). This only allows such models to mechanistically estimate the intrinsic demographic performance of a population, but strongly limits their application when projecting population trends in the long term. For Apennine brown bears, the deterministic model tells us that the current vital rates of the core population are not compatible with a more than negligible numerical increase. In practical terms, by neglecting the additional threats caused by demographic stochasticity, the further loss of genetic diversity and deleterious genetic effects, the risk of disease outbreaks and the effect of environmental variation, this bear population is likely to remain small and exposed to a relatively high risk of extinction, if the average survival or reproductive rates do not increase.

The complementary contribution of the individual-based stochastic population model was to add realism to the estimated population growth rate by contemplating the net effect of individual and environmental variability in demographic rates. The role of stochasticity in extinction processes, as defined by the so-called small-population paradigm (Caughley 1994), is to randomly amplify temporal fluctuations in population size and vital rates (Beissinger & McCullough 2002). While stochastic processes always increase the variance and unpredictability of population changes over time (Engen et al. 1998), in small populations they can trigger extinction according to an inverse exponential relationship between population size and extinction probability (Ovaskainen & Meerson 2010). Stochastic individual-based models have often been chosen for their flexibility in incorporating structural knowledge about the life-history of the study species, and for their ability to reveal the underlying demographic mechanisms of population rate of change. They were employed to estimate extinction probability and minimum viable population sizes for several terrestrial mammals, including grizzly bears (Ursus arctos horribilis) in North America (Knight & Eberhardt 1985), tigers (Panthera tigris) in Nepal (Kenney et al. 1995) and brown bears in the Spanish Cantabrian mountains (Wiegand et al. 1998). As expected, the stochastic model produced a lower estimated population growth rate for Apennine brown bears than the one provided by the deterministic version. Although small in absolute terms, this difference marks the distinction between a marginally positive (deterministic) and a slightly negative (stochastic) population trend. In both cases, it is unlikely that, under the current conditions, the bear population in the core of its distribution would act as a source facilitating a range expansion. This is consistent with accumulated evidence of a very slow population expansion, if any, beyond the historical core range, with no stable reproductive nuclei in the peripheral portions of the range (Ciucci et al. 2017).

Sensitivity of population growth rate to changes in vital rates differed between the deterministic and stochastic models. Whereas sensitivity analysis through the stochastic population model emphasised the relative importance of cub survival, when compared to that of adult females (Table II), population growth rate through the deterministic model was essentially a function of adult female survival, as expected for species with long generation times (Gaillard et al. 2005). As a matter of fact, the loss of one reproductive female required on average almost 12 years for its replacement. As a consequence, cub survival was 10 times less important in affecting population growth rate than adult female survival according to the deterministic model, but only 2.5 times less when considering demographic stochasticity.

The comparison between different management scenarios showed that conservation actions that would translate into a reduction in adult females’ mortality are those with the highest probability of rapidly generating an increase in population growth rate and, expectedly, range expansion (Table IV). Alternative management scenarios, especially those contemplating measures aimed at increasing food availability, all exhibited a marginal effect on population growth rate, or at best a slow turnaround in the numerical decline of the population (Table IV). Our projections suggest that they should not be disregarded per se as ineffective or unworthy, but they should be seen more as complementary conservation interventions to the main effort of effectively reducing mortality, rather than stand-alone strategies for the recovery of the population (see scenario 7, Table IV). We therefore suggest that economic investment in this direction should not be considered a first conservation priority and should be contemplated only in combination with an effective set of actions aimed to reduce human-caused mortality.

Our projections under different management scenarios also emphasise the susceptibility of the Apennine brown bear population to possible catastrophic events (Figure 2). Just one catastrophic event in the next 100 years would be enough to double the probability of population extinction, irrespective of the expected interaction between demography and the negative genetic effects that a sudden reduction in population size could cause. For the scope of our study, we defined a catastrophic event as a 50% increase in the annual mortality rate (Reed et al. 2003), with no specific reference to its underlying cause. In this sense, and given the ecological and socio-economic context in which Apennine brown bears live, human-related mortality causes should be considered the first candidate for the risk of a catastrophic die-off to occur. Even though the number of Apennine brown bears found dead each year has remained constant during the last few decades, there have been years when human-caused bear fatalities far exceeded average levels (e.g. 32 bears were found dead between 1980 and 1985, with 14 bears retrieved dead in 1982 only; Ciucci & Boitani 2008). Moreover, the possibility that disease outbreaks transmitted by livestock and other domestic animal could cause an increase in bear mortality should not be discounted. Although there is no direct evidence of disease outbreaks in other bear populations (Fey et al. 2015), the spatial proximity between livestock, stray dogs and bears in the PNALM (Ciucci & Boitani 2008), and the evidence of several bears being positive for canine distemper virus (CDV) and Brucella spp. (Marsilio et al. 1997; Di Francesco et al. 2015), should be considered a strong enough premise for the risk of disease-related mortality. Such risk should also be evaluated in the light of the very small Apennine brown bear population size. When only about 14 females of reproductive age are estimated to live in the population (Tosoni et al. 2017b), even a few of these females dying of disease could have serious demographic consequences on population performance and persistence.

When interpreting projections of our viability analysis, it should be considered that the risk of inbreeding depression, which we did not contemplate in our population projections, is likely relevant and expected to further decrease the chances of future persistence of the Apennine brown bear population (Lorenzini et al. 2004; Ciucci & Boitani 2008). Therefore, we caution that persistence probabilities for scenarios predicting a decrease in population size actually underestimate the additional threat represented by a further loss of genetic variation in the population. Even under these liberal conditions, our findings suggest that there is little chance for Apennine brown bears to significantly expand their range beyond the PNALM in the near future, despite the fact that this has been recognised as a fundamental goal of any conservation strategy for this small bear population (Anonymous 2011; Ciucci et al. 2017). Promoting such geographic expansion also means enhancing our understanding of the spatial and demographic dynamics between the core and the periphery of the species distribution. At present, we know that about half of the bear cubs born each year in the core population are missing after 1 year, even though it is not clear if this is due to local mortality or to a high emigration rate (Gervasi et al. 2017). Although only one ascertained case of sexually selected infanticide (SSI) has been reported for the Apennine brown bear population during the last 10 years (P. Ciucci, personal observation), the possibility that human-related mortality could induce frequent social disruption and induce males’ predatory behaviour towards cubs should not be a priori disregarded. The importance of SSI in the demography of brown bear populations has been empirically shown in Scandinavia (Swenson et al. 2001), and theoretically put in evidence for the endangered Pyrenean population (Chapron et al. 2009). Therefore, identifying the relative magnitude of local cub mortality vs. emigration will contribute to clarifying the spatial and demographic structure of this bear population, and it will favour proper allocation of conservation resources to those actions that maximise the chances of population expansion. For the same reasons, it will be necessary to clarify what mortality rates and causes await young bears eventually leaving the core for the more peripheral portions of the range.

To enhance our chances of success in the conservation of Apennine brown bears, it is crucial to conduct a monitoring programme (sensu Elzinga et al. 2001) designed to reveal in due time the direction and extent of the effects of conservation interventions on the bear population. Towards this aim, the findings from our demographic projections are useful to inform such monitoring in three respects: (i) population size should be estimated at time intervals which are consistent with the most likely generation time of the population (i.e. 11.4 years). As future estimates of population size correspond to a given probability of persistence (Figure 1), this in turn should be used to evaluate the success of conservation interventions. More specifically, it is crucial that no further reduction in population size, compared to the current numbers, occurs during the next few decades, as this would correspond to a rapid increase in the extinction probability (p > 0.2; Figure 1(c)); (ii) besides the population trend, monitoring efforts should be focussed on the most relevant vital rates as the inherent causes of population change. Our analysis suggests that particular effort should be employed to monitor survival of the adult females, their reproductive performance, and cubs’ survival. Accordingly, compared to current rates, a combined 10% reduction in adult female mortality and a 10% increase in recruitment, through an increase in either reproductive rates or cub survival, represents a realistic conservation and monitoring objective to pursue for the next decade (see scenario 4, Table IV), as this would further reduce extinction risk (Figure 3) while allowing for population expansion; and (iii) monitoring and demographic analyses should be expanded to the more peripheral portions of the range, to assess the spatial structure on a wider scale and eventually integrate source–sink dynamics into conservation planning. Following the fate of juvenile bears dispersing from the core population is a priority to start shedding light on these processes.

Acknowledgements

We are grateful to E. De Matthaeis and E. Duprè for their valuable support. We thank Guillaume Chapron and an anonymous referee for their useful comments and suggestions on a previous version of this manuscript.

Disclosure statement

No potential conflict of interest was reported by the authors.

Additional information

Funding

This work was funded by the Italian Ministry of the Environment and the Unione Zoologica Italiana (UZI).