Integrating snake distribution, abundance and expert-derived behavioural traits predicts snakebite risk

1. Despite important implications for human health, distribution, abundance and be - haviour of most medically relevant snakes remain poorly understood. Such data deficiencies hamper efforts to characterise the causal pathways of snakebite en - venoming and to prioritise management options in the areas at greatest risk. 2. We estimated the spatial patterns of abundance of seven medically relevant snake species from Sri Lanka, a snakebite hotspot, and combined them with indices of species' relative abundance, aggressiveness and envenoming severity obtained from an expert opinion survey, to test whether these fundamental ecological traits could explain spatial patterns of snakebite and envenoming incidence. The spatial intensity of snake occurrence records in relation to independent environ - mental factors (fundamental niches and land cover) was analysed with point pro - cess models. Then, with the estimated patterns of abundance, we tested which species' abundances added together, with and without weightings for aggressive - ness, envenoming severity and relative abundance, best correlate with per capita geographic incidence patterns of snakebite and envenoming. 3. We found that weighting abundance patterns by species' traits increased correla - tion with incidence. The best performing combination had three species weighted by aggressiveness and abundance, with a correlation of r = 0.47 ( p < 0.01) with snakebite incidence. An envenoming severity and relative abundance-weighted

Regional variability of snakebite and envenoming burden within and among nations have been attributed to human factors like poverty and occupation (e.g.agriculture; Harrison et al., 2009).However, despite the logical relevance, propensity to bite, frequency of venom injection (Pucca et al., 2020), distribution and abundance of venomous snakes in determining snakebite and envenoming burden has received less attention (Murray et al., 2020).Incorporating these factors into snakebite burden mapping could improve the extent to which development of preventive and therapeutic interventions can target areas or periods of highest risk and need.
Previous analyses show that there are (a) socio-economic (e.g.poverty, % of people working in agriculture; Harrison et al., 2009), (b) environmental (e.g.rainfall, temperature and vegetation greenness; Ediriweera et al., 2016;Molesworth et al., 2003) and (c) ecological (e.g.potential venomous snake abundance and richness; Bravo-Vega et al., 2019;Yañez-Arenas et al., 2014, 2016) factors driving geographical snakebite patterns.Abundance and richness of venomous snakes are relevant explanations of snakebite risk because snakebites are human-wildlife contacts (Glaudas, 2021), hence the number of encounters should be a function of human and snake abundance (Bravo-Vega et al., 2019).
Whether human-snake encounters result in bites and bites result in envenoming depends on traits like snake behaviour, venom injection and its toxicity to humans.Some species exhibit a greater propensity to strike or deliver venom after a bite (Healy et al., 2019;Pucca et al., 2020).The chance of fatal envenoming is also affected by venom composition and its systemic effects on the human body (Healy et al., 2019;Pucca et al., 2020).A logical pathway to better understand the relationship between snake species and snakebites is therefore to map the spatial variation of bite risk as a function of species' abundances weighted by their aggressiveness and the severity of the envenoming illness they produce.For instance, highly abundant but unaggressive species should cause less disease than less abundant but highly aggressive and venomous snakes.Such knowledge is still difficult to find in primary scientific literature, and have to be inferred from surrogates and/or expert knowledge.
Field-collected abundance data are missing for most snakes (Luiselli et al., 2020), but estimating abundance patterns across geographic ranges is possible, in theory, by characterising the land cover, climate or topography of the locations where species occur (Martínez-Meyer et al., 2013).Recent work shows that many organisms are more abundant in locations with optimal environmental conditions (Osorio-Olvera et al., 2020).Also, snakes are ectotherms, hence the effect of climate on the geography of population sizes should be stronger than for endotherms (Kearney et al., 2018).
Measures of aggressiveness and envenoming severity for some medically important snakes are available in scholarly publications (Gerardo et al., 2019;Tasoulis & Isbister, 2017) but for many they remain obscure.However, experienced field herpetologists know species' propensity to bite when encountered (Gibbons & Dorcas, 2002), or how abundant species are in relation to each other (Steen, 2010).
Likewise, medical practitioners treating snakebite victims know the severity of the illness produced by various snake venoms.Therefore, exploring the contribution of these traits to explain snakebite and envenoming patterns constitutes another important step towards incorporating causality in snakebite risk analyses.Furthermore, such data are invaluable as snakebite research is still lacking basic data (Longbottom et al., 2018).
To this end, we estimated the spatial patterns of abundance of seven venomous snakes of medical relevance.Using expert-derived data, we explore its contribution towards risk in the island nation of Sri Lanka, a global snakebite hotspot (Kasturiratne et al., 2005), with the highest resolution estimates of snakebite and envenoming incidence yet available (Ediriweera et al., 2016) to explore: (a) can occurrence-based distribution models of snakes broadly predict snakebite incidence patterns?;(b) can expert-derived data improve these predictions?; and (c) are there specific combinations of species and traits that better explain snakebite and envenoming incidence?combination of two species was the most strongly associated with envenoming incidence (r = 0.46, p = 0).

