Ecological Niche Modelling of King Cobra, Ophiophagus hannah (Cantor, 1836) in Nepal

The fragile ecosystem of greater Himalaya is home to diverse flora and fauna which are vulnerable to climate change impact. This study aimed to determine the suitable habitat of King Cobra Ophiophagus hannah (Cantor, 1836) in the current scenario and near-future scenario of the year 2040-2060 (RCP 2.6, RCP 4.5, and RCP 8.5). Geographic coordinates of its occurrence were obtained from published literature and environmental layers were obtained from worldclim.org and processed using ArcGIS and software R. The model was developed using MaxENT algorithms. The model was validated using the area under the curve (AUC) and True skill statistics (TSS), which showed that the model is very good (AUC =0.86) and (TSS=0.6). The results showed that altitude is a prime factor affecting the distribution of King Cobra in Nepal with a percent contribution of 31% followed by actual evapotranspiration 14.4% and least by Mean diurnal range (Bio2nep) 7%. The model predicted an area of 41,214 Km2 as suitable habitat for King Cobra in Nepal. The Chitwan National Park possesses a maximum suitable area (726.89 km2) followed by Chitwan-Buffer Zone (601.68 sq. Km2) and least by Sagarmatha National Park (1.73 km2) and Shey-Phoksundo-buffer zone (1.73 Km2). Our data indicate urban areas being the most suitable sites followed by open forest covers. The projection indicates contraction of 8% and 4% in a suitable area for RCP 4.5 and RCP 2.6 respectively, whereas RCP 8.5 showed expansion of 4%. Having probably suitable habitat in urban areas research proposes frequent awareness camping about conservation and protection of vulnerable King Cobra in Nepal.

King Cobra dwell close to water sources in the edges of the forest, mangrove swamps, pristine rainforests, plantations, agricultural field and open scrublands in the vicinity of water bodies and prefers to build a nest in the area with leaf litter height of average 10-17cm (Rao et al., 2013). Cantor (1836) found King Cobraliving in hollow trees and observed on branches of such trees in Calcutta, India. King Cobra inhabits tropical and subtropical forests but prefers moist temperate forests (Shah, 2000;Marshall et al., 2020). King Cobra is reported from almost all Terai and Mid-hill districts of Nepal (Pandey et al., 2018;Thapa et al., 2019;Baral and Koirala, 2020;Devkota et al., 2020b;Rai, 2020;). Mating, nesting, and juveniles of King Cobra are reported from an altitude range of 1000 to 1500 m .
An inspiring number of King Cobra related short notes, checklists on distribution, habit, and habitat are published from different parts of the country in the last decade. Published anecdotal notes on distribution, habit, and habitat cannot represent the entire landscape of Nepal and related biologically meaningful data (Baral and Koirala, 2020;Devkota et al., 2020;Rai, 2020). With the increasing number of a report on injured specimens, killing, and nesting sites of King Cobra from Kaski, Palpa, Parbat, Kathmandu valley, Makawanpur, Sindhuli and Sarlahi (Baral et al., 2018Thapa et al., 2019;Bhattarai et al., 2020;Devkota et al., 2021) further studies on the distribution of King Cobra in Nepal is crucial for its conservation. Furthermore, intentional human killing, increasing human population, habitat fragmentation due to land-use change are major threats to King Cobra in Nepal (Devkota et al., 2021).
Nepal is vulnerable to climate change impact which will adversely affect the habitat of endangered flora and fauna in the future (Gaire et al., 2014;Bhattacharjee et al., 2017). Herpetofauna is affected by climate change (Araújo et al., 2006) and studies on its impact are limited in Nepal. Taking into account, climate change consequences, conservation threats, and insufficient literature (Schleich and Kästle, 2002;Baral et al., 2018;Bhattarai et al., 2018;Pandey et al., 2018;Devkota et al., 2021)this study used the MaxENT algorithm to determine the distribution of King Cobra in Nepal. MaxENT is a Species Distribution Model (SDM) widely used for predicting habitat or suitable areas of flora and fauna in particular geographic space and time (Phillips and Schapire, 2004). MaxENT has been used in some Asian regions for predicting the habitat of mammals like pangolin in Nepal (Suwal et al., 2020), Golden Langur in Bhutan (Thinley et al., 2019), Himalayan Musk Deer (Lamsal et al., 2018), and Bats in Nepal (Thapa et al., 2021).
According to Siqueira et al. (2014), information on biogeographical variation in species endemic richness is critical for understanding and conservation of biological diversity, and to develop rigorous conservation plans for a region. Understanding the distribution of King Cobra and factors governing the same would help in its conservation in Nepal (Vetaas and Grytnes, 2002).
The objectives of the study are to a) determine the suitable habitat of King Cobra; Ophiophagus hannah (Cantor, 1836) in the current scenario in Nepal b) investigate the impact of climate change on King Cobra's near-future Representative Concentration Pathway (RCP) 2.6, 4.5 and 8.5 scenarios in about the year 2040-2060 in Nepal. The result of the study would help in developing appropriate mitigation strategies using deterrents or barriers to reduce human King Cobra Conflict. The ecological idea about the species would help local government, civil organizations, and wildlife officials to plan conservation strategies.

