Snakebites are associated with poverty, weather fluctuations, and El Niño.

Snakebites are environmental and occupational health hazards that mainly affect rural populations worldwide. The ectothermic nature of snakes raises the issue of how climate change's impact on snake ecology could influence the incidence of snakebites in humans in ways that echo the increased predation pressure of snakes on their prey. We thus ask whether snakebites reported in Costa Rica from 2005 to 2013 were associated with meteorological fluctuations. We emphasize El Niño Southern Oscillation (ENSO), a climatic phenomenon associated with cycles of other neglected tropical diseases (NTDs) in the region and elsewhere. We ask how spatial heterogeneity in snakebites and poverty are associated, given the importance of the latter for NTDs. We found that periodicity in snakebites reflects snake reproductive phenology and is associated with ENSO. Snakebites are more likely to occur at high temperatures and may be significantly reduced after the rainy season. Nevertheless, snakebites cluster in Costa Rican areas with the heaviest rainfall, increase with poverty indicators, and decrease with altitude. Altogether, our results suggest that snakebites might vary as a result of climate change.

INTRODUCTION

Although the technology for antivenom manufacturing is now available in the public domain and, in principle, most deaths and sequelae attributable to snakebite envenoming are avoidable (, ), snakebite envenoming remains one of the leading causes of deaths due to neglected tropical diseases (NTDs) (). With an estimated 2.5 million people affected by snakebites annually (of whom 400,000 will develop permanent sequelae and 85,000 will die) (), without proper estimation of snakebite burden in terms of disability-adjusted life years (), and with the exclusion of snakebites from packages for the control of NTDs (, ), it can be argued that snakebite envenomation is a neglected tropical disease. Moreover, it is likely that the actual magnitude of snakebite envenoming is much higher than hospital-based estimates, as judged by community- and household-based surveys (). From an ecological point of view, snakebites mainly occur in agricultural niche construction by humans (, ), where snakebites reflect the interaction between humans and snakes that emerges from the invasion, disruption, or destruction of snakes’ natural habitats by humans (). This interaction closely resembles competition, because humans and venomous snakes are harmful to each other (); yet, in an evolutionary context, venoms evolved to improve the fitness of snakes as predators of species other than humans (). As top predators, snakes play a major role in ecosystem regulation (), and some species of medical importance () have been resilient to global changes in biodiversity.

Snakes are ectothermic organisms whose distribution, movement, foraging patterns (), and life history strategies () change as a function of weather fluctuations. With warming, increased predation by some snakes has been reported (), supporting the notion that climate variability could influence snakebites, a related antagonistic interaction. Independently of changes in snake biology, the increased vulnerability of destitute populations exposed to snakes because of a lack of appropriate shelter in the face of extreme weather events () can increase the occurrence of snakebites, underscoring the relevance of socioeconomic factors to understanding changes in venomous snakebite patterns associated with weather, climate change, and variability, thus raising the need for understanding the interplay of both factors given the observations and expectations of more frequent extreme weather events with climate change (). Costa Rica (CR) is an ideal setting for studying snakebites because most of the population affected by snakebites seeks treatment in health posts (), where free treatment is provided and snakebite reporting is mandatory. This wealth of high-quality data contrasts with data records from other regions of the world where underreporting is associated with poorly developed health systems (, ). CR is located in an area where El Niño Southern Oscillation (ENSO) has marked effects on weather patterns (). Most snakebites in CR are attributable to the terciopelo Bothrops asper (Fig. 1A), and terciopelos and their bites are widespread across the neotropics (). All these factors make CR an ideal setting for determining whether snakebite dynamics are influenced by changing weather patterns and for understanding how poverty might influence the spatial distribution of snakebites given the high degree of socioeconomic inequity across the country ().

Fig. 1

Snakes and snakebites in CR.

(A) The terciopelo B. asper. (B) Average annual snakebite incidence, by canton, from 2005 to 2013. County color indicates snakebite incidence rate, county boundary color indicates relative risk, and a marking described in the map legend indicates the primary cluster.

RESULTS

A total of 6424 snakebites were reported in CR from 2005 to 2013. The 9-year average incidence rate was 15.24 per 100,000, ranging from 10.63 per 100,000 to 22.98 per 100,000 when the entire population was assumed to be at risk. Nevertheless, those figures underestimate the incidence rate in the at-risk population (mainly rural), with the average rate jumping to 41.27 per 100,000, ranging from 30.53 per 100,000 to 58.94 per 100,000 with a steadily decreasing at-risk population (fig. S1). Average snakebite incidence rates for each canton during the study period are shown in Fig. 1B. The highest incidence of snakebites occurred in southern CR and in the northern canton of La Cruz, bordering Nicaragua. Snakebites usually occurred in suburban or rural regions, as revealed by a nonlinear negative correlation with population density (fig. S2). Figure 1B shows that snakebites are spatiotemporally clustered in the Caribbean coast and southern CR, a pattern that shows low variability through time (fig. S3).

