Hunga火山喷发后热液喷口动物与共生菌群对种群衰退的基因组响应对比及其进化动态研究

本综述深入探讨了2022年Hunga海底火山喷发对西南太平洋Lau海盆热液喷口动物-微生物共生体系种群瓶颈效应的基因组响应。研究发现宿主种群全基因组多样性变化微弱但有效种群大小（Ne）长期下降，而共生菌群则表现出基因组变异显著减少（包括特定栖息地菌株的丢失）和塔吉马D值（Tajima's D）负向偏移。环境样本中检测到的游离态共生菌提示其可能成为宿主重新招募共生体的储备库。该研究为自然扰动后动物-微生物种群的进化动力学提供了新见解，对评估深海生态系统恢复力具有重要意义。

　　

引言

遗传多样性对物种长期存续至关重要，它使种群能够适应环境变化（Kardos等2021；Thompson等2023）。种群遗传变异受自然选择类型和强度、遗传漂变、基因流水平及有效种群大小（N_e，即 actively interbreeding 的个体数量）等多种相互作用因素影响（Charlesworth 2006；Reed 2008；Walsh和Blows 2009；Assis等2016；Ellegren和Galtier 2016）。由自然灾害或人类活动引起的瓶颈效应是常见的进化情景，会强烈降低种群的有效大小和相关遗传变异。根据其严重程度，遗传瓶颈可通过限制自然选择效能和增加遗传漂变影响，显著降低受影响种群的适应性和适应潜力，最终增加灭绝风险（Whitlock 2000；Markert等2010；Gamblin等2024）。或者，瓶颈事件后的选择也可能有利于具有恢复力特性的基因型或谱系，最终促进对扰动条件的快速适应（Coleman和Wernberg 2020综述）。由于瓶颈事件后进化轨迹的不确定性，关于有效种群大小和遗传多样性的信息日益被纳入受威胁物种的管理计划（Willi等2022；Ottewell和Byrne 2022）。

深海热液喷口是神秘的海洋栖息地，拥有独特的动物和微生物群落（van Dover 2000）。这些生态系统中的许多无脊椎动物物种与化学合成细菌存在营养共生关系，这些细菌利用化学还原剂产生的能量将无机碳转化为有机物（Sogin等2020，2021）。值得注意的是，许多喷口相关类群在国际自然保护联盟红色名录中被列为易危或濒危物种（https://www.iucnredlist.org）。这部分是由于热液喷口的岛屿状性质以及喷口特有物种的分布受限。因此，从深海热液喷口商业开采矿产资源的计划引发了对喷口相关物种和生态系统恢复力的担忧（Levin, Baco,等2016；van Dover等2018；Orcutt等2020；Brunner等2022），凸显了需要更好地理解这些类群对环境扰动的进化响应。

2022年，海底Hunga火山（原Hunga Tonga–Hunga Ha?apai）的喷发在西南太平洋Lau弧后盆地的深海热液喷口种群中导致了大规模死亡事件（Beinart等2024），为研究喷口相关物种遗传瓶颈的影响及其在自然灾害后的恢复过程提供了机会。Lau海盆喷口群落主要由不同共生的腹足类（Alviniconcha boucheti, A. kojimai, A. strummeri, Ifremeria nautilei）和一种双壳类（Bathymodiolus septemdierum）组成。它们共同代表了该区域热液生态系统的基础动物区系（Podowski等2009，2010）。所有这些类群都携带硫氧化弯曲菌或γ变形菌共生体，这些共生体在动物幼虫阶段从环境中获得，随后寄居在动物鳃组织的细胞内或细胞外（Suzuki, Kojima, Sasaki,等2006；Suzuki, Kojima, Watanabe,等2006；Duperron 2010；Ikuta等2021）。Lau海盆内的宿主-共生体关联是物种特异性的，每个宿主类群只携带一种或两种共生体物种，尽管共生体菌株组成可能随当地环境条件而变化（Breusing等2022；Breusing, Xiao,等2023）。

