中西部黑腿蜱种群扩张的遗传与景观连通性研究：对莱姆病传播的启示

《Ecology and Evolution》：Genetic and Landscape Connectivity of Blacklegged Ticks During Range Expansion in Select States of the Midwestern USA

本研究通过精细空间种群基因组学分析，揭示了美国中西部黑腿蜱（Ixodes scapularis）在范围扩张过程中的遗传多样性与景观连通性模式。研究发现威斯康星州北部种群作为扩散源，存在东西向基因流动走廊；而密歇根州种群与周边地区存在显著遗传分化。景观遗传模型（IBR）表明夏季土壤湿度高、人类干扰低的区域更利于基因流动。该研究为蜱媒疾病防控提供了重要的空间遗传学依据。

　　

ABSTRACT

自20世纪70年代以来，美国中西部地区经历了黑腿蜱（Ixodes scapularis）的范围扩张，这种蜱虫是莱姆病病原体伯氏疏螺旋体（Borrelia burgdorferi）的主要传播媒介，导致莱姆病发病率上升。公共卫生监测表明，威斯康星州西北部是这些蜱虫的避难所，并向密歇根等邻近州扩散。然而，重新出现和入侵的过程仍不清楚。为了改善蜱虫管理，我们研究了中西部（东部北中部）地区的环境变量是否限制了蜱的扩散，以及是否可以识别连通性走廊。通过开发精细尺度的空间种群基因组数据，我们的分析揭示了威斯康星州遗传多样性丰富的种群，其中北部种群促进了该州最近的扩张。我们确定了几个促进威斯康星州蜱虫扩散的东西向基因流动走廊。沿威斯康星州密西西比河存在一个独立的扩散网络，向南延伸。相比之下，密歇根州的种群与威斯康星州和印第安纳州的种群表现出显著的遗传分化，遗传多样性低，州内基因流动高。我们还发现了连接密歇根半岛区域的高景观连通性，以及密歇根湖南部附近陆地的显著基因流动。地理隔离，以及夏季土壤湿度低和人类干扰高的景观，被发现在整个区域限制了基因流动，尽管这些影响较小。通过认识到景观连通性影响了不同遗传种群的扩散，并针对看似不太有利于蜱虫的景观进行针对性干预，可以加强该地区黑腿蜱的管理。

1 Introduction

自20世纪70年代以来，美国中西部受到北方宿主搜寻型黑腿蜱（Ixodes scapularis）生态型重新出现和范围扩张的影响，增加了地区莱姆病发病率。黑腿蜱，特别是北方生态型，是导致莱姆病的细菌伯氏疏螺旋体（Borrelia burgdorferi sensu stricto）以及几种其他重要人类疾病的主要传播媒介。疾病发病率的增加通常伴随着蜱虫的到来和定殖。

监测数据显示，黑腿蜱从威斯康星州西北部的重新出现始于20世纪70年代和80年代，这是由于重新造林的努力和可能的气候变暖。逐渐地，入侵进程遍及威斯康星州西部，然后继续向东推进，导致在20世纪90年代和21世纪初全州范围内的入侵。自20世纪90年代以来，邻近的密歇根州经历了类似的入侵模式，蜱虫首先在靠近西部的地区（即威斯康星州和印第安纳州）扩张，并逐渐入侵密歇根州中部和东部县。尽管监测数据记录了逐县的范围扩张，有助于病媒控制，但它们可能无法准确描述连通性模式和扩散机制。

一个关键的挑战是描述和识别祖先结构、基因流动以及介导黑腿蜱在中西部景观上扩散的景观决定因素，这将为了解范围扩张动态和日益增加的蜱传疾病风险提供见解。种群和景观基因组学为理解遗传变异与景观之间的关系提供了一个有用的框架。由世代扩散种群驱动，基因组将记录中性变异的模式，主要是基因流动和遗传漂变的产物，考虑到突变罕见且选择针对基因组的微小区域。如果跨区域采样的种群之间的遗传分化程度在空间上相关，则可以揭示范围扩张的趋势和扩散程度。此外，可以研究种群间的中性遗传分化，以识别随时间推移促进或阻碍基因流动的地理和景观因素，例如景观覆盖类型和配置、气候条件和栖息地质量。如果功能性障碍阻碍了病媒的传播，则可以在某些区域实施管理计划，以潜在地限制疾病传播。这些功能性障碍可能存在于未入侵区域或得到增强，以防止疾病出现。

