Spatial Analysis and Seasonal Variation of Snakebites of the Genera Bothrops and Crotalus in the State of São Paulo

Abstract:

Snakebites are classified as a neglected tropical disease and are associated with poverty and climatic oscillations. Accidents caused by snakes can lead to death and cause serious sequelae. This study aims to analyze the spatial distribution patterns of snakebites caused by the genera Bothrops and Crotalus in the State of Sao Paulo, between 2013 and 2022. Snakebite data were gathered from the National System of Notifiable Diseases (SINAN), cartographic and demographic data from the Brazilian Institute of Geography and Statistics (IBGE), and climatic data from the WorldClim platform. All data were organized according to the four seasons (spring, summer, fall, and winter). Spatial data analyses were conducted using ArcGIS Pro 3.4 and GeoDa 1.22, employing univariate and bivariate spatial autocorrelation techniques based on Global and Local Moran’s I indices. The results revealed spatial clustering patterns for both Bothrops and Crotalus in all seasons. The main clusters for Bothrops were in the southern and northwestern regions, while Crotalus clusters were concentrated in the central, northwestern, and northeastern regions. A positive spatial autocorrelation between precipitation and Bothrops incidence rates was observed in three seasons. Non-parametric statistical tests also indicated significant seasonal differences in incidence rates for both snake genera.

Keywords:
Bothrops; Crotalus; snakebite epidemiology; neglected tropical disease; spatial statistics

1. Introduction

Snakebite accidents involving humans are classified by the World Health Organization as neglected tropical diseases (WHO 2023a). Nevertheless, these accidents represent a public health problem because they can cause sequelae and may even result in death (Baldassin et al. 2021). Globally, it is estimated that between 4.5 and 5.4 million people are bitten by snakes each year (WHO 2023b). Of this total, approximately 1.8 to 2.7 million develop clinical illness, and between 81,000 and 138,000 die from complications related to envenomation (WHO 2023b). Brazil is among the countries most impacted by snakebite envenomation, both in absolute number of cases and in the territorial extent of risk areas (Silva 2024).

In Brazil, there are about 435 species of snakes, with 38 belonging to the Viperidae family and 34 to the Elapidae family, the two families responsible for most snakebite accidents in the country (Guedes et al. 2023). The Viperidae family includes the genera Bothrops (e.g., jararaca, jararacuçu, urutu), Crotalus (e.g., rattlesnake), and Lachesis (e.g., bushmaster). The Elapidae family is represented by the genus Micrurus (e.g., coral snake) (Pinho and Pereira 2001; Guedes et al. 2023).

A recent study that described the epidemiology of accidents caused by venomous animals in the State of São Paulo between 2010 and 2022 showed a total of 26,280 snakebite accidents, with 67 deaths, making snakes the fourth leading cause of envenomation in the state, after scorpions, spiders, and bees (Eloy 2023). However, the lethality of snakebite accidents was the highest when compared to other venomous animals.

Globally, snakebite accidents are associated with poverty and climatic fluctuations (Chaves et al. 2015). In Brazil, epidemiological studies indicate that low educational attainment and work in rural areas are among the main determinants of snakebite accidents (Baldassin et al. 2020). One of the first studies to address the spatial distribution and association of accidents caused by snakes from the Elapidae and Viperidae families in Brazilian municipalities was conducted for the period 2007-2015 (Santos 2018). The results indicated that human factors contributed the most to the occurrence of these accidents, followed by spatial effects and environmental factors (Santos 2018).

Bivariate spatial autocorrelation analyses between snakebite accidents by snake genus and socioeconomic and environmental factors in the state of Paraná, Brazil, showed significant spatial correlations between Bothrops incidence rates and environmental and socioeconomic factors (Kono et al. 2024). Regarding the seasonal incidence of snakebites in Brazilian regions for the period 2001-2012, a scientific investigation showed that the highest rates occurred between November and April in all regions, with more pronounced seasonal differences observed in the South region (Chippaux 2015).

More recent investigations have applied predictive modeling techniques to study snakebite incidence rates, but without fully capturing the spatial complexity of the variables involved. For example, the Generalized Least Squares (GLS) model with spatial correlation structure in the residuals was used to assess the association between Bothrops snakebite incidence rates and potential explanatory factors in Brazilian municipalities (Silva et al. 2024). Among the findings for the Southeast region, the authors identified a positive association between Bothrops incidence rates and the presence of native forest, historical forest loss, and the high climatic suitability of the species Bothrops jararaca.

In the specific context of the State of São Paulo, other approaches have deepened the analysis of snakebite accidents from different spatial and temporal perspectives. A recent study addressed the epidemiological profile, seasonality, vegetation cover, and the spatial distribution pattern of snakebite accidents in the state of São Paulo for the period 2010 to 2022 (Freitas et al. 2025). The results corroborate some of the findings by Chippaux (2015), indicating that snakebite incidence is higher between November and April. Additionally, the spatial autocorrelation analysis using the Getis-Ord (Gi*) statistic identified clusters of high incidence rates caused by Bothrops, Crotalus, and Micrurus in several Regional Health Departments.

