2024年通过一项密集的短期调查，揭示了日本北海道中部地区蜱虫及其传播病原体的空间分布情况以及相关的生态因素

《PLOS One》：Spatial distribution of ticks and tick-borne pathogens in central Hokkaido, Japan and associated ecological factors revealed by intensive short-term survey in 2024

　　摘要：蜱传病原体通过蜱虫叮咬传播，导致人类和家养动物感染疾病。为了预测蜱传疾病的发生，需要高分辨率地了解感染风险分布及其生态驱动因素。我们的目标是绘制日本北海道中部蜱虫和病原体感染蜱虫的空间分布图。2024年5月13日至6月26日期间，在171个地点持续采样收集成年蜱虫和若虫蜱虫，并对蜱传病原体进行检测，随后进行生态位建模。共收集到9,908只蜱虫，主要在硬蜱属（Ixodes spp.）蜱虫中检测到了地方性蜱传病原体（即蜱传脑炎病毒、礼佐病毒、北地纳伊罗病毒、莱姆病群螺旋体及复发性发热群螺旋体）。根据蜱虫采集和病原体检测数据，预测了蜱虫和病原体的潜在适宜栖息地。模型具有可接受的预测性能（留一法交叉验证的中位AUC=0.83，中位TSS=0.59）。硬蜱（Ixodes persulcatus）和卵蜱（Ixodes ovatus）被确定为决定所有病原体分布的主要蜱虫种类。此外，不同病原体和蜱虫物种的适宜栖息地存在差异。在考虑的环境变量中，积雪深度对蜱虫和病原体分布差异有显著影响。本研究的结果扩展了我们对蜱传病原体感染空间风险分布及其生态背景的理解。

引言：与硬蜱（蜱螨目，Ixodidae）相关的蜱传病原体通常由蜱虫和野生动物在野外维持，并通过蜱虫叮咬感染人类和家养动物，从而引发疾病[1–3]。硬蜱是蜱传病原体的宿主，能够长期携带这些病原体并将其传播给动物[4,5]；病原体可以在蜱虫的多个发育阶段存活（即跨阶段传播），部分病原体还可以通过卵传播（即经卵传播）。未感染的蜱虫通过吸食受感染动物的血液或同时吸食同一宿主的受感染蜱虫而感染病原体[4,5,6]。在蜱传病原体的自然生命周期中，野生动物不仅作为病原体的放大宿主促进其传播，也是蜱虫的宿主[7]。因此，了解病原体、蜱虫和野生动物之间的关联对于阐明蜱传病原体在自然界的循环及预测蜱传疾病的发生至关重要。每种蜱传病原体都有其特定的媒介蜱虫种类和偏好动物宿主种类，这些宿主能够维持、放大并传播病原体。例如，蜱传脑炎病毒（TBEV）主要由硬蜱属（Ixodes ricinus）蜱虫传播，小型啮齿动物是TBEV循环的重要宿主[1]。虽然特定的媒介蜱虫和动物宿主种类与蜱传病原体的生命周期直接相关，但某些蜱虫与非宿主动物之间的关联也不容忽视[7]。蜱虫-宿主关联对蜱传病原体的传播尤为重要；具有长距离扩散能力的野生动物（如鸟类）可以将受感染的蜱虫带到其非地方性分布区域[8,9]。此外，由于蜱虫和野生动物的生存依赖于环境资源，蜱传病原体的循环间接受到地形、气候和地貌等环境因素的影响[10,11]。这些病原体、蜱虫、野生动物和环境之间的直接和间接相互作用影响着蜱传疾病的发生。