最早关于黑腿蜱遗传分化的研究集中在大陆尺度的多样性模式，并发现北方和南方蜱虫之间存在深刻的遗传分歧。随着莱姆病疫情的恶化，最近的工作重点集中在基因流动和美国东北部蜱虫的传播上。遗传变异和多样性模式被用来阐明东北部的扩散路线，并发现了由近程和远程扩散产生的定向基因流动的证据。已经清楚的是，东北部蜱虫种群的遗传结构遵循一种梯度，这是由于距离隔离（IBD）模式导致附近州之间蜱虫的基因流动受限。东北部的空间入侵模型比来自现有隐性避难所的局部重新出现的替代模型更受支持。虽然现有文献先前已经确定了环境对黑腿蜱栖息地适宜性和丰度的影响，但尚不清楚这些因素是否影响整个景观的遗传连通性。由于北方生态型在东北部和中西部地区都占主导地位，最近的比较基因组学研究开始更详细地调查中西部黑腿蜱的种群遗传学。虽然中西部蜱虫在遗传上与东北部蜱虫相似，但中西部最近的范围扩张事件被认为是独立于东北部的。尽管有大量的研究工作，这些研究中有限的空间采样密度不允许测试中西部地区的详细范围扩张假设和景观连通性。

监测数据表明，威斯康星州西北部的局部避难所可能播种了中西部的扩张；然而，这一假设尚未经过遗传学验证。中西部来自避难所的扩散动态应该在遗传数据中显而易见，这取决于通过宿主网络从避难所扩散的黑腿蜱的数量。生态变异也可能限制该地区蜱虫的定殖、持久性和扩散，例如有利森林的可用性、土壤条件和气候 envelope。当前的分布生态位模型表明，中西部上游的大部分景观适合黑腿蜱的定殖和范围扩张。尽管如此，可以合理地预期这些景观将为该地区的黑腿蜱创造一些空间异质性的基因流动障碍，尽管有限。为了进一步了解扩散、定殖、遗传和环境之间的相互作用，我们在中西部进行了更全面的黑腿蜱空间采样。利用大量个体的种群基因组样本，我们着手检验以下假设：（1）正如监测数据所暗示的那样，威斯康星州黑腿蜱的扩张来自威斯康星州西北部的局部避难所；（2）正如生态研究所暗示的那样，密歇根州最近的扩张最初源于威斯康星州的蜱虫种群。此外，由于蜱虫在中西部地区的快速扩张，我们检验了假设（3）连通性主要由地理距离驱动，而不是主要由景观变异解释。

2 Materials and Methods

2.1 Sample Collection

在2021年至2023年期间，在整个美国中西部地区（在一些文献中称为北中部地区）收集了宿主搜寻型黑腿蜱（若虫和成虫）。采样工作集中在记录了最早重新出现和范围扩张的威斯康星州，以及最近被入侵的密歇根州。同时，我们将采样工作扩展到爱荷华州和威斯康星州边境密西西比河流域的其他中西部蜱虫种群，以及印第安纳州和俄亥俄州更南部的区域。这些种群可能代表了过去未冰川地区的黑腿蜱避难所，或者尽管在1800年代经历了森林砍伐但仍持续存在的避难所。在每个地点通过拖拽法采集多个个体。我们在威斯康星州的21个地点、密歇根州的15个地点以及其他邻近州的6个地点采集了样本。我们将所有42个采样点纳入我们的研究，经过质量控制后，总共517个个体作为最终样本集。这42个种群按地理区域分组以便讨论，即威斯康星州（WI）、密歇根州上半岛（UP MI）、密歇根州下半岛（LP MI）、威斯康星州-爱荷华州边境（WI-IA）、印第安纳州（IN）和俄亥俄州（OH）。