Synthesis and applications.
We show that snakebite risk is explained by abundance, aggressiveness and envenoming severity of the snake species most frequently involved in envenoming cases.Incorporating causality via ecological information of key snake species is critical for snakebite risk mapping, helping to tailor preventive measures for dominant snake species and deploying the necessary antivenom therapies.

K E Y W O R D S
abundance, aggressiveness, fundamental niches, neglected tropical diseases, occurrence data, point process models, snakebite envenoming

| Overview
Estimating snake abundance patterns and testing its association with incidence comprised: (a) identifying study species, collecting occurrence records from private and public repositories (GBIF and VertNet), testing for spatial autocorrelation and filtering (1-10 km) and measuring the effect of observation bias on the environmental conditions of occurrence records; (b) selecting and reducing climatic dimensions via an approximation to snake species' fundamental niches (FNs), developing land cover data and its derivatives (proportion of agricultural land, distance to forests and tree cover) and compiling topographic data (digital elevation model and topographic wetness index); (c) fitting, validating, correcting observational bias (Warton et al., 2013) and selecting point process models (PPMs); (d) collecting expert data on snake species' relative abundance, aggressiveness and envenoming severity to use as weights; and (e) testing the ability of different combinations of estimated snake abundances with and without weights to predict the estimates of snakebite and envenoming incidence.The various ecological and statistical concepts used in these analyses are described in greater detail below and in Figure 1.
The statistical method used were PPMs, which comprise regression tools for the analysis of point intensity (average number per unit) in a Cartesian plane as a function of environmental factors.PPMs have higher levels of control, customisation, statistical rigour, diagnostics and geostatistical extensions to model spatial dependence (Baddeley & Turner, 2005), many of which are lacking in mainstream methods for species distribution modelling (e.g.Maxent; Renner et al., 2015).The types of PPMs include Poisson-assume spatial independence between points (Baddeley et al., 2016)-Area-i nteraction-assume moderate levels of spatial dependence between points and estimate it with a random effect for point pairs and log-Gaussian Cox processes can model highly clustered points via a spatially-varying random effect (Baddeley & van Lieshout, 1995).

| Study species
Seven medically relevant, truly venomous, snake species in the Viperidae and Elapidae families in Sri Lanka were included in the study: spectacled (Indian) cobra Naja naja, common krait Bungarus caeruleus, Sri Lankan krait B. ceylonicus, Russel's viper Daboia russelii, saw scaled viper Echis carinatus, Sri Lankan green tree pit viper Trimeresurus trigonocephalus and the cryptic species complex of hump-nosed pit vipers (Hypnale hypnale, H. nepa and H. zara, hereafter Hypnale spp.).Occurrence records from published work, local and international museum databases, public repositories (GBIF and VertNet) and personal observations of local herpetologists were compiled (step 1 of Figure 1).

