Abstract: Objective　To analyze the changes in the suitability of Oncomelania hupensis snail habitats in China and to identify climatic factors that significantly affect the suitability in China from 1990 to 2020, so as to provide insights into formulation of schistosomiasis control measures. Methods　All schistosomiasis⁃endemic villages recorded in the Information Management System for Parasitic Disease Control of Chinese Center for Disease Control and Prevention were selected as study areas, and national O. hupensis snail data in China from 2011 to 2021 were collected, while national O. hupensis snail data in China from 1990 to 1993 were collected from the Annual Report of National Schistosomiasis Epidemiological Observation Points in China. The temperature and cumulative precipitation during the spring active period (between April and June), summer dormancy period (between July and August), autumn active period (between September and October), and hibernation period (between December and February of the next year), as well as the maximum temperature of the warmest month, minimum temperature of the coldest month, temperature annual range, annual precipitation, precipitation of the wettest month, precipitation of the driest month, precipitation seasonality, precipitation of the wettest quarter, precipitation of the driest quarter, precipitation of the warmest quarter, and precipitation of the coldest quarter were calculated using historical monthly meteorological data from the WorldClim database during the same period of O. hupensis snail surveys, and annual mean temperature, isothermality, mean diurnal temperature range, temperature seasonality, mean temperature of the wettest quarter, mean temperature of the driest quarter, mean temperature of the warmest quarter, and mean temperature of the coldest quarter were calculated using daily mean surface air temperature (2 m height) data from the ECMWF Reanalysis version 5 (ERA5) hourly data database. Four topographic variables, namely elevation, slope, terrain and distance to water bodies, were obtained from the Resource and Environmental Science Data Platform of the Chinese Academy of Sciences, and the normalized difference vegetation index (NDVI) was obtained from the Geographic Remote Sensing Ecological Network Platform. Soil characteristic data were obtained from the harmonized world soil database, and the long⁃term series population density variable was estimated by previous studies based on Global Human Settlement Layer⁃Population and WorldPop databases. Data were processed to a 1 km2 resolution via resampling. An ecological niche model for O. hupensis snails was constructed using machine learning models and O. hupensis snail breeding data in China from 2011 to 2021, and the suitability of O. hupensis snail habitats was estimated and compared during the study periods between 2018 and 2020 and between 1990 and 1993, to analyze the changes in the suitability of O. hupensis snail habitats. O. hupensis snail data from both periods were calculated separately and analyzed using five models, including neural networks, support vector machines (SVM), naive Bayes, classification and regression tree (CART), and extreme gradient boosting (XGBoost) models, and the sensitivity, specificity, and area under the receiver operating characteristic (ROC) curve (AUC) of models were calculated. The Shapley Additive Explanations (SHAPs) algorithm was employed to calculate the mean contribution of each variable to the predictions of models. Partial dependence plots (PDPs) were constructed using the pdp package to analyze the response of O. hupensis snail habitat suitability to individual variables. Based on historical climatic data and current non⁃climatic data, a counterfactual scenario map was generated using models trained with 2011—2015 datasets. The counterfactual simulated results were compared with the actual O. hupensis snail habitat suitability from 2018 to 2020 to quantitatively assess the effect of climate changes on the suitability of O. hupensis snail habitats from 1990 to 2020.  Results　The areas of O. hupensis snail habitats decreased from 3 726.641 km2 in 2011 to 3 562.876 km2 in 2015 in China, and then rebounded to 3 692.687 km2 in 2021. The XGBoost model exhibited an AUC of 0.983, a sensitivity of 0.993 and a specificity of 0.915, which was selected for subsequent analyses. Through SHAPs calculation combined with PDPs analysis, the four most significant climatic factors affecting the suitability of O. hupensis snail habitats in China included mean temperature of the coldest quarter [1.489 ℃, 95% confidence interval (CI): (1.105 ℃, 1.875 ℃)], precipitation of the wettest quarter [0.706 mm, 95% CI: (0.275 mm, 1.138 mm)], mean diurnal temperature range [0.610 ℃, 95% CI: (0.423 ℃, 0.797 ℃)], and precipitation of the warmest quarter [0.388 mm, 95% CI: (0.199 mm, 0.577 mm)]. Modeling results indicated that the areas with a reduced suitability of O. hupensis snail habitats accounted for 12.201% [95% CI: (11.588%, 12.830%)] of total areas along the Yangtze River basin from 1990 to 2020, while the areas of expanded O. hupensis snail suitable habitats in eastern coastal and high⁃latitude regions accounted for 10.009% [95% CI: (9.380%, 10.637%)] of total areas of the Yangtze River basin. Substantial expansions of O. hupensis snail suitable habitats were found along the lower reaches of the Yangtze River, with a 68.233% [95% CI: (67.463%, 69.002%)] increase in Zhejiang Province and 57.648% [95% CI: (56.878%, 58.417%)] increase in Jiangsu Province, respectively. The middle reaches of the Yangtze River experienced an overall reduction in the suitability of O. hupensis snail habitats, with a 45.784% [95% CI: (45.015%, 46.554%)] reduction in Anhui Province and 33.307% [95% CI: (32.538%, 34.077%)] reduction in Hubei Province, respectively. The areas with an increased suitability of O. hupensis snail habitats accounted for 2.461% to 6.166% of the total provincial areas in two provinces of Sichuan and Yunnan in the upper reaches of the Yangtze River, and in Guangdong Province, Guangxi Zhuang Autonomous Region and Fujian Province, respectively, while the areas with a decreased suitability accounted for 0.890% to 11.891% of the total provincial areas, respectively. In addition, counterfactual analysis revealed the areas with an increased suitability of O. hupensis snail habitats due to climate changes along the Yangtze River basin accounted for 20.738% [95% CI: (19.968%, 21.507%)] of the total basin areas, and the area with a reduced suitability accounted for 4.678% [95% CI: (3.909%, 5.448%)]. Conclusions　The areas with a decreased and increased suitability for O. hupensis snail habitats were roughly equal in size in China from 1990 to 2020, with the areas with a decreased suitability mainly concentrated in the middle reaches of the Yangtze River region, and the areas with an increased suitability mainly located in the eastern coastal areas and high⁃latitude regions. Climate change is the primary factor driving the increased suitability of O. hupensis snail habitats in China.