2.2 Genomic Data Processing

黑腿蜱要么在提取前用剃须刀片切碎，要么用液氮用研杵捣碎。然后使用DNeasy Blood and Tissue kits（Qiagen LLC）提取基因组DNA。所有提取的蜱虫都进行了基因型测序（GBS）程序，使用限制性内切酶（PstI/MspI）的双酶切来降低基因组复杂性。黑腿蜱样本单独进行条形码编码，并使用Illumina DNA Prep Kit（Illumina Inc.）准备进行测序。Illumina双末端测序在Novaseq6000和X Plus平台上进行。结果经过交叉检查以确保一致性。DNA输入量低的密歇根州若虫蜱进行了额外的全基因组扩增（REPLI-g）步骤。先前的工作表明这种方法对种群基因组分析是无偏的。为了确保这一点，我们使用两个对照样本（经过和未经过基因组扩增测序）测试了单核苷酸多态性（SNP）呼叫的基因型一致性。

原始序列数据根据条形码使用STACKS 2.64中的PROCESS_RADTAGS进行解复用并分配给个体ID，其中禁用了限制性内切酶检查以保留最多的读数。使用CUTADAPT 3.5去除适配器，同时使用TRIMMOMATIC 0.39去除低质量碱基和短读数。使用FASTQC v0.11.9进行读取质量检查。然后将双末端读段与黑腿蜱参考基因组（NCBI: GCF_016920785.2）使用BWA-MEM 0.7.17-r1188进行比对。比对后的读段使用SAMTOOLS 1.16.1转换为二进制文件格式并进行排序。排序后的比对读段然后使用GATK 4.4.0.0调用SNP，基本上遵循Germline Short Variant Discovery管道中详细说明的最佳实践。通过移除缺失率高的样本（99%）来最小化变异的缺失，使用PLINK2 v2.00a3 SSE4.2测量。使用VQSR知情的硬过滤器按照推荐方案进一步过滤得到的变异。

最后，使用VCFTOOLS 0.1.17进一步处理变异调用格式（VCF）文件，仅保留至少在80%样本中存在的双等位基因SNP。平均变异深度（DP）为29（最小DP≈6，最大DP≈129，标准差SD≈8）。由于连锁不平衡可能会使一些下游分析产生偏差，我们使用PLINK2移除相关的SNP（窗口大小50 kb，步长10，R²阈值0.1）。大多数分析需要零缺失数据，因此使用BEAGLE v5.4估算剩余的缺失变异。我们采用样本的初始主成分分析来识别和移除一个遗传异常值（> 6×标准差）的样本。支持性脚本已在Dryad上提供。我们保留了单例等位基因（约占总变异的27%）用于下游分析，因为它们提供了关于入侵过程的信息，并且在我们研究中质量高（平均DP覆盖率≈29，最小DP≈6，最大DP≈124，SD≈8）。由于分子生态学研究中通常过滤次要等位基因，我们将我们的结果与移除单例（次要等位基因计数；MAC = 2）的变异集进行了比较以进行演示。

2.3 Genetic Distance Statistics and Isolation by Distance

为了表征中西部地区42个种群之间的遗传差异，我们使用R 4.3.1中的ADEGENET 2.1.10、GRAPH4LG 1.8.0和MMOD 1.3.3计算并比较了几个成对遗传距离。具有空间有限基因流动的种群可能表现出IBD模式，即成对遗传距离与地理距离正相关。选择“1 – 共享等位基因比例”（也称为Dps）和线性化Fst进行进一步分析，因为它们是景观遗传学中的标准指标，并且与地理距离表现出高度相关性。我们使用R包VEGAN 2.6-4中的Mantel检验测试了IBD，使用大圆距离（使用R包SF 1.0-16中的st_distance函数计算）作为地理距离的度量，并使用R包STATS 4.3.1中的LM函数生成R²值。计算了遗传距离的多变量Mantel相关图，以获得每个距离类别的IBD强度，并估计遗传相似性减弱和IBD开始的距离阈值。

