Bayesian spatial modelling and landscape correlates of snakebite risk in Fuzhou, China, 2021-2025

Abstract

Objective: To evaluate whether combining Bayesian spatial modelling with landscape-ecological analysis could identify high-risk townships and account for their spatial clustering in Fuzhou, Fujian Province, China.

Methods: In this retrospective study, patients with clinically confirmed venomous snakebite treated at a single provincial tertiary referral centre in Fuzhou between January 2021 and December 2025 were included in descriptive demographic, temporal, and taxonomic analyses. In addition, patients with concordant township-level residence and bite location were included in spatial analyses. Five-year cumulative incidence per 100 000 population was calculated for each of 141 townships. A Besag-York-Mollié 2 (BYM2) model was used to obtain smoothed relative risks, and Global Moran's I and Getis-Ord Gi* statistics assessed spatial clustering and hotspots. Associations between landscape metrics and township-level case counts were examined using negative binomial regression with log population included as an offset.

Results: A total of 741 patients were eligible for this retrospective analysis, and 560 of them were included in spatial analyses. Viperidae accounted for 589 of 741 cases (79.49%); patients aged 40-79 years accounted for 85.02% of cases, with a modest male predominance (1.28:1). Bayesian-smoothed incidence showed significant spatial clustering (Moran's I=0.322, P=0.001), and Getis-Ord Gi* analysis identified 13 hotspot townships, with smoothed relative risks ranging from 3.40 to 22.16. Forest connectivity (COHESION IRR=3.37, 95% CI 1.72-6.59) was positively, and forest patch isolation (ENN_MN IRR=0.59, 95% CI 0.49-0.72) inversely, associated with incidence; cropland shape complexity (LSI IRR=1.48, 95% CI 1.22-1.80) was also positively associated. These landscape covariates accounted for the spatial structure identified by BYM2, with no significant residual autocorrelation in model residuals (Moran's I=-0.110, P=0.968).

Conclusions: Landscape configuration, particularly forest connectivity and cropland edge complexity, is statistically associated with the spatial distribution of snakebite risk in Fuzhou and appears statistically sufficient to account for the clustering identified by Bayesian spatial smoothing. These findings, while ecological and correlational rather than causal, offer a replicable framework for prioritising township-level snakebite prevention resources.

Keywords

Snakebite / Spatial epidemiology / Bayesian spatial modelling / Landscape ecology / Risk assessment

References

[1] World Health Organization. Snakebite envenoming: fact sheet. [Online]. Available from: https://www.who.int/news-room/fact-sheets/detail/snakebite-envenoming. [Accessed on 6 February 2026].

[2] Williams SS, Wijesinghe CA, Jayamanne SF, Buckley NA, Dawson AH, Lalloo DG, et al. Delayed psychological morbidity associated with snakebite envenoming. PLoS Negl Trop Dis 2011; 5(8): e1255.

[3] Harrison RA, Hargreaves A, Wagstaff SC, Faragher B, Lalloo DG. Snake envenoming: A disease of poverty. PLoS Negl Trop Dis 2009; 3(12): e569.

[4] Kasturiratne A, Lalloo DG, de Silva HJ. Chronic health effects and cost of snakebite. Toxicon X 2021;9-10:100074.

[5] World Health Organization. Snakebite envenoming: A strategy for prevention and control. Geneva: WHO; 2019.

[6] Lai R, Yan S, Wang S, Yang S, Yan Z, Lan P, et al. The Chinese guideline for management of snakebites. World J Emerg Med 2024; 15(5): 333-355.

[7] National Health Commission of the People's Republic of China. Diagnosis and treatment guidelines for common animal injuries (2021 edition). 2021. [Online]. Available from: https://www.nhc.gov.cn/wjw/c100175/202108/126e735653fe4ea78912fcc7d2e7dace.shtml. [Accessed on 25 April 2026].

[8] Huang Z, Lv Xd, Huang DL, Wang F, Wei WX, Pang GD, et al. Retrospective analysis of clinical and epidemiological characteristics of Chinese cobra bites in a South China Emergency Department. Sci Rep 2025; 15: 21960.

[9] Liang M, Liang P, Luo W, Wu Y, Peng M, Yu P. Epidemiological analysis of 952 venomous snake bite in Wuzhou City of Guangxi Zhuang Autonomous Region. Zhonghua Wei Zhong Bing Ji Jiu Yi Xue 2022; 34(1): 85-90.