北海道是日本最北端的岛屿县，其气候和生物特征与本州岛不同，并与欧亚大陆东北部共享蜱虫和蜱传病原体。在该地区，五种与硬蜱相关的蜱传病原体已导致或可能导致人类疾病：蜱传脑炎病毒（TBEV）、礼佐病毒（YEZV）、北地纳伊罗病毒（BJNV）、莱姆病群螺旋体（LDB）和复发性发热群螺旋体（RFB，即Borrelia miyamotoi）[1,12–19]。蜱传脑炎、莱姆病和复发性发热在日本属于法定报告疾病，2024年北海道分别报告了2例、18例和11例病例[20]。自2021年首次报告以来，北海道已发现22例由礼佐病毒引起的新型蜱传病毒感染[13,15,16,21,22]。日本尚未报告北地纳伊罗病毒感染的人类病例，但该病毒在北海道的硬蜱中广泛存在[23]。虽然北海道几乎所有蜱虫叮咬病例均由硬蜱（Ixodes persulcatus）和卵蜱（Ixodes ovatus）引起[24,25]，但也有其他蜱虫种类存在，例如传播严重发热伴血小板减少综合征的Haemaphysalis longicornis[27]。尽管先前的研究为这些病原体和硬蜱提供了有价值的见解，但其地理覆盖范围仍不完整[12,23,28–30]。因此，仍缺乏对其空间分布的高分辨率了解，这阻碍了对该地区蜱传疾病生态学和有效风险评估的理解。为填补这一知识空白，需要跨蜱虫/病原体种类进行空间上全面且可比较的分析。本研究旨在探讨北海道中部（日本人口最密集的地区之一）蜱虫和蜱传病原体分布的环境决定因素。基于1.5个月内171个地点持续采样获得的蜱虫和病原体数据，以及五种病原体的检测结果，我们使用景观、气候、地形和野生哺乳动物分布作为解释变量，进行了生态位建模（ENM）[31]，以预测蜱虫和蜱传病原体的潜在适宜栖息地。此外，通过量化各物种的空间聚集度和蜱虫与病原体的生态位相似性，研究了蜱虫/病原体物种之间的分布差异。我们的发现有助于识别与蜱虫和蜱传病原体空间分布相关的生态因素，并为研究区域内的蜱传疾病预防措施提供参考。

