On the elephant trails: habitat suitability and connectivity for Asian elephants in eastern Indian landscape

Study area

Odisha state is situated in the eastern part of India (81°24′–87° 29′E and 17°48′–22°34′N), covering an area of 155,707 km² (Fig. 1). Affected by the Southwest and Northeast monsoons, the state bears an annual mean temperature of 15 °C–38 °C and annual precipitation of 960–1,870 mm. During the monsoons, the humidity level is high (c. 80%) for four months (July to October). Summer (March to June) and winter (November to February) are the other two seasons. The state has 51,345 km² of forest cover, represents ∼33% of its geographical area (India, 2017). The forest types of Odisha are mainly categorized into tropical moist deciduous, tropical dry deciduous, and a few areas of semi-evergreen forests (Champion & Seth, 1968). Odisha can be divided into three major physiographical regions: the middle mountainous region, the western rolling uplands, and the eastern coastal plains (Sinha, 1971). The mountainous region constitutes approximately three-quarters of the state. Generally, the topography is undulating, except for the high plateaus and coastal plain areas. The mountainous region of the state is one of the most mineral-rich areas of the country. Around 60% of the state’s population is dependent on agriculture, and small-scale cultivation is the predominant livelihood strategy of the people (Panda & Padhi, 2021). Rice is the primary crop in the state, grown during the monsoon season (July to October) and harvested between October and December (Hoda et al., 2021). The state has a population density of 269 per km² and 41.9 million people; primarily (80%) live in rural areas. Also, the state has the third-largest tribal population in the country, accounting for 22.8% of the total population of the state (Dungdung & Pattanaik, 2020).

The state has 19 Wildlife Sanctuaries, one National Park, two tiger reserves, and one biosphere reserve, constituting 4.73% of the state’s geographical area. Additionally, the state has also declared three elephant reserves, namely, Mayurbhanj Elephant Reserve, Mahanadi Elephant Reserve, and Sambalpur Elephant Reserve. The Odisha state has 30 administrative districts and 50 forest divisions.

Asian elephant presence data, spatial filtering and spatial shift of conflict events

We collected presence data (n = 6,352) of the Asian elephants through direct observations and newspaper reports between January 2011 and January 2022. The direct observations included direct sightings, indirect signs of the species (i.e., footprint, dung and feeding signs) through opportunistic surveys of Asian elephants or other large mammals, and human-elephant conflict incidences (attack on humans and livestock, raid crops and damaging houses or other infrastructures) in different protected areas (PA) and forest divisions. The opportunistic direct observations of Asian elephants were obtained from various sources, including incidental observation, transect walk and forest road survey. Data were collected without following a systematic approach across state forest divisions in protected and unprotected areas of forest patches. During transect or forest road surveys, two observers who were familiar with Asian elephant signs have recorded the location. Our survey yielded 2,432 direct observation and 195 conflict locations covering nine protected areas and 22 forest divisions across the state.

Our opportunistic survey may be prone to geographic bias as it may not have sampled the full range of environmental conditions in which the species occur (Phillips et al., 2009; Yackulic et al., 2013; Guillera-Arroita et al., 2015). Therefore, to understand the presence of Asian elephants across the state, we also collected information on Asian elephants daily from regional Odia newspapers from 2017 to 2021. We selected two newspapers due to their wide reach throughout the state and availability of online publishing of local news. In addition to elephant presence information, we also collected geographical information such as the name of the reported village, Gram Panchayat, Tehsil, and district from the newspapers. We projected geospatial coordinates of elephant incidents referring the centroid of the reported village as an approximate location of the elephant presence. We derived the digitized village boundary map from the Survey of India. The average area of villages was 2.92 km², and the mean distance between centroid and village boundary was 0.44 ± 0.25 SD km, which indicates a maximum error of <1 km in the projected presence locations. Conflict records involving house damage, crop raid or human injury and death by Asian elephants were also collected at the same time from newspaper and location were manually filtered to avoid duplicate records. Given the absence of other sympatric large-sized species, its presence and signs are easily detected and identified; consequently, the presence data is highly accurate. We yielded 3,725 presence records of elephants from newspapers covering 43 forest divisions and assumed that all potential elephant habitats were adequately covered. We pooled the elephant locations from direct observations and newspapers from 2011 to 2022 to create a comprehensive dataset.

