Identifying priority core habitats and corridors for effective conservation of brown bears in Iran

Iran lies at the southernmost range limit of brown bears globally. Therefore, understanding the habitat associations and patterns of population connectivity for brown bears in Iran is relevant for the species’ conservation. We applied species distribution modeling to predict habitat suitability and connectivity modeling to identify population core areas and corridors. Our results showed that forest density, topographical roughness, NDVI and human footprint were the most influential variables in predicting brown bear distribution. The most crucial core areas and corridor networks for brown bear are concentrated in the Alborz and Zagros Mountains. These two core areas were predicted to be fragmented into a total of fifteen isolated patches if dispersal of brown bear across the landscape is limited to 50,000 cost units, and aggregates into two isolated habitat patches if the species is capable of dispersing 400,000 cost units. We found low overlap between corridors, and core habitats with protected areas, suggesting that the existing protected area network may not be adequate for the conservation of brown bear in Iran. Our results suggest that effective conservation of brown bears in Iran requires protection of both core habitats and the corridors between them, especially outside Iran’s network of protected areas.

Large carnivore populations increasingly face threats from habitat loss, fragmentation and isolation worldwide [1][2][3][4][5] . Additionally, large carnivores often reside in and disperse across landscapes characterized by heterogeneous patterns of habitat suitability and mortality risk [6][7][8][9][10] . Existing protected areas often fail to support viable populations of large carnivores in many parts of the world as a result of their large area requirements, low densities and high dispersal abilities [6][7][8][9][10][11] . Thus, conservation planning for large carnivores requires assessing the efficacy of existing protected areas, prioritizing establishment of new protected areas in strategic locations, and protection of dispersing individuals as they move between these networks of protected areas [12][13][14][15][16] .
One of the approaches to ensure regional large carnivore long-term viability is based on establishing and protecting large core habitat patches, and a network of connectivity linkages, and low-risk areas among them 15,17,18 . Thus, spatially connected networks of core habitat patches often emerges as a priority to reach landscapes of coexistence between humans and large carnivores 8,10,19-22 on national and even international population level approach 23,24 .
Conservation prioritization depends on accurate assessment of the importance of areas as habitat core areas supporting local populations and also prediction of patterns of connectivity among these core population areas. Habitat suitability modeling and connectivity analyses are foundational to these two critical objectives. In recent years several studies have used multi-scale habitat modeling as the foundation for conservation prioritization 11,22,[25][26][27][28][29] , concluding that multi-scale optimization greatly improves prediction of habitat quality 30 .
Connectivity analyses are often based on resistance surfaces extracted from habitat models [25][26][27] . While recent studies have shown that habitat selection is often a poor proxy for resistance to movement and dispersal [31][32][33] , and that genetic and movement data should be used when they are available 34 . However, when such data are not extant, as in the case for brown bear in Iran, habitat quality can be used as the basis of predicting landscape resistance to movement 32,33 .
Data collection and spatial filtering. Brown bear occurrence points (n = 208) were obtained from a variety of sources including opportunistic direct observations, signs, mortality records and camera trap photos which were gathered by game rangers of Iranian Department of Environment (DoE) and researches during 2005-2020. Prior to analysis the quality, reliability and precision of all data were evaluated, and only locations with high certainty of identification and location were included. Given the challenges of acquiring data in the Iranian context, we were not able to collect more presence points of brown bears in Iran as this species has not been surveyed systematically in the country.
To reduce spatial autocorrelation in the occurrence points, which can bias model predictions, we spatially filtered points prior to analysis 26,61 . Based on the mean size of male brown bear home ranges in Turkey (83 km 2 ) 62 , circles of 5 km-radius were centered on each presence point to exclude proximal points using the Spatially Rarify Occurrence Data tool in the SDMtoolbox 63 . After this spatial filtering we retained a total of 184 presence points in the final dataset used for habitat modeling.
Environmental variables for habitat modeling. We selected topographic, climatic, land cover and human footprint variables for habitat modeling of the brown bear in Iran based on previously published habitat relationships models for the species 25 . A digital elevation model (DEM) from the 30 m Shuttle Radar Topography Mission (SRTM, downloaded from http://earth explo rer.usgs.gov), was used to calculate slope (using Surface Tool in Spatial Analyst Tools) and surface roughness variables (Geomorphometry and Gradient Metrics toolkit) 64 in ArcGIS 10.3. The national land-cover map of Iran 65 was used to create three land use variables: orchard, grassland and forest. A circular focal mean moving window with 2.5 km radius was used to create density maps for each land cover class (orchard, grassland and forest). The 16-day composite MODIS data (MOD13A1 V6 map at 500-m cell size; http://earth explo rer.usgs.gov) was used to calculate the mean normalized difference vegetation index (NDVI) values of the year 2010. Due to the importance of water resources for wildlife 53 , distance from rivers was calculated. Human footprint was represented by several variables, given that distinguishing the effects of different kinds of human perturbation is important (i.e. human population density, human infrastructure and road network) 66 . Out of 19 bioclimatic variables (http://world clim.org) 67 , we selected six variables that we believed to be most relevant to predicting the distributions of the species: annual mean temperature (BIO1), max temperature of the warmest month (BIO5), min temperature of the coldest month (BIO6), annual precipitation (BIO12), precipitation of warmest month (BIO13) and precipitation of driest month (BIO14).
Multi-scale species distribution modeling. In addition to the 184 occurrence records, we created 1000 pseudo-absence points randomly placed across the study area outside of the 5 km-radius circles 25 . We used two methods to reduce multicollinearity among variables: (1) the MaxentVariableSelection package 68 in R 69 was used to exclude variables by setting a contribution threshold of 1%, regularization multiplier of 1-5 with increments of 0.5 and inter-correlation of 0.7. Variables with the highest area under the curve (AUC) of receiver operating characteristic (ROC) and the lowest Akaike Information Criterion (AIC) were chosen; and (2) the Variance Inflation Factor (VIF) of the dataset was checked using r-package usdm 70 to exclude variables with VIF > 3 71 .
We used random forest (RF) 72 to predict habitat suitability for the brown bear, using with a multi-scale approach 73,74 . Multiple-scale modeling with random forest has been shown to outperform other approaches for ecological prediction 17,75 . Six scales (1, 2, 4, 8, 16 and 32 km) 25 were calculated for each variable using variable radius focal mean analysis. These selected scales span the range of brown bear habitat requirements, from resources within core habitats to the extent of reported home ranges of the species 48 , and correspond to the scale range previously shown to be relevant to large carnivore's 11,75 and for brown bears in particular 25 .
For each variable at each scale, we ran univariate RF and the performance of each scale of habitat suitability model was evaluated using AUC and True Statistic Skill (TSS). The scale with the highest performance was chosen for including in a final multivariate optimized model 75,76 . Variable contributions for each model were calculated, and we produced response curves of presence points to the retained predictor variables. All of these analyses were carried out using the Random forest r-package 73 . Transforming habitat suitability to landscape resistance. To estimate landscape resistance 77 , we converted the habitat suitability maps to resistance maps using the Eq. 1 78 : where R represents the cost resistance value assigned to each pixel and HS represents the predicted habitat suitability derived from the suitability models described above 32,78 . We rescaled the resistance values to a range between 1 and 10 by linear interpolation, such that minimum resistance (Rmin) was 1 when HS was 1, and maximum resistance (Rmax) was 10 when HS was 0 78 .

Connectivity modeling.
We used the universal corridor network simulator (UNICOR) 79 to create two sets of connectivity predictions including (1) resistant kernels 38 and (2) factorial least-cost paths 37 . The factorial least-cost path analysis implemented in the UNICOR simulator applies Dijkstra's algorithm to resolve the singlesource shortest path issue from every mapped species occurrence location on a landscape to every other occurrence location 79 . The analysis produces the sum of predicted least-cost paths from each source point to each destination point. The resistant kernel algorithm calculates the cumulative resistance cost-weighted dispersal kernel www.nature.com/scientificreports/ around each source point up to a user-defined dispersal threshold, providing the rate of organism movement through every pixel in the landscape as a function of the density and number of source points, the dispersal ability of the species, and the resistance of the landscape 38 . The cumulative resistant kernel surface reflects the spatial incidence function of the expected rate of movement of each species through each pixel in the landscape 80,81 .
To account for uncertainties regarding brown bear dispersal abilities 22,40 , we used five distance thresholds in the resistant kernel analyses: 50,000, 100,000, 200,000, 300,000 and 400,000 cost units, which represent movement abilities of 50, 100, 200, 300 and 400 km, respectively, through optimum low resistance habitat. The longest dispersal distances recorded for brown bear are 90 km for a female and 467 km for a male in Norway 82 , therefore, this range of modeled dispersal distances brackets the expected dispersal capability of the species.
We calculated the factorial least-cost path network without a dispersal threshold 77 to provide a broad-scale assessment of the regional pattern of potential linkage and to map potential long-distance corridors. The buffered least-cost paths were then combined through summation 37 to produce maps of connectivity among all pairs of presence points.
The resistant kernel connectivity maps were used to identify brown bear core areas 9,15 . We defined core habitat patches as contiguous patches with resistant kernel values > 25% of the highest recorded for the species 29,40 . To evaluate the effectiveness of the current network of protected areas in providing connectivity for this species in Iran, we quantified the extent and percentage of predicted core areas and corridors for this species within the current network.
Conservation prioritization of core habitats. We prioritized core habitat patches based on probability of connectivity (dPC) 83 and integral index of connectivity (dIIC) 84 for all identified core habitats across the five dispersal distance scenarios (i.e. 50, 100, 200, 300 and 400 km) in Conefor 2.6 85 . The dPC and dIIC indices are frequently used as connectivity measures in conservation prioritization studies 42,43,86,87 . The dIIC index considers both habitat amount and habitat reachability across the habitat network, and linkages as dispersal events between patches 84,86 . dIIC also quantifies the loss of connectivity if a patch is removed from the habitat network and can be decomposed into dIICflux (dIICf), dIICconnector (dIICc), dIICintra (dIICi). dPC also considers both habitat amount and habitat reachability, and the probability of dispersal between patches, and can be decomposed into dPCflux (dPCf), dPCconnector (dPCc) and dPCintra (dPCi) 84,86 . dPCintra and dIICintra measure intrapatch connectivity, while flux fraction of a particular node (dPCflux and dIICflux) reflects both patch attributes (e.g., area of suitable habitats) and its position within the landscape, while connector fraction (dIICconnector and dPCconnector) depend only on the topological position of a patch in the landscape 88 . Connector fraction quantifies the importance of the node as a stepping-stone for dispersal, i.e. facilitating dispersal between distant nodes 86 . We used a distance-probability value of 0.5 and 0.05 for minimum and maximum dispersal distances, respectively, as recommended by Saura and Torne 84 .
To quantify differences in the predicted extent and configuration of habitat, we calculated a suite of fragmentation metrics on predicted core habitat patches using FRAGSTATS software 30,45 . To conduct the FRAGSTATS analysis, we first converted the UNICOR resistant kernel outputs into patches by applying a cutoff value 31 . For each species, any values above 25th percentile of the highest dispersal scenario were reclassified as 1, representing habitat patches of high connectivity. Everything else was reclassified as 0. Then, we calculated four class level metrics using FRAGSTATS v4.2.1 45 including: (1) percentage of the landscape (PLAND), which quantifies the habitat patches of high connectivity as a percentage of the study area; (2) area-weighted mean radius of gyration (GYRATE_AM) or correlation length, which provides a measurement of the extensiveness of habitat patches of high connectivity; (3) largest patch index (LPI), which represents the percentage of the landscape comprised by the largest habitat patch of high connectivity; and (4) number of isolated patches (NP), which provides a measure of the degree of fragmentation. These metrics have been used frequently in past connectivity research 31,41,44,88 , and were shown through simulation to be strong predictors of functional landscape connectivity and gene flow 80 .
Based on the 16 km scale model, forest density, roughness, NDVI and human footprint were the most influential variables predicting potential habitat for bears in Iran. Excluding human footprint, all these variables emerged also as the most influential predictors in the other spatial scales considered (Table S3). Roughness, forest density, and NDVI (both natural and artificial vegetation; i.e., forests and orchards) showed a positive relationship with the probability of bear occurrence increased (Fig. S1). Plain areas in the central portion of Iran, www.nature.com/scientificreports/ along with the northern and southern parts of the country, had the highest resistance, and extensive areas of Zagros and Alborz had the lowest resistance for the brown bear at the scale 16 km (Fig. 2). Other scales showed similar patterns (Fig. S2).
Brown bear core habitat and connectivity network. Our connectivity simulation modeling for the brown bear revealed that most core habitat and connectivity areas are concentrated in the northern and northeastern (Alborz Mountains) and northwestern to southwestern parts of the study area (Zagros Mountains) (Fig. 3). Overall, we identified 15 core areas (107,233.67 km 2 ) at a dispersal distance of 50,000 cost units (corresponding to limit of female dispersal ability), of which two were particularly larger than 30,000 km 2 (Core 1 and 2) with a total area of 75,911.57 km 2 located in Zagros and Alborz Mountains. Thirty-one percent (31%) of Core 2 is covered by protected areas. In contrast, only 13.8% of Core 1 is covered by protected areas. Depending on assumed brown bear dispersal distances, between 17.77 and 22.10% of predicted bear core areas overlap with protected areas (Table 2). Also, the density of roads inside predicted core habitats for the same dispersal distance of 50 km is 239.63 m/km 2 (Table 2) (Table S4).
The strongest predicted connectivity for brown bear was between the northern and western bear population core areas, specifically between core areas 1, 2, 3 and 7 (Fig. 3). Based on our calculation, approximately 21% of the entire landscape was favorable to brown bear movements, but with varying degrees of strength. For this species, the lowest degree of connectivity was estimated for northeastern, northwestern and southeastern parts of the brown bear Iran distribution (Fig. 4). Based on the density of least-cost paths, the southern and northern nuclei were predicted to be the most isolated. The 29.43% of this corridor network falls within protected areas, but most predicted corridor paths are bisected multiple times by roads ( Fig. 4 and Table 3).
Identification of top-ranked brown bear core habitats. The least important predicted habitat core changed with simulated dispersal ability. For example, at a dispersal distance of 50 km, Cores 13, 14 and 15 had the smallest contribution to overall habitat connectivity, while for distances of 100 km Cores 14, 15 and 12 were replaced with Cores 9, 6 and 8. This result revealed that with increasing dispersal distance, there were important changes in the relative importance of different patches and an overall increase in connectivity importance across patches. Based on dPCc, Core 3 was the most important stepping stone among other patches at dispersal ability 50-300 km. At dispersal ability of 50 and 100 km, Core 1 and 5 had the next most important contributions as stepping stones. Core 2 had large contribution as a stepping stone at dispersal ability 100 km only. Over dispersal distances of 300 and 400 km, there was a decreasing trend in importance of core patches as stepping stones in the connectivity network of the brown bear. Based on dIICc, Core 1 and 2 were the most important patches at dispersal distance 50-200 km.
The contribution of core habitats to landscape connectivity revealed a different pattern of ranking according to the dPC index and dispersal distance scenarios (Fig. 5). From the 50 km dispersal distance to 400 km,  5 and 6). Based on dPCf, from 50 to 300 km core 3 was more important than core 2, but this pattern was not observed in dIICf (Figs. 5 and 6). From 50 to 300 km, values of dPCc for Core 3 were higher than Core 1 and 2, but this trend was different for dIICc. The number of isolated patches decreases with higher dispersal ability, whereas LPI, PLAND and CL rise with higher dispersal ability (Table 4). For this species, 9-13.5% of the landscape is occupied by connected habitat patches depending on dispersal ability.

