靶向组装技术高效恢复潮间带泥滩氨单加氧酶基因多样性及其生态意义

本综述系统评价了靶向组装（targeted assembly）技术在恢复微生物功能基因（如amo）方面的显著优势。相比传统全组装方法（single/multi-sample assembly），该策略以更低的计算资源（computational resource）和时间成本（running time），获得了更高的遗传多样性（genetic diversity）、更完整的操纵子结构（operon）及更少的嵌合序列（chimeras），为研究氮循环（nitrogen cycling）和生态系统功能（ecosystem functioning）提供了高效技术路径。

　　

微生物群落在地球生物圈中介导多种生物地球化学循环过程，这些过程通过它们携带的功能基因维持生态系统的多功能性和稳定性。尽管研究不断进展，在复杂环境中获取高质量的微生物功能基因序列数据集仍然具有技术和计算资源上的挑战。本研究以编码氨单加氧酶的amo基因家族为例，旨在通过靶向组装从鸟枪法宏基因组中恢复重要的微生物功能基因。

与传统组装方法（如单样本组装和多样本组装）相比，靶向组装恢复了更高的amo基因多样性，同时所需计算资源和运行时间显著减少。此外，通过靶向组装恢复的amo基因嵌合体更少，同时恢复了更多的amo操纵子。不仅观察到了常见的amoABC亚基，还发现了较少见的亚基，如amoX和amoE。值得注意的是，通过靶向组装恢复的古菌amoA亚基代表了氨单加氧酶的最多“超级进化枝”，包括NT-α、NT-γ、NP-γ、NP-ζ和NP-η，证明了靶向组装相对于传统方法的优势。恢复的amo基因集合也观察到了可比的空间模式，如分类群-面积关系和距离-衰减关系。本研究展示了一条从鸟枪法宏基因组中以最小计算资源和运行时间恢复微生物功能基因的有效途径。

微生物群落通过介导必需元素的多种生物地球化学循环并凭借其携带的功能基因维持生态系统稳定性和多功能性，在地球生物圈中发挥关键作用。然而，从这些复杂群落中恢复关键功能基因仍然具有挑战性。不同技术存在各自的优势和局限性。在本研究中，以amo基因家族为例，我们表明靶向组装能够以最小的计算资源和运行时间，从鸟枪法宏基因组中准确快速地恢复高质量amo序列。与传统的全组装方法相比，通过靶向组装恢复的amo序列具有更多的操纵子、更高的（系统发育）遗传多样性和更少的嵌合体。本研究为恢复微生物功能基因提供了一条有效的替代途径，特别是在计算资源有限的情况下。

微生物通过介导必需元素的多种生物地球化学循环并维持生态系统稳定性，在地球生物圈中发挥关键作用。一个这样的例子是氮循环，通过该循环，不同形式的氮（例如铵、硝酸盐、亚硝酸盐、一氧化氮、一氧化二氮和氮气）在地球各圈层之间转化和循环。不同氮形式在不同氧化还原状态之间的转化由相应微生物功能基因编码的酶/蛋白质执行。例如，氨单加氧酶是由amo基因编码的一种关键氮转化酶，它催化氨氧化为羟胺，随后进一步氧化为亚硝酸盐。因此，为了研究微生物在介导各种生物地球化学循环和生态系统功能中的作用，一个关键问题是恢复生态系统中的相应功能基因。

在过去的二十年中，高通量分析技术迅速发展并常规用于研究复杂生态系统中微生物群落的组成和功能潜力。已经开发了多种主要方法，包括高密度功能基因阵列（例如GeoChip）、特定标记基因的靶向扩增子测序和总DNA的鸟枪法宏基因组测序。其中，功能基因阵列通过将环境DNA与已知序列设计的探针杂交来捕获微生物谱，因此被称为封闭格式技术。高通量测序方法，包括扩增子和鸟枪法宏基因组测序，具有发现靶向基因新序列变异的潜力。尽管这些方法差异很大，但它们在过去被广泛用于研究微生物介导的生物地球化学过程和生态系统功能。

