Applying species distribution models in public health research by predicting snakebite risk using venomous snakes' habitat suitability as an indicating factor.

Snakebite envenoming is an important public health problem in Iran, despite its risk not being quantified. This study aims to use venomous snakes' habitat suitability as an indicator of snakebite risk, to identify high-priority areas for snakebite management across the country. Thus, an ensemble approach using five distribution modelling methods: Generalized Boosted Models, Generalized Additive Models, Maximum Entropy Modelling, Generalized Linear Models, and Random Forest was applied to produce a spatial snakebite risk model for Iran. To achieve this, four venomous snakes' habitat suitability (Macrovipera lebetinus, Echis carinatus, Pseudocerastes persicus and Naja oxiana) were modelled and then multiplied. These medically important snakes are responsible for the most snakebite incidents in Iran. Multiplying habitat suitability models of the four snakes showed that the northeast of Iran (west of Khorasan-e-Razavi province) has the highest snakebite risk in the country. In addition, villages that were at risk of envenoming from the four snakes were identified. Results revealed that 51,112 villages are at risk of envenoming from M. lebetinus, 30,339 from E. carinatus, 51,657 from P. persicus and 12,124 from N. oxiana. Precipitation seasonality was identified as the most important variable influencing distribution of the P. persicus, E. carinatus and M. lebetinus in Iran. Precipitation of the driest quarter was the most important predictor of suitable habitats of the N. oxiana. Since climatic variables play an important role in shaping the distribution of the four venomous snakes in Iran, thus their distribution may alter with changing climate. This paper demonstrates application of species distribution modelling in public health research and identified potential snakebite risk areas in Iran by using venomous snakes' habitat suitability models as an indicating factor. Results of this study can be used in snakebite and human-snake conflict management in Iran. We recommend increasing public awareness of snakebite envenoming and education of local people in areas which identified with the highest snakebite risk.

Introduction

Snakebite envenoming is known as an important public health problem and the cause of medical emergencies around the globe[1-14]. On Earth, between 421,000 and 1.2 million people are envenomed by venomous snakes annually and around 125,000 deaths per year are attributable to snakebite envenoming[1,6,14,15]. Snakebite envenoming is mostly described as a neglected public health issue in the tropics[2,6,16,17], however, it is also an important challenge for public health in temperate areas like Iran[18].

Iran is home to 81 snake species[14,19] of which 25 are venomous (Nine species are sea snakes and 16 species are terrestrial snakes). Macrovipera lebetinus (Linnaeus, 1758), Echis carinatus (Schneider, 1801), Pseudocerastes persicus (Duméril, Bibron & Duméril, 1854) and Naja oxiana (Eichwald, 1831) are widespread in Iran and are responsible for the most snakebite incidents in the country[18,20]. A study reported 53,787 cases of snake bites between 2002 and 2011 in Iran[18]. Despite considerable research into the phylogeny, taxonomy, morphology and ecology of venomous snakes in Iran[21-28] snakebite envenoming has received less attention[18,29]. In fact, snakebite is an important uninvestigated public health problem and conservation challenge in Iran[30-32]. Thus, more effort should be made to identify areas with high snakebite risk and reduce envenoming risk from snakes.

Species Distribution Models (SDMs) have found an important application in biodiversity research[33-35]. They are employed in studying habitat suitability[36-38], identifying environmental drivers of species distribution[39-44] and predicting impacts of climate change on biodiversity[45-50]. Species distribution models are successfully used to identify suitable habitats of species even in areas with no distribution records[51-55]. Thus, these models can be used to identify suitable habitats of venomous snakes as proxies of snakebite risk[12,56-58] in data poor regions like Iran.

The main goal of this paper was to apply SDMs and produce a spatial risk model for snakebite in Iran. Five distribution modelling methods[35] and distribution data of M. lebetinus, E. carinatus, P. persicus and N. oxiana were used to produce their habitat suitability models. Then, the five habitat suitability models of each species were combined by ensemble approach[59] and finally the four species ensemble models were multiplied to identify potential snakebite risk. We also determined a number of villages whose population are at risk of envenoming by these four snakes in Iran. Since it believed that E. carinatus is responsible for the most mortality and morbidity cases of snake bite in Iran[18] we hypothesized that more villages are at risk of envenoming form E. carinatus compared to M. lebetinus, P. persicus and N. oxiana.

Results

Model performance

All models developed in this study performed well based on the three model performance evaluation metrics, AUC, TSS and Boyce index (see plates a in Figs. 1, 2, 3, 4).