Despite the findings provided by Freitas et al. (2025), other studies that also analyzed snakebite accidents by snake genus in the state of São Paulo have reported partially divergent results. Considering the four snake genera occurring in the State of São Paulo between 2013 and 2022, the lowest absolute and relative numbers of snakebite accidents were recorded for the genera Lachesis and Micrurus, which also did not present a spatial pattern significantly different from random at the municipal level (Sotocorno et al. 2024). Thus, there is divergence among similar research results produced by different authors, contrasting with part of the findings by Freitas et al. (2025).

Despite the advances identified in the literature, three major aspects remain uninvestigated in the context of the State of São Paulo: (a) the distribution of spatial patterns and spatial outliers of snakebite incidence rates for Bothrops and Crotalus across the different seasons; (b) the degree of spatial association between snakebite incidence rates and climatic factors such as temperature and precipitation; and (c) the significance level of differences among seasonal snakebite incidence rates in the state. It is expected that the results of this research may support the formulation of public policies aimed at more effective distribution of antivenom, as well as the implementation of educational and preventive measures directed at the most vulnerable populations.

The central hypothesis of this study is that the spatial associations between snakebite incidence rates and climatic factors differ from random. Based on the foregoing, the present study aims to analyze the spatial distribution pattern of snakebite accidents caused by Bothrops and Crotalus snakes in association with climatic factors in the State of São Paulo, during the period from 2013 to 2022. The research questions addressed are: “What is the spatial distribution pattern of Bothrops and Crotalus snakebite incidence rates in each season of the year?”; “Is there spatial association between Bothrops or Crotalus snakebite incidence rates and temperature and precipitation in each season?”; and “Are the differences in Bothrops or Crotalus incidence rates across seasons statistically significant?”.

2. Materials and methods

2.1 Characterization of the study area

The State of São Paulo is composed of 645 municipalities, covering a territorial area of 248,219.485 km². With a population of 44,411,238 inhabitants and a population density of 178.92 inhabitants per km², its Human Development Index (HDI) is 0.806, making it the second highest in the country (IBGE 2021; IBGE 2022). The municipalities of the state are grouped into 17 Regional Health Departments (DRS), which serve as administrative divisions of the State Health Department (Figure 1) (São Paulo, 2025). These departments are responsible for coordinating the actions of the department in their respective regions, promoting intersectoral collaboration among municipalities and other civil society entities.