Alviniconcha、Ifremeria和Bathymodiolus共生体在喷发前在东Lau扩张中心（ELSC）的热液喷口丰富，尽管A. boucheti、I. nautilei和B. septemdierum在最北端的Kilo Moana喷口场在2009年至2015年间变得不活跃时已经经历了一定程度的种群损失（Reysenbach和Seewald 2015）。火山影响进一步减少了最靠近Hunga喷发的喷口场所有三种软体动物属的种群，其中Alviniconcha全生体显示出比Ifremeria或Bathymodiolus全生体更强的种群下降（Beinart等2024）。所有三个属在Tow Cam喷口场完全或几乎完全灭绝，该地经历了高达150厘米的显著火山灰沉积（Beinart等2024）。Alviniconcha在Tahi Moana喷口场也完全灭绝，该地同样有大量火山灰沉积（12-23厘米）。在ABE喷口场，火山灰沉积厚度为8-15厘米，只有小片Alviniconcha存活。相比之下，稀疏的Ifremeria和Bathymodiolus种群在Tahi Moana和ABE仍然存在。基于通过视频调查对喷发前后种群的半定量评估，所有三个属在Tu'i Malila喷口场维持了看似正常的种群，该地没有火山灰沉积（Beinart等2024）。

在本研究中，我们比较了喷发前后的Alviniconcha、Ifremeria和Bathymodiolus共生体种群，以评估宿主和共生体种群遗传多样性和有效种群大小的初始变化。作为重新定殖潜力的一个决定因素，我们进一步评估了细菌共生体的游离形式是否仍能在环境中检测到。

材料与方法

2.1 样品收集、DNA提取和宏基因组测序

通过远程操作车辆从东Lau扩张中心（ELSC）的五个喷口场（Kilo Moana, Tow Cam, Tahi Moana, ABE, Tu'i Malila）在2009年、2016年和2022年R/V Thomas G. Thompson或R/V Falkor考察期间收集蜗牛和贻贝样品。2022年考察在Hunga喷发后约3个月进行。从汤加火山弧（距ELSC站点约700公里）未受影响的热液喷口群落在2016年和2022年考察期间获得了额外样品。本研究这些样品仅用于宿主的转录组重建以提高组装质量和完整性（见下文）。

所有蜗牛通过远程操作车辆Jason II或Ropos收集。在船上，每个动物的含共生体鳃组织被切除并快速冷冻或在-80°C下保存在RNALater（Thermo Fisher Scientific Inc., Waltham, MA, USA）中。DNA在罗德岛大学使用E.Z.N.A Mollusc DNA Kit（Omega Bio-tek, Norcross, GA, USA）提取，并送至Psomagen Inc.（Rockville, MD, USA）进行宏基因组文库制备和测序。文库使用plexWell384（seqWell, Beverly, MA, USA）和Twist 96-Plex（Twist Bioscience, San Francisco, CA, USA）试剂盒制备，并在NovaSeq 6000和NovaSeq X Plus仪器上进行150 bp双端测序，平均深度为每个样品34,112,119总读数。从国家生物技术信息中心（NCBI）下载了额外的序列数据（BioProjects: PRJNA523619, PRJNA855930）。我们之前的工作表明，通过我们开发的变异调用流程可以减少不同文库制备和测序运行的批次效应（Breusing, Xiao,等2023）。因此，我们有信心将这些不同数据集合并进行分析。在进一步分析之前，使用Trimmomatic（Bolger等2014）从原始读数中去除适配器和低质量碱基。每个宿主物种的身份通过形态学确定，并在必要时通过COX1条形码基因测序进一步验证（Folmer等1994）。

2.2 宿主和共生体遗传参考的生成

为了评估宿主和共生体种群遗传多样性的变化，我们首先组装了宿主转录组和共生体泛基因组作为读数映射的参考。

宿主：对于每个Alviniconcha和Ifremeria物种，从ELSC北部、南部和中部以及汤加火山弧的一个代表性喷口场的一个样品中提取总RNA，使用Zymo Direct-zol RNA Miniprep Kit（Zymo Research Inc., Irvine, CA, USA）。RNA提取物送至Psomagen Inc.（Rockville, MD, USA）使用TruSeq Stranded mRNA Kit（Illumina Inc., San Diego, CA, USA）制备文库，并在NovaSeq 6000仪器上进行2×150 bp双端测序，每个文库深度约为100 M读数。RNA-seq读数进行质量过滤、错误校正和组装，如Breusing等（2022）；Breusing, Xiao,等（2023）所述。对于Bathymodiolus septemdierum，使用Breusing, Xiao,等（2023）现有的高质量转录组进行分析（NCBI accession: GJZH00000000）。

