Estimating geographic patterns of ophidism risk in Ecuador

Abstract

The species richness of venomous snakes in Ecuador (~39 species) is among the highest in the world. However, until now no information exists regarding geographic patterns of ophidism. In this study, we present a detailed spatial snakebite risk map which was built by stacking weighted ecological niche models of the 19 snake species responsible for the majority of Ecuador’s envenomation cases. Our weights were based on the proportion of cases reported for each species on local epidemiological studies. Based on our analyses, we identify 184 densely populated rural communities with high snakebite risk that should be monitored by health organizations. We also identified three densely populated rural locations (Palora Metzera, Sangay and Shell) that may require special attention because they had much higher snakebite risk values than the rest.

Keywords:

• Snakebites
• Ecuador
• ophidism
• correlative ecological niche models

1. Introduction

Ophidism is considered a neglected tropical disease, despite its high incidence, and the considerable health problems (serious and permanent functional sequelae in affected persons) and high number of deaths that it causes in human populations [1,2]. Recent global estimates suggest that a minimum of 421,000 envenomings and 20,000 deaths owed to venomous snakes take place each year, and that these numbers could be up to 1,841,000 bites and 94,000 deaths [3]. And, based on the fact that envenomings occur in about one in four people who are bitten, these same authors estimated that between 1.2 and 5.5 million of snakebites occur annually.

Ecuador has a great diversity of venomous snakes (~36 species, which account to 5 and 19% of the diversity of the world and the Americas, respectively) [4], and one of the highest prevalence of snakebites in the continent [5]. However, epidemiologic studies in the country remains scarce [4]. Most of the information consists of isolated hospital reports for certain regions of the country [2]. To date, there has been one comprehensive study of temporal and social aspects of snakebites, but a geographical analysis of snakebite risk is still lacking [2].

Ideally, strategies to deal with snakebites would be based on assessments of true rates of bites and envenomation via community-based epidemiological studies, independent of the vagaries of hospital reporting, but this procedure would be expensive to apply at large scales (e.g. Ecuador) [6]. As a result, correlative ecological niche models have been recently used to infer snakebite risk at regional and continental scales [5,7].

We applied this niche modelling technique to estimate potential snakebite risk (defined here as the probability of being bitten by a venomous snake in a given location) in Ecuador integrating: (1) maps of potential distribution and environmental suitability of snake species of medical importance, and (2) information that approximates the probability of being bitten by a specific species. We also identified vulnerable human communities for snakebites based on a densely populated rural communities map.

2. Methods

2.1. Biological data