Rainfall indices in the State of São Paulo show regional variations. According to the average monthly precipitation between 1970 and 2000, as extracted from the WorldClim platform (https://www.worldclim.org), the southern, southeastern, and southwestern portions of the state record precipitation values ranging from 95 mm to 275 mm, whereas in the northern, northwestern, and parts of the central-western and western-northwestern regions, values range from 89 mm to 135 mm.

According to the Köppen-Geiger climate classification, the State of São Paulo includes types A (tropical) and C (humid temperate) climates. The tropical climate A is subdivided into Af (humid tropical with no dry season) and Aw (humid tropical with dry winter). The temperate climates C are subdivided into Cwa (hot with dry winter), Cwb (temperate with dry winter), Cfa (hot with no dry season), and Cfb (temperate with no dry season) (Martinelli 2010). The region is characterized by rainfall in the summer and a dry climate in the winter.

Figure 1:
Spatial distribution of Regional Health Department - DRS across SP.

2.2 Data source and organization

Cartographic data from the 2018 (IBGE 2025), annual population estimates from 2013 to 2021, and the 2022 census population (IBGE 2022) were used. Data on snakebite accidents were obtained from the Notifiable Diseases Information System (SINAN) (DATASUS 2023), disaggregated by month and by snake genus (Bothrops and Crotalus), for the period from 2013 to 2022. The analysis considered snakebite incidence by municipality of occurrence, excluding residence and notification locations. The data was organized and integrated into the ArcGIS Pro geographic information system, version 3.4.

Temperature and precipitation data were extracted from the WorldClim platform in raster format, with a spatial resolution of 30 arc-seconds (Fick and Hijmans 2017). These data were converted to vector format in ArcGIS Pro, and zonal statistics were used to estimate the mean temperature and mean precipitation for each of the 645 municipalities, considering the four seasons of the year.

The numbers of accidents were summed for each month during the study period. The absolute values were then grouped by season: spring (September, October, and November), summer (December, January, and February), fall (March, April, and May), and winter (June, July, and August). Seasonal accident rates were calculated based on the average population during the period and expressed per 100,000 inhabitants. Considering that snakebite notification is mandatory, null values were replaced with zero, assuming no cases in the corresponding municipalities.

2.3 Spatial autocorrelation analysis

An exploratory spatial data analysis was performed through quantitative thematic mapping, using both incidence rates and absolute values. While incidence rates were mapped using a choropleth approach with the Natural Breaks classification method, absolute values were represented using proportional symbol mapping, based on the direct ratio technique (Dent et al. 2009). Next, a spatial autocorrelation analysis was carried out to determine whether the aggregated data exhibited spatial dependence, that is, whether the variation of a variable could be influenced by its geographic location rather than occurring entirely at random (Camargo et al. 2004).

To analyze the spatial autocorrelation of the aggregated mean incidence rates of snakebite accidents among municipalities for each season, both the Global Moran’s I index (Moran 1950) and the Local Moran’s I (LISA - Local Indicators of Spatial Association) (Anselin 1995) were employed. The global analysis describes the overall spatial pattern of the variable under study, while the local analysis identifies specific spatial clusters in particular areas within the region of interest. Moran’s I values range from -1 to +1, with positive values indicating positive spatial autocorrelation (spatial similarity) and negative values indicating negative spatial autocorrelation (spatial dissimilarity).

The spatial autocorrelation analysis was conducted in two stages. First, only the aggregated incidence rates of Bothrops or Crotalus bites were considered, separately for each season. Then, to assess the strength of the relationship between temperature or precipitation and the incidence rates of accidents caused by Bothrops or Crotalus, two variables were analyzed simultaneously. This second case corresponds to the calculation of the bivariate Moran’s I index, which examines how one variable at a given location influences another variable in neighboring locations (Anselin et al. 1995; Grekousis 2020; Anselin 2024).

For both the univariate and bivariate analyses, the Global Moran’s I index (Moran 1950) was first calculated, followed by the identification of spatial clusters through Anselin’s Local Indicator of Spatial Association - LISA (Anselin et al. 1995; Grekousis 2020; Anselin 2024). In both Moran analyses, the spatial relationship were defined using inverse distance, with a 40-kilometer radius established with the support of Ripley’s K function (Ripley 2004), and 999 permutations were applied.

The bivariate spatial autocorrelation analysis, which determines the strength of the association between two variables, assessed the relationships between the incidence rates of Bothrops or Crotalus and precipitation and temperature, separately for each season. Only the variables that exhibited clustered and positive spatial patterns advanced to the application of Anselin’s Bivariate Local Indicator of Spatial Association (Anselin 1995; Grekousis 2020; Anselin 2024).

The identification of clustered spatial patterns in both univariate and bivariate analyses was based on three main criteria: positive or negative Moran’s I, associated with a p-value ≤ 0.05, and a Z-score ≤ -1.96 or ≥ 1.96. While the bivariate analysis was performed in GeoDa 1.22, the univariate analysis and map layouts were carried out in ArcGIS Pro.

2.4 Seasonal incidence comparison: Friedman and Bonferroni tests

Finally, the Friedman and Bonferroni statistical tests were applied as part of the nonparametric statistical analysis (Siegel and Castelan JR. 2006; Volpato and Barreto 2016). These tests allow for verification of whether seasonal differences in snakebite incidence rates, for each snake genus, are statistically significant at a 95% confidence level (α = 0.05). The Friedman test was first applied, as it is based on medians - appropriate since the incidence rates did not follow a normal distribution and the data were dependent. The Friedman test aims to identify whether statistically significant differences exist within a dataset. When such differences were confirmed, the Bonferroni test was applied to determine, pairwise, between which variables these differences occurred.

3. Results and discussion

3.1 Quantitative thematic mapping of snakebites

To analyze the distribution pattern of snakebite incidence rates caused by Bothrops and Crotalus snakes in association with climatic factors in the State of São Paulo, during the period from 2013 to 2022, the results aim to answer questions related to the type and location of spatial patterns of incidence rates for each snake genus throughout the seasons, as well as the patterns associated with precipitation and temperature variables in each season, and the significance level of the differences in incidence rates among the different seasons.

Between 2013 and 2022, a total of 14,716 snakebite cases caused by Bothrops and Crotalus were recorded in the State of São Paulo. Of these, approximately 83.4% were attributed to Bothrops and 16.6% to Crotalus. During the study period, the average incidence rate of snakebites caused by Bothrops was 94.04 per 100,000 inhabitants, while the average rate for Crotalus was 24.21 per 100,000 inhabitants.

The choropleth mapping using the Jenks classification algorithm revealed distinct characteristics in the spatial distribution of incidence rates caused by Bothrops (Figure 2-a) and Crotalus (Figure 2-b) in the State of São Paulo. Our findings regarding the predominant locations for the two studied genera are in line with ecological niche studies, which have shown an inverse correlation in the geographic distribution of Bothrops and Crotalus, considering climatic and environmental factors (Yañez-Arenas et al. 2014).

The main cluster of Bothrops incidence rates is found in the Regional Health Department (DRS) of Registro, in the southern portion of the state, followed by smaller clusters in the DRS of Sorocaba, Taubaté, São José do Rio Preto, and Araçatuba - the latter located in the northwestern part of the state (Figure 2-a). In contrast, small clusters of Crotalus incidence rates are distributed across several DRS regions of the state, with concentrations in Marília, Bauru, and São José do Rio Preto, as well as in the northeastern parts of Franca, Ribeirão Preto, and São João da Boa Vista (Figure 2-b).

Figure 2:
Spatial distribution of incidence rates caused by Bothrops and Crotalus in São Paulo State (2013 - 2022).

For the same study period, the proportional symbol maps (Figure 3) provide additional insights, highlighting different behavioral patterns between the two snake genera. Accidents caused by Bothrops are distributed throughout the state, with concentrations particularly in the DRS of Registro, Baixada Santista, Greater São Paulo, and Taubaté (Figure 3-a), on the São Paulo coast. In the case of Crotalus envenomation, although cases are also spread across much of the state, their occurrence is considerably lower than that of Bothrops, showing clusters in the DRS of Bauru, Piracicaba, Campinas, and São João da Boa Vista, and an almost complete absence in the DRS of Registro (Figure 3-b).

3.2 Univariate spatial autocorrelation of accident rates

The spatial autocorrelation analysis was based on the null hypothesis (H₀) that the spatial distribution of incidence rates is the result of random chance, with the aim of validating the alternative hypothesis (H₁) that the spatial distribution is not random. Clustered spatial patterns were detected for all seasons of the year, for both Bothrops and Crotalus, with a confidence level greater than 95% (p-value < 0.05), indicating positive spatial autocorrelation (Table 1). For Bothrops, the Moran’s I index was moderate (between 0.4 and 0.6) in three of the four seasons (spring, summer, and fall). In the case of Crotalus, although all Moran’s I values for the seasons were very weak (between 0.0 and 0.2), the detected clusters were highly statistically significant.

Figure 3:
Spatial distribution of accidents caused by Bothrops and Crotalus, in Sao Paulo State (2013 - 2022).

Table 1:
Spatial autocorrelation between municipalities in Sao Paulo State for the incidence of snakebite using the Moran’s global index.

Gender	Season	Moran’s Index	Z-score	Pseudo p-value	Interpretation
Bothrops	Spring	0.428121	24.046791	< 0.0001	Clustered
	Summer	0.404293	22.787114	< 0.0001	Clustered
	Fall	0.435169	24.803348	< 0.0001	Clustered
	Winter	0.375180	20.942193	< 0.0001	Clustered
Crotalus	Spring	0.083019	4.693183	< 0.0001	Clustered
	Summer	0.129168	7.271683	< 0.0001	Clustered
	Fall	0.134222	7.578390	< 0.0001	Clustered
	Winter	0.092670	5.281886	< 0.0001	Clustered

The results of the local spatial autocorrelation analysis, using the Local Moran’s I index to identify the incidence rate patterns of Bothrops and Crotalus in the State of São Paulo, showed that the locations of the main clusters varied distinctly between the two snake genera across the four seasons (Figures 4 and 5). Variations in the size of high-incidence clusters of Bothrops and Crotalus were also observed throughout the year. These results are consistent with those presented in the choropleth maps.

The largest high-incidence clusters of Bothrops were identified in the DRS of Registro and Sorocaba in three seasons (spring, summer, and fall) (Figures 4a, 4b, 4c). Notably, in the DRS of Registro, spatial autocorrelation of accident rates remained high in all seasons analyzed. Our results demonstrate that the highest concentrations of Bothrops accidents in the state occur in the Dense Ombrophilous Forest of the Atlantic Forest biome (MAPBIOMAS Brasil 2025) - a region that also has high rainfall indices (Martinelli 2010). A related study demonstrated a spatial association between Bothrops incidence rates and forest formations in the state of Paraná (Kono et al. 2024).

Furthermore, consistent clusters of high Bothrops incidence rates were observed along the borders of the DRS of São José do Rio Preto and Araçatuba, near the Tietê River, across all seasons. The western and northwestern regions are among the hottest and driest in the state (Martinelli 2010). Although the best climatic conditions for the Bothrops genus are environments with high precipitation and temperature (Brasil 2001), our findings for the northwest of the state partially diverge from those of Silva et al. (2024), who analyzed the spatial distribution of incidence rates at the national scale. The state-level analysis revealed local patterns that were concealed in the national analysis due to the influence of much higher incidence rates in areas of northern Brazil.

The spatial distribution of Bothrops and Crotalus accident rates presented in the univariate spatial autocorrelation maps of this study is consistent with the results of the choropleth mapping. However, the autocorrelation findings here differ from a study also conducted in the state of São Paulo (Freitas et al. 2025). In that study, the authors found clusters of high Bothrops incidence rates only along the Dense Ombrophilous Forest, including the DRS of Greater São Paulo. In the present study, the clusters along the forest are strongly distributed in the DRS of Registro and part of the DRS of Sorocaba - areas of Dense Forest and deforestation activity, respectively - as well as in São José do Rio Preto and Araçatuba.

Figure 4:
Spatial clusters and outliers of Bothrops incidence rate across seasons.

An analysis of high-incidence spatial clusters for Crotalus (Figure 5) showed that the largest areas are concentrated in the DRS of Marília and Bauru, followed by portions in northern São José do Rio Preto and areas in northeastern Franca, Ribeirão Preto, and São João da Boa Vista. The predominance of high-incidence clusters for this snake genus occurs in the Cerrado biome, as Crotalus snakes are typically found in open fields, dry areas, sandy or rocky terrain, and are rarely found along the coast (Brasil 2001). Another important aspect is that the most notable clusters were found in the fall season (Figure 5-c).

With respect to the Crotalus genus, the clustered patterns of high incidence rates found by Freitas et al. (2025) also differ from those in the present study. Additionally, those authors also identified a clustered pattern in the DRS of Campinas, which was not found in this investigation.

Figure 5:
Spatial clusters and outliers of Crotalus incidence rate across seasons.

3.3 Spatial analysis of accident rates associated with climatic factors

The spatial analysis of accident rates in association with seasonal temperature and precipitation was conducted based on the null hypothesis (H₀) that the distribution of this relationship occurs at random, aiming to validate the alternative hypothesis (H₁) that the distribution is not random. Analyses were performed with a confidence level greater than 95% (p < 0.05). With respect to Crotalus, no significant correlations were identified between incidence rates and climatic factors.

The results of the bivariate spatial autocorrelation, which evaluated the distribution patterns between seasonal incidence rates of Bothrops accidents and seasonal means of temperature or precipitation, revealed two very weak positive spatial clusters for precipitation in spring (I = 0.139) and winter (I = 0.159), as well as a weak positive cluster in fall (I = 0.265) (Table 2). On the other hand, no relevant spatial association was identified between Bothrops rates and precipitation in summer (I = 0.042).

The spatial identification of the three main clustered patterns of high Bothrops incidence rates is revealed in the DRS of Registro, Sorocaba, Baixada Santista, and Taubaté (Figure 6). Conversely, the analysis also revealed spatial outliers, that is, areas with high incidence rates associated with low precipitation in the Northwest region of the state, during spring, fall, and winter (Figures 6-a, 6-b, 6-c). This result reinforces the ability of Bothrops snakes to inhabit intensely modified environments (Sazima, 1988; Silva et al. 1988), as observed in the Northwest of the state, where the proportion of native forest is extremely low, but the increase in vegetation cover, especially in areas of permanent preservation, has been recorded in recent decades, according to the Forest Inventory of the State of São Paulo (São Paulo 2020).

Bothrops snakes are found in rural areas and on the outskirts of major urban centers, with a preference for humid environments such as forested areas and cultivated lands, as well as locations favorable to the presence of rodents (Brasil 2001). In contrast, Crotalus snakes inhabit open fields, dry, sandy, and rocky areas, being rare along the coastal region (Brasil 2001). Additional studies have shown that snakebite accidents, both in Costa Rica and Brazil, are more likely in areas with higher temperatures (Chaves et al. 2015; Chippaux 2017), but it is important to note that areas located at latitudes farther from the Equator experience effects distinct from those in milder locations, as in the southernmost regions of Brazil (Chippaux 2017).

Table 2:
Spatial autocorrelation between climatic variables and incidence rate by snakebites in Sao Paulo State using bivariate Moran’s global index.

Climatic variable	Season	Bothrops		Crotalus
		Moran’s index	Interpretation	Moran’s Index	Interpretation
Precipitation	Spring	0.139	Clustered	-0.066	Random
	Summer	0.042	Random	-0.027	Random
	Fall	0.265	Clustered	-0.088	Clustered
	Winter	0.159	Clustered	-0.090	Clustered
Temperature	Spring	-0.075	Random	0.060	Random
	Summer	0.010	Random	0.015	Random
	Fall	-0.034	Random	-0.036	Random
	Winter	0.029	Random	0.075	Random

3.3 Comparison between accident rates and climatic factors

Accident notifications occurred predominantly in the summer, followed by fall, spring, and winter, respectively, for both snake genera (Figure 7-a), which corroborates the findings of Chippaux (2015), showing that the seasonal incidence of snakebites in the Southeast was highest between November and May, that is, predominantly in the summer. Figure 7-a displays the absolute number of accidents by snake genus in each season of the year, with the Y-axis representing the number of accidents and the X-axis separating by genus and season. The same applies to Figure 7-b; however, the incidence of accidents is expressed per 100,000 inhabitants on the Y-axis.

Figure 6:
Spatial clusters and outliers between mean precipitation and incidence rates by Bothrops in the Spring, Fall and Winter seasons.

However, when considering each snake genus, the differences in incidence rates shown so far in the graphs of Figure 7-a do not allow us to determine whether the differences are significant. The analysis was conducted under the null hypothesis (H₀) that there is no statistically significant difference between the seasons, versus the alternative hypothesis that at least one of the seasons differs significantly from the others. Analyses were performed with a confidence level greater than 95% (p < 0.05). Since snakebite incidence rates in the State of São Paulo do not follow a normal distribution, it was necessary to use nonparametric statistics, as recommended in the literature (Siegel and Castellan 2006; Volpato and Barreto 2016). The Friedman and Bonferroni techniques are appropriate for dependent groups larger than two.

Figure 7:
Snakebites in the Sao Paulo State, accidents in (a) and incidence rate in (b) (2013-2022).

The results of the comparative statistical tests using the Friedman and Bonferroni techniques revealed significant differences in the incidence rates of each snake genus across the seasons. For the Bothrops genus, the Friedman test indicated statistically significant differences [χ² (3) = 350.084; p < 0.0001]. The multiple comparison analysis performed using the Bonferroni test showed significant differences in incidence rates between the following seasonal pairs: winter and spring (p < 0.0001); winter and fall (p < 0.0001); winter and summer (p < 0.0001); spring and fall (p = 0.011); and spring and summer (p < 0.001).

The literature has shown that increased hunting activity of Bothrops snakes is more intense with higher precipitation and temperature (Chaves 2015) and in the presence of forest formations (Kono et al. 2024; Silva et al. 2025). Although the spatial association between snakebite incidence rates and temperature in the State of São Paulo was not different from random, the results of the seasonal comparison statistics showed the opposite.

For Crotalus incidence rates across the seasons, the Friedman test also revealed significant differences during the study period [χ² (3) = 125.214; p < 0.0001]. The multiple comparison analysis, again performed using the Bonferroni test, showed statistically significant differences in the following comparisons: winter and spring (p < 0.0001); winter and summer (p < 0.0001); and winter and fall (p < 0.0001).

4. Conclusion

This study analyzed the spatial distribution patterns of snakebite incidence rates caused by Bothrops and Crotalus snakes in the State of São Paulo across the four seasons, during the period from 2013 to 2022. It also examined the spatial association between these incidence rates and climatic variables (temperature and precipitation), as well as tested whether seasonal differences in incidence are statistically significant for each genus, thus achieving the previously established objectives.

The spatial distribution of high-incidence patterns and spatial outliers for Bothrops revealed large clusters in the DRS regions of Registro, Sorocaba, and Taubaté along the Dense Ombrophilous Forest, followed by secondary clusters in the DRS regions of Araçatuba and São José do Rio Preto, near the Tietê River. Spatial outliers of low incidence rates - low values adjacent to high values - were mainly identified in the DRS of São José do Rio Preto.

Regarding spatial pattern identification, the major clusters were positively associated with mean precipitation in spring, fall, and winter, while the secondary clusters were negatively correlated with precipitation during these three seasons. This reveals that Bothrops snakes are capable of inhabiting areas heavily altered by human activity. On the other hand, the spatial analyses of Crotalus incidence rates in association with precipitation were random throughout the entire year.

Additionally, Bothrops and Crotalus incidence rates showed distinct spatial patterns, reflecting differences in habitat preference. Crotalus incidence clustered in drier areas, predominating in regions where the Cerrado biome is more prevalent - mainly in the DRS regions of Bauru, Marília, northern São José do Rio Preto, and São João da Boa Vista. In contrast, the results of the comparative statistical tests indicated significant differences in Bothrops or Crotalus accidents between seasons, with winter resulting in the most significant differences.

These results highlight the importance of spatial autocorrelation as a measure of spatial dependence in the analysis of this phenomenon, considering spatial relationships and demonstrating that the presence of clusters indicates the presence of underlying processes. Future studies should explore spatial associations between snakebite incidence and environmental and anthropogenic factors to support even more targeted public health interventions and prevention strategies.

ACKNOWLEDGEMENT

The development of this research was supported by the Programa Institucional de Bolsas de Iniciação Científica (PIBIC), through a scientific initiation scholarship granted by the Conselho Nacional de Desenvolvimento Científico e Tecnológico (CNPq), and by the Pró-Reitoria de Pesquisa (PROPe) da Unesp (Universidade Estadual Paulista “Júlio de Mesquita Filho”), awarded to the main author at different periods.

REFERENCES

• Anselin, L. 1995. Local indicators of spatial association - LISA. Geographical Analysis, 27(2), pp.93-115.
• Anselin, L. Syabri, I. and Smirnov, O. 2002. Visualizing Multivariate Spatial Correlation with Dynamically Linked Windows. In: L. Anselin, S. Rey. New Tools for Spatial Data Analysis: Proceedings of the Specialist Meeting Santa Barbara, University of California: Center for Spatially Integrated Social Science.
• Anselin, L. 2024. An introduction to spatial data science with GeoDa: exploring spatial data (vol. 1). CRC Press.
• Araújo, S. C. M. Ceron, K. Guedes, T. B. 2022. Use of geospatial analyses to address snakebite hotspots in mid-northern Brazil - A direction to health planning in shortfall biodiversity knowledge areas. Toxicon, 213(15), pp.43-51.
• Baldassin, J. C. S. Francisco, S. R. Silva, R. W. Moura, R. F. and Pombo, A. P. M. M. 2021. Perfil epidemiológico e dinâmica da distribuição dos acidentes ofídicos em humanos no estado de São Paulo. Revista Brasileira de Geografia Médica e da Saúde - Hygeia, 17, pp. 216-226.
• Brasil. FUNDAÇÃO NACIONAL DA SAÚDE (FUNASA), 2001. Manual de diagnóstico e tratamento de acidentes por animais peçonhentos [online]. Available at: <Available at: https://www.icict.fiocruz.br/sites/www.icict.fiocruz.br/files/Manual-de-Diagnostico-e-Tratamento-de-Acidentes-por-Animais-Pe--onhentos.pdf > [Accessed on: March 03, 2025].
  » https://www.icict.fiocruz.br/sites/www.icict.fiocruz.br/files/Manual-de-Diagnostico-e-Tratamento-de-Acidentes-por-Animais-Pe--onhentos.pdf
• Câmara, G. Carvalho, M. S. Cruz, O. G. and Correa, V. 2004. Análise de Dados de Área. In S. Druck, M. S. Carvalho, G. Câmara, A. M. V. Monteiro, (eds). Análise Espacial de Dados Geográficos Brasília: EMBRAPA.
• Chaves, L. F. Chuang, T. W. Sasa, M. et al. 2015. Snakebites Are Associated with Poverty, Weather Fluctuations, and El Niño. Science Advances, 1(8), pp. e1500249-e1500249.
• Chippaux, J. P. 2015. Epidemiology of envenomations by terrestrial venomous animals in Brazil based on case reporting: from obvious facts to contingencies. Journal of Venomous Animals and Toxins including Tropical Diseases, v. 21, n. 13, pp.
• Chippaux, J. P. 2017. Incidence and mortality due to snakebite in the Americas. PLOS Neglected Tropical Diseases, 11(6), pp. 1-39.
• Dent, B. D.; Torgoson, J. F.; Hodler, T. W. 2009. Cartography: Thematic Map Design New York: McGraw-Hill.
• Eloy, L. J. 2023. Acidentes por animais peçonhentos: Série Histórica 2010 - 2021. Boletim Epidemiológico Paulista, 20(219), pp.1-10.
• Fick, S. E. and Hijmans, R. J. 2017. WorldClim 2: new 1-km spatial resolution climate surfaces for global land areas. International Journal of Climatology, 37(12), pp. 4302-4315.
• Freitas, G. D. de. Lacerda, A. B. Azevedo, T. S. de. Oliveira, A. de. Spinola, R. M. F. Dourado, F. S. Wen, F. H. and Chiaravalloti-Neto, F. 2025. Quais características e onde há maior risco de acidente ofídico no estado de São Paulo? Revista Brasileira de Epidemiologia, 28, pp. e250026.
• Guedes, T. B. Entiauspe-Neto, O. M. and Costa, H. C. 2023. Lista de répteis do Brasil: atualização de 2022.Herpetologia Brasileira, 12(1), pp. 56-161.
• Grekousis, G. 2002. Spatial analysis methods and practice: describe - explore - explain through GIS Cambridge: Cambridge University Press.
• IBGE, Instituto Brasileiro de Geografia e Estatística. 2025. Malha Municipal [online] Available at: <Available at: https://www.ibge.gov.br/geociencias/organizacao-do-territorio/malhas-territoriais/15774-malhas.html?edicao=24048&t=acesso-ao-produto > [Accessed on: March 03, 2025].
  » https://www.ibge.gov.br/geociencias/organizacao-do-territorio/malhas-territoriais/15774-malhas.html?edicao=24048&t=acesso-ao-produto
• IBGE, Instituto Brasileiro de Geografia e Estatística. 2023a. IBGE define bioma predominante em cada município brasileiro para fins estatísticos [online] Available at: <Available at: https://agenciadenoticias.ibge.gov.br/agencia-noticias/2012-agencia-de-noticias/noticias/40519-ibge-define-bioma-predominante-em-cada-municipio-brasileiro-para-fins-estatisticos > [Accessed on: March 03, 2025].
  » https://agenciadenoticias.ibge.gov.br/agencia-noticias/2012-agencia-de-noticias/noticias/40519-ibge-define-bioma-predominante-em-cada-municipio-brasileiro-para-fins-estatisticos
• IBGE, Instituto Brasileiro de Geografia e Estatística. 2023b. Censo Demográfico 2022: Resultados Gerais da Amostra [online] Available at: <Available at: https://censo2022.ibge.gov.br/ > [Accessed on: March 03, 2025].
  » https://censo2022.ibge.gov.br/
• IBGE, Instituto Brasileiro de Geografia e Estatística. 2025. Produção Agropecuária [online] Available at: <Available at: https://www.ibge.gov.br/explica/producao-agropecuaria/sp > [Accessed on: March 03, 2025].
  » https://www.ibge.gov.br/explica/producao-agropecuaria/sp
• Kono, I. S. Pandolfi, V. C. F. Marchi, M. N. A. Freitas, N. and Freire, R. L. 2024. Unveiling the secrets of snakes: Analysis of environmental, socioeconomic, and spatial factors associated with snakebite risk in Paraná, Southern Brazil. Toxicon , 237, pp.
• MAPBIOMAS Brasil, 2025. Uso e cobertura do Solo [online] Available at: <Available at: https://mapbiomas.org/ > [Accessed on: March 03, 2025].
  » https://mapbiomas.org/
• Martinelli, M. 2010. Clima do Estado de São Paulo. Confins, 8, pp77-78.
• Mise, Y. Y. Lira-da-Silva, R. M. and Carvalho, F. M. 2016. Agriculture and snakebite in Bahia, Brazil - An ecological study. Annals of Agricultural and Environmental Medicine, 23(3), pp. 416 419.
• Moran, P. A. P. 1950. Notes on continuous stochastic phenomena. Biometrika, 37, pp. 17-23.
• Pinho, F. M. O.; Pereira, I. D. 2001. Ofidismo. Rev. Assoc. Med. Bras., 47(1), pp.24-29
• Ripley, B. D. 2004. Spatial statistics New Jersey: Wiley-Interscience.
• São Paulo, 2025. Departamentos Regionais de Saúde [online] Available at: <Available at: https://www.saude.sp.gov.br/ses/institucional/departamentos-regionais-de-saude/ > [Accessed on: March 03, 2025].
  » https://www.saude.sp.gov.br/ses/institucional/departamentos-regionais-de-saude/
• São Paulo, 2020. Inventário Florestal do Estado de São Paulo: mapeamento da cultura vegetal nativa. [online] Available at: <Available at: https://smastr16.blob.core.windows.net/home/2020/07/inventarioflorestal2020.pdf > [Accessed on: July 14, 2025].
  » https://smastr16.blob.core.windows.net/home/2020/07/inventarioflorestal2020.pdf
• Santos, C. J. C. 2018. Análise espacial da ocorrência de acidentes ofídicos no Brasil e seus determinantes ambientais e socioeconômicos PhD. Pontifícia Universidade Católica de Goiás.
• Sazima, I. 1988. Um estudo de biologia comportamental da jararaca, Bothrops jararaca, com uso de marcas naturais.Mem Inst Butantan, 50(3), pp. 83-99.
• Silva, F. F. B. Moura, T. A. Siqueira-Silva, T. Gutiérrez, J. M. and Martinez, P. A. 2024. Predicting the drivers of Bothrops snakebite incidence across Brazil: a spatial analysis. Toxicon , 250(), pp. 108107.
• Sotocorno, G. P. Haiachi, G. H. Serrano, M. C. M. Ramos, A. P. M. and Pugliesi, E. A. 2024. Análise dos Padrões Espaciais de Acidentes por Serpentes dos Gêneros Lachesis e Micrurus no Estado de São Paulo (2013-2022). In: Encontro Nacional de Ensino, Pesquisa e Extensão, Presidente Prudente, SP. Inteligência artificial: impacto nas profissões.
• Siegel, S. and Castellan, N. J. 2006. Estatística não-paramétrica para as ciências do comportamento Porto Alegre: Artmed.
• Volpato, G. and Barreto, R. 2016. Estatística sem dor Botucatu: Best Writing.
• World Health Organization (WHO), 2023a. Neglected tropical diseases [online] Available at: <Available at: https://www.who.int/news-room/questions-and-answers/item/neglected-tropical-diseases > [Accessed on: March 03, 2025].
  » https://www.who.int/news-room/questions-and-answers/item/neglected-tropical-diseases
• World Health Organization (WHO), 2023b. Snakebite envenoming [online] Available at: <Available at: https://www.who.int/health-topics/snakebite > [Accessed on: March 03, 2025].
  » https://www.who.int/health-topics/snakebite
• Yañez-Arenas, C. Peterson, A.T. Mokondoko, P. Rojas-Soto, O. and Martõânez-Meyer, E. 2014. The Use of Ecological Niche Modeling to Infer Potential Risk Areas of Snakebite in the Mexican State of Veracruz. PLoS ONE, 9(6), pp. e100957.

• DATA AVAILABILITY
  The datasets generated during and/or analyzed during the current study are available from the corresponding author on reasonable request.

Data availability

The datasets generated during and/or analyzed during the current study are available from the corresponding author on reasonable request.

Data citations

IBGE, Instituto Brasileiro de Geografia e Estatística. 2025. Malha Municipal [online] Available at: <Available at: https://www.ibge.gov.br/geociencias/organizacao-do-territorio/malhas-territoriais/15774-malhas.html?edicao=24048&t=acesso-ao-produto > [Accessed on: March 03, 2025].

IBGE, Instituto Brasileiro de Geografia e Estatística. 2023a. IBGE define bioma predominante em cada município brasileiro para fins estatísticos [online] Available at: <Available at: https://agenciadenoticias.ibge.gov.br/agencia-noticias/2012-agencia-de-noticias/noticias/40519-ibge-define-bioma-predominante-em-cada-municipio-brasileiro-para-fins-estatisticos > [Accessed on: March 03, 2025].

IBGE, Instituto Brasileiro de Geografia e Estatística. 2023b. Censo Demográfico 2022: Resultados Gerais da Amostra [online] Available at: <Available at: https://censo2022.ibge.gov.br/ > [Accessed on: March 03, 2025].

IBGE, Instituto Brasileiro de Geografia e Estatística. 2025. Produção Agropecuária [online] Available at: <Available at: https://www.ibge.gov.br/explica/producao-agropecuaria/sp > [Accessed on: March 03, 2025].

Publication Dates

• Publication in this collection
  14 Nov 2025
• Date of issue
  2025

History

• Received
  31 Mar 2025
• Accepted
  28 Aug 2025