2.4 Genomic Summary Statistics

基因组多样性指数通常为了解基因流动和种群动态提供见解。计算了预期杂合度、核苷酸多样性、稀有的私有等位基因丰富度和Tajima's D，以总结先前定义的42个种群的遗传变异模式。前两个指数描述了遗传多样性的原始模式，而私有等位基因丰富度衡量了基因流动影响独特遗传变异信号的程度，全基因组范围的Tajima's D指示了最近的种群扩张（如果为负）或收缩（如果为正）。在移除少于八个样本的42个种群的一个子集后计算了私有等位基因。

2.5 Population Structure

为了识别种群祖先的模式，我们采用了主成分分析（PCA）并使用了两种聚类算法，SNMF（稀疏非负矩阵分解）和CONSTRUCT（连续结构）。PCA是一种非参数降维工具，用于可视化个体之间的遗传差异。使用软件PLINK2估算SNP呼叫中缺失的变异并计算主成分特征向量。聚类算法可以识别最佳拟合的推定祖先种群数量（K）以及每个个体不同基因组祖先的比例。使用R包LEA 3.12.2，通过确定具有最低交叉熵的模型来评估K=1和K=7个可能祖先种群簇之间统计上最拟合的模型，然后使用SNMF函数在该模型下估计个体祖先比例。将个体样本的估计祖先在42个采样地点合并，以显示基因组祖先的空间模式。然而，祖先估计可能会被IBD模式混淆，这可能导致祖先种群数量的高估。因此，采用CONSTRUCT作为控制方法，在校正IBD的同时估计祖先比例。这种贝叶斯方法评估空间（IBD校正）和非空间模型（类似SNMF）的预测准确性，其中执行不同K值（1-7）和（非）空间模型的交叉验证以找到最佳模型。交叉验证后，运行最佳模型进行200,000步马尔可夫链蒙特卡洛（MCMC）过程以稳定祖先估计的后验估计，并在两次独立运行中重复此过程。

2.6 Barriers and Connectivity of Gene Flow

由于复杂的基因流动模式可能由于蜱虫从避难所重新出现或快速范围扩张而产生，我们实施了有效迁移表面的快速估计（FEEMS），它可以揭示累积基因流动的空间模式。简而言之，跨离散空间样本的迁移模式在图的边上被估计为迁移表面，以帮助推断时间平均的区域迁移障碍和连通性。FEEMS在惩罚似然框架下实现了高斯马尔可夫随机场模型。在留一法交叉验证框架中评估了作为平滑正则化参数的不同lambda值，并选择产生最小交叉验证误差的最佳lambda值。

遗传图网络是表征基因流动模式的另一种方法，可用于通过使用包含更多近期基因流动信号的汇总统计量Dps来揭示扩散网络。我们使用GRAPH4LG从遗传数据估计了扩散网络。该软件实现了基于规则的最小生成树修剪方法和快速贪婪模块化优化，以生成带有模块的遗传图。节点大小与连接到每个节点的链接的逆权重的平均值成比例地测量，这是比较地点水平扩散的良好代理。

2.7 Landscape Ecological Data Collection and Processing

为了检查景观因素对基因流动的影响，从公共来源收集了假设与蜱虫生态学相关的景观生态数据，包括来自WORLDCLIM的气候数据、人类发展指数、来自SOILGRIDS的土壤数据、来自PRISM的30年正常太阳辐射（水平）以及来自NLCD 2021 USFS Tree Canopy Cover数据集的森林覆盖。首先，移除了高度相关的变量（阈值>0.7），同时保留了与微气候干燥风险相关的太阳辐射变量。经过变量选择后，最终的变量集包括海拔、最干季度降水、年降水量、人类发展指数、0-5厘米深度15年平均土壤粘土和有机质含量、30年正常太阳辐射（水平）和森林覆盖。所有最终变量都使用R包TERRA 1.7.78转换为相同的地理投影（NAD83）、范围（美国中西部）和分辨率（30弧秒或约1平方公里）下的栅格格式。使用SF 1.0-16将采样位置转换为空间对象。具有缺失栅格像元值的空间位置从其最近邻采样。

