http://orcid.org/0000-0001-5496-0144
Figures
Abstract
Many large carnivore populations are expanding into human-modified landscapes and the subsequent increase in coexistence between humans and large carnivores may intensify various types of conflicts. A proactive management approach is critical to successful mitigation of such conflicts. The Cantabrian Mountains in Northern Spain are home to the last remaining native brown bear (Ursus arctos) population of the Iberian Peninsula, which is also amongst the most severely threatened European populations, with an important core group residing in the province of Asturias. There are indications that this small population is demographically expanding its range. The identification of the potential areas of brown bear range expansion is crucial to facilitate proactive conservation and management strategies towards promoting a further recovery of this small and isolated population. Here, we used a presence-only based maximum entropy (MaxEnt) approach to model habitat suitability and identify the areas in the Asturian portion of the Cantabrian Mountains that are likely to be occupied in the future by this endangered brown bear population following its range expansion. We used different spatial scales to identify brown bear range suitability according to different environmental, topographic, climatic and human impact variables. Our models mainly show that: (1) 4977 km2 are still available as suitable areas for bear range expansion, which represents nearly half of the territory of Asturias; (2) most of the suitable areas in the western part of the province are already occupied (77% of identified areas, 2820 km2), 41.4% of them occurring inside protected areas, which leaves relatively limited good areas for further expansion in this part of the province, although there might be more suitable areas in surrounding provinces; and (3) in the eastern sector of the Asturian Cantabrian Mountains, 62% (2155 km2) of the land was classified as suitable, and this part of the province hosts 44.3% of the total area identified as suitable areas for range expansion. Our results further highlight the importance of increasing: (a) the connectivity between the currently occupied western part of Asturias and the areas of potential range expansion in the eastern parts of the province; and (b) the protection of the eastern sector of the Cantabrian Mountains, where most of the future population expansion may be expected.
Citation: Zarzo-Arias A, Penteriani V, Delgado MdM, Peón Torre P, García-González R, Mateo-Sánchez MC, et al. (2019) Identifying potential areas of expansion for the endangered brown bear (Ursus arctos) population in the Cantabrian Mountains (NW Spain). PLoS ONE 14(1): e0209972. https://doi.org/10.1371/journal.pone.0209972
Editor: Aristeidis Parmakelis, National and Kapodistrian University of Athens, GREECE
Received: March 27, 2018; Accepted: December 15, 2018; Published: January 4, 2019
Copyright: © 2019 Zarzo-Arias et al. This is an open access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.
Data Availability: All relevant data are within the paper and its Supporting Information files.
Funding: The research was financially supported by the Gobierno del Principado de Asturias (with FEDER co-financing) with the project CN-17-023 to VP and M.M.D. VP was financially supported by the Excellence Project CGL2017-82782-P financed by the Spanish Ministry of Economy, Industry and Competitiveness (MINECO), the Agencia Estatal de Investigación (AEI) and the Fondo Europeo de Desarrollo Regional (FEDER, EU). M.M.D. (RYC-2014-16263) and F.D. (RYC-2013-14662) were supported by Ramon & Cajal research contracts from the Spanish Ministry of Economy, Industry and Competitiveness. The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.
Competing interests: The authors have declared that no competing interests exist.
Introduction
As a consequence of the implementation of major conservation and management actions [1,2], many large carnivore populations are expanding into human-modified landscapes [3–7], which may provoke an increase in several types of conflicts, e.g., livestock predation, crop damage and, more rarely, attacks on humans [8–10]. Predicting potential range expansion areas is an important step towards proactive management strategies minimizing conflict, thereby enhancing large carnivore population viability [11–13]. This is particularly important for small and isolated populations that are confined as a result of the expansion of humans and habitat degradation, and for which spatial expansion is therefore essential for their conservation. Different habitat suitability models, like maximum entropy models [14–16], are nowadays commonly used for exploring the availability of favourable habitats and the likely spatial distribution of population expansions, as well as the environmental factors determining them [17–20]. Since reliable absence data is frequently difficult to obtain, these models are usually based on species presence-only occurrence data such as maximum entropy models. In human-dominated landscapes, where fragmentation and loss of good habitat, primarily large patches of continuous forest with little human encroachment [21–23], is continually arising, these models have become very popular in carnivore population studies and conservation [24–27], given the abovementioned difficulty in obtaining absence data [28]. Further, a key ingredient of these models is the spatial scale considered. While large-scale spatial models are useful for understanding broad population patterns and processes related to the distribution of a species, high-resolution models offer information on specific niche requirements of locally adapted populations [29,30]. Thus, combining different spatial scales provides the opportunity to address and improve our knowledge about the relationships between species and the environment by providing more accurate predictions on species distributions [22,31].
The brown bear Ursus arctos is one of the most widespread large carnivores in the world, occupying different countries in North America, Europe and the north of Asia (http://www.iucnredlist.org/details/41688/0). Brown bears were historically persecuted and nearly eliminated from much of Western Europe in the 20th century in order to avoid conflicts and as a result of hunting [5]. The effects of direct persecution were aggravated by other threats like habitat loss and fragmentation due to the expansion of the human population, which conflicts with the large spatial requirements of this species [32].
Northern Spain is home to the last two isolated populations of brown bear in south-western Europe, which have been protected for more than 30 years. The main population, which is estimated to consist of approximately 200 individuals (95% CI: 183–278; [33]) only, inhabits the Cantabrian Mountains (NW Spain) and is divided into two tenuously connected subpopulations [34], with most of the population inhabiting the western sector of the Cantabrian Mountains belonging to the province of Asturias (hereafter referred to as “Asturias”). Previous studies have shown that both subpopulations are increasing in number and range, especially in the western region [33–37]. This current positive trend may cause re-colonisations of areas where bears have disappeared and/or the occupation of new ones [3].
The small number of bears, their complete isolation and the limited connectivity between the two subpopulations make the study of the range expansion of this population particularly interesting. In addition, as its expansion might increase conflicts with human populations and activities (i.e. livestock predations, crop and apiary damages, attacks on humans), the identification of potential new colonization areas is of great importance for supporting conservation management actions. Here, by using maximum entropy models, we aim to identify the potential range expansion areas of Cantabrian brown bears. We have specifically: (1) identified suitable areas for bears at a coarse scale (5 x 5 km); (2) explored, at a finer scale (1 x 1 km), the best areas highlighted in the previous model; and (3) evaluated which environmental variables determine habitat suitability for this population on both scales.
Materials and methods
Study area
Asturias is one of the four regions of north-western Spain still inhabited by brown bears. The region is characterised by an oceanic climate with mild temperatures and high humidity, with annual mean temperatures ranging from 14°C on the coast, to 2–3°C on the highest points of the mountains (http://www.worldclim.org/). Asturias comprises 10,602 km2 with more than a million inhabitants distributed both in big cities and small towns from the coast to the Cantabrian Mountains. Population density is ca. 100 inhabitants/km2 (Instituto Nacional de Estadística http://www.ine.es/jaxiT3/Datos.htm?t=2886) and road density is 47.4 km/100 km2 (http://www.seap.minhap.gob.es/index.html). Elevation ranges from 0 to 2648 m a.s.l. (http://www.sadei.es/). The Cantabrian Mountains are principally covered by forests of chestnut (Castanea sativa), oak (Quercus petraea, Q. robur, Q. pyrenaica and Q. ilex) and beech (Fagus sylvatica), alternating with pastures and brushwood, and subalpine matorral above 1700 m [35,38].
Bear occurrence data
We used a database of brown bear occurrences within Asturias that consists of geolocalized direct observations; indirect signs of the presence of the species (i.e., footprints, hair and scats); and damage records caused by bears to livestock, beehives, crops, and human activities and infrastructures. The database covers observations from 1995 to 2016 and has been compiled by the regional government of Asturias (S1 Fig). Observation data came from several sources, primarily: (a) systematic direct bear observations by regional government field staff (the Patrulla Oso, i.e., the Bear Patrol, and all the guards of the Principado), as well as by the Brown Bear Foundation (FOP, Fundación Oso Pardo), the Asturian Foundation for the Conservation of Wildlife (FAPAS, Fondo para la Protección de los Animales Salvajes) and personal observations of the authors; and (b) camera traps that were randomly located by the FAPAS and Bear Team during the last twenty years, mainly in forested areas where bears are less visible. Any indirect observations, i.e. tracks, signs and damages, were done by trained personnel. We removed observations with obviously erroneous or doubtful spatial locations (e.g., incomplete coordinates or poorly georeferenced observations).
Definition of utilized brown bear range and identification of potential expansion areas
We evaluated the potential areas of bear expansion by using two different spatial scales [31]. We first defined a coarse scale of 5 x 5 km (25 km2), which is approximately the average size of Spanish brown bear home ranges [29,38,39]. For this coarse scale, we opted to model the suitability for distribution range rather than for occurrence [13], to obtain a primary and more general distribution of the favourable habitat for the species. Therefore, we first defined the distribution range of Asturian brown bears as pixels with 3 or more years of presence. We binarily classified the model into suitable and unsuitable areas. Within areas identified as suitable range from the coarse scale model, we then considered a model with a finer scale of 1 x 1 km to enable identification of suitable areas for brown bears within their area. For this scale, we therefore used raw observations as model input so that the model describes suitability for bear presence within areas that are regarded as within their suitable range.
Environmental variables
We used 25 layers related to human infrastructures, vegetation and geomorphology as predictors in our habitat suitability models (S1 Table). Climatic variables came from the WorldClim 2.0 database (http://www.worldclim.org/) and were described by [40]. We used a discrete land cover layer available from the Cartografía Temática Ambiental del Principado de Asturias 1989–1998 (1:50,000), which we converted into percentage area within each pixel. We discarded any land cover class that did not possess at least 1% mean occurrence. We used each of the remaining classes as separate layers for our models, where each layer describes the relative abundance of a particular class. We also utilized the total number of land cover classes per pixel as well as the Shannon index of relative occurrences. We employed a Normalized Difference Vegetation Index (NDVI) as an index of greenness, which was provided by the Instituto de Recursos Naturales y Ordenación del Territorio (INDUROT). We used a digital elevation model from MDT200. From this model we derived elevation and mean aspect of the slopes. We obtained human population density from the Sociedad Asturiana de Estudios Económicos e Industriales (http://www.sadei.es/es/portal.do) and highway, road and footpath density as well as river density from geophysical layers maintained by the Principado of Asturias.
In case of a spatial correlation > 0.7 (Pearson coefficient) between two layers we only retained the layer that we regarded as most relevant for bear biology. For instance, we removed all climatic variables except seasonal variability in precipitation, since they were all correlated with elevation which is more ecologically relevant for bears [41,42].
We retained 19 uncorrelated variables: % shrublands, % gorse, % heath, % fern % forests, % pastures, % planted forest, % planted conifer forest % cliffs, mean value per pixel of the Shannon index of land cover heterogeneity and the normalized difference vegetation index, number of land cover classes, elevation, slope, precipitation seasonality, river, human population, highway, road and footpath density. All the variables were projected to the same reference system (ETRS89 / ETRS-UTM29) and scaled to a 5 x 5 and a 1 x 1 km resolution. For the final modelling, we removed any cells which had less than 50% of its area within the limits of the study region in each scale (most of these cells occurred along the coast) and we included all the variables mentioned above.
Modelling with MaxEnt
We used the software MaxEnt version 3.4.1 [15] called from the statistical environment R [43] version 3.3.3 using the packages dismo version 1.0–12 [44] and ENMeval version 0.2.0 [45]. Because of the limited number of observations in eastern Asturias, which could have been partly due to observer bias, we only trained models on data from east of the A-66 and AS-1 highways (Fig 1). These roads appear to function as partial physical barriers for the dispersal of brown bears between west and east subpopulations [31,34]. With this restriction, we retained 4878 bear locations which were clustered into 99 5 x 5 km pixels of bear range that were used as training locations for the coarse scale model, and 1091 1 x 1 km pixels within the suitable bear range which contained observations that were used for the finer scale model. For each scale, we utilized the centre coordinate of each cell, which was either classed as bear home range (coarse scale) or that contained a bear observation (finer scale), as occurrence data. We used 500 iterations, a convergence threshold of 10−5 and values from all cells in the west section of Asturias to build our first two models. Then, for both scales, we projected the model output onto the whole of Asturias once an optimal model structure had been selected (see below) using values from the eastern section of Asturias for the coarse scale and identified suitable range in the eastern part for the finer scale as background points.
Binary classification of a coarse (5 x 5 km) MaxEnt model broadly identifying suitable range areas, including favourable areas that are part of the current bear range, as well as protected areas (a); and model output from the fine scale (1 x 1 km) MaxEnt model identifying the probability of bear occurrence within the identified suitable bear range in the coarse model (b).
Selecting and evaluating the models
MaxEnt associates the presence data to environmental values using 5 different feature types (linear, quadratic, product, threshold and hinge), which represent different types of parameterizations. These feature types represent different transformations of the covariates, allowing the modelling of potential complex relationships between the variables and preventing overfitting [46]. The mean of each feature has to be close (within some error limits) to the empirical average over the presence locations [47]. The software also controls over-fitting using a regularization parameter which penalizes variables with low contribution to the model. Although machine learning algorithms such as the one used in MaxEnt generally favour more complex model solutions than likelihood based algorithms, overfitting can still be problematic since models essentially parameterize random spatial noise [48]. Hence, to identify an optimal model structure we evaluated candidate models with all types of feature combinations, each run over a set of regularization multipliers ranging from 0 to 10 (S2 Fig (13)). Each model included the same set of 19 uncorrelated environmental variables. We identified the best combination of feature types and regularization multiplier using Akaike’s Information Criterion corrected for small sample sizes (AICc) [49]. We calculated the AICc values from the raw model output where the sum of the log transformed raw values were treated as equivalent to model likelihood [48]. Although potentially at odds with machine learning philosophy, AIC values have been shown to be an efficient and reliable method for identifying optimal levels of model complexity in MaxEnt applications. Following [50], we regarded models within 2 AICc units of each other as having equivalent empirical support. Of these models, we selected the model containing the least number of parameters, and if two models had the same number of parameters we chose the model with the lowest number of feature types. We regard this quantitative method of model selection beneficial, since it provides objective and rigidly defined criteria for evaluating the numerous alternative models associated with any given MaxEnt exercise.
To evaluate model performance we first used the AUC (Area Under a Receiver Operating Characteristic–ROC–Curve) value [51], which indicates how efficiently a model differentiates between occurrence and background locations. AUC values from 0.7 to 1 generally suggest that the model has adequate predictive ability [52]. Second, we calculated three model performance metrics based on cross validation using a checkerboard method to separate our occurrence data into training and testing [45]: AUCtest (ability of testing locations to distinguish between background and presence locations), AUCdiff (difference in the ability to distinguish between presence and background locations between training and test data [48] and ORmin (proportion of test locations with a value below the lowest value of training locations [45].
Evaluating variable contributions
To evaluate the relative contribution of each environmental variable to the models we used a jacknife procedure and a heuristic method provided by MaxEnt. The jacknife test shows the gain in AUC value of each variable when used in isolation and the lack of gain when removed from the full set of variables (S3 Fig). The heuristic method calculates the percent contribution of each variable as the proportional contribution to the model training gain for every iteration of the model fitting process [14].
Presentation of model output
For both models we used the complementary log-log (cloglog) format as model output, as it has an intuitive interpretation and is monotonically related to other potential output formats [15]. This format allows interpreting model output as a probability of occurrence. However, as we were interested in suitable areas for bear range in our coarse scale model, and not necessarily the relative suitability within identified range, we have presented the coarse scale model as a binary classification, which broadly identifies favourable bear range. For this purpose we used the 10 percentile training presence in the cloglog values as a threshold for suitable bear range. This threshold selects the value above which 90% of the training locations are correctly classified, and is one of the most common thresholds used in MaxEnt habitat suitability models [53]. However, we have presented the unclassed cloglog output of our coarse scale model in S4 Fig. For the fine scale we maintained the cloglog values in order to represent habitat suitability for bear occurrence.
Identification of potential bear expansion areas and evaluation of the suitability of used range and expansion areas
For the coarse scale, cells classed as suitable but which were not part of the identified distribution range were regarded as potential expansion areas. To evaluate whether Asturian brown bears have expanded their range incrementally out of an initial core area (S5 Fig), we calculated the NODF nestedess index as an index of spatial nestedness over time [54]. This value describes the extent to which cells included in the range of a given year also form part of the range in subsequent years [14], which can be regarded as a temporal analogue to spatial nestedness [55]. Following Eriksson and Dalerum [14], we calculated the nestedness on annually identified pixels recognized as suitable range. The index value ranges from 100, indicating complete nestedness, to 0, indicating an anti-nested pattern. To evaluate if our observed values differed from random expectations we compared our nestedness index to values derived from 999 null models constrained to retain the original marginal sums. We only conducted this analysis on the coarse scale, as we do not regard absences of observations within pixels in the finer scale to be especially informative.
We evaluated whether suitability at the fine scale differed between the used range and the potential range expansion areas using a linear model. We used the log transformed raw model output as a response variable, and a 3 level factor as a predictor. This factor consisted of the classes obtained from the coarse model “used range” (suitable areas already occupied by bears), “expansion areas west” (unoccupied areas west of highways A-66 and A-S1) and “expansion areas east” (unoccupied areas east of highways A-66 and AS-1). We opted to separate the expansion areas in eastern and western Asturias because we suspect there may be observation bias in the eastern part, which may have underestimated the used range. We also added up to 7th order polynomials of the spatial coordinates to the model predictors to account for spatial autocorrelation. We selected this level of complexity for eliminating spatial autocorrelation by selectively adding polynomial complexity until we could not detect further autocorrelation using Moran’s I values calculated on the residuals [56]. For our data set, a polynomial approach was more efficient in removing spatial autocorrelation than approaches directly defining spatial autocorrelation in the model correlation matrix [57] or approaches using spatial eigenvectors as predictors [58]. To explore pairwise differences between the classes we used least square means with a Tukey correction for multiple comparisons [59,60].
Results
MaxEnt model selection and evaluation
The optimal coarse scale model identified by the AICc values included linear features, a regularization multiplier of 10, and discarded 12 variables of the whole set because of their limited contribution (S2 Table). The fine scale model included linear, quadratic, product and threshold features with a regularization multiplier of 0, and included all the selected variables with the minimum contribution being 0.6. Table 1 shows the different values of the evaluation metrics (see Methods) for the five candidate models for each of the coarse and fine scales with the highest empirical support. The best models at each scale had mean AUC values of 0.782 and 0.7368 respectively, showing adequate predictive ability. However, several models of both coarse and fine scale were regarded as having equal empirical support (Table 1). More complex models generally showed clear signs of over fitting, whereas less complex models lost predictive abilities (S2 Fig). The 10 percentile method yielded a threshold of 0.359 in cloglog values, which we used to defined suitable areas for bear range in the coarse scale.
MaxEnt variable contribution
The two variables that contributed the most to the coarse scale model (5 x 5 km) were altitude and slope (Table 2), which were confirmed both using the jacknife tests and the heuristic evaluations (S3 Fig and S2 Table). Both had a positive effect on the suitability for bear distribution (S6 Fig). Some variables of human impact (i.e., highways and footpaths) contributed modestly to this model with a negative influence on habitat suitability for bears (S5 Fig). Noteworthy is that several variables had no influence on the coarse scale model. These included precipitation seasonality, human population density, NDVI, pastures, cliffs, conifer plantations, rivers, shrubland, land cover classes and roads. The finer scale model (1 x 1 km) was most influenced by percentage of forest cover, precipitation seasonality and human population density (Table 2), with forest cover and precipitation seasonality being positively associated and human population density negatively associated with probability of bear presence (S7 Fig). As with the coarse scale model, the variable contributions were confirmed by both the jacknife tests and the heuristic evaluations (S3 Fig and S2 Table).
The percentage values are based on a heuristic method that estimates the proportional contribution of each variable to the model training gain for every iteration during model fitting.
Predicted range distribution and potential expansion areas
Brown bears occupy a total of 2430 km2 in Asturias. Most bear territory was identified in the south western part of the province, and our coarse scale model confirmed that most suitable areas are in the southern parts of the province (Fig 1A). The NODF nestdeness value characterizing the type of spatial expansion of brown bears in Asturias was marginally less nested than random expectations (NODFobs = 18.42, NODFexp = 18.70 ± 0.136, z = -1.97, p = 0.05). This suggests that the expansion has been not been caused by annual range expansions in which individuals have settled in neighbouring areas, but rather that they have dispersed in different directions and hence moving the utilized range over time.
Our coarse scale model identified 4977 km2 as suitable bear range, which represents close to half of the territory of Asturias (Fig 1A). However, most of the suitable areas in the western part of the province were already occupied (77% of identified areas, 2820 km2), 41.4% of them occurring inside protected areas (Fig 1A). This leaves relatively limited good areas for further expansion in this part of the province, although there might be more suitable areas in surrounding provinces. In the eastern part of the Asturias, however, 62% (2155 km2) of the land was classified as suitable, and this part of the province boasted 44.3% of the total areas identified as suitable bear range. Forty-one percent of all expansion areas (suitable areas outside the current bear range) were identified in protected areas (Fig 1A).
Our finer scale model identified the most favourable areas in the southwestern part of the range (Fig 1B). In line with this observation, there were significant differences in suitability between the used areas and the unused areas in the western and eastern parts of the province (F2,2 = 821.8, p < 0.001), with the utilised areas having significantly higher suitability than the unused areas both in the west (t = 26.7, df = 1, p < 0.001) and in the east of Asturias (t = 26.8, df = 1, p < 0.001). The areas in eastern Asturias had, however, higher suitability than the unused areas in the western part of the province (t = 26.8, df = 1, p < 0.001) (Fig 2).
Discussion
Our study highlights that (a) a large portion of the Cantabrian Mountains belonging to Asturias is potentially suitable as bear range and (b) only most of the suitable western areas of the province have already been occupied. In fact, according to our coarse scale model, more than 75% of the favourable areas of western Asturias have records of brown bear presence. Because dispersal out of this area seems to be limited [34–36], if the population continues to expand as previously projected [33] it may either experience a local density increase or bears will expand into areas relatively far from the core of the population. Anecdotal information suggests that bears have already moved out of the most favourable areas and most of the observations made in cells outside of our classed bear range have been recorded in recent years.
Despite the consistent positive trend in the population size of the western subpopulation of the Cantabrian brown bear, the eastern subpopulation has had a substantially smaller population increase [35]. Although we cannot rule out that at least part of this lesser increase may be caused by observer bias, the stark contrast in population growth between the western and eastern subpopulations suggests high mortality rates in the eastern subpopulation and/or a relatively limited dispersal of bears out of western Asturias. Indeed, some previous studies have suggested a narrow connection between the two subpopulations [36,61]. Additionally, prior studies hypothesised a prevalence of low quality food-items in the eastern sector of the Cantabrian Mountains [39,62].
Our results also highlight the importance of using different scales to model habitat suitability. Indeed, the most important environmental variables differed between the coarser- and finer-scale models. This difference may be related to the relative influence of landscape features depending on the scale. While a broader spatial scale yields a general environmental description for the entire species distribution (even among different populations), a finer spatial scale is more related to local requirements [29]. On the one hand, altitude and slope were the most influential environmental variables in the coarse spatial scale model. This agrees with previous studies on brown bear distribution patterns [22,29,31,38,41,63,64]. Brown bears tend to appear in high and rugged areas, especially females with cubs trying to avoid infanticide. They also appear quite elusive nearby human settlements, in order to escape from disturbances produced by high human activity. Our results point to more inaccessible areas, where human density and activities are scarce, as suitable areas for the bears. Although with less influence, the coarse spatial scale model was also negatively affected by footpaths and highways. The latter has previously been negatively associated with habitat suitability for brown bears [22], suggesting that it may represent a limitation for brown bear dispersal.
On the other hand, forest cover (positive effect) and human population density (negative effect) were the most influential variables in the fine scale model. Indeed, forest cover represents crucial food and shelter for the Cantabrian brown bear, as one of its main food resources in the Cantabrian mountains are acorns [62,65,66]. On the other hand, human density may be related to bear avoidance of human disturbance, as proved in other studies showing that bears are more detectable further from human settlements, where human activity is greater [22,31,67,68].
We expect that the results of our models can also be used as a practical tool for brown bear damage prevention and conflict mitigation. Indeed, knowing the areas into which the population of the Cantabrian brown bear is likely expand in the near future would allow authorities and conservation organizations to focus information campaigns and pre-emptive damage control actions on these areas. Such proactive approaches are important for successful large carnivore conservation and management [69]. For Cantabrian brown bears, damages to apiaries are the main source of conflict, and damage prevention strategies have been shown to be effective to avoid them [70]. We suggest that our maps of potential brown bear range expansion areas should be overlaid with spatially explicit data on apiaries to allow for the identification of high risk areas where conflicts may occur. In addition, since brown bears have disappeared from certain areas some decades ago, local communities are no longer familiar with how to coexist with this large carnivore. Thus, local information campaigns directed at inhabitants of areas of potential bear expansion and based on studies like the present one may represent a crucial strategy to prevent human-wildlife conflicts.
Supporting information
S1 Fig. Brown bear occurrence data and location of the study area in Europe.
https://doi.org/10.1371/journal.pone.0209972.s001
(PDF)
S2 Fig.
Evaluation metrics for 130 candidate models containing different levels of complexity defined by a range of five feature type combinations including linear (L), quadratic (Q), product (P), threshold (T) and hinge (H) features, each evaluated over a range of regularization multipliers ranging from 0 to 10, for (a) the coarse and (b) fine scales of the distribution of the Cantabrian brown bear in Asturias. Evaluation metrics include delta AICc, which is the difference in AICc (Akaikes Information Criterion corrected for small sample sizes, calculated as the sum of the log transformed raw output penalized by the number of model parameters), AUC test, which is the AUC (area Under the receiving operator characteristics Curve) score for the testing data set, AUC diff, which is the difference in AUC scores between the training and testing data sets, and OR min, which is a threshold dependent statistic corresponds to the proportion of testing localities that have MaxEnt output values lower than the value associated with the training locality with the lowest value.
https://doi.org/10.1371/journal.pone.0209972.s002
(PDF)
S3 Fig.
Jacknife evaluations of variable contributions to the (a) coarse and (b) fine scale models. The variables with the highest gain when used in isolation are slope for the coarse scale (a) and forest cover foir the fine scale model (b). These variables therefore seem to have provided the most useful information by themselves for each scale. The variables that decreased the gain most when omitted, and thus possessed the greatest amount of information not present in the other variables, were slope for the coarse scale (a) and population density for the fine scale model (b).
https://doi.org/10.1371/journal.pone.0209972.s003
(PDF)
S4 Fig. Output of the coarse scale model with a 5 x 5 km resolution.
The map presents a clog-log transformation of the raw MaxEnt output, which can be interpreted as a probability of brown bear range occurrence.
https://doi.org/10.1371/journal.pone.0209972.s004
(PDF)
S5 Fig.
Schematic examples of incremental range expansion (a) out of an initial core area as well as (b) a patchy range expansion were no area is occupied two consecutive years, their nestedness values as well as the association matrices used to calculate nestedness.
https://doi.org/10.1371/journal.pone.0209972.s005
(BMP)
S6 Fig. Associations between predicted suitability estimated from the coarse scale model each of the included environmental predictors.
https://doi.org/10.1371/journal.pone.0209972.s006
(PDF)
S7 Fig. Associations between predicted suitability estimated from the fine scale model each of the included environmental predictors.
https://doi.org/10.1371/journal.pone.0209972.s007
(JPG)
S1 Table. Description, source and original format of the 25 environmental variables initially developed for the construction of the models.
Variables marked with * are the ones not correlated and ultimately used in the modelling.
https://doi.org/10.1371/journal.pone.0209972.s008
(PDF)
S2 Table. Variable contribution to the construction of the coarse and fine scale models.
https://doi.org/10.1371/journal.pone.0209972.s009
(PDF)
S3 Table. Centre coordinates of the 5 x 5 km grids classed as bear home range used as bear occurrence data in the coarse scale model.
https://doi.org/10.1371/journal.pone.0209972.s010
(CSV)
S4 Table. Centre coordinates of the 1 x 1 km grids that contained a bear observation used as bear occurrence data in the fine scale model.
https://doi.org/10.1371/journal.pone.0209972.s011
(CSV)
Acknowledgments
We thank the Servicio de Espacios Protegidos y Conservación de la Naturaleza, Dirección General de Biodiversidad of the Principado of Asturias for allowing us to use the brown bear database. Special thanks to the Patrulla Oso members for their continuous help in the field: David Cañedo, Juan Díaz García, Miguel Fernández Otero, Damián Ramos and Fernando Samoano.
References
- 1. Treves A, Karanth KU. Human-carnivore conflict and perspectives on carnivore management worldwide. Conserv Biol. 2003;17(6):1491–1499.
- 2. Linnell JDC, Swenson JE, Anderson R. Predators and people: conservation of large carnivores is possible at high human densities if management policy is favourable. Anim Conserv. 2001;4(4):345–349.
- 3. Chapron G, Kaczensky P, Linnell JDC, von Arx M, Huber D, Andren H, et al. Recovery of large carnivores in Europe’s modern human-dominated landscapes. Science. 2014;346(6216):1517–1519. pmid:25525247
- 4. Pyare S, Cain S, Moody D, Schwartz C, Berger J. Carnivore re-colonisation: reality, possibility and a non-equilibrium century for grizzly bears in the southern Yellowstone ecosystem. Anim Conserv. 2004;7(1):71–77.
- 5. Zedrosser A, Steyaert SMJG, Gossow H, Swenson JE. Brown bear conservation and the ghost of persecution past. Biol Conserv. 2011;144(9):2163–2170.
- 6. Breitenmoser U. Large predators in the Alps: the fall and rise of man’s competitors. Biol Conserv. 1998;83(3):279–289.
- 7. Kaczensky P, Chapron G, von Arx M, Huber D, Andrén H et al. Status, management and distribution of large carnivores—bear, lynx, wolf and wolverine—in Europe. Rep to EU Comm. 2013;272.
- 8. Stahl P, Vandel JM, Herrenschmidt V, Migot P. Predation on livestock by an expanding reintroduced lynx population: Long-term trend and spatial variability. J Appl Ecol. 2001;38(3):674–687.
- 9. Athreya V, Odden M, Linnell JDC, Krishnaswamy J, Karanth U. Big cats in our backyards: persistence of large carnivores in a human dominated landscape in India. PLoS One. 2013;8(3):2–9.
- 10. Penteriani V, Delgado MDM, Pinchera F, Naves J, Fernández-Gil A, Kojola I, et al. Human behaviour can trigger large carnivore attacks in developed countries. Sci Rep. 2016;6(1432):20552.
- 11. Treves A, Naughton-treves L, Mladenoff DJ, Wydeven AP, Treves A, Naughton-treves L, et al. Predicting human-carnivore conflict: a spatial model derived from 25 years of data on wolf predation on livestock. 2004;18:114–125.
- 12. Kaartinen S, Luoto M, Kojola I. Carnivore-livestock conflicts: determinants of wolf (Canis lupus) depredation on sheep farms in Finland. Biodivers Conserv. 2009;18(13):3503–3517.
- 13. Ericksson T.; Dalerum F. Identifying areas for an expanding wolf population. Biol Conserv. 2018;In Press.
- 14. Phillips SJ, Anderson RT, Schapire RE. Maximum entropy modeling of species geographic distributions. Ecol Modell. 2006;190:231–259.
- 15. Phillips SJ, Anderson RP, Dudík M, Schapire RE, Blair ME. Opening the black box: an open-source release of Maxent. Ecography (Cop). 2017;40(7):887–893.
- 16. Merow C, Smith MJ, Silander JA. A practical guide to MaxEnt for modeling species’ distributions: What it does, and why inputs and settings matter. Ecography (Cop). 2013;36(10):1058–1069.
- 17. Jerina K, Debeljak M, Džeroski S, Kobler A, Adamič M. Modeling the brown bear population in Slovenia: a tool in the conservation management of a threatened species. Ecol Modell. 2003;170(2–3):453–469.
- 18. Rondinini C, Stuart S, Boitani L. Habitat suitability models and the shortfall in conservation planning for African vertebrates. Conserv Biol. 2005;19(5):1488–1497.
- 19. Muntifering JR, Dickman AJ, Perlow LM, Hruska T, Ryan PG, Marker LL, et al. Managing the matrix for large carnivores: a novel approach and perspective from cheetah (Acinonyx jubatus) habitat suitability modelling. Anim Conserv. 2006;9(1):103–112.
- 20. Jedrzejewski W, Jedrzejewska B, Zawadzka B, Borowik T, Nowak S, Mysłajek RW. Habitat suitability model for Polish wolves based on long-term national census. Anim Conserv. 2008;11(5):377–390.
- 21. Huck M, Jędrzejewski W, Borowik T, Miłosz-Cielma M, Schmidt K, Jędrzejewska B, et al. Habitat suitability, corridors and dispersal barriers for large carnivores in Poland. Acta Theriol (Warsz). 2010;55(2):177–192.
- 22. Mateo-Sánchez MC, Cushman SA, Saura S. Scale dependence in habitat selection: the case of the endangered brown bear (Ursus arctos) in the Cantabrian Range (NW Spain). Int J Geogr Inf Sci. 2014;28(8):1531–1546.
- 23. Mateo-Sánchez MC, Gastón A, Ciudad C, García-Viñas JI, Cuevas J, López-Leiva C, et al. Seasonal and temporal changes in species use of the landscape: how do they impact the inferences from multi-scale habitat modelling? Landsc Ecol. 2016;31(6): 1261–1276.
- 24. Schadt S, Revilla E, Wiegand T, Knauer F, Kaczensky P, Breitenmoser U, et al. Assessing the suitability of central European landscape for the reintroduction of Eurasian lynx. J Appl Ecol. 2002;39:189–203.
- 25. Kuemmerle T, Radeloff VC, Perzanowski K, Kozlo P, Sipko T, Khoyetskyy P, et al. Predicting potential European bison habitat across its former range. 2011;21(3):830–843. pmid:21639048
- 26. Whiteman A, Passoni G, Rowcliffe JM, Ugarković D, Kusak J, Reljić S, et al. Identifying key denning habitat to conserve brown bear (Ursus arctos) in Croatia. Wildl Res. 2017;44(4):309–315.
- 27. Swanepoel LH, Lindsey P, Somers MJ, van Hoven W, Dalerum F. Extent and fragmentation of suitable leopard habitat in South Africa. Anim Conserv. 2013;16(1):41–50.
- 28. Bassi E, Willis SG, Passilongo D, Mattioli L, Apollonio M. Predicting the spatial distribution of wolf (Canis lupus) breeding areas in a mountainous region of central Italy. PLoS One. 2015;10(6):1–14.
- 29. Martin J, Revilla E, Quenette P, Naves J, Swenson JE, Ame C. Brown bear habitat suitability in the Pyrenees: transferability across sites and linking scales to make the most of scarce data. 2012;49(3):621–631.
- 30. Martin J, Basille M, Van Moorter B, Kindberg J, Allainé D, Swenson JE. Coping with human disturbance: spatial and temporal tactics of the brown bear (Ursus arctos). Can J Zool. 2010;88(9):875–883.
- 31. Mateo-Sánchez MC, Gastón A, Ciudad C, García-Viñas JI, Cuevas J, López-Leiva C, et al. Seasonal and temporal changes in species use of the landscape: how do they impact the inferences from multi-scale habitat modeling? Landsc Ecol. 2016;31(6):1261–1276.
- 32. Swenson JE, Gerstl N, Dahle B, Zedrosser A. Action plan for the conservation of the brown bear (Ursus arctos) in Europe. Nature and environment Council of Europe. 2000;69.
- 33. Pérez T, Naves J, Vázquez JF, Fernández-Gil A, Seijas J, Albornoz J, et al. Estimating the population size of the endangered Cantabrian brown bear through genetic sampling. Wildlife Biol. 2014;20(5):300–309.
- 34. Gonzalez EG, Blanco JC, Ballesteros F, Alcaraz L, Palomero G, Doadrio I. Genetic and demographic recovery of an isolated population of brown bear Ursus arctos L., 1758. PeerJ. 2016;4.
- 35. Martínez Cano I, González Taboada F, Naves J, Fernández-Gil A, Wiegand T. Decline and recovery of a large carnivore: environmental change and long- term trends in an endangered brown bear population. Proc R Soc B. 2016;283(1843):20161832. pmid:27903871
- 36. Pérez T, Naves J, Vázquez JF, Seijas J, Corao A, Albornoz J, et al. Evidence for improved connectivity between Cantabrian brown bear subpopulations. Ursus. 2010;21(1):104–108.
- 37. Palomero G, Ballesteros F, Nores C, Blanco JC, Herrero J, García-Serrano A. Trends in number and distribution of brown bear females with cubs-of-the-year in the Cantabrian Mountains, Spain. Ursus. 2007;18(2):145–157.
- 38. Penteriani V, Delgado MM, López-Bao J V., Vázquez García P, Monrós JS, Vigón Álvarez E, et al. Patterns of movement of three rescued and released female brown bears in the Cantabrian Mountains, northwest Spain. Ursus. 2018;2–7.
- 39. Naves J, Wiegand T, Revilla E, Delibes M. Endangered Species Constrained by Natural and Human Factors: the Case of Brown Bears in Northern Spain. 2003;17(5):1276–1289.
- 40. Hijmans RJ, Cameron SE, Parra JL, Jones PG, Jarvis A. Very high resolution interpolated climate surfaces for global land areas. Int J Climatol. 2005;25(15):1965–1978.
- 41. Kobler A, Adamic M. Identifing brown bear habitat by a combined GIS and machine learning method. Ecol Modell. 2000;135:291.
- 42. Posillico M, Meriggi A, Pagnin E, Lovari S, Russo L. A habitat model for brown bear conservation and land use planning in the central Apennines. Biol Conserv. 2004;118(2):141–150.
- 43.
R Core Team (2018). R: A language and environment for statistical computing. R Foundation for Statistical Computing, Vienna, Austria. URL https://www.R-project.org/.
- 44. Hijmans RJ, Phillips S, Leathwick JR, Elith J. Package “dismo.” Circles. 2017;9(1).
- 45. Muscarella R, Galante PJ, Soley-Guardia M, Boria RA, Kass JM, Uriarte M, et al. ENMeval: An R package for conducting spatially independent evaluations and estimating optimal model complexity for Maxent ecological niche models. Methods Ecol Evol. 2014;5(11):1198–1205.
- 46. Elith J, Phillips SJ, Hastie T, Dudík M, Chee YE, Yates CJ. A statistical explanation of MaxEnt for ecologists. Divers Distrib. 2011;17(1):43–57.
- 47. Phillips SJ, Dudík M. Modeling of species distribution with Maxent: new extensions and a comprehensive evalutation. Ecograpy. 2008;31(December 2007):161–175.
- 48. Warren DL, Seifert SN. Ecological niche modeling in Maxent: the importance of model complexity and the performance of model selection criteria. Ecol Appl. 2011;21(2):335–342. pmid:21563566
- 49. Akaike H. A new look at the statistical model identification. IEEE Trans Automat Contr. 1974;19(6):716–723.
- 50.
Burnham K. P., & Anderson DR. Model selection and multimodel inference. Springer Science & Business Media. 2003.
- 51. Fielding AH, Bell JF. A review of methods for the assessment of prediction errors in conservation presence / absence models. Environ Conserv. 1997;24(1):38–49.
- 52. Araújo MB, Pearson RG, Thuiller W, Erhard M. Validation of species–climate impact models under climate change. Glob Chang Biol. 2005;11:1504–1513.
- 53.
Young N., Carter L., Evangelista P. A MaxEnt model v3. 3.3 e tutorial (ArcGIS v10). Fort Collins, Colorado: Natural Resource Ecology Laboratory at Colorado State University and the National Institute of Invasive Species Science. 2011;1–30.
- 54. Almeida Neto M, Guimaraes P, Guimaraes PR Jr, Loyola RD, Ulrich W A consistent metric for nestedness analysis in ecological systems: reconciling concept and measurement. Oikos 2008; 117:1227–1239
- 55. Ulrich W, Almeida-Neto M, Gotelli NJ. A consumer’s guide to nestedness analysis. Vol. 118, Oikos. 2009:3–17.
- 56. Moran A. A test for the serial independence of residuals. Biometrika. 1950;37(1/2):178–181.
- 57.
Pinheiro JC, Bates DM (2000) "Mixed-Effects Models in S and S-PLUS", Springer. Science & Business Media. 2000.
- 58. Borcard D, Legendre P All-scale spatial analysis of ecological data by means of principal coordinates of neighbour matrices. Ecol Mod. 2002;153:51–68.
- 59. Tukey JW. Comparing individual means in the analysis of variance. Int Biometric Soc. 1949;5(2):99–114.
- 60. Lenth R. V. (2016). Least-squares means: the R package lsmeans. Journal of statistical software, 69(1), 1–33.
- 61. Mateo-Sánchez MC, Cushman SA, Saura S. Connecting endangered brown bear subpopulations in the Cantabrian Range (north-western Spain). Anim Conserv. 2014;17(5):430–440.
- 62. Naves J, Fernández-Gil A, Rodríguez C, Delibes M. Brown bear food habits at the border of its range: a long-term study. J Mammal. 2006;87(5):899–908.
- 63. Penteriani V, Zarzo-Arias A, Delgado MM, Dalerum F, Gurarie E, Peón Torre P, et al. Do brown bear females with cubs settle in risky areas in spatially confined populations? The case of the Cantabrian Mountains, Spain. Submitted.
- 64.
A.P. C. El Oso en León. Univ León, Spain. 1995.
- 65. Rodríguez C, Naves J, Fernández-Gil A, Obeso JR, Delibes M. Long-term trends in food habits of a relict brown bear population in northern Spain: The influence of climate and local factors. Environ Conserv. 2007;34(1):36–44.
- 66. Fernández-Gil A. Comportamiento y conservación de grandes carnívoros en ambientes humanizados. Osos y Lobos en la Cordillera Cantábrica. 2013;296.
- 67. Clevenger AP, Purroy FJ, Pelton M. R. Movement and activity patterns of a european brown bear in the Cantabrian Mountains, Spain. Bears Their Biol Manag. 1990;8:205–211.
- 68. Oberosler V, Groff C, Iemma A, Pedrini P, Rovero F. The influence of human disturbance on occupancy and activity patterns of mammals in the Italian Alps from systematic camera trapping. Mamm Biol. 2017;87:50–61.
- 69. Ripple WJ, Beschta RL. Trophic cascades in Yellowstone: the first 15years after wolf reintroduction. Biol Conserv. 2012;145(1):205–213.
- 70. Naves J., Ordiz A., Fernández-Gil A., Penteriani V., Delgado M.M., López-Bao J. V., Revilla E. DM. Spatiotemporal patterns and environmental correlates in apiaries depredated by brown bears in the Cantabrian Mountains, Spain. Submitted.