传统上，功能基因的扩增子测序用于靶向功能集合的群落水平研究，而鸟枪法宏基因组测序具有解析靶向功能基因相对丰度变化的优势。这种不同的应用主要是由于这些方法的技术差异。具体来说，扩增子测序通常比鸟枪法宏基因组捕获更多特定靶向区域的序列变异，但由于PCR扩增偏差和无法提供整个微生物群落的全面视图，无法像鸟枪法宏基因组那样量化靶向基因的相对丰度。例如，常用的细菌amoA引物组（1F/2R）的覆盖率低至29.61%。相比之下，鸟枪法宏基因组测序可能因脱靶测序而遭受信息损失，因为它通常捕获大量非靶基因，导致恢复的靶基因较少。因此，这些不同的方法通常以互补的方式应用。值得注意的是，许多研究尚未挖掘的一个宝藏是功能基因序列变异所携带的遗传和系统发育信息，这可能为靶向功能基因的生态和进化提供独特见解。对于全长基因、操纵子和基因簇来说尤其如此。为了获得这些必要信息，需要有效恢复全长靶向功能基因甚至操纵子。

以前，已经努力通过靶向组装（即提取特定基因的读段进行组装）从鸟枪法宏基因组中恢复微生物16S rRNA和功能基因。随着宏基因组数据分析的最新进展，我们预计可以更有效地执行此过程，并显著减少计算资源和运行时间的消耗。在本研究中，以amo基因家族为例，我们旨在通过针对特定基因家族的靶向组装方法从鸟枪法宏基因组中恢复amo基因。amo基因家族编码主要存在于氨氧化细菌（AOB）和氨氧化古菌（AOA）中的氨单加氧酶，通过将氨转化为羟胺（进一步氧化为亚硝酸盐）在氮循环中起关键作用。我们比较了靶向组装与常规全组装方法的性能，表明靶向组装在计算资源和运行时间方面显著优于全组装。更重要的是，我们注意到靶向组装比全组装恢复了更高的遗传多样性和更多的操纵子，为生态和进化研究带来了更多信息。同时，我们还表明基于组装的方法比扩增子测序方法恢复了更高的靶向基因的分类和系统发育多样性。本研究展示了一种简单但有效的方法，用于从鸟枪法宏基因组中恢复具有生态重要性的微生物功能基因，并促进微生物生态研究中群落水平和核苷酸水平生物多样性的联系。

用于本研究中分析的潮间带泥滩样本于2021年6月从位于中国青岛的典型潮间带泥滩区域（东经120.75°，北纬36.46°）收集。为了调查区域微生物多样性的背景库，设计了一个嵌套采样方案。在此采样方案中，设置了多个不同半径（1、5、10、20、50、100和200米）的采样区域从一个中心点出发。对于每个增加的采样区域，收集四个样本。在低潮时沉积物暴露于空气中时，共收集了25个沉积物样本。对于每个采样点，收集五个表面土壤核心（深度约15厘米）并均质化。收集的沉积物样本立即在冰上运输到实验室。每个样本保留约200克沉积物，其中100克暂时存储在4°C直至理化性质分析，而另外100克存储在-80°C用于总DNA提取。

使用DNeasy PowerSoil Kit（Qiagen, Hilden, Germany）从0.5克沉积物样本中提取总DNA，遵循制造商的说明。使用NanoDrop ONE Spectrophotometer（NanoDrop Technologies Inc., Wilmington, DE）通过260/280和260/230 nm比率评估DNA质量，高质量DNA存储在-80°C。为了扩增氨氧化细菌（AOB）的amoA基因，使用引物对1F/2R（1F, 5′-GGGGTTTCTACTGGTGGT-3′; 2R, 5’-CCCCTCKGSAAAGCCTTCTTC-3′; K = G或T, S = G或C）。扩增子测序和鸟枪法宏基因组测序在Illumina HiSeq 2500平台（Illumina, Inc., San Diego, CA, USA）上由Magigene Biotechnology Co., Ltd.（广州，中国）进行。