2.8 Isolation by Resistance (IBR) and Environment (IBE) Analyses

虽然直线地理距离通常用作由于空间有限基因流动导致的遗传分化的预测因子（IBD），但景观矩阵可以通过称为阻力隔离（IBR）的过程限制基因流动。为了测试IBR，将景观特征分配成本，并解释为在种群对之间阻碍基因流动的生物学相关因素。除了IBD和IBR之外，环境隔离（IBE）是第三个过程，其中采样点的栖息地特征可能导致站点之间的遗传分化，例如当相似或不同的站点栖息地条件影响站点之间的扩散率时。

使用所有景观变量通过优化阻力值来建模IBR，使用线性化Fst计算的观察到的遗传距离，这是累积景观对遗传多样性影响的首选指标。将景观数据聚合到约10公里分辨率以进行更有效的优化。使用R包RADISH中的最大似然估计进行IBR建模，使用“最大似然种群效应”（MLPE）回归来解释种群效应（n=42）。在这里，使用AIC评估所有景观变量组合中最受支持的模型，并与IBD零模型进行比较，并获得每个建模景观因素的系数。使用lm函数计算模型R²值。接下来，我们采用来自RADISH的最受支持的IBR电导模型并将其转换为阻力表面。然后阻力表面被OMNISCAPE v0.6.2使用以可视化景观知情的基因流动。OMNISCAPE输入文件是由RADISH生成的最佳拟合栅格阻力模型，但重新缩放后最小值为1以符合程序要求。

为了评估IBE，使用了具有随机化的多元矩阵回归（MMRR），在R包ALGATR 1.0.0中实现，同时考虑了整体景观地理（IBD或IBR）的影响。具体来说，在回归框架中，使用地理距离矩阵或从最佳IBR模型生成的阻力矩阵与IBE一起共同估计对遗传距离（响应变量）的影响。我们使用随机排列（n=999）方法来解释种群依赖性并计算IBE模型的统计显著性。所有景观变量都在模型中进行了测试，没有进一步的变量选择，因为选择是在上游数据准备阶段完成的。计算模型R²值以评估回归解释的总变异。

为了控制任何潜在的更深层次历史种群结构的偏差，这可能会使景观模型的结果产生偏差，我们另外使用了线性混合效应模型（LMM）与R包corMLPE 0.0.3和NMLE 3.1-162，添加了一个随机效应来控制K=3对最受支持的IBR模型和IBE模型。每个种群对被分配一个独特的祖先对ID。为了比较模型，我们计算了AIC_C，在R包MUMIN 1.48.4中实现，以选择用于解释的最佳模型。

3 Results

3.1 Genetic Distance, IBD, and Summary Statistics

生物信息学处理产生了约80K个未连锁的SNP，代表了所有地理区域和样本（平均缺失率17%，随后估算，平均深度≈29，最小深度≈6，最大深度≈129，标准差≈8）。我们保留了所有高质量的单例用于分析，因为使用MAC过滤器移除它们不会影响结果的解释。具体来说，对MAC敏感的种群遗传汇总统计量仅略微改变了数值，并且跨种群和区域的相对差异与保留所有等位基因的分析相同。种群之间的遗传距离度量通常较低，表明遗传分化较低，但幅度因统计量而异。为了选择最具代表性的遗传距离指标，我们检查了交叉相关性和与地理距离（IBD）的相关性，使用Mantel检验。Dps（mantel系数=0.21）和线性化Fst（0.28）都为美国中西部的地理隔离提供了中等强度且显著的信号。线性化Fst的IBD模型的R²值为8%。使用所选遗传距离的Mantel相关图显示，IBD模式较弱，但在约300公里尺度上导致遗传差异，表明扩散效应在更大地理尺度上减弱。

