Range dynamics of Walterinnesia morgani (Serpentes, Elapidae) during climatic oscillations in Iran

Reptiles have a crucial part in maintaining global biodiversity and the functioning of dynamic ecosystems, owing to their ecological roles and functions. Nevertheless, these organisms are susceptible to human-induced disruptions and the deterioration of their habitats, leading to their categorization as the third most endangered group of vertebrates on a global scale. Understanding the spatial distribution of reptiles is crucial due to their often specific habitat needs and limited vagility. Morgan’s black cobra ( Walterinnesia morgani ) is a secretive venomous snake species that has thus far received little attention in Iranian scientific literature. The aim of the present study was to determine the existing distribution pattern of the cobra and to speculate on how climatic changes might affect it. Maximum entropy modeling was used to examine a dataset consisting of 16 occurrence records gleaned from field observations and the literature. The niche of the species was predicted using current and future climate change forecasts and bioclimatic and topographical characteristics. The models predicted a future reduction in the wide distribution region of W. morgani in southern and western Iran. It was discovered that climatic factors like temperature range, precipitation dynamics, and river proximity all played a key role in shaping the pattern of distribution. The predicted suitable areas for W. morgani were dependent on water sources; however, future scenarios showed a decline in suitable habitats. This study underscores the importance of conservation efforts in light of the potential implications of climate change on this species. To further understand the range shifts and adaptive strategies of the species, further study of its ecology and dispersal dynamics is required.