共生体：对于每个样品，过滤后的宏基因组读数使用metaSPAdes（Nurk等2017）组装，使用21到121的kmer，然后使用metaBAT2（Kang等2019）、MaxBin2（Wu等2016）和metaWRAP（Uritskiy等2018）将组装结果分箱成宏基因组组装基因组（MAGs）。使用DAS Tool（Sieber等2018）识别高质量、非冗余MAGs，并使用GTDB-Tk（Chaumeil等2019）进行 taxonomic 分配。用于分析的额外MAGs从NCBI BioProjects PRJNA523619和PRJNA855930下载。使用fastANI（Jain等2018）计算每个细菌类群MAGs之间的平均核苷酸同一性（ANIs），以确定它们是否代表同一物种（即ANI > 95%）。每个细菌物种的MAGs，其CheckM（Parks等2015）完整性得分 > 90%且污染得分 < 10%，进一步用于使用Panaroo（Tonkin-Hill等2020）进行泛基因组重建，如Breusing等（2022）；Breusing, Xiao,等（2023）所述。由于来自Kilo Moana和Tow Cam的A. boucheti共生体高质量MAGs稀疏，这些喷口场MAGs的完整性阈值降低至80%。

2.3 种群结构、核苷酸多样性、杂合性和Tajima's D的评估

对于每个宿主-共生体对，我们进行了场点特异性和全海盆分析，以调查与Hunga喷发相关的遗传多样性变化。在场点特异性分析中，我们只关注每个宿主-共生体对在喷发前后采样且样本量足够（即每个场点 > 10个体）的喷口场，并比较不同年份的模式。由于Alviniconcha在最北端站点完全灭绝，因此对于该属物种，这些比较只关注受影响较小的喷口场Tu'i Malila或ABE，这些地方Alviniconcha在喷发前后都存在。相反，受影响较重的喷口场被纳入I. nautilei和B. septemdierum的比较中，它们在喷发后在整个ELSC仍然存在。在全海盆分析中，我们将一个宿主-共生体对的所有采样喷口场合并为“喷发前”和“喷发后”组。对于两组，我们使用相同数量的总样本，并旨在包括每个可用喷口场至少10个样本。

宿主：每个选定样品的宏基因组读数使用Bowtie2（Langmead和Salzberg 2012）以非常敏感的局部模式映射到相应的宿主转录组。生成的BAM文件使用Picard（https://github.com/broadinstitute/picard）和LoFreq（Wilm等2012）进一步处理，以标记光学重复、重排indels和重新校准碱基质量。随后使用Angsd（Korneliussen等2014）估计每个分析组的基因型可能性，考虑Hardy–Weinberg平衡的偏差（Vieira等2013），并使用折叠位点频率谱推断汇总统计。每个位点的最大和最小深度基于全局深度分布设置，而保留位点所需的有数据个体最小数量在场点特异性分析中设置为50%，在全海盆分析中设置为75%。对于种群结构的估计，变异位点经过严格过滤，只包括高质量单核苷酸多态性（SNPs）（-baq 1 -C 50 -minInd (nInd*0.75) -minMapQ 30 -minQ 20 -uniqueOnly 1 -remove_bads 1 -only_proper_pairs 1 -SNP_pval 1e-6 -sb_pval 0.05 -hetbias_pval 0.05 -minMaf 0.01 -skipTriallelic 1）。由于Alviniconcha全海盆分析的一些样品稍后制备，且初步分析表明宿主数据中存在由Twist 96-Plex文库制备和测序引起的批次效应，我们在Alviniconcha宿主物种中应用了两个额外的SNP过滤器以去除此效应（-qscore_pval 0.05 -edge_pval 0.05）。使用Angsd的thetaStat和realSFS程序分别计算这些SNPs的每位点核苷酸多样性和种群间配对F_STs，而观察和预期杂合性（即给定基因座具有两个不同等位基因的个体比例）从Angsd的HWE汇总文件计算。我们进一步推断了Tajima's D，这是一个常见的群体遗传统计量，指示特定基因组区域的突变模式是否遵循中性假设（Tajima 1989）。负值暗示低频多态性过量，表明正选择或近期瓶颈事件后的种群扩张。正值表示稀有等位基因缺失，可能由于平衡选择或种群收缩。Tajima's D值在500 bp滑动窗口中计算，使用宽松的过滤标准以防止统计量偏向正值（-baq 1 -C 50 -minMapQ 30 -minQ 20 -uniqueOnly 1 -remove_bads 1 -only_proper_pairs 1）。