Key words: Oncomelania hupensis, Climate change, Habitat suitability, Machine learning

摘要： 目的　分析1990—2020年我国钉螺生境适宜度变化情况并识别对其影响显著的气候因素，为制定血吸虫病防控措施提供参考。方法　以中国疾病预防控制中心寄生虫病防治信息管理系统记录的全部血吸虫病流行村作为研究区域，收集2011—2021年全国螺情数据；于《全国血吸虫病流行病学观测点年报》收集1990—1993年全国螺情数据。另外，利用与钉螺调查数据相同时间段的WorldClim数据库历史逐月气象数据计算研究区域在钉螺春季活跃期（4—6月）、夏蜇期（7—8月）、秋季活跃期（9—10月）、冬眠期（12月—次年2月）的气温和累积降水量及最暖月最高温、最冷月最低温、年温度范围、年降水量、最湿月降水量、最干月降水量、降水季节性、最湿季降水量、最干季降水量、最暖季降水量和最冷季降水量；基于ERA5小时级数据库中每日2 m气温数据计算年均温、等温性、平均日温差、温度季节性、最湿季均温、最干季均温、最暖季均温和最冷季均温；自中国科学院资源与环境科学数据中心获取海拔、坡度、地形和水体距离4个地形变量；自地理遥感生态网获取归一化植被指数（normalized difference vegetation index，NDVI）；自世界土壤协调数据库获取土壤特质相关数据；基于全球人类居住区人口数据集、WorldPop全球人口数据集等数据库获取长时间序列人口密度变量估计值。以重采样方式将上述数据处理为1 km2分辨率数据。基于机器学习模型及2011—2021年我国钉螺孳生数据构建钉螺生态位模型，估计并对比2018—2020年和1990—1993年2个时段研究区域的钉螺生境适宜度，分析生境适宜度的变化趋势。分别计算2018—2020年和1990—1993年数据，采用神经网络、支持向量机（support vector machine，SVM）、朴素贝叶斯、决策树（classification and regression tree，CART）和极端梯度提升（extreme gradient boosting，XGBoost）5种模型进行建模分析，计算模型灵敏度、特异度和受试者工作特征曲线下面积（area under the curve，AUC），并选择最优模型进行分析。计算沙普利加性解释（Shapley additive explanations，SHAPs），估计各变量对模型预测结果的平均贡献度；采用pdp包构建局部依赖模型（partial dependence plots，PDPs），分析钉螺生境适宜性对单变量的响应。基于历史气候数据和当前非气候数据，采用2011—2015年数据训练所得模型生成反事实情景地图。比较反事实模拟结果与2018—2020年实际钉螺生境适宜度，定量评估气候变化在1990—2020年对钉螺生境适宜度的影响。结果　我国有螺面积从2011年的3 726.641 km2降至2015年的3 562.876 km2，随后又回升至2021年的3 692.687 km2。XGBoost模型的AUC值为0.983，灵敏度和特异度分别为0.993和0.915，选择该模型用于后续分析。通过计算SHAPs并结合PDPs分析发现，影响我国钉螺生境适宜度的4个最重要气候因素依次为最冷季均温［1.489 ℃，95%可信区间（confidence interval，CI） ：（1.105 ℃，1.875 ℃）］、最湿季降水量［0.706 mm，95% CI：（0.275 mm，1.138 mm）］、日温差均值［0.610 ℃，95% CI：（0.423 ℃，0.797 ℃）］和最暖季降水量［0.388 mm，95% CI：（0.199 mm，0.577 mm）］。模型分析结果表明，1990—2020年长江流域生境适宜度下降面积占长江流域总面积的12.201%［95% CI：（11.588%，12.830%）］，东部沿海和高纬度地区钉螺适宜生境扩张区域占长江流域面积的10.009%［95% CI：（9.380%，10.637%）］。长江下游地区钉螺适宜生境增幅较大，浙江省和江苏省适宜度上升面积分别占全省总面积的68.233%［95% CI：（67.463%，69.002%）］和57.648%［95% CI：（56.878%，58.417%）］。长江中游地区钉螺生境适宜度整体降低，其中安徽省和湖北省钉螺适宜度降低的面积分别占全省面积的45.784%［95% CI：（45.015%，46.554%）］和33.307%［95% CI：（32.538%，34.077%）］。长江上游四川省、云南省以及广东省、广西壮族自治区和福建省钉螺生境适宜度上升面积占全省总面积的2.461% ~ 6.166%，下降面积占全省总面积的0.890% ~ 11.891%，变化幅度相对较小。通过反事实分析发现，长江流域因气候变化导致的钉螺生境适宜度增高的区域占流域总面积的20.738%［95% CI：（19.968%，21.507%）］，适宜度下降的区域占4.678%［95% CI：（3.909%，5.448%）］。结论　1990—2020年我国钉螺生境适宜度下降和上升区域面积相当，其中适宜度下降区域主要集中于长江流域中游地区，适宜度上升区域主要位于东部沿海及高纬度地区。气候变化是驱动我国钉螺生境适宜度增高的主要因素。

关键词: 钉螺, 气候变化, 生境适宜度, 机器学习