[10] Shen Y, Hu S, Bao S, Tu M. Clinical characteristics of snakebite patients attending a tertiary hospital in Hangzhou, Southeast China. Trans R Soc Trop Med Hyg 2024; 118(4): 287-292.

[11] Tianxing L, Tingting C, Zhe Z. Epidemiological analysis of snake bites in Yangjiang City of Guangdong Province in the past five years. J Clin Emerg 2020; 21(9): 727-730.

[12] Pintor AFV, Ray N, Longbottom J, Bravo-Vega CA, Yousefi M, Murray KA, et al. Addressing the global snakebite crisis with geo-spatial analyses—Recent advances and future direction. Toxicon X 2021; 11: 100076.

[13] Rue H, Martino S, Chopin N. Approximate Bayesian inference for latent Gaussian models by using integrated nested laplace approximations. J R Stat Soc Ser B 2009; 71(2): 319-392.

[14] Riebler A, Sørbye SH, Simpson D, Rue H. An intuitive Bayesian spatial model for disease mapping that accounts for scaling. Stat Methods Med Res 2016; 25(4): 1145-1165.

[15] Martín G, Yáñez-Arenas C, Rangel-Camacho R, Murray KA, Goldstein E, Iwamura T, et al. Implications of global environmental change for the burden of snakebite. Toxicon X 2021;9-10:100069.

[16] Tiglao PJ, Aoki Y, Huy NT. Human–snake conflict and snakebite envenoming: A planetary health perspective. Trans R Soc Trop Med Hyg 2025; 119(12): e18-e19.

[17] Osinski Francisco E, Fernandes M, Somariva Prophiro J. Land use and land cover factors associated with snakebites: A systematic review. Rev Fac Ciênc Méd Sorocaba 2025; 27: e67141.

[18] Rai A, Chettri M, Dewan S, Khandelwal B, Chettri B. Epidemiological study of snakebite cases in sikkim: Risk modeling with regard to the habitat suitability of common venomous snakes. PLoS Negl Trop Dis 2021; 15(11): e0009800.

[19] Alcântara JA, Bernarde PS, Sachett J, da Silva AM, Valente SF, Peixoto HM, et al. Stepping into a dangerous quagmire: Macroecological determinants of Bothrops envenomings, Brazilian Amazon. PLoS One 2018; 13(12): e0208532.

[20] Jucá TL, Normando LRO, Ibrahim AB, Chapeaurouge A, Monteiro-Moreira ACO, Mackessy SP. Drought, desertification and poverty: A geospatial analysis of snakebite envenoming in the Caatinga Biome of Brazil. Int J Health Plann Manage 2021; 36(5): 1685-1696.

[21] Row JR, Blouin-Demers G, Lougheed SC. Movements and habitat use of eastern foxsnakes (Pantherophis gloydi) in two areas varying in size and fragmentation. J Herpetol 2012; 46(1): 94-99.

[22] Marshall BM, Crane M, Silva I, Strine CT, Jones MD, Hodges CW, et al. No room to roam: King Cobras reduce movement in agriculture. Mov Ecol 2020; 8: 33.

[23] Stephenson MD, Schulte LA, Klaver RW. The relative contributions of habitat area, configuration, and vegetative diversity on snake and lizard presence in agricultural landscapes. Conserv Sci Pract 2024; 6(4): e13100.

[24] Bauder JM, Breininger DR, Bolt MR, Legare ML, Jenkins CL, Rothermel BB, et al. Multi-level, multi-scale habitat selection by a wide-ranging, federally threatened snake. Landsc Ecol 2018; 33(5): 743-763.

[25] Blouin-Demers G, Weatherhead PJ. Habitat use by black rat snakes (Elaphe obsoleta obsoleta) in fragmented forests. Ecology 2001; 82(10): 2882-2896.

[26] Guo CP, Zhong MJ, Wang XY, Yang SN, Tang K, Jia LL, et al. An updated species checklist of amphibians and reptiles in Fujian Province, China. Biodiv Sci 2022; 30(8): 22090.

[27] Zhang X, Liu L, Zhao T, Zhang W, Guan L, Bai M, et al. GLC_FCS10: A global 10 m land-cover dataset with a fine classification system from sentinel-1 and sentinel-2 time-series data in Google earth engine. Earth Syst Sci Data 2025; 17(8): 4039-4062.