原始氨氧化细菌（AOB）amoA扩增子序列使用DADA2流程（v1.34.0）处理，该流程旨在从Illumina序列数据中解析精确生物序列，不涉及序列聚类。不是通过长度去除引物，而是使用DADA2中的“removePrimers”命令从正向和反向读段中识别和去除PCR引物，允许两个错配。引物去除后剩余的不配对读段在进一步由DADA2流程处理之前被排除。包括质量过滤、样本推断、配对读段合并、去重复和嵌合体识别在内的过程使用默认参数进行。使用功能基因数据库NCycDB通过DIAMOND程序（v2.1.12）（选项：-k 1 -e 0.00001 -id 0.3）区分高质量amoA基因和pmoA基因。通过MEGAN6对amoA ASV进行 taxonomic assignment。

原始读段首先使用Trimmomatic进行修剪。使用PEAR（v0.9.11）软件合并原始Illumina配对末端读段，随后使用MEGAHIT（v1.2.9）组装宏基因组数据集的合并配对末端和单末端序列。采用并比较了三种类型的宏基因组组装方法，包括单样本组装、多样本组装和靶向组装。对于单样本和多样本组装，所有读段作为输入提供给MEGAHIT。此外，还采用了两种现有的靶向组装工具，Xander和SAT-assembler。将得到的contig与NCycDB进行比较，使用DIAMOND（选项：-k 1 -e 1 -id 0.3）识别靶向contig。对于靶向组装，首先使用DIAMOND以相同的宽松参数（选项：-k 1 -e 1 -id 0.3）将配对末端读段与NCycDB进行搜索。选择此宽松截止值以捕获更广泛的相关序列，确保包含所有与amo和pmo基因相关的读段。然后将与amo和pmo基因最佳匹配的读段提取作为输入提供给MEGAHIT和SPAdes（v4.2.0）。所有提取的读段用作靶向组装的输入，未使用丰度截止值。对于所有三种组装方法，使用了一系列kmer大小（-k-min 29 -k-max 141）。对于每次组装，保留由各种k-mer大小合并产生的共有contig数据集用于后续分析。

然后使用CD-HIT（v4.6.2）对靶向contig进行冗余去除。使用FragGeneScan（v1.18）（选项：-complete = 1 -train = complete）预测contig中的基因。使用USEARCH（v11）程序的两种方法，包括uchime_ref和uchime_denovo，使用默认参数检测潜在嵌合体。此外，还应用了Mothur（v2.20）中实现的两种方法，即chimera.perseus和chimera.vsearch，使用默认设置。当使用uchime_ref和chimera.vsearch算法时，使用NCycDB中的amo序列作为数据库判断嵌合序列。预测的编码序列还使用DIAMOND程序针对NCycDB功能基因数据库进行搜索，以进行amo和pmo基因的功能基因注释，应用更严格的阈值（选项：-k 1 -e 0.00001 -id 0.3）以确保更高特异性地识别相关基因。然后识别并保留包含amo基因的contig用于进一步分析。

使用MetaQUAST（v5.2.0）评估三种组装方法的组装质量，包括N50、N90、L50和L90等参数。使用EggNOG Database和NCBI Database进行amo操纵子功能注释。通过Kraken2（v2.1.3）程序进行 taxonomic assignments。使用Bowtie2（v2.5.4）将原始读段与contig进行比对。使用Samtools排序并将SAM文件转换为BAM格式。使用R包“vegan”（v2.7-1）计算alpha和beta多样性。DDR的空间周转率计算为对数转换的地理距离和群落相似性[基于1 –（Bray-Curtis距离度量的相异性）]之间的线性最小二乘回归关系。TAR的斜率系数基于观察到的群落丰富度在对数转换空间中使用线性回归估计。分别对细菌amo和古菌amo基因分析DDR和TAR。使用MAFFT（v7.525）对提取的基因进行序列比对。使用对齐的序列通过FastTree（v2.2）推断系统发育关系。将树文件上传到Interactive Tree of Life进行可视化。