| Testing for autocorrelation and observational bias
To obtain the datasets to analyse as Poisson PPMs, we first applied a spatial filter to each species' raw data using distance thresholds of 0-10 km.The result was 20 different datasets for each species where the minimum distances between points were effectively 0.5-10 km.
Then we tested for spatial autocorrelation with K-envelopes with 39 iterations on each spatially filtered dataset.With the K-envelope tests, we identified the filtering distance that decreased spatial autocorrelation (test results are shown in Supporting Information-Autocorrelation tests).
To identify observation bias towards roads and cities and control it in the modelling process (Warton et al., 2013), we used the filtered datasets and a raster layer of minimum distance to roads and cities (https://data.humdata.org/).To do so, we compared the distance to roads and cities of occurrence records with the distance to roads and cities of random samples within a polygon enclosing the occurrences of each species.Then we ran two unpaired Wilcoxon tests, to see whether species occurrences were closer to roads and cities than expected by chance (Kadmon et al., 2004).To see if occurrences biased towards roads and cities were also environmentally biased, we compared random samples close and far from roads and cities.

| Estimating fundamental niches (FNs)
The FN is defined as the full set of abiotic conditions (e.g.temperature, humidity or rainfall) tolerated by a species (Soberón & Peterson, 2005), and is often estimated from geographical occurrence records (latitude, longitude of locations where organisms have been observed).Here we approximated FNs using occurrence records within and outside Sri Lanka by (a) extracting climatic conditions from the WorldClim v2.0 dataset (Fick & Hijmans, 2017), (b) summarising climatic preferences with minimum volume ellipsoids (Qiao et al., 2016) and (

| Land cover data
We produced a layer of five land cover classes (agriculture, forest, degraded forest, urban and tea) of Sri Lanka for the years 2004-2017 using 30 m atmospheric-corrected reflectance products with unsupervised ISODATA clustering with ENVI version 5.4.
Overall accuracy of the classified images was >95% (validated with 600 random samples).To analyse snake occurrences, we upscaled the final land cover layer by majority vote to the 1 km resolution of FNs.

F I G U R E 1 Diagram of the modelling process, from handling occurrence data to spatial association with incidence data
We further derived individual binary land cover classes, distance to closest forest patch and proportion of agricultural land to analyse snake datasets.Other land cover-related variables were as follows: proportion of land covered by trees (obtained from The Earth Observatory, https://earth obser vatory.nasa.gov/)and topographic wetness index generated from SRTM (Jarvis, 2008).All data were projected and analysed in the datum of Sri Lanka (SLD99, EPSG:5235).

| Point process models
The analysed data are collections of geo-referenced points of seven different snake species, hence a marked point process.Data were collected opportunistically over decades, which precludes its analysis as a marked point pattern to estimate interactions between species (Baddeley et al., 2016).Consequently, we compared three types of PPMs for estimating abundance patterns: Poisson, Area-interaction and log-Gaussian Cox process (LGCP; step 3, Figure 1).
To select a suitable model formula, we sought for lower AIC, correct computation of statistical effects and convergence of the likelihood function.We began with a list of 35 model formulas containing at least the FN layer and combinations of individual land cover classes, tree cover, proportion of agricultural land, distance to closest forest patch and distance to roads or cities, per snake species.To control observation bias, it was assumed to be represented by distance to roads and cities only, and to predict spatial trends, we set its value to zero.
Roads and cities could also be favoured environments; however, dedicated studies beyond the scope of our study are still needed.
After selecting the model formula, we fitted the Area-interaction and LGCP models.To select between Poisson, Area-interaction or LGCP, we analysed the spatially smoothed residuals and lurking plots (residuals along x and y coordinates), fitted the K-function to model smoothed residuals, simulating K-and L-envelopes with the predicted intensity patterns for each model type (Baddeley et al., 2016) and measuring the models' capacity to predict spatially partitioned datasets using the partial ROC test (Peterson et al., 2008; step 3.2, Figure 1).The partial ROC test measures the area under a curve formed by the proportion of predicted points versus the proportion of the area of a raster surface at increasing threshold values, and then divided by the random expectation (Supporting Information-Validation statistics, Methods).

| Association of snake species with land cover
We tested if species were more likely to be found in certain land cover classes by comparing the observed and expected frequency of observation per class.We used the 1-km-filtered datasets and a land cover layer of the areas without land-use changes between 2004 and 2017.Observed frequencies were determined by counting the occurrence records per class and species, and the expected frequencies were equivalent to the proportion of area covered by each land cover class.Given that observational bias was significant (methodological overview and PPMs), we divided the number of observed records per land cover class by the average number of records of all species per class to reduce the greater representation of urban areas.The resulting observed and expected frequencies were analysed with a Chi-squared test, estimating p-values with 10,000 Monte Carlo realisations of the samples as implemented in R's chisq.test(Hope, 1968; R-Development-Team, 2019).

| Snake ecological/biological traits and expertderived data
A questionnaire-based survey was conducted at a workshop attended by herpetological and medical experts in Sri Lanka to collect biological and clinical data on the study species.The index of aggressiveness was collated by agreement among three herpetologists, while the index of envenoming severity by three medical experts.
Both indices are ratings on a 0-10 scale (min-max for each trait, Figure 1, steps 1 and 4).Data on species' relative abundance and associations with different types of land cover were agreed upon a different group of nine herpetologists (Supporting Information-Questionnaire).Relative abundances were provided first as the estimates of population densities in 1-km 2 optimal habitat by two field herpetologists.Given the uncertainty, we transformed expert estimates to relative abundances.These indices (Table 3) were used as weights in the analyses of the next section.

| Association of incidence patterns with species combinations
We performed an exploratory analysis of snakebite and envenoming incidence with the estimated patterns of snake abundance.First, we generated all possible additive combinations of two to seven species with and without weighting (1920 total combinations) with each and all species' traits: estimated abundances summed and multiplied by the aggressiveness, envenoming severity and relative abundances.Then we tested the spatial association of snake combinations with snakebite and envenoming incidence measured as the fraction of the population estimated to be at risk (Ediriweera et al., 2016), with the modified t-test to account for spatial autocorrelation, equivalent to a Pearson correlation test (Clifford et al., 1989; Figure 1, step 5).
Adjusting for relative abundance was necessary as the estimated point intensities depend on the number of occurrence records.If these do not reflect relative abundances, point intensities will not reflect such characteristics of the snake assemblage.In addition, correcting for relative abundance with expert opinion partially addresses the issue of detectability differences among snake species.

| Occurrence data
We found significant spatial autocorrelation and observation bias, where most snake species were more likely to be recorded near

| Fundamental niches
The climatic variables selected for each species' FNs are shown in Table 1, along with the estimated centroids (Table 1).Mos species' FNs were characterised with three bioclimatic variables, but for Trimeresurus trigonocephalus, we used the first three principal components (Table 1; Figure 2).Centroids listed in Table 1 are those of the grey, outer ellipses in the leftmost column of Figure 2.These ellipses encompass 99% of the climatic conditions that snake species encounter across their distributions.The smaller orange ellipses inside the grey ellipses represent the climatic conditions that snakes encounter inside Sri Lanka.These ellipses show that the climatic conditions that occur in Sri Lanka are only a portion of what non-endemic species (Bungarus caeruleus, Echis carinatus, Daboia russelii and Naja naja) encounter elsewhere, and are in most cases fully contained inside the grey ellipses.

| Point process models
All species' models, except Echis carinatus, were very similar in all measures of goodness-of-fit (residuals, K-and L-envelopes and K-residuals; Supporting Information-Model fit tests), although Area-interaction models had smaller residuals than Poisson.
Fitted environmental effects of LGCP models were slightly different from Poisson, but intensity trends were identical (1:1 correspondence).Consequently, we only compared Area-interaction and Poisson PPMs, and based on smaller residuals we kept Areainteraction models for all species except Echis carinatus (Poisson; Table 2).
All selected models contained the FN variable (DNC spp -distance to niche centroid; Table 2), with a negative effect indicating that all species' abundance patterns are significantly affected by climate.The variable distance to roads and cities was also present in all final models.Given that we found significant observational bias by using distance to roads and cities as surrogate, it was set to zero for predicting point intensity.
Most models retained tree cover, distance to forests or proportion of agriculture (Table 2).The absence of land cover in all selected models indicates that there is no evidence of preference or avoidance of any specific classes.This coincides with the observed/expected ratios which rarely exceeded the 1 or −1 (preference/avoidance) thresholds (Figure 3, full details below).
The models of Echis carinatus, Hypnale spp., Naja naja and Trimeresurus trigonocephalus contained proportion of agricultural land and distance to closest forest patch.Also, all species' models except Bungarus caeruleus TA B L E 1 Venomous snake species reported in Sri Lanka, where families are Elapidae (E) and Viperidae (V).Global occurrences are GBIF, VertNet and expert-contributed records combined (in parentheses), and which occur inside and outside Sri Lanka.The division for Sri Lankan records are the full data points within Sri Lanka before (Total) and after filtering.The variables column lists the climatic factors selected for species, and the location of the multivariate centroid of the ellipsoids formed by the three variables for each species.For the meaning of the climatic axes, we refer the reader to Fick & Hijmans (2017) 266.4, 25.6, 94.1 266.4, 25.6, 94.1 F I G U R E 2 From left to right.Grey ellipses represent the species-wide niches, while orange represent occupied niches within Sri Lanka.
In the central panel is the predicted point intensity of each species within Sri Lanka (grey areas are those where the models predicted zero potential abundance).The rightmost panels show the Mahalanobis distance to the species-wide centroid of the ellipse in the left-side panels (values of the centroids are shown in Table 1)   (omission = 0.125, p = 0; Supporting Information-Validation statistics).

| Snake associations with land cover classes
Most of the study species occurrences were not randomly distributed among land cover classes (Figure 3).Bungarus caeruleus was more frequently found in urban areas and less so in agriculture, whereas Echis carinatus was more frequently reported in agriculture.Daboia russelii was detected more frequently in degraded forest and less in tea plantations, and Hypnale spp. was reported in tea plantations, less so in agriculture, urban areas and forests.Naja naja was detected more frequently in degraded forest than the rest of the land cover classes.
Bungarus ceylonicus and Trimeresurus trigonocephalus showed marginal association with forest and tea, and avoided urban and agriculture (Figure 3).Despite statistically significant results for most species, some ratios were close to 1 (random expectation threshold), suggesting that better data are needed to distinguish associations from random.

| D ISCUSS I ON
Snakebite is a major health hazard in low-income tropical countries (Harrison et al., 2009).Geographical estimates of venomous snake abundance provide a causal basis for estimating risk to populations, which we estimated for Sri Lanka using occurrence records of seven species.We integrated abundance models with snake aggressiveness and envenoming severity, and found the combinations of species and traits that best explain the estimates of snakebite and envenoming incidence (Ediriweera et al., 2016).In contrast with previous analyses of snake distributional ecology (Bravo-Vega et al., 2019;Yañez-Arenas et al., 2014), our results show that at a higher resolution, risk to local populations is driven by snake abundances, aggressiveness and venom toxicity.While these analyses successfully characterised a key role for snake ecology, links between snake traits and geographic variability of snakebite and envenoming risk were still lacking.
The spatial association of incidence patterns with snake species' traits is because venom transmission after a snakebite is akin to infectious disease transmission which is a mass action process (Bravo-Vega et al., 2019).Snakebite and envenoming illness are driven by species' abundance, aggressiveness and envenoming toxicity, among other factors like the propotion of dry bites (Pucca et al., 2020), activity overlap with humans (Goldstein et al., 2021), and species' tendency to seek shelter in houses (Ariaratnam et al., 2008).
Estimating these factors and their geographic variability should be prioritised to characterise snakebite fundamental epidemiology and further improve analyses.For instance, envenoming severity for D. russelii is known to vary geographically (Wüster, 1998).In this sense, our analyses parse out the explanatory value of aggressiveness and severity for snakebite and envenoming incidence at 3 km resolution.
The unexplained variability in our models may relate to the abovementioned snake factors and further human risk factors absent in the additive combinations, like socio-economic (Harrison et al., 2009), occupational (Mise et al., 2019) and cultural (Kumar, 2006)  Social and cultural factors also explain the epidemiologic relevance of some snake species.For instance, D. russelii bites mainly occur in the lower limbs of young men because they traditionally practice barefoot rice paddy farming (Kularatne, 2003).In contrast, Hypnale spp.bites tend to occur in the upper limbs of women, because Hypnale spp.'s primary habitat is leaf litter (Kularatne et al., ), and women tend to be in charge of home gardens.Finally, the majority of B. caeruleus bites occur among the poor while they are asleep in the floor (Kularatne, ).Including some of these factors in the analyses will help clarify each snake species' roles and find how these affect the geography of snakebite incidence.
In fact, land cover may be used to represent social characteristics at large scales (Ganzeboom et al., ); however, snake species' relationship to land cover is still uncertain, hence we require better data in this respect.In the light of existing limitations, characterising species' traits relevant to snakebite human risk factors is key to of women, because Hypnale spp.'s primary habitat is leaf litter (Kularatne et al., 2011), and women tend to be in charge of home gardens.Finally, the majority of B. caeruleus bites occur among the poor while they are asleep in the floor (Kularatne, 2002).Including some of these factors in the future analyses will help clarify each snake species' roles and find how these affect the geography of snakebite incidence.In fact, land cover may be used to represent social characteristics at large scales (Ganzeboom et al., 1992); however, snake species' relationship to land cover is still uncertain, hence we require better data in this respect.In the light of existing limitations, characterising species' traits relevant to snakebite human risk factors is key to implement adequate mitigation strategies such as mosquito nets to prevent B. caeruleus envenoming, promoting footwear and various degrees of agricultural mechanisation against D. russelii bites in the geographical areas where species are most abundant (Figure 2).
Combining publicly and privately available data with expert knowledge allowed us to overcome some data limitations to characterise snake abundance and biology as the driving factors of snakebite and envenoming incidence.Therefore, we provide a causal basis for the role of environmental conditions on snakebite burden.
Better understanding the human-snake-environment triad provides a novel avenue towards better eco-epidemiological understanding of snakebite and paves a path for anticipating its potential responses to global change: population growth, snake biodiversity loss, climate and land-use change.
c) calculating the Mahalanobis distance to the centroid of the ellipsoid (distance to niche centroid, DNC).Reported FNs were modelled with either a subset of three bioclimatic variables or with the first three components of a principal component analysis.The resulting FNs were then used to represent climate in the PPMs (see below) to estimate the geographic variability of species' abundances (Figure 1, steps 2.1-2.1.1).Resolution of the analysed data are approximately 1 km.Full details are given in Supporting Information-Methods.
Formulas of the models selected for each species.Tr = Proportion covered by trees; DF = Distance to closest forest patch; R = Distance to closest main road; Agr = Proportion of agricultural land; To = Topographic wetness index; DNC spp .= Distance to species' niche centroid.Estimates of variable effects appear in Supporting Information.'Type' indicates whether the model used for each species as Area-interaction (AI) or Poisson PPM 3.3.1 | PPM validationResiduals plots indicate bias was not widespread, and therefore the models adequately represented point intensity trends.Partial ROC tests indicate that models performed better than random on partitioned validation data for all species: Bungarus caeruleus 1.33 (omission = 0.05, F I G U R E 3 Land cover associations of the different snake species.y-axis shows the ratio − 1 of observed/expected presences within each land cover class.Positive ratios mean the species has been y amount of times more frequently observed in that land cover class than expected by chance.Negative ratios mean the opposite.Each species' panel contain the Chi-squared test results for reference Bungarus ceylonicus 1.74 (omission = 0.05, p = 0), Daboia selii 1.13 (omission = 0.025, p = 0.02), Echis carinatus 1.60 (omission = 0.125, p = 0), Hypnale spp.1.16 (omission = 0.15, p = 0), Naja naja 1.14 (omission = 0.1, p = 0.01) and Trimeresurus trigonocephalus 1.52 of Bungarus caeruleus, Daboia russelii and Trimeresurus trigonocephalus weighted by aggressiveness and adjusted for relative abundance had the highest, positive and significant spatial association with snakebite incidence (r = 0.47, df = 163, p = 0).The same additive combination of species without F I G U R E 4 Species' combination that maximised association with snakebite (left) and envenoming (right) incidence data.Insets in the top-right corners show scatterplots of analysed data Scores of aggressiveness, bite severity and overall hazard score of each venomous Sri Lankan snake species, curated by Sri Lankan herpetologists and medical experts.Relative abundances are mean and min-max of estimates adjusting for relative abundance had a lower association with snakebite incidence (r = 0.44, df = 111, p 0), and this was even lower without weighting for aggressiveness (r = 0.41, df = 183, p = 0).
data to exist, are an approximation to incidence variability.Nevertheless, snake species correlated with snakebite and envenoming incidence could be those more frequently involved in human-snake interactions.The different subsets of species for snakebite and envenoming indicate that snake biological factors are the primary drivers of spatial heterogeneity.The absence of important species from both subsets do not detract from their actual con-, N. naja and E. carinatus do not explain the east versus west envenoming pattern as would be expected is because the first two are almost equally abundant on both sides, and the latter is present in small areas, contributing few data points.The remaining, B. ceylonicus, rarely causes severe envenoming illness because it is unaggressive, and T. trigonocephalus is only mildly venomous.For these reasons incidence patterns are broadly driven by a few species at the 3 km scale.
(Ediriweera et al., 2017)d populations.Also, incidence data are the estimates of a national snakebite survey, which despite being the highest resolution It is possible then that abundance hotspots of B. caeruleus and D. russelii represent habitat surrogates for other non-or mildly venomous snakes frequently encountered (de Silva & Aloysius, 1983), but for which we lack occurrence or abundance data.Regarding envenoming, the severity-weighted combination of B. caeruleus and Hypnale spp.had more explanatory power than all other possible combinations and is consistent with the epidemiology of snakebites in Sri Lanka(Ediriweera et al., 2017).Envenoming incidence is higher in eastern Sri Lanka, where the highly venomous B. caeruleus is very abundant.Whereas in the western wet zone, Hypnale spp. is more abundant but causes less severe envenoming than B. caeruleus.The high abundance of both and greater envenoming severity of B. caeruleus than Hypnale spp.(Table3) explain why people seek formal medical treatment more often after a snakebite in dry eastern Sri Lanka(Ediriweera et al., 2017).The fact that D. russelii