方法：
**蜱虫采样**：研究区域位于日本北海道中部（44°00.0′N – 42°63.3′N，142°02.5′E – 140°91.25′E），东西长约124公里，南北长约209公里，以札幌市为中心（人口约200万[https://www.city.sapporo.jp/toukei/jinko/jinko.html，2025年2月5日访问]。札幌市周围的平原地区被三个山区环绕，本研究分别称为北部山区、南部山区和东部山区（图1）。下载：PNG（大图）/ TIFF（原始图）。

**采样方法**：2024年，在随机选定的171个地点（图1右侧面板中的黑色点）使用0.70米×1.00米的法兰绒布标记法收集成年蜱虫和若虫蜱虫。为确认采样点在研究区域内随机分布，使用ArcGIS pro ver. 3.3.2进行了平均最近邻分析[32]。采样时间为5月13日至6月26日，与研究区域的蜱虫物候期一致[24,28]。所有采样均在无降水日的白天进行，且每次采样由第一作者（M.I.）单独完成，每次30分钟，以保持采样工作的稳定性。选择基于时间的采样方法是因为它有助于在不同植被高度和密度的栖息地中标准化采样工作。所有采样地点均为公共可进入区域，无需特殊许可，且本研究未涉及濒危或受保护物种。成年蜱虫和若虫蜱虫的物种及性别通过显微镜根据形态特征进行鉴定[26,33]。成年蜱虫按物种、性别和地点分组，若虫蜱虫按属和地点分组。运输过程中死亡的蜱虫未被使用。每个样本池被转移到tacoTM Pre-loaded Steel Bead Tube（GeneReach Biotechnology Corp., 台中市）中，然后用蒸馏水清洗一次。样本池在-80°C下保存，直至进行核酸提取。

**核酸提取和病原体检测**：使用blackPREP蜱虫DNA/RNA Kit（Analytikjena, 德国耶拿）按照制造商说明提取DNA和RNA，但需用Micro Smash MS-100R（TOMY, 东京）以3,000 rpm转速、30秒的时间和4°C的温度将蜱虫匀浆三次。提取的DNA和RNA在-80°C下保存，以待进一步实验。使用逆转录定量实时聚合酶链反应（RT-qPCR）或qPCR在qTOEWR3 G Real-Time PCR Thermal Cycler（Analytikjena, 德国耶拿）上检测研究区域内的五种地方性蜱传病原体基因组。TBEV RNA的检测方法参照Achazi等人（2011）[34]，使用One Step PrimeScript III RT-qPCR mix（Takara Bio Inc., 滋贺县），引物1：5′–TGGAYTTYAGACAGGAAYCAACACA–3′，引物2：5′–TCCAGAGACTYTGRTCDGTGTGGA–3′，探针：5′–FAM–CCCATCACTCCWGTGTCAC–MGB–NFQ–3′，热循环条件为：50°C 15分钟，95°C 2分钟，随后45个循环，每个循环95°C 15秒，60°C 34秒。YEZV RNA的检测方法同样使用One Step PrimeScript III RT-qPCR mix，引物1：5′–CATACAGGAAGGCCATCTCATT-3′，引物2：5′–AGCCCTTGACACTGCATTT–3′，探针：5′–/56–FAM/ACCTACTAC/ZEN/TGGATGTGGAAGGCAGA/3IABkFQ/–3′，热循环条件为：52°C 5分钟，95°C 10秒，随后40个循环，每个循环95°C 5秒，60°C 30秒。BJNV RNA的检测使用One Step PrimeScript III RT-qPCR mix，引物1：5′–ATGGCTCATGAACGCAAC-3′，引物2：5′–GGTYTTGATGACYTCTGGGTC–3′，探针：5′–/56-FAM/AAAGCTCTGACCGTTTACCCAGCC/3IABkFQ/–3′，热循环条件为：52°C 5分钟，95°C 10秒，随后40个循环，每个循环95°C 5秒，60°C 30秒。LDB和RFB基因组的检测采用多重qPCR方法，参考Barbour等人（2009）[35]和Takano等人（2014）[12]的方法，使用Premix Ex Taq（Probe qPCR）（Takara），引物1：5′–GCTGTAAACGATGCACACTTGGT–3′，引物2：5′–GGCGGCACACTTAACACGTTAG–3′，LDB探针：5′–6FAM–TTCGGTACTAACTTTTAGTTAA–3′，RFB探针：5′–VIC–CGGTACTAACCTTTCGATTA–3′，目标基因为16S rRNA，热循环条件为：95°C 20秒，随后45个循环，每个循环95°C 5秒，60°C 30秒。所有实验均进行技术重复，包括阳性对照和无模板阴性对照。在两次重复实验中均显示指数扩增曲线的样本（Ct值<45）以及在相应阴性对照中无扩增的样本被视为阳性。对于仅在一次重复实验中显示指数扩增的样本，需进行技术重复实验；只有在两次实验中均显示可重复扩增的样本才被严格归类为阳性。使用R中的PooledInfRate包根据样本池的基因组检测结果，估计每种病原体的阳性率（%）及其95%置信区间（95% CI）[36]。LDB也被称为Borrelia burgdorferi sensu lato（s.l.），它由多种蜱虫种类携带[2,3]。在我们的研究区域，LDB基因型与蜱虫种类之间的关联已通过多项研究得到证实[2,5,18,37]。此外，在日本，大多数人类莱姆病病例是由I. persulcatus携带的LDB基因型（即Borrelia garinii和Borrelia afzelii）引起的[5,18,37]。因此，在后续步骤中，从I. persulcatus中检测到的LDB被称为pLDB，并与总LDB分开进行分析。为了使用蜱虫收集数据和病原体检测数据进行后续的空间分析和建模，每个调查地点的数据被转换为存在或不存在的状态：如果收集到至少1只成年蜱虫，则认为蜱虫存在；如果至少有一个阳性样本池，则认为病原体存在；所有其他地点都被视为不存在蜱虫和病原体的地点。为了确保形态鉴定的最高准确性，仅使用成年蜱虫进行这一定义。

为了研究每种蜱虫和病原体的空间聚集情况，使用R语言中的spdep包在正空间自相关的替代假设下进行了Moran’s I检验[38]。Moran’s I的计算方法如下[39]：其中，表示在地点收集到的每种蜱虫成虫的数量，表示病原体的存在（1）/不存在（0），表示调查地点的总数，在本研究中为171个。是一个行标准化的空间权重矩阵的元素，定义为：其中表示地点之间的距离（公里）。控制距离衰减的参数（r）从0.1变化到20，以检查不同空间尺度上的空间自相关性（补充S1图）。Moran’s I接近-1意味着相邻地点彼此不同，而Moran’s I接近1意味着物种呈簇状分布。Moran’s I接近0表示随机模式。p值<0.05被认为具有统计学意义。在全球空间分析之后，使用spdep包计算了局部空间统计量G-star，以分析每种蜱虫和病原体的局部聚集情况[38]。我们使用了在每个地点收集到的每种蜱虫成虫的数量或病原体的存在（1）/不存在（0），以及用于Moran’s I的相同空间权重矩阵。使用ggplot2包[40]绘制Z值，以可视化每种物种的热点。

在生态位模型（ENM）中，包括了20个景观、气候、地形和野生哺乳动物分布的变量（表1和补充S2图）。六个景观变量（即森林面积、草地面积、湿地面积、耕地面积、城市面积和水域面积）是从JAXA提供的高分辨率土地利用覆盖图ver. 23.12中提取的，该图的原始分辨率为10×10米（https://www.eorc.jaxa.jp/ALOS/lulc/data/index_j.htm#japan_v23.12，访问日期2024年7月8日）。这些景观类别是根据JAXA的人造卫星观测结果分类的。五个气候变量（即年最大积雪深度、6月的降水量、6月的平均温度、6月的日照小时数和1月的最低温度）和两个地形变量（即平均海拔和最大坡度角）来自日本国土交通省提供的国家土地数值数据，原始分辨率为1×1公里（https://nlftp.mlit.go.jp/ksj/index.html，访问日期2024年7月8日）。这些气候值是基于日本气象厅30年的气象观测数据计算的，地形数据基于日本地理空间信息局创建的数字高程模型（DEM）。哺乳动物的分布信息来自日本环境省提供的鸟类和哺乳动物分布调查（https://www.biodic.go.jp/youchui/youchui_top.html#mainText和https://www.biodic.go.jp/kiso/do_kiso4_mam_f.html#mainText，访问日期均为2024年7月26日），以及北海道政府提供的2020年浣熊捕获信息（https://www.harp.lg.jp/opendata/dataset/1742.html，访问日期2024年7月26日），还有北海道研究组织提供的梅花鹿狩猎信息图（https://www.hro.or.jp/industrial/research/eeg/development/datamap/deermap.html，访问日期2024年7月26日）。这些哺乳动物地图采用相同的5×5公里网格。对于每种哺乳动物，我们将一个或多个数据源视为存在的网格定义为存在网格，而其他网格则定义为不存在网格（棕熊Ursus arctos yesoensis的分布、貉Nyctereutes procyonoides albus的分布、狐狸Vulpes vulpes schrencki的分布、浣熊Procyon lotor的分布以及梅花鹿Cervus nippon yesoensis的分布）。所有气候和地形变量都采用相同的1×1公里网格。景观变量是从10×10米网格聚合到1×1公里网格的。对于哺乳动物变量，每个5×5公里网格的数据被均匀分配到其内的所有25个1×1公里网格中。尽管5公里网格是哺乳动物分布数据的最高分辨率，但统一缩放到1公里可能会忽略细尺度变化。使用ArcGIS Pro ver. 3.3.2 [31]提取解释变量。

生态位模型使用广义线性模型（GLM）进行，分辨率为1×1公里。我们使用R语言中的stats包[41]来拟合GLM。作为响应变量，使用了每种蜱虫/病原体的存在/不存在状态，假设它遵循二项分布（使用logit链接函数）。每个连续解释变量都通过减去其平均值并除以其标准差来进行标准化。为了避免模型复杂性并保持统计功效，没有考虑环境变量之间的交互项。为了变量选择，我们通过两步过程开发了最佳预测模型：首先使用Akaike信息准则（AIC）进行自动选择，然后进行手动选择。作为识别候选变量集的初始步骤，我们使用R语言中的MuMIn包[42]从所有解释变量组合中选择了AIC最佳的模型。从这个AIC最佳模型开始，按照以下描述进行了手动变量选择，同时考虑了预测准确性和模型复杂性。在单次回归分析中，p值<0.05的解释变量被优先考虑纳入模型。使用接收者操作特征曲线下面积（AUC）[43]和真实技能统计量（TSS）[44]（通过留一法交叉验证（LOOCV）[31]计算）来选择要纳入模型的解释变量。AUC的范围是从0到1，AUC>0.7是预测的可接受阈值，TSS的范围是从-1到1，TSS>0.5是预测的可接受阈值。为了避免多重共线性，使用performance包[45]计算了方差膨胀因子（VIFs）。当VIFs超过5时，逐一排除相关的解释变量，并选择具有最佳预测准确性的变量组合（即更高的AUC和TSS）。我们还使用上述的Moran’s I检验探索了每个模型残差中的空间自相关性。当在GLM的残差中检测到空间自相关性时，我们使用R语言中的mgcv包[46]拟合了一个广义加性模型，以对网格坐标（经度lon和纬度lat）进行空间平滑处理。我们考虑了一组基于薄板回归样条的平滑项：经度（s(lon）和纬度（s(lat）的一维平滑，用于大尺度趋势，以及张量积交互（ti(lon, lat）用于更局部的模式。从非空间模型开始，逐步添加这些项，并为每种物种选择最简单的模型，直到残差不再显示显著的空间自相关性。在选择最佳模型后，我们进行了分层5折交叉验证，以评估最佳模型的预测性能。选择这种方法是为了保持每折中存在/不存在的比例平衡。作为性能指标，我们计算了AUC/TSS，并报告了每个蜱虫/病原体在五个折中的平均值和标准差（SD）。最后，使用每个蜱虫/病原体的最佳模型来预测每个网格的存在概率，然后将其映射到研究区域。预测概率使用最大化TSS的阈值转换为二进制格式，以便以存在/不存在格式绘制潜在分布图。对于可视化，使用了R语言中的ggplot2包[40]。

为了量化每种病原体和每种蜱虫之间的空间重叠，计算了Schoener’s D和Warren’s I。这两个指标都衡量了两种物种相对于其可用性的空间利用程度[47–49]。Schoener’s D定义为，而Warren’s I定义为。在上述方程中，表示ENM预测的物种在网格上的存在概率。Schoener’s D和Warren’s I的范围是从0（两种物种之间没有重叠）到1（完全重叠）。为了计算95%置信区间，我们采用了1000次重复的自助法。

**结果**
日本北海道中部的蜱虫和蜱传病原体的流行情况
平均最近邻分析表明，171个调查地点是随机分布的（z分数=?0.502，p值=0.616）。在171个调查地点中有159个地点收集到了一个或多个蜱虫，总共收集到4,608只成年蜱虫和5,300只若虫蜱虫（表2、图2和补充S1表）。成年蜱虫被鉴定为九个物种：I. ovatus、I. persulcatus、Ixodes pavlovskyi、Ixodes tanuki、Haemaphysalis megaspinosa、H. longicornis、Haemaphysalis flava、Haemaphysalis japonica和Haemaphysalis concinna。I. ovatus和I. persulcatus是最常见的物种，分别占收集到的成年蜱虫的82%和11%。每种I. tanuki和H. concinna只收集到一只成年蜱虫；因此，这些蜱虫物种被排除在空间分析之外。总共9,745只活蜱虫被汇集到1,026个样本中，用于RT-qPCR和qPCR（表2）。TBEV、YEZV、BJNV、LDB和RFB分别在7个（平均Ct值：23.4–38.3）、8个（平均Ct值：24.2–36.2）、11个（平均Ct值：15.6–26.9）、476个（平均Ct值：30.6–42.0）和16个（平均Ct值：26.8–38.4）样本中被检测到（补充S2表）。我们的筛选方法及其额外分子确认（包括测序）的详细理由和结果在补充S1 Text1中提供。所有五种病原体基本上都在Ixodes spp.蜱虫中被检测到。TBEV和YEZV在I. ovatus和I. persulcatus中的流行率相似。BJNV几乎只在I. persulcatus中被检测到，其比例显著更高（1.92，95% CI：0.95–3.50），而I. ovatus的比例为0.03，95% CI：0.00–0.13。LDB在I. ovatus、I. persulcatus和I. pavlovskyi中的流行率相当。由于北海道大多数莱姆病病例是由I. persulcatus携带的LDB基因型引起的[2,18,36]，因此在后续步骤中，从I. persulcatus中检测到的LDB被单独作为pLDB与总LDB分开分析。至于RFB，在I. persulcatus中的估计阳性率（1.27，95% CI：0.52–2.64）和I. pavlovskyi（6.53，95% CI：0.39–27.28）显著高于I. ovatus（0.19，95% CI：0.08–0.37）。

**表2. 收集到的蜱虫和蜱虫中检测到的病原体。**
**图2. 蜱虫数量和病原体存在情况。**
对于用于空间分析的七种蜱虫，每个调查地点收集到的成年蜱虫数量用圆圈的大小表示。每种病原体的存在和不存在地点分别用橙色和蓝色点表示。右下角的卫星图像也显示在图1和图3中。底图图像使用Esri的ArcGIS底图提供。内容是Esri的知识产权，并经许可在此使用。版权所有 ? 2025 Esri及其许可方。所有权利保留。
蜱虫收集和病原体检测的结果被转换为存在/不存在数据（补充S3表）。值得注意的是，尽管收集到的H. flava的数量低于H. longicornis（表2），但H. flava的存在地点数量高于H. longicornis。为了研究蜱虫和病原体分布的空间聚类性，我们使用了Moran’s I检验方法，该方法基于收集到的每种蜱虫成虫的数量以及每种病原体的存在/缺失数据（见表3和补充表S4）。无论距离衰减参数（r）如何，除了H. flava之外，所有蜱虫物种的Moran’s I值均为显著正数。当距离衰减较强时，H. flava的Moran’s I值在统计上并不显著。除了YEZV之外，所有病原体物种的Moran’s I值在强或弱参数设置下均为显著正数。在任何参数设置下，均未检测到YEZV的空间聚类性（见补充表S4）。下载：PNG（较大图像）/ TIFF（原始图像）。

表3. 使用Moran’s I检验研究每种蜱虫/病原体物种的空间聚类性结果。https://doi.org/10.1371/journal.pone.0349386.t003

Moran’s I值和p值是使用R语言中的spdep包计算的[38]，其中参数用于控制距离衰减（分别代表强或弱的距离衰减）。完整结果见补充表S4。为了揭示每种蜱虫和病原体的热点区域，计算并可视化了局部空间统计量G-star（见补充表S3）。分析结果与图2中绘制的原始蜱虫采集和病原体检测数据基本一致。Ixodes spp.蜱虫及其相关病原体（如LDB）在北部和南部山区表现出明显的聚集现象。我们发现，H. flava和YEZV的分布相对分散，但它们在局部仍呈现聚集特征。

通过ENM（生态网络模型）的LOOCV（Leave-One-Out Cross-Validation）方法评估了模型性能和解释变量。对于大多数蜱虫/病原体物种，模型预测的准确性达到了可接受的标准（AUC > 0.7且TSS > 0.5）。然而，H. flava的模型预测性能较差（AUC = 0.71且TSS = 0.41，未达到TSS阈值，见补充表S5）。总体而言，AUC的中位数为0.83（四分位数范围：0.76–0.86），TSS的中位数为0.59（四分位数范围：0.53–0.64）。此外，还进行了分层5折交叉验证以评估最佳模型的预测稳定性。除了H. flava之外，所有物种的AUC和TSS均值均高于可接受阈值（AUC > 0.7且TSS > 0.5），且标准差较小（见补充表S5）。H. flava的AUC和TSS低于可接受阈值。最佳模型中包含的解释变量见图3和补充表S6。年度最大积雪深度（snow）被纳入所有蜱虫物种的最佳模型中，对Ixodes蜱虫具有正效应，而对Haemaphysalis蜱虫具有负效应。海拔高度（elev）在蜱虫和病原体的最佳模型中均被纳入，且具有正效应。在景观变量中，森林面积（forest）的大小具有正效应，被纳入Ixodes ovatus、I. persulcatus、H. flava、YEZV和LDB的最佳模型中；而草地面积（grass）的大小具有负效应，被纳入I. pavlovskyi、H. megaspinosa和pLDB的最佳模型中。与景观、气候和地形变量相比，哺乳动物的分布较少被纳入最佳模型中。对于I. pavlovskyi、H. megaspinosa、H. longicornis、H. flava、H. japonica和pLDB的模型，应用了空间平滑处理（见补充表S6）。最终模型的残差中未检测到显著的空间自相关性（Moran’s I检验，p > 0.05，见补充表S5）。

下载：PNG（较大图像）/ TIFF（原始图像）

图3. 最佳模型中包含的解释变量。最佳模型中每个解释变量的系数大小用点的大小表示。所有连续解释变量均已标准化。各系数的p值表示如下：*** < 0.001；** < 0.01；* < 0.05。GLM的详细输出见补充表S6。https://doi.org/10.1371/journal.pone.0349386.g003

图4. 最佳模型预测的每种蜱虫/病原体物种在每个网格上的存在概率。预测的适宜栖息地以二进制（存在/缺失）格式绘制在图4中。Ixodes蜱虫的适宜栖息地覆盖了整个山区，除了研究区域的南部。I. persulcatus的适宜栖息地完全包含在I. ovatus的适宜栖息地范围内。Haemaphysalis蜱虫的适宜栖息地主要集中在研究区域的南部，而H. flava的适宜栖息地则较为分散。尽管某些北部山区存在多个H. japonica缺失样本，但这些样本仍被预测为高存在概率区域（见图2和图4）。高山地区（参见补充表S2中的海拔高度）通常被预测为病原体的适宜栖息地。LDB的适宜栖息地与I. ovatus和I. persulcatus有大量重叠，而某些I. persulcatus的适宜栖息地（即北部山区的一些网格）则不适合pLDB。南部山区的大片区域被预测为RFB的适宜栖息地。TBEV的适宜栖息地主要分布在北部山区。BJNV的适宜栖息地分布在北部和南部山区。YEZV的适宜栖息地分布较为分散。

下载：PNG（较大图像）/ TIFF（原始图像）

表4. 每种病原体与蜱虫之间的Schoener’s D和Warren’s I值。https://doi.org/10.1371/journal.pone.0349386.t004

讨论
了解蜱传病原体的分布及其决定因素对公共卫生保护至关重要。许多研究试图通过长期收集的人类疾病发生数据[50]、媒介蜱虫记录[51]以及野生动物中的特异性抗体[52]来揭示病原体感染的风险分布。尽管感染病原体的蜱虫的分布对于理解感染潜力也很重要[53]，但由于蜱传病原体（尤其是蜱传病毒）在蜱虫中的低流行率，空间建模面临重大挑战[51]。在本研究中，通过短期密集采集蜱虫并随后进行病原体分子检测，我们获得了足够的蜱虫和蜱虫内病原体的可比数据集，从而能够预测它们的适宜栖息地。有趣的是，蜱传病原体的预测适宜栖息地与其媒介蜱虫物种的栖息地不同；除了LDB与I. ovatus/I. persulcatus之间的Schoener’s D和Warren’s I值接近1之外，其他情况下这些值均不接近1。我们的发现表明，媒介蜱虫物种的分布仅部分反映了感染病原体的蜱虫的分布，应在每个地区全面调查感染地方性病原体的蜱虫的局部分布。通过考虑十八个环境变量，我们的生态网络模型显示，景观、气候和地形变量影响了多种蜱虫和病原体的分布。其中一个气候变量——年度积雪深度对蜱虫的分布有显著影响。积雪深度对Ixodes蜱虫的正面影响可能源于积雪增加了它们的越冬存活率，因为积雪可以保护它们免受极端寒冷事件的影响[54]。此外，先前的日本研究报告指出积雪深度对Haemaphysalis蜱虫的生命周期有负面影响[28,55]；积雪深度减少了作为Haemaphysalis spp.蜱虫重要宿主的梅花鹿的数量[7,56]。相比之下，年度积雪深度未被纳入TBEV、YEZV、BJNV和pLDB的最佳模型中，在TBEV、YEZV和RFB的单变量回归分析中也不显著（数据未显示）。这些发现表明，积雪深度可能与蜱虫和病原体在研究区域内的分布差异有关。在考虑的因素中，野生哺乳动物的分布对蜱虫或蜱传病原体的影响较小。这可能是由于我们的采样数据与环境预测因子在空间和时间上的不匹配所致。研究区域内公开可用的野生哺乳动物分布数据（以5×5公里网格的形式记录存在/缺失情况）可能无法捕捉到哺乳动物密度的细微变化。为了进一步研究哺乳动物物种与蜱虫/病原体物种之间的相互作用，需要最新的高分辨率密度数据。在我们的短期蜱虫采集中，Ixodes spp.是北海道中部采样点中最常见的蜱虫物种，并且与Haemaphysalis蜱虫相比，更频繁地检测到其感染病原体。此外，Schoener’s D和Warren’s I的结果表明，所有病原体的空间分布主要与I. persulcatus和I. ovatus重叠。在北海道，人类蜱虫叮咬主要由成年I. persulcatus引起[24,25]，这也出现在感染YEZV和RFB的人类患者身上[15,57]。我们的发现与先前的研究一致，表明成年I. persulcatus是该地区向人类传播蜱传病原体的主要媒介。此外，我们的结果还表明，野外I. ovatus的数量比I. persulcatus更多，因此感染TBEV或YEZV的I. ovatus也更多。实际上，从北海道蜱虫样本中分离出的所有TBEV都来自I. ovatus[30,58,59]。因此，为了防止疾病的发生和病原体的进一步传播，不仅需要密切监测I. persulcatus，还需要监测I. ovatus。Moran’s I检验和局部空间统计量表明，几乎所有蜱虫/病原体物种的分布都在局部到区域尺度（半径数十公里范围内）呈现空间聚类；而H. flava和YEZV则被认为与其他蜱虫或病原体物种相比空间分散程度更高；这一结果也在生态网络模型预测的潜在分布中得到验证。通常，由于扩散能力有限，物种在空间上呈现聚类。蜱虫和蜱传病原体的扩散几乎完全依赖于野生动物，因此宿主偏好对其分布模式有很大影响。实际上，H. flava从鸟类身上的采集频率高于其他Haemaphysalis物种[8]。与H. flava不同，H. japonica似乎在适宜栖息地之间的扩散受到限制，因为在一些预测为其适宜栖息地的地区并未采集到这种蜱虫（例如北部山区的西南沿海地区）。与H. flava类似，YEZV的散布分布可以通过鸟类频繁引入来解释。Nishino等人（2024）报告称，在从北海道迁徙鸟类中采集的2,534只蜱虫中检测到8株YEZV[9]，而在同一样本中未检测到TBEV（作者个人沟通）。此外，有报道称TBEV的微焦点在数十年内空间稳定[53]，而YEZV的微焦点则在短时间内出现[29]。尽管需要更多关于它们对鸟类偏好的定量研究，但我们的结果支持YEZV的感染风险在空间上比TBEV更为零星。本研究强调了短期密集采集蜱虫并随后进行生态位建模的有效性。然而，我们的研究也存在一些局限性。首先，每次采样点仅进行30分钟的标记采集可能会遗漏蜱虫和病原体的真实存在情况，尽管采样时间（2024年5月13日至6月26日）是根据先前研究中主要蜱虫物种的物候数据设定的[24,28]。其次，虽然RT-qPCR和qPCR检测的Ct值临界值（<45）旨在最大化检测灵敏度，但这也引入了假阳性的可能性。我们无法完全排除假阳性的可能性；在不同病原体物种中，确认性测序的成功率仅部分达到预期。第三，为了有效利用某些病原体物种的少量样本（这些病原体的流行率极低，例如TBEV和YEZV），在建模中选择了LOOCV和分层5折交叉验证方法；这种选择可能被认为在空间预测上过于乐观。在未来的研究中，需要实施空间分块的交叉验证并扩展数据集，同时整合其他相关数据（例如疾病发生情况），以减少空间预测的乐观性并建立更可靠的分布预测[53]。尽管存在这些局限性，我们相信本研究的结果对于蜱虫及其传播病原体的生态学研究至关重要，同时也有助于制定针对研究区域内蜱传疾病的对策。

结论：在2024年的短时间内（大约1.5个月），我们在北海道中部的171个地点进行了蜱虫采集，随后对这些地方性的蜱传病原体进行了筛查和空间分析。研究结果揭示了该地区七种蜱虫和五种蜱传病原体的潜在分布情况。Ixodes persulcatus和I. ovatus被确定为决定病原体分布的主要蜱虫种类。然而，除LDB外，病原体的空间分布与蜱虫的分布并不完全一致。积雪深度似乎是导致蜱虫和病原体适宜栖息地差异的原因之一。鸟类频繁的传播可能使得H. flava和YEZV的分布比其他蜱虫和病原体更为分散。尽管由于病毒流行率较低以及模型验证过程中的固有局限性，这些预测结果需要谨慎解读，但我们的研究结果可以帮助确定需要进一步监测的重点区域。这些发现扩展了我们对蜱虫及其传播病原体的空间分布及其相关生态因素的理解，突显了在短时间内进行密集蜱虫采样的有效性。

支持信息与致谢：我们感谢Kawabata博士、Takano博士和Kobayashi博士提供TBEV、LDB和RFB的qPCR阳性对照样本。同时，我们也感谢Nakao博士、Isoda博士、Itakura博士和Tabata博士对这项研究的支持。