共生体：使用Bowtie2以非常敏感的局部模式将每个选定样品的宏基因组读数映射到相应的共生体泛基因组。与宿主一样，BAM文件进一步处理以去除PCR重复、indel重排和碱基质量重新校准。为了最小化不均匀读数深度对变异检测的影响，读数比对下采样至样品中最低数量的比对读数，使得每个样品泛基因组覆盖率至少为30倍。使用FreeBayes（Garrison和Marth 2012）调用种群结构分析的变异，遵循Breusing等（2022）；Breusing, Xiao,等（2023），排除缺失数据超过25%的样品和变异位点。然而，由于大多数来自Tu'i Malila的A. boucheti共生体样品和大多数来自Tow Cam的A. kojimai共生体样品达到此排除阈值，这些全生体全海盆分析中允许的缺失数据量增加至35%。我们之前的工作（Breusing, Xiao,等2023）表明，应用的SNP过滤器在减少不同文库制备和测序运行对热液喷口共生体群体基因组分析的影响方面是有效的，并且相对于已知生物信号，在最终数据中未观察到强批次效应。基于最终变异数据集的排序结果中可见的喷口特异性聚类用于界定位置特异性、亚种水平的共生体种群或菌株。这些菌株定义通过使用DESMAN的STRONG实现进行菌株反卷积分析进一步补充（Quince等2017，2021）。使用python scikit-allel模块（https://github.com/cggh/scikit-allel）计算喷口特异性共生体种群之间的配对F_STs，而基于双等位基因SNPs，使用Romero Picazo等（2019）的方法为每个喷口种群推断每位点和宿主间核苷酸多样性。每位点核苷酸多样性衡量共生体种群在每个基因组位点的遗传变异，而宿主间核苷酸多样性是指示两个宿主个体的共生体种群遗传差异程度的指标。对于Tajima's D的计算，FreeBayes运行时不设最小等位基因频率阈值，生成的VCF文件经过最小过滤以去除indels周围5 bp范围内的SNP位点。Tajima's D在500 bp窗口中估计，使用允许为单倍体数据计算汇总统计的VCFTools重新实现（https://github.com/jydu/vcftools）。

变异数据的主成分和非度量多维标度图以及核苷酸多样性、Tajima's D和观察杂合性的箱线图和小提琴图在R（R Core Team 2024）中生成。使用vegan包（Oksanen等2024）进行PERMANOVAs以确定喷发对宿主和共生体群体遗传变异的影响。使用配对t检验或Wilcoxon秩检验评估组间核苷酸多样性、Tajima's D和观察杂合性的统计差异以及组内预期和观察杂合性。基于错误发现率程序进行多重检验校正。在宿主中，参数估计统一子采样至每组100个数据点以减少与过采样相关的偏差。

2.4 宿主有效种群大小变化的评估

我们使用基于2022年样品的折叠位点频率谱信息来确定宿主种群在Hunga喷发后是否经历了瓶颈效应。使用Angsd的realSFS程序计算单个喷口场的折叠位点频率谱，使用与Tajima's D相同的过滤标准。然后使用Stairway Plot2（Liu和Fu 2020）推断宿主有效种群大小（N_e）的变化，假设软体动物的突变率为每个核苷酸位点和世代1.645e-9，世代时间为1年（Allio等2017）。此外，我们评估了2022年样本的长期和当代N_e估计。长期有效种群大小根据公式计算：N_e = θ/(4 μ)，使用Watterson's θ估计在基因组窗口中的中位数和上述突变率（μ）为每个位点和世代1.645e-9。当代有效种群大小根据时间等位基因频率变化推断，遵循Jorde和Ryman（2007）的方法，如Pinsky等（2021）所实现，使用每个基因座一个SNP的稀疏SNP数据集和最小次要等位基因频率0.01。当可用时，使用2016年样本作为参考种群。否则，使用2009年样本。