Figure 2 shows the spatial distribution of factors that might be associated with snakebites and their spatial relationship stratified by geographically weighted regression (GWR). We considered elevation (Fig. 2A) as a proxy for temperature, which linearly decreases with altitude () and rainfall (Fig. 2B), both being major environmental factors associated with terciopelo abundance and activity (). Poverty gap index (Fig. 2C) and percentage of destitute housing (Fig. 2D) were analyzed to account for socioeconomic factors that might determine the occurrence of snakebites, as observed for other NTDs in CR (). Results from GWR showed a negative association between average canton elevation and the cumulative number of snakebites (Fig. 2A). Rainfall demonstrated different impacts on snakebites. In wetter regions (Southern Pacific basin), higher precipitation was negatively associated with snakebites (Fig. 2B). In contrast, high precipitation in drier regions (Northern Pacific basin) increased the risk of snakebites (Fig. 2B). Overall, a high poverty gap index was predictive of snakebites (Fig. 2C), and a similar pattern was observed for the percentage of destitute housing, especially in rural areas (Fig. 2D). GWR was highly successful in explaining snakebite spatial patterns because 82% of the deviance was explained by the model (local deviance explanation ranged from 0.47 to 0.88), with the best performance observed in the Pacific basin of CR (fig. S2 and table S1), where snakebite hotspots occur (Fig. 1B and fig. S4).

Fig. 2

Variables spatially associated with snakebites in CR.

(A) Altitude. (B) Rainfall. (C) Poverty gap index. (D) Destitute housing. Coefficients are shown (in the legend of each panel) only when pseudo t values are significant (P < 0.05).

Figure 3 shows temporal patterns of snakebite incidence in CR. Figure 3A depicts the monthly snakebite incidence during the study period; most peaks coincide with the hot and cold phases of ENSO in CR. Seasonal patterns (Fig. 3B) show low monthly variability. Nevertheless, the median snakebite incidence was highest in June and July, with a secondary smaller peak in November. Snakebite incidence was significantly autocorrelated (that is, temporally associated), with snakebites in previous months occurring at lags of 1, 3, 4, 7, 8, and 14 months, as revealed by an autocorrelation function (ACF) (Fig. 3C); had a significantly negative cross-correlation [cross-correlation function (CCF)] with rainfall 11 months before snakebites occurred (Fig. 3D); and had a significantly positive cross-correlation with temperature 1 month before snakebites occurred (Fig. 3E). A cross-wavelet coherency analysis (Fig. 3F) showed that snakebite incidence and sea surface temperature 4 (SST4) were associated nonstationarily (that is, nonconstantly) at interannual scales, with cycles of about 18-month periods in snakebites being synchronous with ENSO activity between 2009 and 2012, a pattern also shown by raw snakebite incidence data (Fig. 3A). On the basis of information from ACFs and CCFs and after a rigorous process of model selection (table S2), we fitted the following model:where N is the logarithm of monthly snakebite incidence, T stands for temperature, R stands for rainfall, t indicates time, and ε ~ N(0, σ2). Parameter estimates are presented in Table 1. Parameters indicate that, on average (exp(μ)), 4.30 snakebites per 100,000 people per month were recorded. Snakebite incidence increased by 24% for every unit increase in temperature (in degrees Celsius) above the average temperature (exp(α)) and decreased by 1% for every 7-mm increase in rain above the average rainfall (exp(β)).

Fig. 3

Temporal snakebite incidence patterns in CR.

(A) Monthly time series from 2005 to 2013; colors indicate the phases of ENSO as explained in the legend (inset). (B) Seasonality in snakebite incidence; monthly box plot shows the log-transformed number of snakebites, and colors indicate the different phases of ENSO. (C) ACF of monthly snakebites. (D) CCF between snakebites and temperature. (E) CCF between snakebites and rainfall. (C to E) Dashed lines represent the 95% confidence interval for correlations expected to arise randomly. (F) Cross-wavelet coherence analysis of snakebites and ENSO. We used SST4 as an ENSO index. In the analysis, a 6-month smoothing window was used. The cross-wavelet coherence scale is from 0 (blue) to 1 (red). Red regions in the plots indicate frequencies and times for which the two series share variability (or power). The cone of influence (in which results are not influenced by the edges of data) and the significantly coherent time-frequency regions (P < 0.05) are indicated by solid black lines.

Table 1

Parameter estimates for the best model explaining the monthly snakebite incidence in Costa Rica.