Several studies indicated that human-elephant conflict incidents are often associated with the proximity of human settlements to forested areas (Sukumar, 1993; Gubbi, 2012; Neupane et al., 2019). Therefore, we also assumed that any conflict-caused elephant in or around the village area was moved in from the nearest forest patch. Subsequently, we spatially simulated and shifted 2,613 conflict events into the nearest forest patch using the “create random points” tool in ArcMap, and the shift distance opted randomly between 5–10 km selected based on home range and daily movement of Asian elephants (Fernando et al., 2008; Baskaran, Kanakasabai & Desai, 2018; Torre et al., 2019). The spatial shift of presence events was intended to model the habitat connectivity for elephants in Odisha rather than spatial trends of conflict. However, our opportunistic surveys and newspaper reports may bias toward easily accessible sub-optimal habitats that have been more extensively reported than other habitats. Therefore, to reduce the sampling bias on model performance and to avoid autocorrelation, we spatially filtered presence locations on the 1 km spatial resolution using the “Spatially Rarefy Occurrence Data for SDMs” tool under SDMtoolbox version 2.4 on ArcMap version 10.8. This procedure resulted in 5,568 spatially independent species presence records from 45 forest divisions (114.06 ± 13.8 SE records/division, 32–492) and 13 protected areas (36.93 ± 8.04 SE records/PA, 2–113) across the state. These presence-only locations were further considered for spatial analysis.

Environmental variables selection