[28] Jiang H, Ku M, Zhou X, Zheng Q, Liu Y, Xu J, et al. CropLayer: A 2 m resolution cropland map of China for 2020 from Mapbox and Google satellite imagery. Earth Syst Sci Data 2025; 17(12): 6703-6729.

[29] Todd BD, Nowakowski AJ, Rose JP, Price SJ. Species Traits explaining sensitivity of snakes to human land use estimated from citizen science data. Biol Conserv 2017; 206: 31-36.

[30] McGarigal K, Cushman SA, Neel MC, Ene E. FRAGSTATS: Spatial pattern analysis program for categorical maps. Amherst, MA: University of Massachusetts; 2002.

[31] Kjoss VA, Litvaitis JA. Community structure of snakes in a human-dominated landscape. Biol Conserv 2001; 98: 285-292.

[32] Palmeirim AF, de Fraga R, Vieira MV, Peres CA. Species-area relationships induced by forest habitat fragmentation apply even to rarely detected organisms. Trop Conserv Sci 2021; 14: 1-7. doi: 10.1177/19400829211017366.

[33] Nordberg E, Ashley J, Hoekstra AA, Kirkpatrick S, Cobb VA. Small nature preserves do not adequately support large-ranging snakes: Movement ecology and site fidelity in a fragmented rural landscape. Glob Ecol Conserv 2021; 28: e01715.

[34] Breininger DR, Mazerolle MJ, Bolt MR, Legare ML, Drese JH, Hines JE. Habitat fragmentation effects on annual survival of the federally protected eastern indigo snake. Anim Conserv 2012; 15(4): 361-368.

[35] Zhao E. Venomous snakes of China. In: Gopalakrishnakone P, Chou LM. (eds.) Snakes of medical importance (Asia-Pacific region). Singapore: National University of Singapore; 1990, p. 243-269.

[36] Zhao EM, Huang MH, Zong Y, Jiang YM, Huang QY, Zhao H, et al. Fauna sinica, reptilia, vol.3: squamata: serpentes. Beijing: Science Press; 1998, p. 1-522.

[37] Devan-Song A, Martelli P, Karraker NE. Reproductive biology and natural history of the white-lipped pit viper (Trimeresurus albolabris Gray, 1842) in Hong Kong. Herpetol Conserv Biol 2017; 12: 41-55.

[38] Mao JJ, Fang SC, Liu LX. The effects of urbanization on human-viper conflicts and the natural history and impacts of removal operations for Protobothrops mucrosquamatus in Northeast Taiwan. Front Amphib Reptile Sci 2025; 3: 1609575.

[39] Yang CK, Yang YJ, Mori A. Ambush site selection by a green bamboo pit viper: Relation to prey abundance and comparison between juveniles and adults. Zool Stud 2024; 63: e55.

[40] Brudvig LA, Damschen EI, Tewksbury JJ, Haddad NM, Levey DJ. Landscape connectivity promotes plant biodiversity spillover into non-target habitats. Proc Natl Acad Sci USA 2009; 106(23): 9328-9332.

[41] Shine R. Reproductive strategies in snakes. Proc R Soc Lond B 2003; 270(1519): 995-1004.

[42] Do MS, Son SJ, Jung JH, Lee SC, Choi G, Nam HK. Effects of environmental factors and intraspecific niche overlap on the body and ecological characteristics of red-tongued pit vipers (Gloydius ussuriensis). Sci Rep 2023; 13: 21310.

[43] Tsai TS, Tu MC. Reproductive cycle of male Chinese green tree vipers, Trimeresurus s. stejnegeri, in northern Taiwan. J Herpetol 2000; 34(3): 424-430.

[44] Ji X, Wang ZW. Geographic variation in reproductive traits and trade-offs between size and number of eggs of the Chinese cobra (Naja atra). Biol J Linn Soc 2005; 85(1): 27-40.

[45] Whitaker N, Shankar PG, Whitaker R. Nesting ecology of the King Cobra (Ophiophagus hannah) in India. Hamadryad 2013; 36(2): 101-107.

[46] Adams WJ, Graf EW, Anderson M. Disruptive coloration and binocular disparity: Breaking camouflage. Proc R Soc B 2019; 286(1896): 20182045.

[47] Frynta D, Štolhoferová I, Elmi HSA, Janovcová M, Rudolfová V, Rexová K, et al. Hooding cobras can get ahead of other snakes in the ability to evoke human fear. Sci Nat 2025; 112: 1.