两种主要方法，包括扩增子和鸟枪法宏基因组测序，通常用于恢复复杂环境中的重要微生物功能基因。通常，扩增子测序比鸟枪法宏基因组恢复更多功能基因变异（例如OTU和ASV），但它面临一个关键问题：恢复的遗传多样性高度受PCR引物覆盖率的限制，这对于功能基因通常很低。相比之下，鸟枪法宏基因组测序具有通过典型组装方法恢复非引物靶向基因的潜力，但恢复的基因数量和丰度通常较低。此外，通过组装方法恢复微生物功能基因存在关键问题。首先，全宏基因组数据集的组装极其消耗资源和时间，限制了从鸟枪法宏基因组中有效恢复靶向功能基因。其次，鸟枪法宏基因组组装可能偏向于生成丰富微生物功能基因的contig，而不是低丰度的功能基因。因此，迫切需要恢复覆盖高遗传多样性的微生物功能基因。

为了克服上述问题，我们测试了通过靶向组装恢复微生物功能基因，并将其与其他方法进行了比较。传统上，鸟枪法宏基因组进行全组装，无论是单样本还是多样本，以生成contig。然后对contig进行基因预测。预测的基因然后与直系同源数据库进行搜索以识别靶向基因，即本研究中的amo基因家族。在这里，使用精心策划的功能基因数据库NCycDB进行功能注释并区分amo和pmo基因。在靶向组装方法中，首先使用宽松参数将鸟枪法宏基因组读段与NCycDB进行搜索。然后提取映射到靶向基因的读段并组装成更长的序列。类似地，通过NCycDB进行基因预测和验证。然后将提取的amo基因进行进一步分析。在本研究中，还将通过鸟枪法宏基因组恢复的细菌amo基因与扩增子测序数据集进行了比较，旨在深入了解这些不同方法在恢复靶向功能基因方面可能存在的差异。

我们首先评估了不同组装方法在恢复靶向功能基因方面的计算效率，重点关注两个最重要的因素，即内存成本和运行时间。正如预期，多样本组装需要最多的计算资源和运行时间，占用720 GB内存并花费945运行小时处理包含25个样本和35.7亿读段的数据集。对于单样本组装，每个样本平均占用约54 GB内存和52运行小时，不包括从contig预测和提取靶向基因的时间消耗。相比之下，从所有样本中提取的读段进行靶向组装仅消耗0.21 GB内存并花费一分钟，加上使用DIAMOND针对NCycDB提取潜在amo读段的25小时。这些结果表明，与传统方法相比，靶向组装大大减少了计算资源和运行时间，使得能够在个人计算机甚至笔记本电脑上对靶向基因家族进行复杂的宏基因组恢复。

其次，我们基于多个指标评估了这些不同方法的组装质量，包括N50、L50、N90和L90。其中，N50和N90分别定义为50%和90%组装的contig/scaffold达到或超过该长度的长度，是用于评估组装连续性的指标。N50和N90值越高表明连续性越好，contig越长。L50和L90通过指示覆盖50%和90%组装长度的contig/scaffold数量来反映组装集中度。L50和L90值越低表示组装越集中。这里的结果表明，不同的组装方法在恢复靶向amo基因时在这些指标上的表现各不相同。对于N90和L90，靶向组装与单样本或多样本组装之间的差异很小，表明靶向组装恢复的功能基因质量与单样本或多样本组装相当。