Discussion
Our results highlight that crucial core habitats for brown bear conservation in Iran are concentrated in the Zagros (core 1) and Alborz (core 2) Mountains with important corridors between them. These two chains of mountains have a high level of biological diversity and endemism 58,89 which makes them two of the most valuable landscapes across the country from a conservation perspective. Local-scale studies revealed low connectivity among patches of brown bear habitat 90 . However, Ashrafzadeh et al. 91 applied circuit theory 36 to model connectivity across the  www.nature.com/scientificreports/ brown bear distribution in Iran and inferred strong connectivity between brown bear core habitats despite clear restrictions such as anthropogenic activity. This is because the circuit theory approach did not incorporate scale dependency of dispersal ability, which has been shown to dominate predictions of connectivity 17,22,80,92 . One of the main strengths of the resistant kernel approach as an alternative connectivity modeling method is its explicit and realistic incorporation of dispersal thresholds 93 . The dispersal ability of a species can affect core habitats and connectivity predictions significantly. Our findings are in accordance with previous studies 22,29,[42][43][44] . Indeed, several past research studies have shown that often dispersal ability is much more influential on connectivity predictions than relative landscape resistance or connectivity algorithm 22,80,92 . Therefore, based on our results, we believe that previous predictions likely overestimate connectivity of brown bear in Iran. In particular, our modeling shows that dispersal limitations that likely pertain to female brown bears result in highly disjointed and fragmented populations. Critically, dispersal of females is required for a population to colonize and reestablish a breeding population in new parts of its former range.
Another important difference between our results and those of Ashrafzadeh et al. 91 is that they reported substantial landscape resistance amongst Alborz and Zagros, and Alborz and Arasbaran 91 , which is inconsistent with the results of our analyses. Our analysis suggests that at realistic dispersal abilities these two main core areas are internally well connected.
These two core habitat parts for the brown bear population (Zagros and Alborz) were predicted to be broken up into a total of fifteen isolated patches if dispersal of brown bear is limited to 50 km (approximately the limit of female dispersal), but with a dispersal ability of 200-400 km (corresponding to the limit of male dispersal) it would result in three to two habitat patches. Our findings implied high sensitivity of the extent and fragmentation of connected habitat as a function of dispersal ability, and also to highly divergent connectivity for male vs female bears. Therefore, considerable effort should be invested in improving understanding of dispersal behavior and functional connectivity of large carnivores such as brown bear to increase the precision of core area delineation and prioritization of corridor paths in Iran.
Our factorial least-cost path analysis identified optimal routes between these areas in order to facilitate connectivity. This spatially-explicit information could guide conservation practitioners to implement landscape  Table 3. The extent and percent of corridors covered by current conservation networks for brown bear in Iran. www.nature.com/scientificreports/ conservation strategies. Similar work using these methods in southern Africa 9,22 and Southeast Asia 15,16,18 applied scenario optimization strategies to evaluate the relative impacts of alternative conservation and development scenarios. Our results provide a strong base for future scenario modeling of this kind in Iran. The information provided here, however, in itself is useful in identifying and prioritizing the most important core areas and corridors for brown bear 9 and identifying where they are most threatened by land use and roads 40 . Based on both AUC and TSS indices, the best habitat suitability model was the one at the 16 km scale, suggesting, as previously seen for this species in Spain 25 , that brown bear habitat selection is dominated by broad-scales    25,90,92,[96][97][98][99] .