We obtained species occurrences through the Global Biodiversity Information Facility (GBIF, http://www.gbif.org) and VertNet (http://www.vertnet.org) websites (accessed to both = March 2017). We eliminated any doubtful (e.g. occurrences that fall in the ocean or in other continents) information and we applied the “Spatially Rarefy Occurrence Data” from “SDMtoolbox” [8] in ArcGIS® 10.3 (©ESRI) to minimize spatial auto-correlation and model overfitting [9]. We used a rarefication buffer of 20 km based on the clustering nature and abundance of presence points; we would have lost many records for some species if we had used a greater distance. Finally, we split presence records randomly into calibration (80%) and evaluation (20%) data-sets (sample sizes for complete, filtered, calibration and validation datasets are presented in Table 1).

Table 1. Summary of ecological niche modelling inputs and evaluations.

Species	T Occ	F Occ	N Calib	N Test	AUC r
Bothriechis schlegelii	416	173	138	35	1.797
Bothriopsis bilineata	59	56	45	11	1.683
Bothriopsis pulchra	15	15	12	3	1.718
Bothriopsis taeniata	41	31	25	6	1.090
Bothrocophias campbelli	9	4	–	–	–
Bothrocophias hyoprora	44	25	20	5	1.537
Bothrocophias microphthalmus	14	3	–	–	–
Bothrops asper	1482	279	223	56	1.746
Bothrops atrox	848	103	82	21	1.081
Bothrops brazili	32	2	–	–	–
Bothrops lojanus	4	3	–	–	–
Bothrops punctatus	23	16	–	–	–
Lachesis acrochorda	4	4	–	–	–
Lachesis muta	85	39	31	8	1.656
Leptomicrurus narduccii	52	27	22	5	1.479
Micrurus ancoralis	43	19	15	4	1.656
Micrurus annellatus	9	0	–	–	–
Micrurus bocourti	3	3	–	–	–
Micrurus dumerilii	47	47	38	9	0.254
Micrurus filiformis	8	0	–	–	–
Micrurus hemprichii	8	0	–	–	–
Micrurus langsdorffi	31	19	15	4	1.378
Micrurus lemniscatus	71	46	37	9	1.986
Micrurus margaritiferus	6	2	–	–	–
Micrurus mipartitus	120	97	78	19	1.929
Micrurus ornatissimus	13	11	9	2	1.918
Micrurus spixii	65	49	39	10	1.745
Micrurus steindachneri	15	8	6	2	1.739
Micrurus surinamensis	33	24	19	5	1.684
Micrurus tschudii	8	8	–	–	–
Porthidium nasutum	105	64	51	13	1.660

Notes:

T Occ = total occurrences, F Occ = filtered occurrences, N Calib = calibration occurrences, N Test = evaluation occurrences, AUC r = AUC ratio of Partial ROC (the absence of a value in this column indicates that the species was not modelled).

We modelled the potential distribution of 19 snake species (Table 1), including only species for which at least five occurrences were available to avoid statistical difficulties associated with small sample sizes [8]. We are confident that excluding these species did not affect our analysis because they have characteristics (e.g. rareness, restricted distributions, low abundance) that would greatly reduce the frequency of bites [5]. Besides, the majority of snakebites in Ecuador are caused by only six species [4,9,10].

2.2. Environmental layers

We used 19 bioclimatic layers (resolution = 30 arc-seconds) from World Clim (http://www.worldclim.org/) [11] to characterize environmental conditions. We masked these layers to the extent of a polygon that represents a hypothetical area of historical accessibility for each species (area M; sensu Soberón and Peterson [12]). We used a vectorial map of terrestrial ecoregions generated by the WWF [13] to determine the M area.

Moreover, we reduced dimensionality and collinearity among bioclimatic layers by applying a principal components analysis (PCA). To do this, we used the “PCARaster” function in the R [14] package ENMGadgets [15]. We retained for the modelling stage the first five components that explained ≥95% of the overall variance.

2.3. Ecological niche modelling

We used Maxent 3.3.3k [16] to estimate environmental suitability and the potential distribution of each species. We chose this program because it is reputed to have better predictive ability than other presence-only data algorithms [17]. We used calibration data across the accessible area (M) for each species with the Maxent bootstrapping/replicated settings and a “cloglog” output. We also carried out ten replicate analyses to explore effects of specific calibration data sets on model outputs. Finally, to obtain the potential distribution maps we thresholded the average environmental suitability models utilizing the “10 percentile training presence” rule to exclude problematic occurrences (e.g. taxonomic misidentifications, presences from sink populations and imprecisions in geo-referencing) [17]. In all the maps, we keep suitability values only within the potential distribution area (Figure 1).

Figure 1. Ecological niche models for the two most dangerous snakes in Ecuador (Bothrops asper and B. atrox). Photo credits: B. asper (Jerry Oldenettel) and B. atrox (Andreas Schlüter). Source: Authors.

2.4. Model evaluation

We assessed predictive performance with the evaluation data-set via the partial ROC (receiver operating characteristic) approach. This technique is based on the traditional ROC [18], but takes into account the area of coverage of the commission error axis by model predictions, and gives preference to omission over commission error in evaluating model strength [19]. We calculated the AUC (area under the curve) ratio for each species using the software Tool for Partial-ROC [20] with an allowed omission of 10% of validation data. AUC ratio values above one indicate that models outperform the null expectation [19].

2.5. Snakebite risk maps

Following Yañez-Arenas et al. [5], we built a snakebite risk map that approximates the probability of being bitten by a venomous snake at a given cell. We estimated snakebite risk by stacking ecological niche models and weighing each species by a number related to the proportion of bites caused by it. Specifically, we summed the environmental suitability values of the cells that were classified as potential distribution for each species. Then, we multiplied these values by the reported proportion of snakebites caused for each species (Table 2) in epidemiological studies carried out in Ecuador [21–23].

Table 2. Species values that weights their probability to cause human envenomation.

Species	Weights based on incidence data
Bothrops asper	80.0
Bothrops atrox	58.0
Bothrops bilineatus	36.0
Bothrops taeniatus	5.0
Bothriechis schlegelii	5.0
Lachesis muta	5.0
Bothrocophias hyoprora	5.0
Porthidium nasutum	5.0
Bothrops pulchra	5.0
Micrurus mipartitus	0.1
Micrurus spixii	0.1
Micrurus ancoralis	0.1
Micrurus lemniscatus	0.1
Micrurus surinamensis	0.1
Micrurus dumerilii	0.1
Micrurus ornatissimus	0.1
Micrurus steindachneri	0.1
Micrurus langsdorffi	0.1
Micrurus narduccii	0.1

Finally, we identified vulnerable rural communities (rural localities with high population density and high snakebite risk) based on the “estimates of rural population density to 2015” global raster [24], and a shapefile of localities obtained from “Geodescargas” (http://www.geoportaligm.gob.ec/; elaborated on January 2013, accessed on June 2015).

3. Results

In all, we gathered 1187 spatially filtered presence records, varying among species from 1 (Micrurus annellatus, M. hemprichii) to 1482 (Bothrops asper). AUC ratio values for almost all species were >1 (with the exception of Micrurus dumerilii), indicating that models outperform null expectations (Table 1).

Modelled snakebite risk showed a broad geographic pattern, wherein the central region (Andes) of Ecuador has low to medium snakebite risk, while the western (Coast) and eastern (Amazonia) regions have in general medium to high risk (Figure 2). Because our snakebite risk index was produced giving greater weight to B. asper and B. atrox, the regions of the map of greatest risk (Figure 2) coincide with the areas of high environmental suitability for these species (Figure 1).

Figure 2. Snakebite risk map for Ecuador. Source: Authors.

We identified 187 rural communities (Table S1) that we consider vulnerable because they are densely populated and have high environmental suitability for the snake species (B. asper and B. atrox) responsible for the majority of accidents in the country. Of these communities, 75% are located in the western region of Ecuador, 14% in the central region, and 11% in the eastern region. In the latter, our model estimates an extensive area with high snakebite risk values, however, there are very few densely populated rural communities (Figure 3). Finally, three locations in the eastern region stand out as having a much higher estimated risk than the rest: “Palora (Metzera)”, “Sangay” and “Shell”.

Figure 3. Categories of snakebite risk for densely populated rural communities in Ecuador. Source: Authors.

4. Discussion

Our study estimated snakebite risk at fine resolution (~1 km²) based on detailed environmental suitability/distribution models for 19 species. A good predictive performance was observed in potential distribution models for almost all venomous snakes, including those of medical importance such as Bothrops asper, B. atrox, B. brazili, Bothrocophias microphthalmus, Bothriopsis bilineata, B. taeniata, Lachesis muta [4,22,23], except Micrurus dumerilii.

This research represents the first of its type in Ecuador, in terms of geographic dimension and scope. Previous studies were carried out at local scales (e.g. Smalligan [22], Theakston et al. [23], Praba-Egge et al. [25]) and the more comprehensive ones do not provide enough spatial details (e.g. González-Andrade and Chippaux [2]).

We predicted that high snakebite risk is associated with tropical and sub-tropical lowlands of Ecuador; a result that is consistent with previous findings [21]. Also, González-Andrade and Chippaux [2] reported that ~56% of snakebites take place in the western region, followed by the eastern (~33%) and central region (~11%), which gives partial support to our analysis of vulnerable human communities. We found more vulnerable communities in the western region followed by the central region and finally the eastern region. Although in the latter our model estimates an extensive area with high snakebite risk, there are very few densely populated rural communities located in this region compared with the other two. Of the 565 densely populated rural communities that we identified, 54% are located in the western region, 41% in the central region and 5% in the eastern region.

A major limitation of the present study is that we were unable to obtain detailed (e.g. at location level) incidence data to validate our snakebite risk map, and our model does not incorporate some socio-demographic factors (e.g. human activities, type of housing, population’s age) that are known to be causal of snakebites [1–3]. However, previous evidence has shown that ecological niche modelling is usually successful at inferring geographic patterns of ophidism at continental [5], country [26] and regional scales [7]. This is particularly useful to infer potential snakebite risk areas when hospital data are scarce or biased. For instance, in many rural communities around the world people prefer to attend traditional doctors instead of a modern clinic [27]. We recognized three rural communities in which health care providers should pay special attention: “Palora (Metzera)”, “Sangay” and “Shell”.

This is a preliminary analysis that represents a significant progress towards understanding geographic patterns of ophidism in Ecuador. Nonetheless, much still remains to be done: First, it is necessary to continue digitalizing occurrences of venomous snakes. Second, it is crucial to generate a snakebite cases reporting system that concentrates data per locality for the whole country. Third, it is important to perform studies to evaluate which factors contribute the most to snakebite risk. And four, anthropogenic climate change should be considered to anticipate strategies for delivering/distribution antidotes and to make people aware and conscious in regions where snakebite risk might increase in the future.

Author contribution

Conceived idea and designed the analyses: Carlos Yañez-Arenas. Analysed the data: Carlos Yañez-Arenas, Luis Díaz-Gamboa, Carlos Patrón-Rivero, Kevin López-Reyes. Wrote the paper: Carlos Yañez-Arenas, Xavier Chiappa-Carrara.

Disclosure statement

No potential conflict of interest was reported by the authors.

Funding

This work was supported by Consejo Nacional de Ciencia y Tecnología [grant number 64636].

Supplemental data

Supplemental data for this article can be accessed here [https://doi.org/10.1080/23766808.2018.1454762].

Acknowledgements

We thank Karla Rodríguez Medina and Eduardo Garza for their review, comments and suggestions to this manuscript.

Associate Editor: Yntze van der Hoek.