全基因组比较揭示细粒棘球绦虫（Echinococcus granulosus sensu stricto）复杂种群动态与寄生适应性的研究意义

本文通过对 111 个细粒棘球绦虫（E. granulosus s.s.）样本全基因组测序分析，揭示其种群结构、基因流等特征，发现 G1 和 G3 基因型存在线粒体 - 核基因组不一致现象，还确定了平衡选择下的关键基因，为包虫病防控提供新靶点。

　　

引言

囊性包虫病（Cystic echinococcosis，CE）是一种由细粒棘球绦虫（Echinococcus granulosus sensu lato）复合体引起的全球性人畜共患寄生虫病，对人类和家畜健康危害极大，每年造成大量的伤残调整生命年损失和经济负担。细粒棘球绦虫狭义种（E. granulosus sensu stricto，E. granulosus s.s.；基因型 G1 和 G3）是导致全球约 88% 人类 CE 病例的主要病原体 。

以往对该寄生虫的研究多依赖线粒体 DNA 来分类基因型和了解遗传多样性，但由于其复杂的生殖生物学特性，线粒体序列可能无法准确推断种群间的进化关系。例如，G1 和 G3 基因型的线粒体基因组存在固定差异，但在三个核基因座上却缺乏遗传差异，这凸显了进行全基因组比较研究的必要性。随着基因组资源的发展，虽为细粒棘球绦虫的研究提供了便利，但全基因组变异数据的缺乏仍限制了对其生殖生物学、传播历史和近期适应性的深入理解。本研究对来自中国关键流行地区的 * E. granulosus s.s.* 进行了全基因组调查，旨在揭示其种群结构、线粒体 - 核不一致性、基因流和选择特征，为寄生虫的传播途径、生殖策略和适应性遗传特征提供新见解，并确定平衡选择下的关键基因，为新的防控干预措施提供潜在靶点。

结果