Extent of corridors (km 2 ) Extent of protected corridors (km 2 ) % of protected corridors
Proper design of protected areas and protected area networks should incorporate spatially explicit prioritization of both core areas and connectivity networks among them, which will ensure that protected areas functionally protect focal populations and enable dispersal among them, which facilitates demographic rescue and gene flow which are crucial for long-term species conservation 21,100 . In our case, the intersection of existing Iranian protected areas with our identified core habitats was relatively low. Only 29.43% of the predicted core area network of brown bear overlapped with protected areas. Therefore, our results can guide to conservation practitioners to establish new protected areas in strategic locations most important within the full conservation network 101 . Protecting a habitat large enough to sustain the life history traits of brown bear as an umbrella species will likely also protect the long-term persistence of many other species 99,102 .
Among the identified core habitats, the highest overlap between core habitats and protected areas was observed for core 1 and core 2. These two core habitats have also been documented to have a high potential for supporting other large carnivores of conservation concern, such as Persian leopards (Panthera pardus saxicolor) 86 . However, the coverage of protected areas could be improved for other identified core habitats which our analysis shows are relatively unprotected. Similarly, Moqanaki and Cushman 41 and Khosravi et al. 42 , predicted that the distribution of protected areas in Iran was an important factor determining the occurrence and dispersal of Asiatic cheetah and sympatric carnivores.
The strategic expansion of protected areas to provide stepping stones and augmentation to the key core areas should be accompanied by the establishment of conservation actions in the linkages between protected core areas to increase functional landscape connectivity for carnivores 9 . The most robust functional corridors for the brown bear were predicted between core 1, 2, 3 and 7 (Fig. 5). We recommend prioritization of establishing new protected areas along these key connectivity routes, land use zonation and management actions along these routes to reduce mortality risk and human-wildlife conflict.
Most protected area networks in developing countries such as Iran are fragmented by roads, and road collisions are a serious threat for carnivores 41,42 . In addition to reducing dispersal success, roads also can cause high rates of direct mortality (wildlife vehicle collisions) which can be serious threat for brown bears 103 . This problem will be aggravated when roads bisect corridors [104][105][106] . Our results identified vulnerable parts of the connectivity network in Core 1 and 2 where roads intersected strong movement corridors 34,40 . Vulnerability of these locations is related to the potential for brown bear vehicle collisions, and concentration of human activities and access in these areas increasing risk of poaching, a common threat for bears in human-modified landscapes 98,100 . Our findings in this regard are similar to those of Moqanaki and Cushman 41 and Khosravi et al. 42 , who both identified multiple instances where primary and secondary roads cross the predicted corridor paths between these core patches. Our results should be combined with those of Moqanaki and Cushman 41 and Khosravi et al. 42 to prioritize road segments for mitigating mortality risk and improving connectivity, such as through fencing and overpass structures to funnel dispersing animals across the road in optimal locations for network connectivity.
Given uncertainty in the functional dispersal distances of existing brown bear populations, plus the high sensitivity of populations to this parameter, it is important to evaluate a range of dispersal distance scenarios, identify dispersal-scale thresholds, and develop multi-scale conservation recommendations. That is one of the strengths of the approach we presented in this paper. By evaluating connectivity and prioritizing core areas and corridors across a range of dispersal distances, we identify the key nodes and linkages across the brown bear population in Iran and quantify its sensitivity to a scale of dispersal behavior. The dispersal distances we modeled bracket the expected functional responses of brown bears. Female brown bears are known to be more risk averse, more philopatric and less mobile than males. The lower end of our dispersal distance simulation suggests that female brown bears are unlikely to disperse distances necessary to recolonize most extant habitat patches. Therefore, for population recovery and range recolonization relocation of female brown bears to areas where male brown bears are known to disperse to may be an important conservation strategy.
When designating vast protected areas is politically intractable 107 , an alternative may be to develop networks of interconnected protected areas. Managing and maintaining functional corridors may be more feasible in some contexts than establishing new protected areas 108 . These connected protected area networks could aid gene flow 21 and prevent isolation of small populations 92 . Due to tendency of bears to avoid humans, connecting habitats for Table 4. FRAGSTATS results for the brown bear. The metrics include: number of individual core patches (NP) largest patch index (LPI), percentage of landscape in connected habitat (PLAND) and correlation length of core habitats (CL) for the brown bear in five levels of dispersal ability (50,000, 100,000, 200,000, 300,000, 400,000). The core habitats were defined as contiguous units with resistant kernel values > 25% of the highest resistance kernel for the species.  27 . However, we also urge caution in the tempting idea that maintaining connectivity can mitigate for habitat loss. In most cases, it cannot. Therefore, carnivore conservation should focus on maintaining core habitat quality and extent as the primary focus and then establishing corridor networks linking protected core areas, with conflict management and efforts to foster human-wildlife coexistence across the landscape. It should also be noted that findings of connectivity studies based on presence points and habitat models are limited by uncertainties 32,33 . Specifically, our analysis uses a relatively modest sample of presence-only data, which could limit the power of our predictions. However, the high performance of our models using bootstrapped and cross-validated model assessment suggests that our habitat predictions are robust. Habitat quality, however, is not the same as connectivity 32,33,75 . It would be better to fit connectivity models with movement 76 or gene flow 108 data instead of habitat models, as dispersal is often related to different factors at different scales than home range habitat selection. Therefore, satellite tracking 76,88 and landscape genetic 31,108 studies are also necessary to make more reliable predictions and validate findings carried out from connectivity prioritizations made on the basis of habitat selection. To conclude, habitat suitability and connectivity models should be considered the first step towards building a nation-wide strategy for corridor improvement 101 . Optimizing habitat protection and connectivity 15,18 should be a core component of efforts to plan the establishment of new protected areas, particularly for large, vulnerable and highly mobile species such as carnivores. Furthermore, to facilitate humanbrown bear coexistence outside protected areas DoE should consider implementing some approaches such as reduction of human-induced mortality, safeguard habitat connectivity; mitigate road effects 109