2.5 热液环境中游离共生体的检测

为了确定Alviniconcha、Ifremeria和Bathymodiolus共生体在环境中的潜在存在，我们对使用悬浮颗粒物玫瑰花采样器（SuPR）（Mclane Research Laboratories Inc., Falmouth, MA USA；Breier等2009）或通用流体获取器（UFO）（National Deep Submergence Facility, Falmouth, MA, USA）和微生物定殖装置收集的深水和水生物膜样品进行了16S V4高变区测序。在每个喷口场，15.71–47.51升水通过2–4个0.22 μM Express Plus膜过滤器（MilliporeSigma, Burlington, MA, USA）泵送，位置靠近动物群落（动物组合上方约10–20厘米；以下称为热液流体样品）和远离热液活动（距离约10秒至100秒米；以下称为海水样品）。每个样品过滤器前组装一个10 μM定制尼龙网预过滤器（Sefar Inc., Depew, NY, USA），以防止收集可能含有共生体的脱落动物细胞。由于SuPR设备在Tu'i Malila部署期间以及ABE和Kilo Moana每次潜水中失败一次，因此在这些地点额外使用UFO获取了10.5升深水样品。为了考虑使用此设备的交叉污染，部署了空过滤器但不泵送，用作阴性对照。对于微生物生物膜的收集，我们使用了加重的PVC饵笼，两端由窗纱网封闭，内含垂直放置、开口的2 mL聚丙烯管，填充破碎矿物基质（安山岩、玄武岩、玻璃珠或Alviniconcha壳）用于微生物定殖。三个这样的笼子在每个喷口场的第一次潜水中部署，并在大约2周后回收。在船上，流体、海水和生物膜样品保存在RNALater（Thermo Fisher Scientific Inc., Waltham, MA, USA）中，并在-80°C冷冻。使用酚:氯仿协议（Data S1）从过滤器样品中提取DNA，SuPR样品使用一半过滤器，UFO样品使用整个过滤器，而生物膜DNA使用E.Z.N.A. Soil DNA Kit（Omega Bio-Tek, Norcross, GA, USA）提取，包括在蛋白酶K消化前额外的45秒珠击步骤。对于SuPR和UFO样品，添加3 Strain Tagged Genomic DNA Even Mix（ATCC, Manassas, VA, USA）作为 spike-in 对照，并制备无DNA提取物的 spike-in 稀释液以评估实验室样品处理过程中的交叉污染。利用这些 spike-ins 的绝对微生物丰度数据将在另一项研究中报告（Hauer等准备中）。最终提取物送至阿贡国家实验室（Lemont, IL, USA）使用515F/806R引物对（Apprill等2015；Parada等2016）进行扩增子文库制备。所有文库在NextSeq2000仪器上进行2×150 bp双端测序，平均深度为每个样品1,455,454总读数。使用FastQC（https://github.com/s-andrews/FastQC）对原始读数进行质量检查，然后使用FastP（Chen等2018）和Trimmomatic（Bolger等2014）分别去除polyG尾和适配器。然后使用Usearch（Edgar 2010）和Vsearch（Rognes等2016）工具合并（-fastq_minmergelen 230 -fastq_maxmergelen 270 -fastq_maxdiffs 10 -fastq_pctid 80）、定向和过滤（–fastq_qmax 42 –fastq_minlen 230 –fastq_truncqual 20 –fastq_maxee 1）读数。过滤后的序列被去噪并分解为扩增子序列变异（ASVs），并通过Vsearch –search_exact命令评估其丰度。使用R包phyloseq（McMurdie和Holmes 2013）和metagMisc（https://github.com/vmikk/metagMisc）去除单例以及具有少于10个平均总读数计数和低于5%流行率的虚假ASVs。ASV序列导入QIIME2（https://qiime2.org）进行分类，使用在515F/806R引物区域训练并考虑海水和海洋生物膜特定分类权重的Naive Bayes分类器（Kaehler等2019）。来自γ变形菌和弯曲菌门的潜在共生体ASVs与已验证的Alviniconcha、Ifremeria和Bathymodiolus 16S rRNA共生体序列参考数据库进行blast。与参考序列有超过99%同一性的ASVs被分类为游离共生体或其近亲；然而，我们承认这里扩增的短16S rRNA基因片段不足以明确确定物种或菌株水平身份，从而确认共生体存在。所有样品的读数计数通过相应ASV在 spike-in 对照样品中的最大读数量进行校正。由于共生体序列在ABE的UFO对照过滤器中丰富，且ABE和Kilo Moana的UFO样品质量无法验证，因此这些喷口场的UFO样品不再用于计算总共生体读数计数。ABE和Kilo Moana的所有其他样品保留用于分析。对于Tu'i Malila，UFO样品中的读数计数通过减去UFO阴性对照中的读数计数进行校正。校正后读数计数少于20,000的样品被排除在进一步分析之外。