Parameter		Estimate ± SE
μ^	Average logarithm of monthly snakebites	1.46 ± 0.06
ϕ^_1	Autoregressive component	0.24 ± 0.09
ϕ^_7	Seasonal autoregressive with a 7-month lag	0.32 ± 0.10
α^	Temperature with 1-month lag	0.21 ± 0.07
β^	Rainfall with 11-month lag	−0.0014 ± 0.0004
σ^_ε	SD of the error	0.44

DISCUSSION

Our analysis shows that snakebites are associated with changes in temperature and rainfall across time, and that unusually high numbers of snakebites occur during the cold and hot phases of ENSO. Spatially, we found that snakebites cluster in the most humid lowland areas of CR, where terciopelos are distributed, and are more frequent in the poorest areas of regions with similar weather patterns. This combination of patterns highlights the fact that snakebites follow meteorological changes, and these patterns most likely reflect the impact of meteorological fluctuations on snake biology. For example, the 7-month autoregressive lag in the number of snakebites in our model coincides with the reproductive phenology of terciopelos (); this can be translated as a lag in the recruitment of juvenile terciopelos, which more frequently bite humans engaged in agricultural activities (). Moreover, terciopelos are known to increase their foraging activity with rising temperatures (, ), which might be reflected in the increased number of snakebites in higher-than-average temperatures associated with the hot phase of ENSO () and the increased number of snakebites at low altitudes, where temperatures and snake abundance are higher. Spatially, the positive association of snakebites with rainfall in the Southern Pacific basin of CR might reflect the increased abundance of terciopelos in humid environments, where the abundance of prey and resources allows for the maintenance of higher snake densities (). Temporally, the delayed negative association with rainfall might represent cascading effects of productivity in the environment (). Basically, a decrease in rainfall associated with ENSO () can increase primary productivity in tropical areas (), which is followed by an increase in the number of organisms that can serve as prey for snakes (for example, rodents) (). Rodents and other prey that benefit from increased productivity can increase the likelihood of a delayed and transient interaction between humans and snakes in rural areas, because snakes will move into contact with humans as the density of their prey crashes once rainfall and productivity return to normal conditions (, , , ). This duality can explain the positive association of snakebites with the hot and cold phases of ENSO. During the hot phase, the incidence of snakebites might rise as a result of increased snake activity at higher temperatures, whereas during the cold phase, snakebites might increase because of abrupt changes in the abundance of prey (, ). During the cold phase, a reduction in seed productivity and its effects on prey availability can force snake movement into foraging areas (), where they come into contact with humans ().

Our results highlight the fact that snakebites occur more frequently in poor settings, following a common pattern for other tropical diseases in the region (, –) and reflecting the general vulnerability of impoverished human populations to the adverse effects of climate change () and neglected diseases (). The latter is a pattern that might be extrapolated to other areas where snakebites are a major health problem. For example, the incidence of snakebites in the northern county of La Cruz, where terciopelos are scarce, is high (). This might highlight the occurrence of snakebites by snakes other than terciopelos (for example, rattlesnakes, which are common in arid areas and occasionally bite humans) () or might even reflect temporal antivenom shortages in neighboring Nicaragua, hence prompting across-the-border movement of patients seeking medical attention.

Finally, our findings highlight the need for increased research on the ecoepidemiology of snakebites, an NTD that, as our results robustly suggest, should be included in the list of diseases or health hazards that are sensitive to environmental changes.

MATERIALS AND METHODS

Data

Using records from the Costa Rican Ministry of Health, we built a database of snakebites in CR, stratified by month and by canton (the geopolitical unit below the province level), from January 2005 to December 2013. Data were collected by health posts and clinics in CR that are under the administration of the Costa Rican Social Security Trust. These public centers provide free treatment of snakebites, and more than 99% of snakebites are diagnosed in health centers based on clinical manifestations. Snakebites inflicted by species of the family Viperidae are characterized, at the clinical level, by local manifestations (that is, pain, edema, hemorrhage, and tissue damage); systemic manifestations are predominantly characterized by coagulation disturbances and bleeding. In the case of coral snakes (family Elapidae), envenoming is characterized by neurotoxic manifestations such as local paresthesias and flaccid paralysis of various muscles ().

For spatial analysis, we considered for each canton (i) its average elevation; (ii) its annual average precipitation; (iii) its poverty gap index, which quantifies the percentage of houses with income below the poverty line (); and (iv) its percentage of destitute housing (that is, lacking services and made of inadequate materials) (). The average elevation for each canton was estimated using a 90-m × 90-m Global Land Survey Digital Elevation Model (www.landcover.org). Elevation was chosen as a covariate in our models because it is an accurate proxy for temperature, which is associated with snake activity. Annual average precipitation for each canton was calculated using a tropical rainfall–measuring mission. Satellite images (3B42 version 7) were acquired from the Goddard Earth Science Data and Information Services Center (daac.gsfc.nasa.gov). The annual average precipitation in CR throughout 2005–2013 was summarized by canton.