第三，我们研究了这些方法恢复的amo基因的生物学特征，强调基因计数、恢复准确性和不同组装方法之间的序列重叠等属性。将恢复的序列与NCycDB进行搜索以区分amo和pmo基因家族。结果，不同组装方法恢复的靶向功能基因数量各不相同，特别是对于pmoB基因。例如，单样本组装恢复了42个amoA基因，而多样本组装恢复了46个，靶向组装识别出35个。恢复的amoC基因数量分别为76、62和71。具体来说，我们进一步检查了与氨氧化细菌（AOB）相关的amoA基因在扩增子测序、单样本组装、多样本组装和靶向组装中的检测率，分别恢复了117、8、15和15条序列。我们还应用SPAdes进行靶向组装，以比较其与MEGAHIT的性能。两个组装器恢复的基因数量总体相当，除了pmoB，其中SPAdes恢复的序列明显更多。此外，还比较测试了两种先前开发的靶向组装器，Xander和SAT-Assembler。值得注意的是，这些不同方法之间的比较不应过度关注，因为先前的方法主要关注单个基因，而当前的方法可扩展到操纵子/簇，但具有以减少资源和运行时间恢复靶向基因的相同思想。由于Xander提供的参考标记基因集仅包括amoA序列，我们将比较分析限制在amoA基因。具体来说，Xander总共产生了18条amoA序列，其中10条分配给氨氧化古菌（AOA），8条分配给氨氧化细菌（AOB）。同样，SAT-Assembler恢复了15条amoA（5条AOA，10条AOB）和10条pmoA序列。此外，将恢复的序列与NCycDB进行搜索以评估准确性。在靶向组装中观察到高比例具有高序列相似性的组装基因，表明与参考序列更匹配。在身份阈值0.3时，单样本和多样本组装之间恢复序列的映射比率低至38.0%，可能是由于组装策略的差异。相比之下，靶向组装与单样本组装之间以及靶向组装与多样本组装之间的比率显著高，分别为92.4%和94.1%。这些发现表明，靶向组装不仅在恢复靶向功能基因计数方面表现出与传统方法相当的性能，而且提供了更高的精度， resulting gene sequences具有高度重叠。

除了性能之外，与功能基因宏基因组组装相关的另一个潜在问题是嵌合体的产生。尽管由于缺乏PCR扩增，嵌合体在鸟枪法宏基因组中很少被关注，但高度相似序列的组装也可能产生类似嵌合体的序列。这里，使用多种嵌合体检测算法检测潜在嵌合序列。如果其与模型比对的分数超过阈值，则查询被预测为嵌合体，或者如果不确定性高则分类为未知。参考数据库由用户提供（uchime_ref和chimera.vsearch）或从提供的序列从头构建（uchime_denovo和chimera.perseus）。结果，当使用uchime_ref模式时，单样本和多样本组装都包含四个嵌合序列以及一个和三个未分类序列。相比之下，靶向组装仅产生非嵌合序列。当使用uchime_denovo模式时，在单样本组装、多样本组装和靶向组装中检测到的嵌合序列数量分别为4、9和0。单样本和多样本组装被识别出相对较高比例的不可分类序列，分别为33.3%和30.6%，而靶向组装导致低得多的百分比，为12.1%。单样本组装、多样本组装和靶向组装中非嵌合序列的比例分别为65.8%、68.7%和87.9%。当使用chimera.vsearch时，在单样本组装中检测到五个嵌合序列（0.12%），在多样本组装中检测到两个（0.05%），在靶向组装中未检测到。相应地，非嵌合序列的比例分别为99.88%、99.95%和100%。当使用chimera.perseus时，在三种组装方法中均未检测到嵌合序列。这些发现表明，靶向组装减少了来自不同源的错误组装读段的可能性。