Introduction
Reptiles, although they account for one-third of global terrestrial vertebrate diversity, have lagged behind other groups, such as birds and mammals, in terms of ecological studies (Biber et al. 2023).Due to the fact that reptiles are typically characterized by specific habitat needs and limited vagility (Guedes et al. 2018), understanding their spatial distribution is crucial.In this regard, revealing species distribution patterns improves our knowledge of the relationship between landscape and species.Climatic oscillations, like topographic and vegetation cover patterns, have contributed to the establishment of the distribution of species in different biomes (Bobrowski et al. 2018;Şahin et al. 2022a).
Although reptile species are widely distributed throughout the world, they are at risk of extinction due to habitat loss and degradation, pollution, invasive species, diseases, and climate change (Gibbons et al. 2000).Anthropogenic impacts-not only habitat use but also misinformation about their venomosity-have disproportionately threatened snake species (Saha et al. 2018).Therefore, evaluating the conservation status of these quite secretive organisms and determining their distribution patterns are important.
The herpetofauna of Iran comprises 81 species of snakes belonging to 34 genera and 7 Families (Rajabizadeh 2018).
The genus Walterinnesia Lataste, 1887 has two species known commonly as desert black snakes or black desert cobras including the type species Walterinnesia aegyptia Lataste, 1887 and Walterinnesia morgani (Mocquard, 1905) (Fig. 1).Both species occur only in the Middle East.The type locality of W. morgani is Khuzestan Province, southwestern Iran (Uetz et al. 2023).Populations of the species in Iran were previously considered as W. aegyptia by Farzanpay (1989), Latifi (1985Latifi ( , 1991Latifi ( , 2000)), and Leviton et al. (1992).Nilson and Rastegar-Pouyani (2007), based on the lower scale row counts around the neck and a banded juvenile pattern (Fig. 2), recognized the eastern populations of the genus (from Turkey and Saudi Arabia to Iran) as W. morgani; however, Rajabizadeh (2018) considered the Iranian populations of the genus as Walterinnesia aegyptia morgani.
Walterinnesia morgani, Morgan's black cobra or black desert cobra, is a venomous snake species that inhabits the Arabian Peninsula (Kuwait and Saudi Arabia), the extreme south of Turkey, Syria, and the majority of Iraq, as well as western and southern Iran (Ilam, Kermanshah, Khuzestan, Bushehr, Fars, and Hormozgan Provinces) (Joger 1984;Gasperetti 1988;Latifi 2000;Uğurtaş et al. 2001;Sindaco et al. 2006;Safaei-Mahroo et al. 2015).In a specific concept, the Zagros Mountains display a wide biodiversity pattern, as well as the Central Iranian Plateau and the northern Persian Gulf, in one of the well-known biodiversity hotspots (Irano-Anatolian Biodiversity Hotspot) in the Palearctic Realm (Gholamifard 2011;Mittermeier et al. 2011).Walterinnesia morgani is one of the representative charismatic species of this region, where its shiny black color draws attention and, as most of the sightings of this species lead to its death, is very important from the aspect of conservation.
Black desert cobras display nocturnal activity patterns (Uğurtaş et al. 2001).During the hottest parts of the day, they seek refuge in burrows or beneath rocks.The seasonal activity of desert cobras may be influenced by local climate conditions (Baran et al. 2006).In some regions, for instance, they may be less active during the hottest period of summer and more active during the cooler months (Baran et al. 2021).Nonetheless, they are capable of being active even in extreme temperatures and are acclimated to harsh desert environments (Baran et al. 2021).
In spite of the known distribution of Morgan's black cobra, W. morgani, (IUCN status is Least Concern (LC)) in the Central and Southern Zagros Mountains, and its medical importance (Abid et al. 2020), the national literature on the species is scarce.The main reason for this is probably the difficulty of studying this mysterious and nocturnally active venomous snake.Therefore, the aim of this study was to reveal the current potential distribution pattern and speculate on the effects of climate oscillations on the potential distribution pattern of the black desert cobra in Iran.

Study area and input data
This study was conducted within the borders of the Islamic Republic of Iran (25.08-39.77°N,44.04-63.3°E;Fig. 1).A total of 16 occurrence records of W. morgani were obtained from the literature (Farzanpay 1989;Latifi 2000;Nilson and Rastegar-Pouyani 2007;Fathinia et al. 2010;Gholamifard and Rastegar-Pouyani 2012;Gholamifard et al. 2012) and our field observations (Suppl.material 1).In cases in which the locality information was not georeferenced, an online geographic system application (i.e., Google Earth Pro) was used to ascertain the most precise location possible.All records were georeferenced using the WGS84 coordinate system and checked with (QGIS Development Team 2023).

Ecological niche modelling (ENM)
Bioclimatic variables and one topographic layer (elevation) were downloaded from the CHELSA database (https://chelsa-climate.org/) at a spatial resolution of 30 arc-second raster grids (Karger et al. 2017;Brun et al. 2022).Additionally, one other topographic variable (distance to river) was obtained from the study of (Gavashelishvili et al. 2018) (Suppl.material 2).Each layer was clipped for the study area in QGIS.Pearson correlations between variables were calculated in R v4.3 (R Core Team 2023) and highly correlated variables were eliminated (r ≥ |0.8|).Six climatic and two topographic variables were retained and used in model construction based on the Pearson correlation coefficient: mean diurnal air temperature range (Bio_2), annual range of air temperature (Bio_7), daily mean air temperatures of the driest quarter (Bio_9), precipitation amount of the driest month (Bio_14), precipitation seasonality (Bio_15), mean monthly precipitation amount of the driest quarter (Bio_17), distance to river (d_river), elevation (elev).The contributions of the variables are given in Table 1.
All used variables were employed to predict the species niche under recent  and future (2071-2100) climate change projections (GFDL-ESM4, IP-SL-CM6A-LR, MPI-ESM1-2-HR, MRI-ESM2-0, and UKESM1-0-LL) with the lowest and the highest limits of the shared socioeconomic pathways (SSPs) from the Coupled Model Intercomparison Project Phase 6 (CMIP6) as these models have the best performance for Eurasia (Eyring et al. 2016;Sun et al. 2022) (Suppl.material 3).Using the R package spThin, to reduce the effects of spatial autocorrelation (Boria et al. 2014) we utilized occurrence data that were separated by more than 2 km (Aiello-Lammens et al. 2015).One of the main issues with ENM is that it may lead to model overfitting if a large number of predictor variables are used in conjunction with a small sample size (Fielding and Bell 1997).Therefore, the Akaike Information Criterion corrected for small sample sizes (AICc) was applied to the potential distribution pattern (Hurvich and Tsai 1989).
Maximum entropy modeling was utilized since it is a robust method that can be applied to presence and pseudo-absence data.Based on presence and pseudo-absence data, this algorithm can predict the presence of a species with a probability between 0 and 1 (Phillips et al. 2009).A total of 2000 background points for W. morgani were sampled at random throughout the survey area.The potential habitat suitability was modeled by implementing MaxEnt 3.4.1 in R using the 'kuenm' package (Phillips et al. 2017;Cobos et al. 2019).To create the models for this cobra species, 80% of the occurrence data was used for generating the candidate models, and the remaining 20% for independent presence as test data.
In order to optimize the model complexity for the cobra species, 31 combinations of MaxEnt's 5 feature classes [hinge (h), threshold (t), product (p), quadratic (q), and linear (l)] along with 17 regulation multiplier values (0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1, 2, 3, 4, 5, 6, 8, 10) were evaluated.This combination of options enabled us to find the best-fitting model representing our data by generating diverse candidate models (Muscarella et al. 2014;Cobos et al. 2019).Following that, the Akaike Information Criterion corrected for small sample sizes (AICc) (with the lowest values) as well as the AUC values (with the highest values) were both used to determine the best models (Hurvich and Tsai 1989).Significance tests were performed using partial ROC (Peterson et al. 2008) and predictive power with a 5% omission rate (Anderson et al. 2003).Lastly, all model inputs were transformed into binary predictions by using the minimum training presence as the threshold to distinguish suitable areas from unsuitable ones (Pearson et al. 2007;Şahin et al. 2022a).

Results
The regions anticipated to be suitable for W. morgani exhibited a significant area under the curve (AUC) value of 0.955±0.014.Furthermore, the recent historical and future bioclimatic predictions (2071-2100) under the SSP126 and SSP585 scenarios had the lowest values for Akaike Information Criterion corrected (AICc) and the delta AICc at 0.169 and 0, respectively.(Figs 2, 3A-J; Table 2).The fact that the AUC data were very close to 1 showed that the potential distribution area revealed by the locality data obtained from the distribution area of the desert cobra displayed a much better performance than a random prediction.For W. morgani, 1054 statistically significant models that had different regularization multipliers and feature classes were evaluated.Linear, quadratic and product (LQP) features along with a regularization multiplier of 0.8, and with the lowest delta AICc, was the best model.The most important variables were d_river (31.9%),Bio_9 (27.7%) and Bio_15 (14.7%), respectively.
According to the analysis, the potential distribution area of W. morgani covers a water source-dependent pattern in western Iran.Additionally, as can be seen in future climatic conditions, the potential distribution range would be narrower than under recent historical conditions.Habitat loss is expected in each future climatic scenario, but its level varies depending on the ssp level (Table 3).For instance, scenarios in ssp 126 level forecast the species range change percentile to be between -9.814 and -18.028%, however, the ssp 585 scenarios may result in a much great-er loss of between 49.337% and 71.262% of the species' suitable habitats in the distribution range (Suppl.material 4: fig.S1A-J; Table 3).As a result of all these scenarios, it is expected that the possible distribution of W. morgani will contract in southwestern Iran in the future.

Discussion
Numerous biotic and abiotic factors have significant impacts on the distribution of species (Pearson et al. 2007).To understand the responses of W. morgani to climatic oscillations during recent history and under future possible scenarios, we employed ENM as a useful analytical instrument to assess the range of the species.Climate changes affect all aspects of biodiversity, from organisms (organismal diversity) to biomes (ecological diversity) (Gaston and Spicer 2004), and pose a significant threat to the integrity of ecosystems (Bellard et al. 2012).Additionally, it is speculated that vertebrate species will face serious problems of adaptation related to alterations in their climatic niches in the near future (Quintero and Wiens 2013).It should be noted that the distribution of reptiles can be significantly affected by direct and/or indirect anthropogenic activities, such as human population size and human-driven climatic changes (Bickford et al. 2010).
Based on these results, most of the suitable predicted areas were slightly wider than the present potential distribution of W. morgani.This might be the possible effect of "d_river" as one of the highest contributors to shape the distribution pattern via the water requirement of the species.Additionally, daily and annual temperature cycles have also contributed to the potential distribution pattern as well as the seasonal factors, especially the driest seasonal precipitation dynamics.This overall bioclimatic and topographic pattern is observed in many herptile species in and nearby the study area (Naumov et al. 2020;Kurnaz and Eroğlu 2021;Şahin et al. 2022a, b;Vaissi 2022;Kurnaz 2023).
On the other hand, our results show that there will be a future decrease in suitable habitats for W. morgani, but the level of the decreasing trend varies depending on the SSP levels and different future scenario sets (Suppl.material 4: fig.S1A-J; Table 3).This trend is compatible with many ecological niche modeling studies that have been applied to herptiles in the study area and nearby regions (Kurnaz and Şahin 2021;Bozkurt 2022).Additionally, this contraction pattern has also been predicted for future distribution trends in many lizard species (Vaissi 2022).Nonetheless, the genus Lacerta and two leopard geckos exhibit the opposite trend in the Anatolian Peninsula and Iran, respectively (Hosseinian Yousefkhani and Nabizadeh 2022;Gül et al. 2023).
When compared to other large vertebrate groups, the ability of reptiles to migrate is quite restricted (Hickling et al. 2006).Moreover, our knowledge about the ecology and dispersal dynamics of W. morgani has received little or no attention.Based on Figs 2, 3, the current research assumed an infinite capacity for dispersal across species and made predictions about its range changes through 2100.Despite contractions and expansions, the potential distribution range of the species under each ssp126 scenario remains less changed (Fig. 3A, C, E, G, I; Table 3); however, under ssp585 scenarios, the species will be contract in southwestern Iran (Fig. 3B, D, F, H, J; Table 3).In the pessimistic scenario (ssp585) for the IPSL-CM6A-LR climate projection, in particular, it is expected that there would be 71.26%habitat loss with only 0.51% habitat gain and as a result, there would be a 70.74% negative range change (Suppl.material 4: fig.S1D; Table 3).Additional factors, such as geographic barriers, may influence the migration rate of a species (Morena-Rueda et al. 2012).In our case, distance to the river and elevation parameters, as topographic variables, make significant contributions to the species distribution.Even though the Zagros Mountains can be assessed as a remarkable geographic barrier, the main limitation factor is the "distance to the rivers" for W. morgani.Because it seems that the cobra tends not to occur very distant from water resources.Moreover, even though the potential distribution map displays suitable habitats in central or eastern regions of Iran, it should be noted that this could be the consequence of bioclimatic and topographic similarities in distinct areas that the species is biogeographically incapable of inhabiting.A similar pattern to ours was observed in several studies which evaluate narrower regions (Heidari 2021;Bozkurt 2022;Şahin 2022b).
On the other hand, the model algorithm that was used in this study did not take into account some parameters, such as parasitism, disease, habitat loss, and fragmentation that might affect the realistic distribution of the species (Todd et al. 2010).In the meantime, a recent study that focused on snake bite risk in Iran showed that the northern and western parts of Iran would have greater risk than the rest of the country (Yousefi et al. 2023).However, our results pointed out that southwestern Iran could be a potentially suitable region for Morgan's black cobra.No common potential distribution model was reported for W. morgani, neither country-based nor across-range, by including up-to-date topographic variables.Thus, our results suggest that W. morgani should be monitored in the near future, with more occurrence records collected in order to better understand its distribution and to increase awareness of human-cobra conflict as a result of anthropogenic activities.To conclude, it is desirable to predict the potential distribution of this species in Iran from the aspect of medical treatment of possible bites of the species and to prepare its antivenom.

Table 1 .
Percentage contribution of the environmental layers used in species distribution modeling of Walterinnesia morgani.

Table 2 .
Summary for selecting the best model for species distribution maps of W. morgani via 'kuenm' package.
AICc: A corrected AIC score, used for a small sample size by increasing the cost for each parameter.wAICc: The model weight is the relative likelihood for each model, divided by the total relative likelihood for all models that were considered.Delta AICc: The difference between the model with the lowest score (the "best" model) and the AICc score for each model.AUC: Area under the curve is a measure of the accuracy of the model.Mean AUC ratio ≥ 1.00, p < 0.05 means predictions are significantly better than a random model.