1. 全基因组多态性检测：研究选取了 111 个中国的 * E. granulosus s.s.* 样本，经定量 PCR（qPCR）鉴定，这些样本的宿主基因组污染率均低于 5%。基于地理分布，样本被分为新疆（XJ，n = 33）、青海和甘南（QG，n = 42）、西藏（XZ，n = 36）三组。全基因组测序（WGS）结果显示，所有样本的测序深度最低为 35.16 倍，平均为 298.55 倍。变异检测共识别出 715,521 个高质量单核苷酸多态性（SNPs），SNP 密度在基因间区域较高，外显子区域较低，且与染色体结构相关。例如，Chr3 基因丰富、重复序列少，多态性密度最低；Chr6 则相反，SNP 密度最高。此外，还鉴定出了群体特异性 SNPs，XJ 组样本虽数量较少，但独特 SNPs 数量最多。
2. G1 和 G3 基因型的线粒体 - 核不一致性：线粒体基因组系统发育树显示，G1 和 G3 基因型存在明显的遗传分化，但 G1 基因型内部无地理结构。G1 样本的线粒体单倍型网络分析也表明，其缺乏明显的地理分化，且新疆样本似乎代表了更古老的单倍型，可能是 G1 基因型后续传播的起源地。相比之下，核基因组的系统发育树和主成分分析（PCA）显示，尽管 G1 和 G3 基因型的线粒体基因组存在明显分化，但 G3 基因型的核基因组与地理相关的 G1 基因型样本聚类紧密，证实了两者之间的线粒体 - 核不一致性，这是典型的深度线粒体分化（DMD）现象，可归因于 “分化逆转” 现象，即它们经历了古代的隔离分化，随后又发生了频繁的核基因流，导致基因组同质化，而线粒体基因组保留了原始遗传变异。
3. 种群结构：基于核基因组的系统发育分析显示出明显的地理结构，XJ 组形成了一个独特的单系分支，QG 和 XZ 组虽未形成明显的单系分支，但存在地理聚类。PCA 分析进一步证实，XJ 组与 XZ、QG 组在 PC1 上明显分离，而 XZ 和 QG 组的 95% 置信区间有部分重叠。ADMIXTURE 分析表明，当K = 2 时，XJ 组与其他两组的遗传共享极少，应被视为一个遗传上不同的种群；当K = 3 时，QG 和 XZ 组的主要遗传成分开始分化，但尚未形成完全遗传独立的种群。遗传重组导致连锁不平衡（LD）衰减，不同种群的 LD 衰减趋势相似，但速率与重组率和遗传混合程度相关。例如，QG 组可能由于经历了更多的遗传混合，重组率较高，LD 衰减更快；XJ 组则因遗传混合有限，重组率较低，LD 衰减较慢。此外，成对 Wright’s F统计量（F_ST）显示，XJ 和 XZ 种群之间的遗传分化最大，XJ 种群的核苷酸多样性（π）和有效种群大小（N_e）最高，XZ 种群最低。
4. 雌雄同体细粒棘球绦虫的异交繁殖：研究发现，细粒棘球绦虫的基因组 LD 衰减速度比雌雄异体的日本血吸虫快，其平均重组率与其他异交物种相当，且染色体低基因密度区域更易发生重组，引入更多杂合位点。部分样本还存在意外的杂合性，多数样本的交配偏好类似随机交配生物，但也有部分样本近交系数（F值）较高，表明近亲交配或自体受精也会发生，这可能与宿主体内虫体密度低或局部寄生扩散有关。这些结果表明，雌雄同体的细粒棘球绦虫具有异交繁殖能力。
5. 影响样本中低次要等位基因频率（Low-MAF）位点的因素：不同种群样本中 Low-MAF 位点数量存在明显差异，XJ 种群样本的 Low-MAF 位点数量较多，XZ 种群最少。研究发现，部分样本为多倍体或具有混合遗传背景，但它们并非导致种群间 Low-MAF 位点数量差异的主要因素。进一步分析表明，有效种群大小和遗传多样性是影响 Low-MAF 位点数量的主要因素，XJ 种群较大的有效种群大小和较高的核苷酸多样性使其杂合性和 Low-MAF 位点增加，而 XZ 种群则相反。
6. 细粒棘球绦虫种群的种群历史：通过多种检测方法确定，种群间的基因流是导致遗传混合的主要因素。D统计量（ABBA - BABA 测试）等分析表明，从 XJ 到 QG 的基因流最为频繁且强度最强，其次是从 XJ 到 XZ 以及 QG 和 XZ 之间的基因流，而从 QG 和 XZ 到 XJ 的基因流相对较弱。SMC++ 计算的有效种群大小历史显示，三个种群的N_e在约 10 万年前逐渐下降，随后在最近几千年经历瓶颈期后扩张。XJ 种群最早经历瓶颈事件（约 5 - 8 千年前），其次是 QG 种群（约 4 - 6 千年前），最后是 XZ 种群（约 1 - 2 千年前），这反映了细粒棘球绦虫种群扩散过程中种群建立的顺序。
7. 种群间缺乏环境适应性进化：利用跨种群方法（F_ST和 XP - nS_L）比较不同气候环境下的种群，发现虽有基因在相关分析中被检测到，但未发现与这些基因相关的连锁选择位点。绝对分歧指数（D_xy）在基因富集区域和重复序列区域的差异不明显，表明种群间的分歧可能主要与近期基因流减少或祖先分歧的单倍型分布有关，而非不同环境压力下的定向选择。
8. 基因组平衡选择表明对宿主内部环境的适应：将三个种群视为一个整体，采用种群内方法（Tajima’s D和nS_L）进行分析。结果显示，在 Tajima’s D和nS_L分析中，重叠于前 1% 窗口的 674 个基因在京都基因与基因组百科全书（KEGG）通路中未显著富集，但在基因本体（GO）富集分析中，许多基因在与膜相关的 GO 术语（GO:0016020）中显著富集，表明它们对寄生生活方式的适应。其中，ECG_09652 基因（编码钠 / 胆汁酸协同转运蛋白，SBAT）参与宿主肠道胆汁酸的摄取，对细粒棘球绦虫的发育至关重要。该基因的启动子区域、外显子和内含子存在强连锁选择位点，这些位点将样本分为三种单倍型。该基因所在的 Chr9 区域处于强平衡选择之下，且平衡选择并非由地理环境或隔离引起，而是可能源于对宿主内部环境的适应。此外，还鉴定出 11 个表面蛋白编码基因也受到平衡选择，这些基因可能是免疫保护的潜在靶点。

讨论

