Determinants and patterns of habitat use by the brown bear Ursus arctos in the French Pyrenees revealed by occupancy modelling

Abstract

The Pyrenean brown bear Ursus arctos population in the mountains between France and Spain is one of the smallest and most threatened populations of large carnivores in Europe. We assessed trends in brown bear habitat use in the Pyrenees and investigated the underlying environmental and anthropogenic drivers. Using detection/non-detection data collected during 2008–2014 through non-invasive methods, we developed dynamic occupancy models, accounting for local colonization and extinction processes. We found two non-connected core areas of occupancy, one in the west and the other in the centre of the Pyrenees, with a significant decrease in habitat use overall during 2008–2014. We also found a negative correlation between human density and bear occupancy, in agreement with previous studies on brown bear habitat suitability. Our results confirm the Critically Endangered status of the Pyrenean population of brown bears.

Introduction

Over the last 5 decades populations of large carnivores have been recovering in Europe following the implementation of conservation policies (Chapron et al., 2014). In parallel, conflicts have arisen between people and wildlife (Treves & Karanth, 2003). The main sources of conflict are damage to livestock, and competition with local hunters for prey species (Ericsson & Heberlein, 2003; Gunther et al., 2004; Piédallu et al., 2016a). For these conflicts to be resolved, or at least mitigated, the expectations of all stakeholders need to be considered, and management decisions should be based on reliable ecological data (Redpath et al., 2013).

The brown bear Ursus arctos is widely distributed throughout Europe, with multiple populations of various sizes and ranges (Swenson et al., 2011), including a large population in Sweden (Kindberg et al., 2011) and a smaller one in the Italian Apennines (Gervasi et al., 2012). One of the smallest and most threatened populations is in the Pyrenees mountains, between south-western France and north-eastern Spain, and is categorized as Critically Endangered based on IUCN Red List criteria (Huber, 2007). Only five individuals were detected in 1995, and individuals were translocated from Slovenia in 1996–1997, 2006 and 2016 as a necessary measure to ensure the survival of the population. The population remains small, and is threatened by demographic stochasticity and inbreeding (Chapron et al., 2009; Swenson et al., 2011).

The distribution of a wild population is a key element in determining its conservation status (IUCN, 2001). However, this is difficult to assess in the case of elusive species with large home ranges (Gittleman & Harvey, 1982), such as the brown bear. To infer the distribution of large carnivores, monitoring often relies on tracks and indirect observations coupled with DNA analyses to identify the species (e.g. Taberlet et al., 1997; McDonald, 2004; Bellemain et al., 2005). The distribution of the brown bear in France remains poorly studied. In Spain, Martin et al. (2012) conducted habitat suitability analyses at a coarse scale on the Cantabrian population and applied their findings in a study of the Pyrenees population, using presence-only methods. Here, we build on these results to address two issues in standard species distribution models.

Firstly, when dealing with free-ranging populations, species detectability is probably < 1, which can lead to false negatives where animals are present but not detected during the survey (Kéry, 2011). Falsely assuming perfect detection can lead to an underestimation of the species’ distribution (Lahoz-Monfort et al., 2014). Site-occupancy models were developed specifically to distinguish between non-detection and absence by modelling the imperfect, possibly heterogeneous observation process (MacKenzie et al., 2002). Secondly, standard species distribution models assume that a species always occupies the most favourable area, and that this area can be reached by dispersal, based on the ecological niche concept (Leibold, 1995). However, natural barriers or dispersal limitations (such as small population size) may prevent a species from reaching a favourable area (Araújo & Guisan, 2006). To address this, static occupancy models were extended to account for colonization and extinction processes: so-called dynamic or multi-season occupancy models (MacKenzie et al., 2003). Although static occupancy models have often been used to study large carnivores (e.g. Carroll et al., 2003; Carroll & Miquelle, 2006; Bayne et al., 2008; Hines et al., 2010), there have been only few applications of dynamic occupancy models (Molinari-Jobin et al., 2012; Miller et al., 2013).