预期杂合度（H_e）、核苷酸多样性（π）和稀有的私有等位基因丰富度在各州之间显示出空间变异，威斯康星州、印第安纳州和俄亥俄州表现出比密歇根州更大的基因组多样性和稀有的私有等位基因丰富度。Tajima's D在各州之间的空间趋势不太明显，但所有种群都具有负值，表明种群正在扩张。请注意，当样本量较小时，Tajima's D估计不太可靠，因此，在少于八个样本的地点，数值应谨慎解释。

3.2 Population Structure

通过PCA估计了个体样本的基因组多样性，并使用前两个主成分进行可视化。PCs的总体趋势显示，威斯康星州几乎所有地点都聚集在图表的一个区域，一些个体样本从主集群分散地偏离。沿PC1具有 divergent 样本的威斯康星州地点在地理上并不聚集。密歇根州样本沿PC1分组，其中一些与威斯康星州的 divergent 个体重叠。所有来自威斯康星州-爱荷华州边境（WI-IA）密西西比河流域、印第安纳州（IN）和俄亥俄州（OH）的个体样本似乎聚集在一起。它们与大多数威斯康星州和密歇根州样本不同。俄亥俄州包含一些在PC空间中与威斯康星州和密歇根州重叠的样本。

通过使用SNMF聚类方法定量地进一步探索了遗传结构。最佳模型估计为K=2，次优模型为K=3。我们检查了K=2到K=7的条形图，以评估不同推定祖先的可解释性。来自相同地理区域的个体在祖先比例上显示出相似性，但严重混合。通过将个体汇总到样本点并将其祖先表示为饼图，进一步检查了祖先比例。虽然K=3在视觉上最可解释，但无论K值如何，威斯康星州都由最多样化的祖先定义。密歇根州的同质性在各种K值中都很突出。同时，分析一致地将印第安纳州、爱荷华州和威斯康星州西南部的种群与威斯康星州和密歇根州的所有其他地点区分开来。所有种群和区域都表现出不同程度的混合。

最后，为了通过考虑空间自相关来证实SNMF结果，我们进行了CONSTRUCT分析。结果显示，从K=3到7，空间和非空间模型在预测准确性方面表现相似。为了进一步选择K=3和7之间的最佳模型，我们检查了不同K值下祖先层的贡献。在K=3时，经过空间校正的模型所有祖先都有意义地贡献了比例（1%），同时保持K的简洁性。CONSTRUCT的空间模型K=3，校正了空间自相关，产生了与SNMF的K=3中发现的混合模式相似的视觉结果。因此，SNMF和CONSTRUCT之间结果的一致性证实了SNMF K=3图的稳健性。

3.3 Gene Flow and Landscape Isolation Modeling

为了检查我们样本点上时间平均和异质性的IBD基因流动模式，我们生成了一个有效迁移表面。威斯康星州在遗传障碍（红色调区域）方面存在空间异质性，而密歇根州在遗传连通性（蓝色调区域）方面相对同质。特别是对于威斯康星州，障碍包围了大多数北部威斯康星州种群。在威斯康星州南部和中部，有几个东西向的基因流动走廊，可能连接到伊利诺伊州。可以发现连接伊利诺伊州、印第安纳州和俄亥俄州的长走廊，但白色表示缺乏基因流动的确定性。我们发现了将密歇根州与中西部其他地区隔开的连续障碍的证据。具有其他lambda值的迁移表面都同意密歇根州与其他地区的分离。除了全范围的连通性和障碍之外，检查特定地点和区域是否充当与其他地点的遗传连通性枢纽也很有帮助。基于Dps的遗传图，揭示了更近期的基因流动，表明在威斯康星州西北部和北部有两个枢纽，另一个枢纽在密歇根州东部。印第安纳州和俄亥俄州与威斯康星州的连接比与密歇根州的连接更紧密。最后，威斯康星州-爱荷华州边境的密西西比河地区和印第安纳州种群构成了它们自己的扩散网络。

先前的基因流动模型对基因流动背后的机制是不可知的。为了评估多变量景观在采样地点种群之间调节基因流动的影响，我们使用RADISH跨遗传距离度量线性化Fst测量了遗传隔离景观阻力（IBR）。所有可能景观组合中最优的IBR模型表明，土壤有机质和粘土含量、最干季度降水和