amo操纵子由多个亚基组成，包括公认的amoC、amoA和amoB。这些亚基的翻译产物共同形成复合氨氧化酶，在将氨氧化为亚硝酸盐的过程中起关键作用，是氨氧化过程中的关键酶系统之一。在复杂环境中，amo操纵子预计比常用的靶向amoA亚基提供更多信息，以了解氨氧化微生物的生态和进化关系。在这里，努力筛选包含amoA、amoB和amoC亚基完整集的contig。结果表明，靶向组装产生了最完整的amo操纵子，总共17个，而单样本组装和多样本组装产生的数量分别为3和6。受到amoX和amoE识别的鼓舞，我们进一步将其他亚基（amoD、amoE和amoX）纳入参考数据库以检查更完整amo操纵子的出现。值得注意的是，从靶向组装数据集中成功恢复了三个包含amoABCX基因完整集的操纵子。先前的研究已经识别出多个新亚基，例如amoXYZ，除了已知的amoCAB亚基。Nitrososphaeraceae科的出现导致amo亚基基因除amoA和amoX外分散在基因组中，amoA和amoX通常连锁。本研究中通过靶向组装和单样本组装获得的contig序列也证明了amoX亚基在amo复合体中的存在。在N. europaea中，检测到基因组中有两个几乎相同的amo操纵子拷贝，由amoC、amoA和amoB（amoCAB）组成，下游跟着两个开放阅读框，即Orf4和Orf5。其中，Orf5被命名为amoD，而Orf4（也称为amoE）被描述为Orf5的精确基因复制，存在于所有β-AOB中。这两个基因，amoD和amoE，具有高度保守的序列，并且在amo操纵子中类似定位，表明它们可能编码在氨氧化中起重要作用的蛋白质。在本研究中，通过靶向组装也识别了amoE基因，为氨氧化研究提供了额外线索。

在基因组学上，生物过程通常由共转录的相邻基因编码，这些基因被称为操纵子。通常，操纵子相关基因通常属于具有相似生物学功能的进化基因家族。在本研究中，在一些序列中也观察到了操纵子相关基因。具体来说，在本研究中通过多样本组装获得的amo操纵子之前或之后很少恢复基因家族。相比之下，在单样本组装尤其是靶向组装恢复的amo操纵子附近或内部发现了诸如LlaDCHIA和ATPase等基因。整合以上信息，结果表明靶向组装在获取amo操纵子方面优于全组装。

然后我们从不同角度研究了不同方法恢复的amo基因家族的多样性，包括系统发育多样性、分类组成和生态模式。

首先，我们探索了不同恢复方法获得的amo基因的系统发育多样性。特别关注古菌amoA基因的谱系。先前的研究构建了古菌amoA基因的高分辨率系统发育树，并定义了多级古菌amoA分类。系统发育包括四个基础谱系，包括NC（Ca. Nitrosocaldales）、NS（Nitrososphaerales）、NT（Ca. Nitrosotaleales）和NP（Nitrosopumilales）。这里，通过比较先前定义的amoA系统发育枝，我们发现通过这些方法获得的潮间带古菌amoA基因被分类为NT（Ca. Nitrosotaleales）和NP（Nitrosopumilales）谱系，并且NP谱系的数量通常大于NT谱系。进一步细分显示，单样本组装包含最少的“超级进化枝”，包括NP-α和NT-γ。多组装识别了四个“超级进化枝”，包括NT-α、NP-γ、NP-η和NP-ζ，而靶向组装识别了五个，包括NT-α、NT-γ、NP-γ、NP-ζ和NP-η。这一观察结果突出了靶向宏基因组组装的优势，特别是其发现全组装未识别序列变异的能力。

此外，构建了一个系统发育树以进一步区分本研究中恢复的amoA和pmoA基因。从NCycDB获取amoA和pmoA基因的参考序列。在单样本组装中，查询和参考序列未能有效聚类，并且amoA和pmoA序列之间的区别不清晰，因为它们混合在一起。相比之下，多样本和靶向组装都导致查询和参考序列更准确的聚类。然而，一部分amoA序列仍然错误分类在pmoA进化枝内。重要的是，靶向组装比多样本组装产生更少的错误分类amoA序列，表明靶向组装可能提供更准确的amoA基因恢复和更好的amoA与pmoA序列区分。

其次，我们研究了两种典型的空间缩放模式，包括DDR和