We identified environmental or anthropogenic drivers and trends in habitat use by brown bears in the French Pyrenees. To do so we fitted dynamic occupancy models to detection/non-detection data obtained through a multi-source systematic monitoring protocol during 2008–2014.

Study area and bear population

The study was conducted on the French side of the Pyrenees at the border between north-eastern Spain and south-western France (Fig. 1). The bears here mostly descend from individuals that were translocated from Slovenia to the Pyrenees in 1996–1997 (two females and one male) and 2006 (four females and one male), although the mother of one belonged to the remnant of the original Pyrenean bear population, which was thought to comprise five individuals in 1995. Field observations indicate that two population cores exist on the French side of the Pyrenees: two male bears have been accounted for in the western core, with the rest of the population found in the central core (Fig. 1).

Fig. 1 The counties of the French Pyrenees, where brown bear Ursus arctos monitoring was conducted. Mountain subsections are delineated in grey.

Methods

Data collection and monitoring

Data were gathered during 2008–2014 by members of the national Brown Bear Network (228 professional members from government agencies and 135 unaffiliated amateur members) under the supervision of the French Game and Wildlife Agency. A systematic monitoring protocol was followed using fixed transects along which the agents looked for bear sign such as hair, scats, claw marks or paw prints. The study area was divided into subsections based on ridge lines and valleys, and the mean area of a subsection was 95 km². The home range of brown bears is c. 85–200 km² for males and and c. 50–100 km² for females (Huber & Roth, 1993; Preatoni et al., 2005); however, the upper limits of these home ranges include the rut during May–June, whereas our analyses were focused on July–November, when home ranges are smaller (Preatoni et al., 2005). Each of the 84 subsections investigated included one transect. Each transect was visited at least once per month during July–November. The length of each transect was proportional to the area of the corresponding subsection, with 0.2 km surveyed per km² of subsection. Signs and observations were validated by experts from the French Game and Wildlife Agency, thus minimizing the risk of false positives resulting from species misidentification (Molinari-Jobin et al., 2012).

Model building and selection

To estimate the probability of bear presence in all the subsections, we built a dynamic occupancy model (MacKenzie et al., 2003) that was parametrized with the probabilities of colonization, γ (the probability of a subsection that was unoccupied the previous year becoming occupied); extinction, ε (the probability for a subsection that was occupied the previous year becoming unoccupied); and initial occupancy, ψ ₁ (the probability of a subsection being occupied in the first year of the study); along with the species detection probability p (the probability of detecting the species in a subsection when present). We used years as primary occasions, between which colonization and extinction probabilities could be estimated, and the months July–November as secondary occasions, during which we considered the occupancy status of a subsection to remain unchanged (the closure assumption). By focusing on the July–November period we excluded the breeding season (April–June), during which male bears in particular are known to increase their movement range while they look for females (Clevenger et al., 1990). Despite this precaution movements in and out of the subsections may still occur and, assuming these movements are random, occupancy should be interpreted as habitat use rather than the proportion of area occupied by the species (MacKenzie & Nichols, 2004). More precisely, the usage made of various habitat components within the home range is usually referred to as third-order selection (Johnson, 1980).