For temporal analysis, we built temperature and rainfall time series with station data and used SST4 to quantify the impact of ENSO on snakebites. We built monthly time series for rainfall and temperature using freely available data (http://ncdc.noaa.gov/cdo-web/). For rainfall, we used data from the Limón and Liberia weather stations, whereas for temperature, we used data from the Limón, Liberia, Palmares, and Alajuela weather stations. For the analysis, all climatic covariates were demeaned. SST4 (aka the El Niño 4 index), which is based on temperature records of an area delimited by the Pacific Ocean (from 5°N to 5°S and from 160°E to 150°W), was used as an index for ENSO, which was obtained from the United States National Oceanic and Atmospheric Administration Climate Prediction Center (http://cpc.ncep.noaa.gov/data/indices/ersst3b.nino.mth.81-10.ascii). ENSO phases were designated on the basis of the time series for the last 30-year anomalies of SST4. When these anomalies were below −2 or above 2, ENSO was considered as going through its cold and hot phases, respectively; ENSO was considered as going through its normal phase when values fell in the interval [−2,2]. On the basis of results from epidemiological studies of snakebites in CR, we assumed that only rural populations were at risk (), and we used county and country data on rural populations () as denominators to estimate the snakebite incidence rate in each county for spatial analysis and in each country for temporal analysis. However, for urban cantons (that is, cantons with more than 80% of their population living in urban centers), we considered all individuals to be at risk of snakebites given the uncertainty of at-risk populations in urban areas, where conservative estimates can be obtained by assuming that everybody is at risk because denominators were either zero or so small that they led to unrealistically high snakebite incidence rates in urban cantons.

Spatial analysis

Spatial clusters of snakebites were detected with the spatial scan statistic (), which relies on creating an infinite number of flexible windows of different sizes and locations. The upper limit of window size was set as reaching 50% of the at-risk population. We assumed that snakebites followed a Poisson distribution, and we considered elliptic scanning windows because CR is longer than it is wide. P values were obtained with 999 Monte Carlo randomizations. To explore the nonstationary impact of spatially explicit climatic and socioeconomic variables on snakebite incidence in CR, we used GWR modeling (). GWR models can be described by the following general equation:where y denotes the cumulative incidence rate of snakebites from 2005 to 2013 at canton i; u and v are centroid coordinates for canton i; x indicates a set of k independent variables (where j = 1,…,k); and β0 and β are coefficient estimates for each location. Snakebites were assumed to follow a Poisson distribution in the GWR model, and population size was used as an offset in the model to account for heterogeneities in population density across the cantons. An adaptive bisquare function was selected as the kernel function in GWR, using a bandwidth selected through minimization of the Akaike information criterion (AIC) (). Globally and locally (that is, for the entire CR and its cantons), we generated estimates of percentages of deviance explained to evaluate the performance of GWR (). For statistical inference, we used pseudo t values ().

Time series analysis

Seasonality of snakebite incidence, defined here as the number of snakebites for 100,000 people at risk, was assessed through box diagrams for each month (). We examined the autocorrelation structure of the snakebite incidence time series by inspecting ACF and partial ACF (PACF) and built a seasonal autoregressive (SAR) autonomous model (), which was used to prewhiten the time series of temperature, rainfall, and SST4. The prewhitened residuals and the residuals from the best autonomous model were used to estimate CCFs. ACF is a graphical summary of the correlation of a time series with itself through several time lags, whereas PACF is a similar function built by only considering consecutive time lags. Prewhitening is a process that “filters” out any correlation structure in an ancillary time series that is similar to that of a focal time series. CCFs are functions depicting the correlation between two time series at different time lags. We fitted a full nonautonomous model that had the autocorrelation structure of the best autonomous model and considered climatic covariates at the significant lags identified in CCFs. This full model was simplified by the process of backward elimination, where covariate exclusion is guided by minimization of the AIC (). The AIC () is a tradeoff function between the number of parameters and goodness of fit. We used cross-wavelet analysis to study the nonstationary (that is, changing through time) association () between SST4 and snakebites to test the hypothesis that interannual variability in snakebite incidence is associated with ENSO. Briefly, a cross-wavelet analysis is composed of a cross-wavelet coherence analysis that determines whether two time series share similar oscillations through time ().

Model diagnostics

Results from the GWR model showed that when standardized residuals were checked, no trend or specific patterns were observed, and Moran’s I test (z-score = −1.32, P = 0.19) indicated no spatial autocorrelation among residuals. Results from the SAR model showed that assumptions of temporal independence and normality were followed by the residuals, ensuring a sound inference.