Study Area
This study was conducted in Nepal (26.347 0 to 30.473 0 N and 80.060 0 to 88.204 0' E) where the elevation ranges from 59m to the top of the world i.e. Mount Everest at 8848.86m(Department of Survey, 2020; Kathmandu Post, 2020). According to the available data from 1980 to 2010, the annual normal precipitation is lowest with 271.4mm in Jomsom and highest (5514.7mm) in Lumle, and the mean normal temperature for Nepal is 10.9 0 C in Jomsom and highest (24.6 0 C) in Janakpur (Department of Hydrology and Meteorology, n.d.).
With an area of 147,516 km 2 , it encompasses various geographical realms and it is a part of one of the 36 global biodiversity hotspots in the eastern Himalaya (Ministry of Forest andSoil Conservation, 2014, 2009;Myers et al., 2000). A total of 25% of the country's land is designated as a Protected Area (PA), with 14 National Parks, one wildlife reserve, six conservation areas, and one hunting reserve (Central Bureau of Statistics, 2020; Ministry of Forest and Soil Conservation, 2009).

Data Collection
A total of 131 occurrence points of King Cobra in Nepal was collected from published literature (Baral et al., 2018;Pandey et al., 2018;Thapa et al., 2019). MaxENT required environmental layers for predicting suitable areas. A total of 19 bioclimatic variables of spatial resolution 30 arc-second (~1km 2 ) available in worldclim.org was obtained for modeling King Cobra for the current scenario. For the future projection of climate suitability of King Cobra in Nepal, future data from worldclim.org was obtained for three Representative Concentration Pathway (RCP) 2.6, 4.5, and 8.5 for the year 2040-2060 from the general climatic model for Interdisciplinary Research on Climate version 6 (MIROC6). The choice of the MIROC6 was motivated by the work on the impact of climate change on future climate and habitat distribution of Himalayan Musk Deer (Lamsal et al., 2018). The digital elevation model (DEM) obtained from worldclim.org was used as altitude (WorldClim, 2020). The land use land cover map was obtained from the Copernicus website and classified as prescribed in the product user manual (PUM) (Buchhorn et al., 2020). Actual evapotranspiration was obtained from Consortium for Spatial Information (CGIAR-CSI) portal, which is an initiative of the Consultative Group for International Agriculture Research (CGIAR) (Trabucco and Zomer, 2018).

Ecological Niche Modeling
Environmental layers were clipped for Nepal using the spatial analyst tool of ArcGIS. The layers were resampled to the same spatial extent and pixel using resample package in software R before running MaxENT. After resampling, data were converted into ASCII format which is readable by MaxENT algorithms using species distribution toolbox (SDMtoolbox) in ArcGIS 10.4. (Brown, 2014;Chen et al., 2020;Elith and Leathwick, 2009). As bioclimatic variables are derived from monthly temperature and precipitation (WorldClim, 2020), these variables have a high correlation among them, which can misinterpret the model output (Yoo et al., 2014). Among 19 bioclimatic variables, five were obtained after the multicollinearity test (Table 1) with a Pearson correlation coefficient (R) value less than 0.8. Altitude, land use land cover, and actual evapotranspiration were other input variables for modeling the habitat of King Cobra in Nepal.

MaxENT settings
After the multicollinearity test, Ecological Niche Model evaluation (ENMeval) was performed in software R using package ENMeval ver. 0.3.1. This was to determine the potentially available background point and features to run the MaxENT model. ENMeval result suggested random multiplier value 2.5, 10,000 background points (pseudo-absence records) and 5,000 iterations. Cross-Validation was selected as a replicate type with 10 replication of the model. The features suggested for the model were Linear, Quadratic, Hinge, Product, and Threshold (LQHPT).

Model Validation and Projection
Model performance was tested using Receiver Operating Characteristics (ROC) curve and the Area under the Curve (AUC) and True Skill Statistics (TSS) (Phillips et al., 2006). The AUC tests the model by providing a threshold-independent assessment of the model using a web-based program that ranges from 0.5-1 (Radosavljevic and Anderson, 2014). If the AUC value is greater than < 0.5 it means the model performs better than random estimation. An AUC value between 0.9 and 1.0 shows excellent model performance; 0.8-0.9 = good; 0.7-0.8 = average; 0.6-0.7 = poor and 0.5-0.6 = insufficient (Phillips et al., 2006). MaxENT output was imported into ArcGIS and made into a binary map (suitable and unsuitable area) and projected using distribution changes between binary SDM toolbox functions.

Model performance and contribution of environmental variables in the present scenario
The model accuracy of King Cobra for the current situation was found good; AUC= 0.858, TSS = 0.6and SD = 0.053. A total of eight variables were used to model the suitable habitat of the Ophiophagus hannah in Nepal. After the multicollinearity test (Pearson correlation coefficient) five bioclimatic variables were selected (Precipitation of coldest quarter, precipitation of the driest month, Temperature Annual Range, Isothermality (BIO2/BIO7) (*100) and Mean Diurnal Range (Mean of monthly (max temp -min temp)). The selection was based on the general thump rule that Pearson's R-value should not exceed 0.8 (Thinley et al., 2020). Two topographical variables (Digital Elevation Model and Land Use Land Cover) and evapotranspiration were selected as predictor variables to run the model in MaxENT. Except forbio7nep, all other variables contributed to better performance of the model as shown by Jackknife of AUC (figure 1). Variables affecting the distribution of King Cobra in Nepal were classified into two categories using topographic and bioclimatic variables based on static and dynamic properties of variables. Land use land cover is assumed to be constant for all RCP scenarios for this study. In the current model, the distribution of King Cobra in Nepal is affected mostly by altitude (31%) followed by actual evapotranspiration (14%) and least by temperature annual range (bio5nep -bio6nep) (2 %) (

Figure 2. Suitable habitat of Ophiophagus hannah in Nepal
Using maximum test sensitivity plus specificity Cloglog threshold of 0.5664, a total of 41,214 km 2 of the suitable area was predicted for King Cobra for the current scenario in Nepal. The King Cobra has connected suitable habitats to thrive in Province 1, 2, Bagmati, Gandaki, and Lumbini. Karnali and Sudur Paschim have the least area with 0.6% (252 km 2 ) and 1%(398 km 2 ) dispersed suitable habitats respectively. For easy interpret the role of variables is represented with curves in Figure 3. Plots represent a MaxENT model for each corresponding variable. The mean diurnal range (Bio2nep) had a better estimation of the suitable habitat at 21 0 C and it picked till 33 0 C to remain constant on the increment of temperature with sudden fall at 26 0 C. Isothermality (Bio3nep) above 42.5 had a better linear contribution in estimating suitable habitat till 57.1 and the habitat suitability remains constant on further increment. Precipitation of driest month (Bio14nep) has peak effect at 3.2 mm (millimeter) while it had linear descending slope on both decreasing and increasing values. Precipitation of quarter (Bio19nep) below 60.2 mmhad better contributions in predicting the habitat. Altitude below 2315 m had contributed better (habitat suitability >0.5) than random for estimation of habitat. Evapotranspiration above 8,156 mm per year contributed well.
Builtup area and deciduous broadleaf forest with a canopy cover of less than 70% are the most suitable land cover as per the model. The data in primary literature were collected from the opportunistic observation of the reference. We overlaid the present predicted suitable habitat of King Cobra over the Protected Areas of Nepal for cross-validation of the claim.
We observed that 89% (36,559 Km 2 ) of the predicted suitable habitats are out of the protected area or in the Buffer Zone. Part of Chitwan National Park, Southern area of Annapurna Conservation area, Parsa National Park, Koshi Tappu Wildlife Reserve, and Shivapuri -Nagarjun National Park among others are found inside the predicted habitat area (Table3).

Future potential habitat distribution of King Cobra
This study is the first of its kind which evaluates the impact of climate change on the habitat of King Cobra in Nepal and Eastern Himalaya. The study evaluates the climate suitability for three IPCC scenarios (RCP 2.6, RCP 4.5 & RCP 8.5) focusing on Nepal Himalaya. The great Himalayas harbors diverse flora and fauna due to its peculiar biogeography and variable climates (Kumar and Chopra, 2009). Due to the changing climate, the flora and fauna of the Himalayas have been under greater threat of extinction and King Cobra is no exception. The future climate suitability of King Cobra in Nepal for three IPPC scenarios was modeled using MaxENT.The model of an average value of bioclimatic variables for the year 2040-2060 in three scenes of RCP 2.6, RCP 4.5, and RCP 8.5 was generated. There was the largest expansion (7045 Km 2 ) suitable habitat area for the King Cobra area at RCP 8.5, while the suitable habitat area has huge range contraction at RCP 4.5 (Table 4). The current and future habitat suitability of King Cobra in Nepal was predicted using occurrence-only data with environment layers (bioclimatic-topographic) using MaxENT software (Phillips and Schapire, 2004). The predicted distribution of King Cobrausing MaxENT modeling is the habitat probability of habitat niche with different climate scenarios. Similar studies have been done with crops, invasive species, and threatened fauna (Lamsal et al., 2018;Chhogyel et al., 2020).
In this study, five climatic factors, two topographic factors, and evapotranspiration factors were identified as potential factors affecting the distribution of King Cobra in Nepal. MaxENT showed that bioclimatic variables like mean diurnal temperature range, temperature seasonality, and the annual mean temperature affected the distribution of King Cobra in Nepal. The results are similar to a study conducted by Jankowsky (1973)who concluded reptiles as a solar exothermic animals with complex physiological and behavioral mechanisms for maintaining body temperature.
The 31% percent contribution showed that altitude had the highest impact on the distribution of King Cobra in Nepal. Altitudinal gradients and the physical environment are prime factors determining spatial and time-related distribution, amplitude, and richness patterns of organisms (Palmer, 1994). The MaxENT predicted that King Cobra in Nepal can live up to 2315m in the current scenario.
Rapid changes in altitude and associated habitat particularly in the mountainous landscape can provide insight into the relationship between environmental heterogeneity and the local or regional species diversity (Bateman et al., 2010). A strong association between altitude, changes in climate, and vegetation was observed (Elkan and Cooper, 1980;Körner, 2007), thus species assemblages can shift rapidly over relatively short distances. These zones of rapid transition can result in areas of increased species diversity (Heaney, 2001). Nevertheless, there is conflicting evidence and debate regarding the determinants of the diversity along such gradients (Rowe, 2009).
The other climatic variables affecting the habitat of the King Cobra are the precipitation of the warmest quarter (21.1 %), precipitation of the driest month (6.8 %), and precipitation of quarter (6.2 %) from the study area. The result conforms with the report of (Sutton et al., 2017) which states that the habitat of snakes might be associated with other climatic factors.
Open habitats within forests are important because they provide access to direct sunlight and thermal mosaic used for thermoregulation (Owen, 1989;Pianka et al., 2003). The MaxENT revealed that King Cobra is mostly associated with urban areas followed by open forested areas. The results concur with the finding of Marshall et al. (2020a) who reported the habitat of the King Cobra is mostly associated with vegetated semi-natural areas. The result also concurs with the finding which states King Cobra dwell nearby agricultural areas (Marshall et al., 2020). The output of the model contradicts a study from Agumbe in India which observed that the King Cobra has a preference for closed-canopy forest (Shankar et al., 2011). A study of King Cobra using radiotelemetry in Thailand (Marshall et al., 2018) supports that the King Cobra dwell around settlements and agricultural lands. Marshall et al. (2018) warned that King Cobra is killed in human settlement areas or urban areas due to different anthropogenic causes (road mortality, pollution, fish traps, and direct persecution). Recent studies published from Nepal show land-use change and habitat fragmentation as the main factors contributing to the human-caused deaths of King Cobra (Devkota et al., 2021).
King Cobrahavehome ranges from 19 ha to 710 ha and is frequently observed in the forest edges Marshall et al., 2019Marshall et al., , 2020. The result of modeled habitat demands the requirement of an awareness campaign to conserve the species, as the King Cobra's suitable habitat is around the human settlement area. Deforestation of forests and habitat fragmentation are the major threats to biodiversity (Harrison et al., 2009). King Cobra, the largest snake is susceptible to be killed or be a victim of human excitement and fear due to its size and other characters. Mitigation of road kills (maintaining slow speed in the evening near herpetofauna habitats), leaving the snakes in their habitat of sighting, and admiring snakes for their ecological roles are potential sources of arguments for the conservation of King Cobra.
The model predicts expansion at RCP 8.5 while there is a reduction in habitat range in RCP 4.5 with an approximate3 0 C rise in temperature. The expansion in habitat range can be due to an increase in temperature of 5 0 C since the King Cobra is solar exothermic species (Jankowsky, 1973). Global warming on a huge scale would challenge the survival of most flora and fauna. Reptiles lack capabilities to shift distribution in with climate change (Araújo and Pearson, 2005). Thus, expansion in suitable habitat does not certify that King Cobra will live in that space. As the scenario threatens the existence of the biotic environment, the extrapolation of data for higher RCP scenarios demands further validation and precautions. In a biodiversity-rich region like Nepal, preventing further loss will remain a challenge as the nation face climate crisis vulnerabilities (Ministry of Forest and Soil Conservation, 2014). Thus, aggressive climate change mitigation measures, ambitious conservation efforts, and food system transformation can support the survival of biodiversity (Sodhi et al., 2010).

CONCLUSION
MaxENT model showed the affinity of King Cobra's habitat towards human settlement. Enlightened human practices can help in the conservation of the King Cobra and its habitat. Habitat suitability in urban and protected areas demands the requirement of a species conservation plan and an immediate awareness campaign. Altitude played a major role to define habitat in the model. A further ecological survey with randomization can validate the predicted habitat. Suitable habitat in the year 2040-2060 for three scenes RCP 2.6, RCP 4.5, and RCP 8.5 remained unchanged for 24%, 20%, and 25% of current habitat respectively. Higher RCP situation in future threatening existence of many species may not be favorable for apex predator such as King Cobra. Results acquired in the modeling need further validation and precaution as the result is based on presence data.