We relied on previous habitat suitability studies of brown bears in Europe to select candidate environmental and anthropogenic covariates for our analysis (Mertzanis et al., 2008; Martin et al., 2010, 2012). We considered five environmental and anthropogenic covariates for each subsection (Table 1; Supplementary Fig. S1). Roughness was calculated as the mean of the absolute differences between the altitude of a subsection and that of its contiguous subsections (Wilson et al., 2007). We used the BD ALTI database (250 m resolution; IGN, 2016a) to calculate the mean altitude of each subsection. Forest cover and shrub cover covariates were extracted from the CORINE Land Cover database (CLC, 2012). Road length was extracted from the ROUTE 500 database (IGN, 2016b). Human density was obtained from the NASA Socioeconomic Data and Applications Center (CIESIN & CIAT, 2005). The maximum correlation between these covariates was 0.51 in absolute value. We used Akaike's information criterion (AIC; Burnham & Anderson, 2002) for covariate selection, and the difference in AIC (ΔAIC) to compare model support with reference to the model best supported by the data. To account for model selection uncertainty we resorted to model averaging, considering all models with ΔAIC < 2. Given the large number of covariate combinations, we used a multistage approach to model selection (Dugger et al., 2011; MacKenzie et al., 2012; Lee & Bond, 2015). We began by selecting the best model structure, focusing on time-varying covariates only, namely year and survey. We considered eight models in total, with either no effect (.) or a year effect on colonization and extinction, and either no effect (.) or a survey effect (where a survey refers to a month, hence a survey-specific covariate) on detection probability (Table 2). As the sampling effort was homogeneous over the study period, we did not consider a year effect on detection. Then, based on previous bear occupancy studies (Nielsen et al., 2006, 2010; Mertzanis et al., 2008; Martin et al., 2010, 2012) and bear biology, we considered specific combinations of the environmental or anthropogenic effects on each of the parameters (ψ ₁, γ, ε and p; Table 1). We tested possible negative effects of the covariates human density and road length on initial occupancy ψ ₁, as a previous study showed that bears avoided human-caused disturbances (Naves et al., 2003; Mertzanis et al., 2008; Martin et al., 2010). Roughness, shrub cover and forest cover were all positively associated with bear presence in previous studies, albeit performed at different scales (Naves et al., 2003; Apps et al., 2004; Nellemann et al., 2007; Martin et al., 2010, 2012). For colonization, we studied possible effects of roughness and human density, which were the most commonly significant covariates in previous studies of bear distribution (Martin et al., 2010). For extinction, we considered the possible effect of the two anthropogenic covariates, human density and road length. We tested the possible effect of human density, roughness and forest cover on detection, as all three could potentially influence the accessibility of bear sign to observers. To account for the variability in subsection size we also included area as a covariate on detection in all models without submitting it to selection. In total, we fitted all possible 8,192 models.

Table 1 Definition of the environmental variables used for the occupancy analysis of brown bears Ursus arctos in the French Pyrenees, and the parameters for which an effect of the covariate was tested. Blank cells indicate the effect was not tested; +/− indicates the predicted sign of the effect of the covariate based on previous studies (see text for references).

Table 2 Model selection with time-varying covariates, with models ranked by the Akaike information criterion (AIC).

¹ψ ₁, initial occupancy probability; γ, colonization probability; ε, extinction probability; p, detection probability; year, year effect on the parameter, which relates to changes between primary occasions (i.e. from 1 year to another in our case); survey, survey effect on the parameter, which relates to the secondary occasions repeated within a year

²ΔAIC, difference between the AIC of the current model and that of the model with the lowest AIC

Given the uncertainty in the selection of the best set of covariates, we resorted to model-averaging for parameter estimation and inference (Burnham & Anderson, 2002). Effect sizes were examined to determine the magnitude of a covariate effect (Nakagawa & Cuthill, 2007). The covariates were standardized prior to the analyses.

To assess any trend in habitat use over the years we estimated the occupancy status of each subsection for each year, and then tested a linear effect of year on the binary occupancy variable using a conditional autoregressive correlation model and an adjacency matrix between the subsections to specify the correlation matrix (Rousset & Ferdy, 2014). A likelihood ratio test was performed to assess the significance of this temporal trend. We applied this procedure to all models with ΔAIC < 2.

We compiled maps for initial occupancy, detection, colonization and extinction by calculating the probability for a given subsection using the model-averaged parameter estimates and the value of the covariates for the subsection.

Analyses were performed in R v. 3.3.0 (R Development Core Team, 2013), using the unmarked (Fiske & Chandler, 2011), spaMM (Rousset & Ferdy, 2014), rgdal, AICcmodavg, classInt, RColorBrewer and spdep packages. The data and R codes are available at http://dx.doi.org/10.5281/zenodo.572637.