1. 细粒棘球绦虫的异交繁殖：本研究从多个基因组层面为细粒棘球绦虫的异交繁殖提供了证据，包括 G1 和 G3 基因型的线粒体 - 核不一致性、基因组的快速 LD 衰减、高重组率和重组结构分布、样本的意外杂合性以及种群间的基因流等。同时，部分样本的高近交系数也表明近亲交配或自体受精的存在，这与细粒棘球绦虫的生活史特征相关。
2. 线粒体 - 核不一致性：细粒棘球绦虫 G1 和 G3 基因型的线粒体 - 核不一致性可能由多种原因导致。“分化逆转” 现象，即地理隔离后的基因流导致核基因组同质化，而线粒体基因组因单亲遗传保留原始变异，是一种可能的解释；“幽灵渗入” 假说认为 G3 的线粒体基因组可能是灭绝基因型的残余并通过基因渗入获得，但这需要进一步验证。此外，核基因组比线粒体基因组具有更多的系统发育信息位点，能检测到更细微的分化，这也解释了为何线粒体基因组地理结构较弱，而核基因组能更好地反映种群间的关系。
3. 细粒棘球绦虫在中国的传播历史：细粒棘球绦虫的传播与畜牧业密切相关，早期可能随中东地区牲畜迁移传播。本研究的遗传数据支持中国细粒棘球绦虫种群最初由中东地区随绵羊迁移引入，新疆是首个引入地。随后，通过 QG 地区作为中间步骤，传播至 XZ 地区，这与藏羊的种群遗传分析和考古记录相吻合。历史上的贸易路线、民族交流和游牧活动等进一步促进了基因流和寄生虫的传播。
4. 种群间的遗传分化：对青藏高原种群的研究未发现明显的与高原环境适应相关的选择信号，基因富集区域的D_xy值与重复区域无显著差异。这可能是由于青藏高原种群的殖民历史相对较短，且寄生虫在寄生生活阶段生活条件稳定，因此目前种群间的遗传分化主要源于地理隔离和近期基因流的减少，而非环境压力下的适应性进化。
5. 细粒棘球绦虫对寄生生活方式的适应性进化：寄生生物依赖膜相关蛋白进行多种重要功能，本研究中受选择的基因在膜相关 GO 术语中显著富集，表明细粒棘球绦虫对寄生生活的适应性。ECG_09652 基因编码的 SBAT 在平衡选择下保持高遗传稳定性，对调节胆汁酸摄取至关重要，可能是治疗 CE 的潜在药物靶点。其他受平衡选择的表面蛋白编码基因则可能是免疫保护的潜在靶点。

结论

本研究通过对细粒棘球绦虫种群的全基因组测序分析，揭示了 G1 和 G3 基因型以及中国三个不同种群间复杂的遗传动态。历史迁移和地理隔离在塑造其遗传景观中起重要作用。不同地理种群的基因组缺乏环境适应性进化的迹象，但存在明显的平衡选择特征，这可能与对寄生生活方式的适应有关。平衡选择下的基因可为药物研发和疫苗设计提供潜在靶点，为包虫病的防控提供了新的思路和方向。

材料与方法

1. 样本采集和 DNA 提取：从新疆、西藏、青海和甘肃的屠宰场采集感染 CE 的牦牛和绵羊的肝脏和肺脏，提取原头蚴并清洗，保存于 - 80°C。使用 Qiagen Blood and Tissue DNA Extraction Kit 提取 DNA，通过 PCR 扩增和测序进行基因分型，再用 qPCR 检测宿主基因组污染比例，选取寄生虫基因组占比超过 95% 的样本进行测序。
2. 序列比对、变异检测和注释：利用 fastp 去除测序读段的接头和低质量碱基，使用 BWA 将合格的双端读段比对到参考基因组，用 samtools 对 BAM 文件进行下采样以标准化测序深度。使用 FreeBayes 检测变异，通过 VCFtools 筛选变异位点，并用 SnpEff 进行功能注释。
3. 遗传结构分析：通过 VCFtools 生成稀疏 SNP 子集，将其转换为 fasta 序列后用 phylip 构建邻接（NJ）树，同时构建线粒体 NJ 树和 G1 基因型线粒体单倍型网络。使用 PLINK 进行 PCA 分析，ADMIXTURE 进行混合分析。此外，还计算了核苷酸多样性（π）、相对分歧指数（F_ST）、绝对差异指数（D_xy）、连锁不平衡衰减、重组率、近交系数（F值）、Low-MAF 位点等，并分析了样本的倍性和混合情况。
4. 种群历史推断：使用 Dsuite 计算D统计量（ABBA - BABA 测试）和 f4-ratio 评估种群间基因流，用 Fbranch 计算i_b指标。利用 SMC++ 基于未分型高质量 SNP 数据集估计有效种群大小历史，通过 vcf2smc 转换数据格式，用 estimate 和 CV 子命令进行计算和交叉验证，根据突变率调整时间并计算当前N_e。
5. 选择信号检测：运用 Selscan、VCFtools 等工具，基于全基因组 SNP 数据集，采用种群内（nS_L和 Tajima’s D）和跨种群（F_ST和 XP - nS_L）统计方法检测选择信号。确定候选选择窗口和基因后，进行 GO 和 KEGG 富集分析，利用相关工具生成 LD 热图、绘制单倍型衰减图，并计算B₂统计量以确定平衡选择区域。