An extensive array of non-correlated environmental variables was selected (r<—0.6—) (Riley, De Gloria & Elliot, 1999; Rood, Ganie & Nijman, 2010; Thapa, Kelly & Pradhan, 2019; Sharma et al., 2020; Wilson et al., 2021; de la Torre et al., 2021). Variables included average annual precipitation (bio12), distance from nearest protected areas, elevation, terrain ruggedness index, and normalized difference vegetation index (Fig. S1; Table S1). Annual precipitation (mean of annual precipitation for the year 2000-2017) was used to represent the vegetation growth and to distinguish coastal areas. This information was acquired from the WorldClim database version 2.0 (http://www.worldclim.org; Fick & Hijmans, 2017). Protected area boundary (includes tiger reserves, National Parks, and wildlife sanctuaries) was obtained from the GeoNode online portal (https://geo.ejatlas.org) and cross-validated with PA map of state forest department. Elephants utilize many of these protected areas throughout the year as a refuge or halt sites while migrating. Therefore, we used the “Euclidean Distance” tool in ArcMap version 10.8 to calculate the shortest distance from the protected area. Elevation data (SRTM 250 m) was downloaded from USGS earth explorer (http://earthexplorer.usgs.gov). Terrain ruggedness index was computed using “Terrain Ruggedness Index (TRI)” tool in QGIS version 3.22.1 Białowieża from elevation data (Riley, De Gloria & Elliot, 1999). Normalized difference vegetation index (NDVI) was downloaded from the Terra-NDVI AppEEARS NASA earth data portal (http://appeears.earthdatacloud.nasa.gov) from 2011 to January 2022. NDVI data was acquired from the middle of January, April, July and October of every year from 2011 to 2021 to represent local seasons. Mean NDVI was pooled from all 45 NDVI layers using the “Raster Calculator” tool in ArcMap version 10.8, further considered one of the predictor variables. All the predictor variables were resampled to 1 km spatial resolution.

To evaluate multicollinearity among variables, we checked the pairwise Pearson correlation coefficient (r) between covariates and considered r =—0.6— as the threshold for variable selection (Khosravi, Hemami & Cushman, 2018). Furthermore, we tested the variance inflation factor (VIF) using package “usdm” version 1.1 in the R studio to sidestep the raster collinearity. We included variables for which the VIF<10, indicating least collinearity between variables (Guisan, Thuiller & Zimmermann, 2017).

Habitat suitability modelling

The habitat suitability model for the Asian elephant in Odisha was generated by ensemble modelling approach using biomod2 version 3.5.1 package in R studio version 1.4.1717 (Thuiller et al., 2021). The precise and more accurate ensemble model was produced with combining generalized linear model (GLM), generalized additive model (GAM), gradient boosted machine (GBM), maximum entropy (MaxEnt), multiple adaptive regression splines (MARS) and random forest (RF) algorithms (Araújo & New, 2007; Marmion et al., 2009). We also estimated the contribution of predictor variables to the ensemble model. Our response variable, i.e., binomial matrix of presence and pseudo-absence (one set of random selection) of the species, was modelled against the predictor variables (Phillips et al., 2009). Considering the optimum ratio, the number of pseudo-absences was kept double of independent presence records (n = 11, 136). We examined the predictive ability of mentioned algorithms using 10-fold cross-validation; 80% of response data was selected randomly for training and 20% for testing the prediction (Guisan, Thuiller & Zimmermann, 2017; Thuiller et al., 2019).

Model performance was tested using receiver operating characteristics (ROC; area under curve, AUC), the ratio of true positives among positive predictions against false positives among negative predictions. AUC = 0.5 indicates that the model’s performance is not better than random, while AUC closer to 1.0 indicates a better performance. Also, Cohen’s kappa (Heidke skill score; KAPPA) and true skill statistic (TSS) were measured for evaluating relative accuracy. The TSS values vary from −1 to +1, where +1 signifies perfect performance, and 0 or less than suggests random prediction. Further, we computed the weighted pool average of models according to their predictive accuracy as independently assessed (Thuiller et al., 2009). Model execution was dominantly executed using binomial logit and Bernoulli distributions. Predicted spatial probabilities were pooled, and the average of the better-fitted models was plotted with predictors for the response assessment. Models, which have >0.79 AUC (around ∼50% top fitted models), were selected for the ensemble (Jarnevich et al., 2018). The ensemble output depicts a relative gradient of habitat suitability for elephants across the study area and, hence, the probability of habitat suitability was rescaled to 0–100.

The habitat suitability was further classified into five equal intervals, i.e., <20% for “least potential,” 21–40% for “moderately low potential,” 41–60% for “moderate potential,” 61–80% for “moderately high potential” and >80% for “highly potential habitat”. Division-wise habitat suitability for Asian elephants was classified, as forest divisions are used as the protection and management unit by the forest department.

Identifying core habitat patches and connectivity modelling

To understand the landscape connectivity of Asian elephants, we selected 19 core habitat nodes, comprising 11 protected areas and eight non-protected areas (Table S2). We identified core habitat nodes, where elephants’ presence has been constantly recorded for last ten years after reviewing previous literature, evaluating empirical observations and consulting expert knowledge. Core habitat nodes were defined as forest patches larger than 100 km², annual home range of elephants (i.e., ranging from 100 to 1,000 km²). Protected area boundaries were obtained from GeoNode online portal, while unprotected core habitat patches were delineated from the Survey of India Topo-sheets (http://onlinemaps.surveyofindia.gov.in). The minimum and conservative distance between nodes were kept at 10 km based on the daily movement of elephant, i.e., 3.8−4.1 km (Chan et al., 2022).

Elephants have high cognitive abilities, retaining valuable information about past locations and strategically navigating to meet their needs (Sukumar, 1993; Sukumar, 2006; Vasudev et al., 2021). Therefore, we opted for a combination of least-cost path and circuit theory tools within the Linkage Mapper program (McRae & Kavanagh, 2011) in ArcMap version 10.8 to identify connectivity between core habitats. As the analysis required a resistance layer, a permeability map was developed for Asian elephants to generate the connectivity models. The inverse of ensemble habitat suitability (/conductance) was generated as a resistance map. This approach assumes that habitat suitability proportionally relates to the animal movement (Mateo-Sánchez et al., 2015). We first defined the least-cost routes connecting all the core habitat patches pairwise to identify the possible link of spatial connectivity. The least-cost model assumes that animals disperse according to previous knowledge of the surroundings. The model estimates the shortest and most cost-efficient distance between the core habitat patches, accounting for resistance in the movement (Adriaensen et al., 2003). Based on the core habitat nodes and resistance surface, the linkage model was run in order to map the least-cost linkages between adjoining core areas. The least-cost path is the route of maximum efficiency between two areas as a function of the distance travelled (length) and the cost traversed (accumulated cost) (Etherington & Holland, 2013). The connectivity linkages were evaluated under the threshold buffer of 50 km from modelled movement pathways, as the continuation of largely unsuitable habitats restricts the elephant movements. We computed the Euclidean distances (EucD) and lengths of Least-cost pathways between nodes. Furthermore, we used the ratio of euclidean distance to the length of the least-cost path (EucD: LCP) to compare the linkage quality and tortuosity. The higher ratio of the former metric indicates difficulties in moving between the core habitat pairs relative to how close they are. Later metrics describe the average resistance encountered for the optimal linking path.

Further, we used Circuitscape through Linkage Mapper to identify important areas for habitat connectivity. Unlike least-cost path models, multiple dispersals can be assessed using the circuit-based approach, which is theoretically based on random walk theory (McRae et al., 2008). As elephants are reported to utilize a wide variety of land use cover, we targeted to model all possible connectivity linkages among the core habitat nodes. For each pair of core habitats, circuitscape calculates multiple connecting paths from an electrical resistance surface. Movement often occurs in areas with the least resistance or higher conductivity. Effective resistance decreases as the number of pathways between core habitats increases. We also calculated current flow centrality by “Centrality Mapper” in the programme Linkage mapper to identify the importance of individual core habitats and linkages to maintain the connectivity within the landscape (Carroll, McRae & Brookes, 2012). Centrality scores of nodes and pathways were classified using the equal break classification, as least, moderate and highly centralized (high centrality is a high priority and so on). If the number of pathways reduces, effective resistance around pathways increases, and movement restricts to pinch points in the landscape. Higher current density within the pinch points indicates the linkages which animals are more likely to use, and no alternative links are available in the surrounding landscape. We mapped pinch points by “Pinch point mapper” in Linkage mapper to find the critical bottlenecks in connectivity; ultimately, these areas require greater conservation attention for the potential connectivity on the larger landscape. Finally, the map of movement corridors was reclassified into three categories, least, moderate and high movement-prone areas, to represent the most essential areas for Asian elephant movements.

Habitat suitability modelling

Of 60 replicated models, 30 were chosen as better-fit models (AUC>0.79). Ten models under RF algorithm, eight models under GBM, three under MARS, three under GAM, and six under MaxEnt algorithm were selected for further spatial ensembling. TSS score of 0.42−0.53 suggested good accuracy of all models (Allouche, Tsoar & Kadmon, 2006). Ensemble model displayed pooled AUC 0.94, KAPPA 0.73 and TSS 0.71 (Table S3). Model-wise AUC, KAPPA, TSS and their weightage for ensembling are given (Table S3). In the ensemble model, GLM had the lowest performance (of 76%) in predicting habitat suitability for Asian elephants in the study area, while the rest algorithms had higher contribution (79–81%).

For the ensemble model, NDVI was the most significant predictor to characterize the habitat suitability or use probability of the Asian elephant (variable importance/contribution 42%; Table S4, Fig. 2). In increase of NDVI, the probability of habitat use by Asian elephant increased. As the terrain ruggedness increased from plain to undulating terrain, the likelihood of habitat use by Asian elephants also increased (19%). The optimal range of elevation for the habitat use of Asian elephants was c. 50–500 m (17%). As the distance to protected areas increased, the likelihood of habitat use of the Asian elephant decreased (13%). The Annual Precipitation (bio12) had little positive contribution to the ensemble model (8%).

The ensemble model predicted that 14.6% (22,442.73 km²) of the state was estimated as a highly potential for Asian elephant use, followed by 13.3% (20,464.37 km²) moderately high potential, 11.34% (17,436.44 km²) moderately potential, 15.14% (23,273.95 km²) moderately low potential areas and the remainder as the least potential (Fig. 3). Out of 50 administrative forest divisions, 46 of which recorded the Asian elephant’s high presence probability (p_suitability >0.6). Focusing on high and moderately high potential areas (Thapa et al., 2018) (Table S5), 40% of forest divisions had more than 500 km² of high predicted probability area of habitat use for Asian elephants (p_suitability >0.6) (Table S5). Of the total area of the PAs, 75.34% (5,481.84 km²) were predicted to be potential habitats for Asian elephants (where p_suitability >0.6; Table S6).

Core habitat and connectivity

We considered 19 core habitat nodes (average area 644.88 ±143.41 (SE) km², range 115.47 to 2,750 km²) and identified 58 potential corridors for the Asian elephant to disperse across the landscape (Fig. 4; Table 1 & Table 2). The average length of corridors (least-cost pathway) was 126 ± 73.38 (SD) km (13.55 to 310.84 km). The mean euclidean distance between core habitat patches was 95.74 ± 59.7 (SD) km (10 to 271.95 km).

Node centrality analyses gave insights into the highly centralized nodes, which have the most significant connections with other nodes. Centrality values for all habitat cores varied considerably, with a mean value of 40.06 ± 13.63 (SD) (18 to 77.07) for the entire network. The highest centralized node (>57.38 centrality score) was only Satkosia (10). The moderate centralized nodes were (>37.69 and <57.38 centrality score) Similipal (1), Jada (2), Khalasuni (7), Anantpur (8), Kapilash (9), Subarnagir (11), Sunabeda (12), Karlapat (14), Kothagarh (15) and Chandaka (18), and the remaining eight found showing least centrality (<37.69 centrality score) (Table 1; Fig. 4).

The mean path centrality score was 10.18 ± 5.2 (SD) (3.49 to 26.57) for the entire network. The path centrality analysis suggested that four paths are highly centralized (>15.39 centrality score), 25 paths are moderately (>7.69 and <15.39) and rest 29 paths are found to be least centralized (<7.69 centrality score) (Table 2; Fig. 4).

Pinch point mapping provided insights on the forest areas connecting Similipal-Jada-Siddhamath-Khalasuni-Anantpur-Satkosia (1-2-3-7-8-10) and Satkosia-Subarnagir-Karlapat-Kothagarh (10-11-14-15) (Fig. 4), as these paths are surrounded by degraded land, densely populated human settlements, and agricultural fields. Rest corridors have broad areas facilitating the species’ movement throughout the state.