Results

Multistage model selection

We found no year or survey effects on any of the parameters ψ ₁, γ, ε or p (Table 2). The ΔAIC of the next two best models (with a year effect on extinction and a survey effect on detection) was > 2, and therefore we used the model with constant parameters as the basic structure for the next step. Despite model uncertainty in the results of the selection procedure on environmental and anthropogenic covariates, some covariates were always included in models with ΔAIC < 2 (Table 3): roughness on detection and colonization probabilities, and human density on extinction and initial occupancy probabilities.

Table 3 Model selection, with environmental and anthropogenic covariates (Table 1). We report only models with ΔAIC < 2. The area of subsections was used in all models in the detection probability, and was not subject to the covariate selection procedure.

The effect of covariates on parameters

We refined the patterns found in the covariate selection step by examining the effect sizes (on the logit scale; Fig. 2). The effect of roughness on detection probability and that of human density on both extinction and initial occupancy probabilities were confirmed, but the colonization probability was not associated with any covariates. We investigated further the links between the covariates and initial occupancy, colonization, extinction and detection probabilities by inspecting the form of these relationships (after back-transformation; Fig. 3). An increase in roughness was associated with an increase in the detection probability, whereas it was more difficult to detect bears (when present) in large subsections (Fig. 3c, d). Initial occupancy was strongly negatively correlated with human density (Fig. 3a), with the least populated areas being more likely to be used by bears. Extinction was also negatively correlated with human density (Fig. 3b).

Fig. 2 Model-averaged parameter estimates (on the logit scale) and confidence intervals of the (standardized) covariates effects (see definitions in Table 1) on (a) initial occupancy, (b) colonization, (c) extinction and (d) detection probabilities. An effect does not appear if the corresponding covariate was not considered in the selection procedure. AREA was used in all models in the detection probability and was not subject to the covariate selection procedure.

Fig. 3 Relationships between model-averaged parameter estimates and the most important standardized covariates (based on the effect sizes in Fig. 2). The parameters of interest are on the y-axes: (a) initial occupancy, (b) extinction and (c–d) detection probabilities. The colonization probability is not displayed because of small effect sizes (see Supplementary Fig. S2 for the relationships between parameters and all covariates). For each relationship, we used the mean values of non-focal covariates. The dashes on the x-axes indicate the observed covariate values.

Distribution maps

The initial occupancy map (Fig. 4b) indicates two population cores (western and central), with the central core extending into south-east Ariège, south-west Aude and Pyrénées Orientales (Fig. 1). The extinction probability in the east of the central core was high (Fig. 4d), which is consistent with the disappearance of bears from that area (Camarra et al., 2012), whereas the colonization probability in the same subsections was close to zero (Fig. 4c). Detection was higher in the central core than in the western core (Fig. 4a). The colonization map indicates that the western population core is more likely to expand to the east, whereas the central one is more likely to expand to the west (Fig. 4c), which is confirmed by yearly occupancy maps (Fig. 5). Occupancy in the west of the western population core (south-west of the Pyrénées Atlantiques) also decreased, while remaining constant in the east of that core (south-east of the Hautes Pyrénées). Overall, we detected a decrease in habitat use by the brown bear during 2008–2014 (P < 0.01 for all models in Table 3), with no new areas being colonized.

Fig. 4 Maps of the model parameters in the mountain subsections of the French Pyrenees, using the model-averaged parameter estimates and mean values of covariates: (a) detection probability, (b) initial occupancy probability, (c) colonization probability, and (d) extinction probability.

Fig. 5 Yearly occupancy probability (ψ _t) of brown bears in mountain subsections of the French Pyrenees from t = 2008 to t = 2014, obtained using the model-averaged parameter estimates, mean values of the covariates, and the formula ψ _t + 1 = (1 − ψ _t)γ  + ψ _t(1 − ε) (MacKenzie et al., 2002).

