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

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.