Figure 1

Variables importance (a), models performance (b) and habitat suitability model (c) of Echis carinatus in Iran. Map was generated using QGIS 3.4.1 (https://www.qgis.org).

Figure 2

Variables importance (a), models performance (b) and habitat suitability (c) of Macrovipera lebetinus in Iran. Map was generated using QGIS 3.4.1 (https://www.qgis.org).

Figure 3

Variables importance (a), models performance (b) and habitat suitability (c) of Pseudocerastes persicus in Iran. Map was generated using QGIS 3.4.1 (https://www.qgis.org).

Figure 4

Variables importance (a), models performance (b) and habitat suitability (c) of Naja oxiana in Iran. Map was generated using QGIS 3.4.1 (https://www.qgis.org).

Habitat suitability of Echis carinatus

Based on ensemble model, southern part of Iran, north of Persian Gulf and vast areas in central parts of the country are identified to have highest suitability for E. carinatus (Fig. 1). Precipitation seasonality and NDVI were the most important determinant of habitat suitability of the species across the country.

Habitat suitability of Macrovipera lebetinus

The most suitable habitats of M. lebetinus are located in the Zagros Mountains, Alborz Mountains, Kopat-Dagh Mountains as well as in some isolated mountains in central Iran (Fig. 2). Precipitation seasonality, precipitation of the driest quarter, and human footprint were the most important predictors of suitable habitats for the species.

Habitat suitability of Pseudocerastes persicus

Central, southwest and northeast of Iran have highest habitat suitability for P. persicus. While, northern parts of the country are not suitable for this species (Fig. 3). Results showed that precipitation seasonality and NDVI were the most important determinant of habitat suitability of the species across Iran.

Habitat suitability of Naja oxiana

Naja oxiana’s the most suitable habitats are located in northeast Iran around Kopat-Dagh Mountains as well as some isolated patches in eastern parts of the country (Fig. 4). Precipitation of the driest quarter was the most important predictor of suitable habitats of the species.

Snakebite risk model

The four venomous snakes’ ensemble habitat suitability models were combined to develop a snakebite risk model for Iran (Fig. 5). Results showed that Khorasan-e-Razavi, east of Semnan, north of Khorasan-e-Jonobi and south of Khorasan-e-Shomali provinces have highest snakebite risk in Iran. West of Khorasan-e-Razavi province has high suitability for the four venomous snakes.

Figure 5

Snakebite risk model in Iran. The snakebite risk model was developed based on combined habitat suitability models’ of Macrovipera lebetinus, Echis carinatus, Pseudocerastes persicus and Naja oxiana. Map was generated using QGIS 3.4.1 (https://www.qgis.org).

Villages at risk of envenoming

Number of villages that are at the risk of envenoming by each of the four snakes (Table 1, Fig. 6) were determined. Results revealed that 51,112 villages are at risk of envenoming from M. lebetinus, 30,339 from E. carinatus, 51,657 from P. persicus and 12,124 from N. oxiana. Area of envenoming risk by each species was estimated (Table 1), M. lebetinus and N. oxiana are identified with largest (362,558 km2) and smallest (121,803 km2) area, respectively.

Table 1

Area (km2) and number of villages that are at the snakebite risk from Macrovipera lebetinus, Echis carinatus, Pseudocerastes persicus and Naja oxiana in Iran.

Species	Villages at high risk	Villages at moderate risk	Area
Macrovipera lebetinus	9824	41,288	362,558
Echis carinatus	5917	24,422	321,089
Pseudocerastes persicus	2150	49,507	655,415
Naja oxiana	3091	9033	121,803

Results are based on the ensemble models.

Figure 6

Villages at risk of envenoming from Macrovipera lebetinus, Echis carinatus, Pseudocerastes persicus and Naja oxiana. Villages at high risk are shown with red circle and villages at moderate risk are shown with yellow circle. Maps were generated using QGIS 3.4.1 (https://www.qgis.org).

Discussion

With this research the first snakebite envenoming risk model was produced at fine resolution (~ 1 km2) in Iran by modeling and multiplying habitat suitability of four medically important venomous snakes which are responsible for the most snakebite incidents in the country[18]. Northeastern parts of Iran were identified to have highest snakebite risk in the country. Our results showed 51,657 villages are at the risk of envenoming from P. persicus which is not in concordance with our hypothesis that states that the highest number of villages are at risk of envenoming form E. carinatus. Results showed that thousands of villages are located in suitable ranges of the four venomous snakes. These villages and villagers are at risk of envenoming from the local species. The snakebite risk model shows which parts of Iran are at risk of envenoming from two, three or even four of snakes. All provinces of Iran, except those in the northwest of the country, are at risk of envenoming from at least two venomous snakes. This highlights importance of snakebite envenoming as public health problem in Iran.

Under climate change some venomous species may expand their distribution ranges[57,58,60], thus, envenoming risk will likely vary[12,17,61]. For instance, Nori et al.[12] modelled distribution of five venomous snakes in Argentina for 2030 and 2080. They found that the snakes’ suitable climate spaces will increase in human populated areas of the country. In another study, Zacarias and Loyola[58] modelled current and future distributions of 13 snakes in Mozambique and showed that venomous snake distribution will change under climate change. They concluded patterns of snakebite risk may change due to climatic changes[58]. Results of current research revealed that climatic variables play an important role in shaping the distribution of four venomous snakes in Iran, thus their distribution may alter with changing climate. It is predicted that suitable habitats of E. carinatus will increase in Iran[62]. This species is an important source of snakebite in the country[18] and its distributional range will likely increase under climate change[62] putting more populations and human settlements at risk of envenoming by this snake until 2070.

Venomous snake populations are declining and many of them are listed by the IUCN Red List as Vulnerable, Endangered or Critically Endangered[63]. Conservation of snakes especially venomous snakes is a big challenge[63] as it is not easy to convince people to conserve venomous snakes which are a significant cause of human mortality and morbidity[6,14,15,56]. It is necessary to identify areas with high risk of snakebite envenoming and prioritize those areas for snakebite risk management in each country. In this study potential risk areas were identified by using venomous snakes’ habitat suitability as an indicator of snakebite risk.

One possible reason of venomous snakes’ conflict with humans is the expansion of agricultural lands. For example, most of the envenoming events occur during agricultural activities in Khorasa-e-Razavi and Sistan and Baluchistan provinces (A. Khani personal observations). According to the Ministry of Health and Medical Education, rural people were victims of around 70% of snake bites in 2011 in Iran[18]. Results of this study can be used to reduce snakebite and venomous snake conflicts with local people, farmers and shepherds in the country. Snakebite risk can be reduced through community education[64-66]. We encourage education of local people about snakebite prevention measures in areas with highest snakebite risk. There are simple solutions to prevent snakebite envenoming in areas with high risk, like villages in west of Khorasan-e-Razavi province. For example, using bed nets and protecting feet, ankles and lower legs by wearing boots can significantly reduce snakebite envenoming[65,67]. We also suggest that areas where snakebite envenoming risk is high should be monitored to determine envenoming events and villagers in these areas must always have access to antivenom supplies[14].

Species distribution models are becoming important tools in public health research[12,17,68-77]. We encourage public hearth researchers to apply SDMs in developing snakebite risk map using venomous snakes’ habitat suitability as an indicator, especially in data poor regions of the world[51-53,78]. Our approach has the potential for practical application in other countries with high snakebite risk.

Materials and methods

Occurrence data

Distribution records of the M. lebetinus, E. carinatus, P. persicus and N. oxiana (Fig. 7) were collected through long-term fieldwork and opportunistic observations from 2005 to 2019, and from published papers (Appendix S1). These four snakes were selected because they are responsible for the most snakebite incidents in Iran[18] and have the widest distribution range across the country[19,79]. By combining presence records from the two sources 89 distribution records were obtained for M. lebetinus, 68 records for E. carinatus, 54 records for P. persicus and 37 records for N. oxiana (Fig. 8 and Table S1).

Figure 7

Photos of the four medically important venomous snake (a) Saw-scaled Viper (Echis carinatus), (b) Central Asian Cobra (Naja oxiana), (c) Levantine Viper (Macrovipera lebetinus), and (d) Persian Horned Viper (Pseudocerastes persicus) modelled in this research to map snakebite risk in Iran. Naja oxiana is listed as Data Deficient, Pseudocerastes persicus as Least Concern and Echis carinatus and Macrovipera lebetinus were not listed in the IUCN Red List[80]. Photos by Masoud Yousefi.

Figure 8

Distribution records of Macrovipera lebetinus, Echis carinatus, Naja oxiana, and Pseudocerastes persicus in a topographic overview of Iran. Provinces numbers should read as follows; (1) West Azerbaijan, (2) East Azerbaijan, (3) Ardabil, (4) Kurdistan, (5) Zanjan, (6) Gilan, (7) Kermanshah, (8) Hamedan, (9) Qazvin, (10) Alborz, (11) Mazandaran, (12) Ilam, (13) Lorestan, (14) Markazi, (15) Qom, (16) Tehran, (17) Khuzestan, (18) Chahar Mahaal and Bakhtiari, (19) Isfahan, (20) Semnan, (21) Golestan, (22) Khorasan-e-Shomali (Northern Khorasan), (23) Khorasan-e-Razavi, (24) Yazd, (25) Khorasan-e-Jonobi (Southern Khorasan), (26) Kohgiluyeh and Boyer-Ahmad, (27) Fars, (28) Kerman, (29) Sistan and Baluchestan, (30) Bushehr, (31) Hormozgan. Map was generated using QGIS 3.4.1 (https://www.qgis.org).

Environmental data

Seven environmental variables (Table 2) related to climate, topography, vegetation, and human footprint were used to develop the four snakes’ habitat suitably models[12,17,56,81-84]. Climatic variables were downloaded from the WorldClim database at 30-s spatial resolution[85]. Normalized Difference Vegetation Index (NDVI)[86] was considered as an indicator of resource availability directly for snakes or indirectly for their prey. Snakebite risk is associated with human population density and activities[56], thus, human footprint index was included in models[83,84]. Human footprint index was produced by combining data on the extent of built environments, population density, electric infrastructure, crop lands, pasture lands, roads, railways, and navigable waterways[83,84]. Topographic heterogeneity was used as topography variable by measuring the standard deviation of elevation values in area grid cells of 1 km2 from a 90 m resolution in the Raster package[87]. Elevation layer was obtained from the Shuttle Radar Topography Mission (SRTM) elevation model[88]. A variance inflation factor (VIF)[89] was calculated for the variables using the ‘usdm’ package[90] in the R environment (v.3.4.3)[91]. to ensure low collinearity among the variables (Table 2).

Table 2

Environmental variables were used in developing habitat suitability of venomous snakes in Iran.

Variable	Description (abbreviation) and references	VIF
Climate	Temperature seasonality (bio4)[85]	1.438
	Annual precipitation (bio12)[85]	2.871
	Precipitation seasonality (bio15)[85]	2.389
	Precipitation of driest quarter (bio17)[85]	3.049
Topography	Topographic heterogeneity (SD)[92]	1.315
Vegetation	Normalized Difference Vegetation Index (NDVI)[86]	3.33
Anthropogenic	Human Footprint (HFP)[83,84]	1.466

Snakebite risk mapping

An ensemble approach[35,59] was applied to model habitat suitability of the M. lebetinus, E. carinatus, P. persicus and N. oxiana, using five methods: generalized boosted models (GBM)[92], generalized additive models (GAM)[93], maximum entropy modelling (MaxEnt)[94], generalized linear models (GLM)[95], random forest (RF)[96] in the biomod2 package[97]. Since these methods need background data points we generated a randomly drawn sample of 10,000 background points (e.g., pseudo-absence points) from the extent of the study area using the PresenceAbsence package[98]. We also estimated contributions of environmental variables in shaping each venomous snake species distribution using the biomod2 package[97]. We calibrated the models using 80% of records drawn randomly and used as training data, and evaluated their performance using the remaining 20% of the data (test dataset). Several metrics were introduced for assessing performance of habitat suitability models[34,99-101]. But it is not certain that which metrics are best for determining habitat suitability models’ performance[34]. Thus, we applied three well known metrics to assess the produced habitat suitability models’ performances. The true skills statistic (TSS), area under the receiver operating characteristic curve (AUC), and the Boyce index[34,99-101]. According to literature[34,101], Boyce index and TSS values range from − 1 to + 1, where + 1 indicates perfect performance and value of zero meaning random predictions. AUC values range from 0 to 1, a value of 0.5 indicates that the performance of the model is not better than random, while values closer to 1.0 indicate better model performance[34,99,100]. Then areas associated with high snakebite risk in Iran were identified by multiplying habitat suitability models of the four species.

To quantify snakebite risk in Iran in more detail number of villages that are at risk of envenoming from the four snakes determined and area of snakebite risk calculated in the Raster package[87]. For this continuous habitat suitability models were converted to suitable/unsuitable maps using maximum test sensitivity with a specificity threshold[102], and then overlaying 185,000 (population in these villages range from less than 50 to 5000 individuals) villages with each snake model.

Supplementary Information.