Discussion

Environmental and anthropogenic effects on model parameters

Human density had a strong negative effect on initial occupancy probability, with the least densely populated areas being the most likely to be used by bears. This result confirms previous analyses suggesting that bears tend to live far from areas with intense human activity (Martin et al., 2010; Long et al., 2011). Several factors, such as the habituation of bears to human presence (Wheat & Wilmers, 2016) or the need for female bears to shield themselves from sexual conflict (Steyaert et al., 2016), may mitigate this effect, but the small size of the Pyrenean brown bear population limits the immediate relevance of these factors, as bears tend to disperse further at low densities, lowering the rate of encounter between individuals and the risk of sexually selected infanticide (Stoen et al., 2006). These results confirm that anthropogenic effects supersede natural elements when it comes to habitat selection by brown bears (Nellemann et al., 2007).

Contrary to our expectation, human density was negatively correlated with the probability of extinction. A possible explanation is the influence of demographic stochasticity in small populations (Gabriel & Bürger, 1992), which gives more weight to extinction events (i.e. in a small population an extinction event may lead the entire population to collapse, which is not necessarily the case in a large population). Human density was lower in the south-east of Ariège and the south-west of Aude and Pyrénées Orientales (Supplementary Fig. S1) than it was in the other areas with high occupancy probability (Fig. 4b), and there were several local extinction events in 2010 and 2011 in these places with lower human density (Camarra et al., 2012).

We found a positive correlation between the detection probability and roughness. A rougher terrain funnelling pathways of bears and people may explain this pattern. The same funnelling effect may explain why signs of bears were easier to detect (when the species was present) in small subsections than in large ones. Overall, species detection was imperfect and estimated to be < 0.6, thus confirming the need to correct for it to avoid underestimating occupancy.

Brown bear habitat use in the French Pyrenees

The occupancy maps of bears in the Pyrenees indicate the existence of two independent population cores, one in the west and another in the centre of the mountain range (Figs 4b & 5). The two cores remained unconnected during the timespan of the study. The dynamics of occupancy over the study period (Fig. 5) showed that habitat use decreased significantly overall. In particular, the extinction of the species in the eastern part of the central core is consistent with the absence of bear tracks in south-east Ariège and south-west Aude and Pyrénées Orientales (Fig. 1) since 2011 (Camarra et al., 2012). These results demonstrate the usefulness of dynamic occupancy models in highlighting trends in habitat use that cannot be identified by static species distribution models (MacKenzie et al., 2003).

That we found many subsections with a high occupancy probability in the western core despite the fact that only 2–3 bears were estimated to live there during 2008–2014 (Piédallu et al., 2016b) suggests a violation of the closure assumption between our secondary occasions (July–November), because there were not enough bears in the population core to occupy all subsections at the same time. We estimated habitat use by brown bears instead of actual occupancy. For species that attack livestock, presence does not have to be permanent to be a source of conflict, and therefore habitat use remains a relevant indicator in the case of large carnivores, which are often characterized by relatively large home ranges (Gittleman & Harvey, 1982) and their use of large areas without actually occupying much land at any given time.

Implications for human–wildlife conflict mitigation

We anticipate that our results will contribute to the scientific evidence that is required for conflict mitigation (Redpath et al., 2013). Attacks on livestock are one of the main causes of negative attitudes towards carnivore presence in general (Kaczensky et al., 2004; Sponarski et al., 2013) and towards brown bears in the Pyrenees in particular (Piédallu et al., 2016a). There is an interest in mapping the areas that are more likely to host bears now and in the future, and the likely hotspots for attacks (Miller, 2015). This could be combined with mapping attitudes towards brown bears (Piédallu et al., 2016a) to identify areas that combine positive attitudes towards bear presence and low risk of attack, which could be primary targets of future management decisions. This could be the first step towards the development of socio-ecological models designed to mitigate human–wildlife conflicts (Aswani, 2011; Dupont et al., 2011; Estoque & Murayama, 2014).