Top predators such as leopard Panthera spp. are often associated with high biodiversity, so the protection of their habitats is one of the most effective ways to conserve biodiversity globally. In this paper, we use ecological niche factor analysis (ENFA), a presence-only environmental habitat-envelope based method to create habitat suitability maps for Persian leopard in Golestan National Park (GNP), Iran. The Persian leopard Panthera pardus saxicolor is an endangered subspecies on the IUCN Red List of Threatened Species. During a one-year field study in 2009, we recorded 120 leopard locations and related these to 14 environmental variables. Our analysis shows that the Persian leopard in this area lives within a very narrow range of conditions and therefore may require rather specific habitat protection and management in this area.
Knowledge of the distribution of threatened species and their habitat requirements is an important element of conservation biology (Engler et al. 2004). Management for threatened species, ecosystem restoration, species reintroductions, population viability analyses and resolving conflicts between humans and wildlife often rely on habitat suitability modelling (LeLay et al. 2001, Hirzel et al. 2001). Multivariate models are commonly used to define habitat suitability and combined with geographical information systems (GIS) allow researchers to create potential distribution maps (Guisan & Zimmermann 2000).
Presence-only modelling techniques are increasingly being used to study the distribution of many different organisms (Robertson et al. 2001, Hirzel et al. 2001). One of these alternative techniques is ecological niche factor analysis (ENFA; Hirzel et al. 2002). Because the results of ENFA analysis are straightforward and often easily interpreted (Rood et al. 2010), and because the software is readily available, we chose to use ENFA for assessing the habitat of the Persian leopard Panthera pardus saxicolor. The canonical depiction of the species' niche relative to its environment allows one to evaluate which part of the available habitat is occupied and to assess to which extent the available habitat is utilised (Titeux et al. 2006, Braunisch et al. 2008). Note that ENFA is more suited to determine a species potential distribution rather than its realised distribution (Jimenez-Valverde et al. 2008).According to the latest assessment of leopard status for the 2008 IUCN Red List of Threatened Species, the Persian leopard should be classified as ‘endangered’ under the category C2a(i) (Khorozyan et al. 2005, Khorozyan 2008). Iran has been the stronghold of the Persian leopard population in the Middle East where the range of approximately 885,300 km2 provides home for an estimated 550-850 individuals (Kiabi et al. 2002, Khorozyan 2008).
The most urgent threat to the Persian leopard is the ever increasing fragmentation of the population into the patchy network of distant and often too small subpopulations. According to Ghoddousi et al. (2008), not a single subpopulation across the entire range is believed to contain more than 100 adult individuals and only a handful of protected areas (all concentrated in Iran) are large enough to maintain viable subpopulations of Persian leopards. Prey reduction because of poaching, infrastructure development, disturbance and habitat loss (collection of edible plants and mushrooms, mining, road construction, deforestation, wild fire and livestock grazing) are the principal factors of the fragmentation, leaving vast tracts of habitats unsuitable for resident Persian leopard subpopulations (Fumagalli 2007, Khorozyan 2008, Ghoddousi et al. 2008). In some areas, Persian leopards regularly attack domestic livestock because of reductions in natural prey, and then clash with rural people who try to eliminate predators by poisoning prey remains (Farhadinia et al. 2007). Previous research on the Persian leopard in Iran has focused on the species' population size and status (Kiabi et al. 2002, Farhadinia et al. 2009), reproductive ecology (Farhadinia et al. 2009), genetic diversity (Farhadinia 2009), territorial marking (Ghoddousi et al. 2008), habitat modelling (Omidi 2008), status analysis and distribution mapping in Iran (Sanei 2007).
Despite the need for efficient habitat protection efforts, more ecological studies on the Persian leopard in Golestan National Park (GNP) are needed. In fact, our study is the first attempt to determine Persian leopard habitat suitability in GNP to date; no field study has been carried out in this region. Hence, in the present study, we used ENFA to: 1) develop predictive habitat suitability maps for Persian leopards, 2) identify the environmental variables important in describing the habitat for this subspecies and 3) quantify the extent and location of potential Persian leopard habitat available for conservation action in GNP.
Material and methods
Study area
The GNP is located in the northeastern part of Iran near the border with Turkmenistan and covers an area of about 919 km2. It is considered one of Iran's most important national parks because of its natural assets such as its verdant, virgin forests. The park is located east of the Caspian Sea between 37°24‘N and 55°58‘E (37.403°N and 55.976°E; Madjnoonian et al. 1999; Fig. 1). The vegetation of the park can be divided into two zones: the Hyrcanian forest in east Alborz (the western section of the park with a high humidity) and the Iran-Turanian vegetation (the eastern section of the park where it is dry; Javanshir 1976). Its maximum and minimum elevations are the summit of Divarkaji and the Tangrah at approximately 2,411 m and 450 m altitudes, respectively. The park includes mountainous areas, hills, fields and plains. The mountainous areas of the park are mostly located in the northern and western parts, with altitude reducing gradually towards the steppe-covered east. The average annual precipitation is 400 mm and annual average temperature is 11.9°C from April to October and 10.5°C from December to March.
GNP is the most important protected area of Iran with great habitats for Persian leopards because of its unique natural situation and the well-chosen location (Agili 2005). One of the unresolved problems of the park at the moment is the Asian highway, which runs through the park and divides it into a northern and a southern part, possibly hampering leopard movement (see Fig. 1).
Leopard distribution data
During a one-year field study in 2009, we recorded 120 detailed point locations of the Persian leopard. These distribution data come from field observations and interviews with biologists and park managers from regional environmental agencies and all were verified by visits to the locations at which Persian leopards were reported. Prey species distribution data were also obtained by means of field observations using GPS and map-based interviews with rangers and staff from the environmental office in GNP. Sampling for evidence of Persian leopard presence was done using tracks, scrapes and scats. We only recorded spots carrying certain signs of the Persian leopard being present: each spot was visited, and based on the certainty of the remains or signs of this animal the records were refined. The Persian leopard presence data are shown in Figure 2.
Environmental data
ENFA needs two types of data to calculate habitat suitability: a map of locations where the species has been detected and a set of quantitative raster maps describing the environment as used by the species under investigation. Independent eco-geographical variables (EGV) quantitatively describe relevant characteristics for each grid cell. These may be topographical features (e.g. altitude and slope), ecological data (e.g. frequency of forests and nitrate concentration), or human infrastructures (e.g. distance to the nearest town and road density). A function of the EGV is then calibrated so as to classify as correctly as possible the cells as suitable or unsuitable for the species. The details of the function and its calibration depend on the analysis (Hirzel et al. 2001). Using expert knowledge and an extensive literature review, we selected 14 variables that might act as direct or surrogate determinants of the current distribution of the Persian leopard in GNP (Table 1): four accounted for environmental traits (habitat structure and geomorphology) and 10 for human impacts and main prey species. Altitude, slope and aspect were calculated from a digital elevation model of 10 m pixel width. The normalised difference vegetation index (NDVI) variable that accounted for habitat structure was derived from Landsat TM images (acquired in August 2007) originally in 30×30 mega pixels. These we then resampled to 10×10 mega pixels. Studies such as Salman Mahini (2007) and Boelman et al. (2011) indicate a positive relationship between NDVI and habitat structure as it shows density of vegetation on the ground. Variables that might potentially account for the human impacts on the Persian leopard territories were distance to road (Asian highway and dirt roads), distance to villages, distance to springs, distance to agricultural lands and distance to rivers. The major prey species of the Persian leopard in the area were goitered gazelle Gazella subgutturosa, red deer Cervus elaphus, wild boar Sus scrofa, wild goat Capra aegagrus and wild sheep Ovis orientalis. We used the presence points for these prey acquired through interview with game guards in the field.
Topographical data (i.e. altitude, slope and aspect) were quantitative and used directly in the ENFA model, but habitat structure, human impact and prey factors were qualitative and were transformed into frequency and distance variables before calculation (Hirzel et al. 2002). Distance variables expressed the distance between the focal cell and the closest cell belonging to a given category. Frequency variables described the proportion of cells from a given category within a circle of a 150 m radius around the focal cell. This radius was applied within Biomapper during our analyses. After the preparation of environmental variables and conversion of ‘Presence’ points to raster grids (10×10 mega pixels), the data were normalised through a Box-Cox transformation (Hirzel et al. 2002). The common resolution of the maps used for the analysis was set at 10 m. This resolution represented a trade-off between accuracy and computation time (Chefaoui et al. 2005).
Statistical analysis
Habitat modelling
The ENFA summarises several EGVs in a few uncorrelated factors retaining most of the information similar to principal component analysis (PCA). The ENFA concept has been incorporated into Biomapper 4 (Hirzel et al. 2007) and follows the procedures outlined by Hirzel et al. (2002). The outputs of the ENFA included eigenvalues and factor scores. The first factor, marginality, describes the distance of the species optimum from the ecological conditions in the study area (i.e. the direction in which the species niche differs most from the available conditions in the study area; Hirzel et al. 2002, Santos et al. 2006). The coefficients of the score matrix as related to the marginality factors indicate the correlation between each EGV and the marginality factor. The greater the absolute value of the coefficient the higher this EGV contributes to the marginality. A low value (close to 0) indicates that the species tend to live in average conditions throughout the study area, whereas a high value (close to 1) indicates a tendency for the species to live in extreme habitats. A positive value means that the species ‘prefers’ the high values of this EGV, while a negative value means that the species ‘prefers’ the low values. The subsequent factors are called specialisation factors and are sorted by decreasing amounts of variance accounted for. These factors describe how specialised the species is in comparison to the available range of habitats in the study area (Hirzel et al. 2002, Santos et al. 2006). Therefore, only a few of the first factors explain the major part of the whole information. Specialisation ranges from 1 to infinity and thus is difficult to interpret. For this reason, it is easier to use the tolerance factor, which measures the choosiness of the species to the available range of EGVs. Tolerance is defined as the inverse of specialisation (1/S) and ranges from 0 to 1, indicating either specialist species (stenoic) who tend to live in a very narrow range of conditions or species that inhabit any of the conditions in the study area (eurioic). With the factor scores computed, a habitat suitability map was created using the harmonic algorithm. The number of significant factors included in the habitat suitability map was decided according to the comparison of the eigenvalues to MacArthur's broken-stick distribution, which is the expected distribution when breaking a stick randomly. The eigenvalues that are larger than expected according to the broken stick distribution may be considered ‘significant’ (Hirzel et al. 2002, Hirzel et al. 2007).
Model validation and accuracy
Before using the ENFA results or habitat suitability map (HSM), we needed to evaluate their accuracy in describing the actual spatial response of the species (Santos et al. 2006). The habitat suitability map was evaluated for predictive accuracy by jackknife cross validation (Boyce et al. 2002).
The species locations were randomly partitioned into k mutually exclusive but identically sized sets. Each k minus 1 partition was used to compute a habitat suitability model and the left-out partition was used to validate it on independent data. This process was repeated k times, each time by leaving out a different partition. This process resulted in k different habitat suitability maps and the comparison of these maps and how they fluctuated, providing an assessment of their predictive power. The number of partitions used was 10. Each map was reclassified into i bins, where each bin i covered some proportion of the total study area (Ai) and contained some proportion of the validation points (Ni; validation points were the observations left out during the cross-validation process). We used the default number of bins (i.e. four bins). The area-adjusted frequency for each bin was computed as Fi = Ni/Ai. The expected Fi was 1 for all bins if the model fitted no better than random. If the model was good, low values of habitat suitability should have a low F (< 1) and high values a high F (> 1) with a monotonic increase between. The monotonic nature of the curve was measured using a Spearman rank correlation on the Fi in a moving window, termed the continuous Boyce Index (Boyce et al. 2002, Hirzel et al. 2004, Santos et al. 2006, Edgaonkar 2008). The Boyce Index measures the correlation between habitat suitability values and the area-adjusted frequency of presence points in the habitat map. The continuous Boyce Index varies from −1 for an inverse model to 0 for a random model to 1 for a perfect model (Boyce et al. 2002, Hirzel et al. 2006). Finally, the habitat suitability map was reclassified into two classes of suitable and unsuitable, and an approximation of the Persian leopard individuals was assessed in the suitable habitat class.
Results
According to MacArthur's broken-stick distribution (Hirzel et al. 2002), the 14 environmental variables considered were reduced to 10 factors (Table 2) that explained 94% of the variance. The percentages explained by the specialisation factor can be seen in Table 2. The presence of Persian leopard was positively associated with NDVI, distance to agricultural lands, aspect and presence of wild boar, wild sheep and red deer on the factor explaining marginality (see Table 2), indicating a preference for these variables. Altitude, distance to villages and frequency of goitered gazelle had the higher coefficient of the first factor showing that the distribution of the species was specially restricted by these variables (see Table 2).
According to the ecological model, Persian leopard presented a tendency to occupy very particular conditions compared to the whole of the GNP (marginality score = 0.755 and tolerance factor = 0.642). The 10 factors retained (out of the 14 computed) accounted for 94% of the total sum of eigenvalues (that is, 100% of the marginality and 94% of the specialisation). The marginality factor alone accounted for 23% of the total specialisation, a rather large value, which means that Persian leopard displays a very restricted range in which they utilise habitats quite different from the average conditions present in the study area. A suitability map was built from these 10 factors for the whole case study (Fig. 3). It is interesting to see that the northern areas with steeper slopes and higher elevations present lower quality habitat for the leopards. The habitat suitability index (HSI) ranged between 0 and 100, with 0 indicating the least suitable habitat and 100 the most suitable. The predictive accuracy of the model was good as the area-adjusted frequency cross validation exhibited values below and above 1 for the low and high suitability bins, respectively. Also, the mean Spearman rank correlation was 0.8. Using the F curve, the habitat suitability map was reclassified into two classes, ‘suitable’ and ‘unsuitable’ (Fig. 4). The HS values > 2 are considered as suitable class. These areas suitable for the Persian leopard are located in the west, central and southeast of the GNP. Also, the unsuitable class was assigned to HSI values ranging from 0 to 2. Our results reveal that the suitable habitat for the Persian leopard comprises an area of 189 km2 with an average altitude of 1,234 m and a slope of 32%. Furthermore, the marginality and specialisation scores show that the Persian leopard lives in very particular conditions in this area. Referring to previous studies about panther home-range size by Kiabi et al. (2002), if we choose 15 km2 as the lowest estimate of the animal's home range, then there would be around 13 Persian leopards in the area. For the highest estimate of the range size, i.e. 78 km2, the number of Persian leopards would be only around two. Based on the field observations, habitat assessment and the interviews with experienced wildlife guards, we suggest that the actual number of the panthers is somewhere in between two and 13.
Discussion
Here we present the first study on Persian leopard's habitat suitability in GNP. Our analysis directly provides two key measurements regarding the niche of the focal species (i.e. marginality and specialisation factors). Furthermore, interpretation of the factors in terms of the EGVs turns out to be very consistent with the published literature as particularly relevant for Persian leopard ecology (Hirzel et al. 2002).
According to the habitat suitability analyses carried out and the environmental niche descriptions, we calculated that the habitat suitable for the Persian leopard comprises an area of 189 km2, averaging 1,234 m in altitude and a slope of 32% encompassing the western, central and southeastern parts of the GNP. Because spatial autocorrelation of the presence data can cause biases in model predictions (Jimenez-Valverde et al. 2008), we checked for autocorrelation among input parameters in the model, but it was not significant.
The presence of Persian leopards had a weak negative association with the distance to villages. Leopards are known to be bold and are commonly found in the proximity of human settlements, where they prey upon livestock (Odden & Wegge 2005). These results are in agreement with previously reported data on Persian leopard habitat use (Omidi 2008, Edgaonkar 2008, Ghoddousi et al. 2008, Farhadinia et al. 2007). In our study, wild goat had the strongest positive association with Persian leopard occurrence. One of the main concerns for the conservation of Persian leopards is habitat fragmentation due to human land use (Acevedo et al. 2007). We feel that this is also an issue in the GNP. Among human interferences, the Asian highway may presumably have the largest negative impact on the Persian leopard's habitat and cause additional mortality because it runs through the park which likely restricts the movements of the Persian Leopard. Part of the prey items for Persian leopard in the GNP are found in the southern section which is divided by this highway. So, the Persian leopard is forced to cross the Asian highway to obtain food, causing collisions with the traffic. In our study area, road accidents and casualties are frequent, suggesting that the negative impact of the main road of the GNP might be significant. As suggested by marginality factor coefficients, watercourses (i.e. springs and rivers) are important habitat parameters for leopard and may restrict its presence. The occurrence of Persian leopard had a low correlation with other parameters such as aspect, NDVI, agricultural lands, wild boar, wild sheep and red deer. Overall, the habitat suitability map represents an overlap between the best habitat for Persian leopard, its preys' habitat (goitered gazelle, red deer, wild boar, wild goat and wild sheep) and leopard observations in this region during the last years.
As Hirzel et al. (2002) suggested, ENFA is a purely descriptive method and cannot extract causal relations. Nonetheless, it provides (at worst) important clues about preferential conditions, and remains a powerful tool to draw potential habitat maps. In this respect, a limitation of the software is that it does not yet provide confidence intervals for distribution maps. Increasingly, conservation managers are demanding risk analyses that incorporate uncertainties in model predictions. These could clearly be obtained through the boot strapping of presence data. Though not yet implemented in Biomapper, this procedure will certainly provide an important and useful extension. A second limitation, less easy to deal with, is that ENFA only handles linear dependencies within the species niche. Multiplicative or non-linear interactions cannot be accommodated in the present approach, except through transformations or non-linear combinations of the original ecogeographical variables. A third limitation is that some EGVs may turn out to be constant in specialisation or in linear combination with other EGVs, which makes species covariance matrix (Rs) singular. This is likely to happen with coarsely measured data or small species data sets. Whenever this happens, Biomapper identifies the constant or correlated EGVs and removes (one of) them from the analysis. An alternative approach would obviously consist in improving the field sample either by increasing the presence data set or by measuring EGVs on a finer scale. Finally, a last important point to emphasise is that our approach characterises ecological niches relative to a reference area. Marginality and specialisation are thus bound to depend on the geographic limits of the study area. Some species may turn out to occur at the very edge of their distribution, and may thus appear quite specialised in the reference set, however widespread they might be otherwise. The niche description methodology derived by Chefaoui et al. (2005), based on such premises, is used in this study to describe the realised niches of panther in GNP. However, ENFA results usually overestimate habitat suitability (Zaniewski et al. 2002, Engler et al. 2004, Acevedo et al. 2007).
In conclusion, understanding the habitat selection processes of wide-ranging large carnivores such as the Persian leopard and thus the likely consequences of habitat transformation is essential for developing conservation strategies to improve long-term viability. Our analyses demonstrate that ENFA is a useful tool to explore the characteristics of the Persian leopard's niche as well as to produce habitat suitability maps that can aid conservation management. Our results indicate that Persian leopard distribution in our study area may be especially constrained by human impacts. Our analysis suggests that the Asian highway fragments Persian leopard habitat and causes death due to road accidents.
In order to protect Persian leopards and their habitats, we suggest several conservation approaches in this region. Firstly, preventing habitat destruction and road widening are important for securing viability of the species. Secondly, habitat connectivity and wildlife corridors such as overpasses and underpasses between isolated habitats can be used to reduce the negative effects of roads and road mortality on wildlife populations and especially the Persian leopard. Thirdly, measures such as ecological compensation need to be carried out to improve the livelihood of local residents, which will decrease the influence of human activity and improve the quality of Persian leopard habitats. Finally, the viability of Persian leopard populations will be enhanced by preserving its preferred types of habitat. More research is needed to understand spatial connectivity and population dynamics of the Persian leopard.
Acknowledgements
we express our gratitude to Editor-in-Chief Ilse Storch and especially Associate Editor John P. Ball for their careful reviewing, editing and useful comments and suggestions on a previous version of our manuscript. We also thank the Iranian Department of the Environment (DoE), particularly the Golestan Provincial Office of DoE, which provided logistical and financial support for field surveys.
References
Figure 3.
Habitat suitability map for the Persian leopard in GNP, as computed from ENFA. The scale shows the habitat suitability values (0 the least suitable, and 100 the most suitable).
Figure 4.
Reclassified habitat suitability map for the Persian leopard in GNP, as computed from ENFA.
Table 1.
Variables used in the spatial modelling of Persian leopard habitat in Golestan National Park, Iran.
Table 2.
Coefficients of the variables generated in ENFA, and percentages explained by marginality (MF) and specialisation factors (SF1 - SF9). EGVs are sorted by decreasing absolute value of coefficients on the marginality factor. The first column shows 100 percent of marginality. On marginality factor, Positive values (+) indicates that leopards are found in locations with higher than average cell values. Negative values (-) indicate that leopards are found in locations with lower than average cell values. Signs of coefficients have no meaning on the specialisation factors.