结果

3.1 宿主-共生体种群遗传结构

每个Alviniconcha和Ifremeria物种的三到五个鳃RNA样品的共同组装产生了36,164–57,396个转录本（63.91–82.73 Mb），代表了宿主转录组的75.80%–81.40%。针对这些参考转录本和先前重建的B. septemdierum转录组（Breusing, Xiao,等2023）的读数映射恢复了5,705–28,463个高质量SNPs用于群体遗传分析。与之前的研究一致（Breusing等2022；Breusing, Xiao,等2023），基于这些SNPs的主成分图显示Alviniconcha、Ifremeria和Bathymodiolus宿主物种中没有地理模式，并且同一喷口场喷发前后的种群之间没有观察到遗传差异。NMDS图通常与PCA图一致，但由于二维高应力值而提供较弱的分辨率。同样，配对F_ST值显示所有宿主物种内的种群分化水平总体较低，Bathymodiolus表现出最少的分化量（Alviniconcha: 0.0210–0.0656；Ifremeria: 0.0162–0.0505；Bathymodiolus: 0.0092–0.0141）。与PCA、NMDS和F_ST结果一致，PERMANOVAs表明喷发对任何宿主物种的种群结构没有显著影响。

在共生体中，79–172个高质量MAGs用于为每个细菌物种创建约2.87–4.61 Mb泛基因组，作为群体遗传分析的参考。针对这些泛基因组进行读数映射后，获得了2843–4417个高质量变异。与之前的报告（Breusing等2022）一致，并且与宿主模式相反，基于这些变异的主坐标和NMDS图以及配对F_ST计算表明Alviniconcha共生体种群在喷口地点之间存在强烈细分（F_ST: 0.1487–0.8617）。Ifremeria（F_ST: 0.0144–0.2504）和Bathymodiolus（F_ST: 0.0507–0.2122）共生体种群的细分存在，但比Alviniconcha共生体弱，其中Bathymodiolus Ca. Thiodubiliella共生体显示出最少的结构。使用DESMAN的菌株反卷积分析结果通常支持排序和F_ST为基础的分析结果，表明来自北部（Kilo Moana, Tow Cam）、中部（Tahi Moana, ABE）和南部（Tu'i Malila）喷口场的共生体种群代表不同的菌株，尽管使用此方法对Ifremeria和Bathymodiolus共生体的菌株水平分辨率明显较弱。来自最北端喷口场的Alviniconcha和Ifremeria共生体菌株在喷发后消失，而Bathymodiolus共生体菌株没有变化。评估喷发对群体基因组变异影响的PERMANOVAs对所有共生体物种均显著，但生物学相关性在Alviniconcha共生体中明显更高（R2 = 6.3%–22.5%） than those of Ifremeria（R2 = 2.0%） or Bathymodiolus（R2 = 2.2%）。然而，由于 dispersion 检验对Gamma1和Ifr-SOX共生体显著，我们不能完全排除这些物种的结果受到喷发前后分析组间 spread 方差的影响。

3.2 遗传多样性变化

除少数例外，在任何宿主物种的全海盆或喷口