原名：Thousands of previously unknown phages discovered in whole-community human gut metagenomes

译名：基于全群落的人类肠道宏基因组研究发现数千个未知噬菌体

期刊：Microbiome

IF：11.607

发表时间：2021.3

通讯作者：Eugene V. Koonin

通讯作者单位：马里兰州贝塞斯达国家生物技术信息中心

大部分的人类相关病毒位于远端胃肠道，由带尾双链DNA噬菌体（dsDNA噬菌体）组成，在最近的超级病毒分类学中，被归类为Uroviricota门下的Caudoviricetes。通过对模型系统的实验和未培养病毒群落的测序研究，噬菌体、细菌和其人类宿主之间的三元相互作用正在以越来越快的速度得到阐明。比较个体体内和体间的人类肠道病毒组揭示了显著的纵向稳定性和常驻噬菌体的高多样性。虽然人类肠道提供了丰富的噬菌体基因组多样性来源，但到目前为止，病毒的探索度远低于整个群落（包括病毒、细菌和古菌）。公共全群落宏基因组数据的快速增长为识别潜伏在宏基因组中的许多新的噬菌体基因组提供了契机。

根据末端基因组的排列，许多完整的噬菌体基因组组装成一个环形的contig（即一个有直接末端重复的contig）。因此，环形完整性可以作为一个特征来确定病毒体和宏基因组基因组中假定的完整噬菌体基因组的存在。然而，dsDNA噬菌体基因组规模相对较小（平均~50kb）和肠道中估计的低病毒与微生物的比率（1：10）共同揭示整个群落的噬菌体宏基因组库中只有相对较少的噬菌体DNA。此外，类似大小的质粒也可以组装成环形contig。最近开发的一种计算方法旨在通过特别关注环形噬菌体基因组的组装以及基于基因内容的质粒自动识别来解决这个问题。质粒和噬菌体之间共享的遗传库，例如，由大肠杆菌噬菌体P1和质粒编码的parABSpartitioning system，可以混淆它们基于注释的自动识别功能，因而需要手动管理。尽管存在有这些挑战，但我们迫切需要通过识别和分类噬菌体以进行基于噬菌体参考序列的分析来减少人类肠道中病毒“暗物质”的数量。

全体病毒圈最近被组织在一个全面的、统一的框架内，这个框架使用蛋白质域编码病毒标志基因来推断主要类群病毒之间的进化联系。随后，这个框架被国际病毒分类委员会（ICTV）作为全面、多层次病毒分类而批准。特别地，dsDNA病毒在其主要衣壳蛋白中具有HK97fold或doublejelly-roll fold，以及参与衣壳成熟的特别的ATP酶，因此似乎具有独立的起源，证明它们分离为两个域（最高的病毒分类单元等级）。带尾的dsDNA噬菌体和其HK97主要衣壳蛋白组成了Caudovirales目，在Caudoviricetes纲、和Uroviricota门之下（同时Uroviricota门也包括动物的远缘疱疹病毒），并进一步分为9个科。随着现在正式承认的仅从序列数据分类病毒的能力，未培养噬菌体基因组的系统发育分析可以用来鉴别分类新的类群。

在此研究中，我们描述了通过对5742个全群落人类肠道宏基因组分析发现的3738个完全组装的噬菌体基因组。以丰度、分类学和基因组组成作为选择基因组进行进一步审查的标准，详细分析了三个感染Bacteroidetes门细菌的新噬菌体科。所有这些候选科，被称为“Quimbyviridae”、“Flandersviridae”和“Gratiaviridae”由感染Bacteroidetes门细菌的噬菌体组成，前两个在人类肠道病毒组中广泛分布且具有高丰度。这些家族和其他尚未分类的噬菌体编码参与细胞对氧化应激反应的酶，这暗示噬菌体在厌氧菌对氧的耐受性中起作用。此外，比较基因组分析暴露了每个科中某些属所特有的遗传cassettes，因此这些遗传cassettes似乎是涉及噬菌体-宿主相互作用的相对较近期获取的基因。将这里识别的所有噬菌体基因组添加到公共数据库中，将大大扩大已知的噬菌体多样性，并增加人类肠道病毒的分类学内容。

1 人肠道宏基因组体噬菌体基因组的鉴定

        从人类粪便样本测序得到的5742个全群落宏基因组从NCBI Assembly数据库下载（2019年8月8日访问）。为了限制搜索可能完整的基因组的空间，从这些组件中提取了95663个环形contig（在contig尾部有50-200bp直接重叠）。接下来，来自CDD数据库的304个噬菌体特异性蛋白质排列和117个自定义排列使用hmmpress（v.3.2.1）转换为隐藏马尔可夫模型（HMMs）。Prodigal （v.2.6.3）预测了95663个contigs中的蛋白质，并使用hmmsearch对304个噬菌体特异HMMs进行搜索，放宽e值截止为<0.05。在校正重新分配的编码后，选择至少一次命中的contigs（n=4,907）进行第二轮搜索，如下所示。使用tRNA-scan-SE（v.2.0）搜索所有的contigs是否存在tRNA。在212个contigs中，发现了一种琥珀色停止密码子抑制剂tRNA。鉴于这种重新分配在人类肠道噬菌体中最常见，这些contigs的ORF被重新预测一遍，其中琥珀色停止密码子被设定为谷氨酰胺密码子。重新翻译的contigs被添加到数据库，所有contigs进行第二次HMM profile搜索，使用更严格的e值截止（<0.01）。当长度超过3kbp且具有至少有一个ORF与e值阈值以下的衣壳、门户或大型终止酶亚基蛋白分布相匹配时，contigs被归类为噬菌体。噬菌体分类用Seeker 和ViralVerify 交叉检查。在这两种工具都将contig归类为非噬菌体的情况下，人工检查蛋白质注释，其中四个病毒可信度模糊的contigs被舍弃。

2 GenBank中噬菌体基因组的收集

        从NCBI分类数据库收集与所有原核病毒对应的分类码，用于从非冗余核苷酸数据库提取超过3kbp的序列（2019年9月访问）。使用entrez direct命令行工具中的“efetch”功能检索下载每个基因组序列的蛋白质预测。缺乏蛋白质预测的基因组序列丢弃。

3 噬菌体基因组的去重和注释

        使用dRep（v. 2.6.2）及其相关依赖的Mash和FastANI，对GenBank集合和人类肠道噬菌体基因组分别以95%的平均核苷酸相似度和80%的基因组长度去重，所有其他设置都为默认设置。从这些contigs中收集蛋白质，并以95%的氨基酸身份相似度和33%蛋白质长度使用mmclust进行聚类。将代表性蛋白序列合成成一个BLAST数据库，并与CDD数据库中的多序列排列（MSAs）相比较，其中使用PSI-BLAST来比较，e值截止为0.01。如果代表性蛋白质序列产生显著结果，则使用最佳hit对蛋白质簇的代表性序列和所有组成成员进行注释。

4 系统发育重建

        大终止酶亚基（TerL）、衣壳或门户蛋白用以前的方法来排列。宏基因组噬菌体的标记蛋白与GenBank噬菌体的标记蛋白结合到一个数据库中，先使用mmclust聚集到50%的氨基酸相似度聚类里。这些聚类用MUSCLE来排列，然后聚类排列彼此之间用HHsearch（v. 3.0）来比较。聚类之间的相似度分数用公式-ln(S_A,B/min(S_A,A,S_B,B))转换为距离，其中S_A,B是profiles A和B之间的相似度。然后利用估计簇距离构造基于非加权对群方法的算术平均树状图（UPGMA）。树的tip（深度<1.5）被用来指导树的叶子上的簇用HHalign进行的成对排列，以形成更大的蛋白质簇。所产生的排列被过滤，以去除间隙大于50%且均匀性低于0.1的部位。使用FastTree，利用氨基酸演化和伽马分布位点率的Whelan-Goldman模型建立近似的最大似然树。对进化树的检查发现了353个几乎相同的PhiX-174序列，这些序列从随后的分析中作为测序试剂中的污染去除。

5 噬菌体基因组分析

        利用Vcontact2（v.0.9.19）构建了噬菌体基因组的基因共享网络，并对已去重的GenBank噬菌体基因组数据库进行默认搜索设置。这些结果被导入到Cytoscape（v.3.8）中进行可视化。

 所选噬菌群的ORF（见正文）还通过HHblits搜索以30%相似度聚类的Uniport数据库和以70%相似度聚类的PDB数据库(可在“http://www.user.gwdg.de/~compbiol/data/hhsuite/databases/hhsuite_dbs”访问，02/2020访问)。使用具有默认设置的DGRScan来检测编码预测逆转录酶(RT)的基因组是否存在与产生多样性的逆转录元件对应的重复序列。为了识别10kb以RT为中心的窗口（DGRScan的默认窗口）之外的重复，模板重复被用作BLASTn queries来查询编码基因组，用以下参数：-dust no -perc_identity 75 -qcov_hsp_perc 50-ungapped - word_size 4。

6 宏基因组中噬菌体基因组的部分丰度

        NCBI GenBank数据库中的去重噬菌体基因组与去重的肠道噬菌体结合到一个数据库中，并使用Bowtie 2索引进行reads招募。使用SRA-toolkit（v.2.10）从NCBI SRA下载了1256种人类肠道病毒组，并用fastp（v.0.20.1）进行质量过滤。质量过滤的病毒组reads用Bowtie 2使用默认设置招募到噬菌体数据库（除了“--no-unal --maxins 1000000” 设定）。按照前面描述的方法计算了每个病毒组中每个噬菌体基因组的长度标准化分数丰度。

7 来自CRISP-spacer匹配的宿主预测

        CRISPR-spacer数据库是根据以前对CRISPR-Cas系统的调查编制的。每个spacer都被用作对噬菌体基因组的BLASTN query，使用单词大小为8和禁用低复杂度过滤。如果spacer长度的95%以上与噬菌体序列的95%相同，则推断出噬菌体宿主预测。

8 抗CRISPR蛋白的预测

        抗CRISPR蛋白（Acrs）的鉴定如之前描述的一样进行。简而言之，每个蛋白都由Acr预测模型从0到1进行评分，其中较高的评分对应于该蛋白成为ACR的可能性更高。然后将蛋白质以50%氨基酸相似度聚类，如果它们符合以下标准，则被视为候选蛋白：1）平均得分为0.9或以上，2）存在于5个或更少基因的聚类中，3）至少一个聚类编码包含HTH域的蛋白，4）簇对任何PDB或CDD数据库序列不产生HHpred probability大于0.9的hit。

1 全群落人类肠道宏基因组体新型噬菌体基因组的鉴定

        我们收集了5,742个全群落组装的宏基因组，以寻找完整的噬菌体基因组。为了将搜索空间限制为完整基因组，我们仅搜索有末端直接重复（50-200bp）的环形contigs（n= 95,663），以匹配已知噬菌体标记profile（即终止酶大亚基、主衣壳蛋白或门户蛋白）的开放读取框（ORF）。总共有3,738个密码编码至少一个通过e值和长度截止标准的ORF（方法）。

        约95%平均核苷酸相似度下去重将噬菌体标记蛋白匹配的contigs数量减少到1,886个。664个contigs的子集编码所有三个标记蛋白，531个contigs编码三个标记中的两个，其余的691个contigs具有一个的可检测的标记蛋白。假定的噬菌体结构的中位数长度为44.9kb，这与最近对dsDNA噬菌体基因组中位数大小的估计一致。为了排除任何污染的contigs（例如，整合有噬菌体的质粒），结果分别用ViralVerify和Seeker进行评估，这是两种经过训练来区分噬菌体基因组和其他序列的生物信息工具。这些工具将几乎所有选择的contigs归类为具有不同可信度的噬菌体，但经过人工检查，发现这些工具存在一些假阴性分类（附加文件2）。虽然我们不能排除一些非噬菌体的contigs被错误地保留的可能性，但结果共同表明，匹配标记的环形contigs的组合主要由完整的噬菌体基因组组成。

        为了确定噬菌体的宿主范围，我们使用了一个来自原核基因组的CRISPspacer的数据库来查询宏基因组噬菌体是否有潜在的匹配物。总共发现553（29%）的去重噬菌体基因组是至少一个允许宿主预测的CRISPR-Cas系统的目标。最常见的预测宿主是Firmicutes（321噬菌体)，其次是Bacteroidetes（143）、Actinobacteria(43)、Proteobacteria（41）和Verrucomicrobia（4）。

        许多噬菌体已被发现编码抗CRISPR蛋白（Acrs）以表达对CRISPR-Cas的防御。考虑到它们在反防御中的功能，Acrs进化迅速，并与实验验证过的的Acrs表现出有限的序列相似性，这使得根据序列反推具有挑战性。然而，一种基于机器学习的方法最近被开发出来，它利用基因组内容来识别候选Arcs。该方法的应用表明，41个噬菌体，其中16个被其推断宿主的CRISPR-Cas系统发现，编码至少一个候选Acr。得分最高的Acrs属于双歧杆菌BifidobacteriumCRISPR-Cas系统靶向的四个噬菌体。这四个噬菌体在≥97%的核苷酸相似度和≥90%长度水平上与培养的双歧杆菌分离株中未被研究的原噬菌体相同，这样，通过CRISPRspacer-protospacer匹配确认了它们的宿主取向性分配。在这些噬菌体中，两个候选Arc编码基因位于大终止酶亚基和整合酶之间。Acr编码基因的定位表明他们不仅表达在最初进入宿主细胞和溶原性阶段，也表达在在过渡到分裂程序阶段防止溶原噬菌体基因组被CRISPR-Cas剪切，如实验证明李斯特菌Listeria感染的噬菌体一样。酶的转录通常由含有HTH结构域的蛋白调节，这些被称为Arc相关蛋白（Acas）。事实上，在这里确定的双歧杆菌Bifidobacterium噬菌体中，一个短的HTH结构域编码ORF位于Arcs的下游，预测到它可以在整个噬菌体生命周期中调节这两个基因的表达。虽然这些未被特征的双歧杆菌Bifidobacterium噬菌体具有Acr基因座的特征，但在本研究中鉴定的绝大多数噬菌体还没有任何可检测的Acrs，它们被CRISPR-CAS攻击。其中一些噬菌体可能编码不同的Arc，而导致我们使用的训练的先前的Arc的方法无法检测到，而另一些可能使用替代的anti-CRISPR策略。

2 肠道噬菌体的分类学分解确定了以前未知的假定家族

        为每个噬菌体基因组中编码的大终止酶亚基(TerL)、主衣壳蛋白(MCP)和门户蛋白构建了系统发育树，使用一种迭代的方法来构建底层的排列。这些树是与从NCBIGenBank数据库中提取的噬菌体基因组中提取的参考蛋白质一起构建的。反映用于识别噬菌体contigs的蛋白质，1480（78%）基因组分配给Uroviricota，360（19%）分配给Phixviricota，46（2%）分配给Loebvirae。Phixviricota包括用作测序试剂的大肠杆菌噬菌体phiX174；然而本分析中检测到的360个Phixviricota噬菌体不包括任何与phiX174密切相关的序列。剩下的分析集中在属于Uroviricota噬菌体的分类分解上，因为这些contigs是迄今为止被恢复的基因组的最大部分。

        Uroviricota门被组织成一个单一的纲（Caudoviricetes）和目（Caudovirales），但是一个新的目被提出来，它包括crAss-like噬菌体，人类肠道病毒中常见的和显然最丰富的噬菌体。我们的资料恢复了141个噬菌体基因组（从601个基因组中去重），这些基因组与crAss-like噬菌体有着系统发育关系，这是另一项单独研究的主题（Yutinet al., 在发表当中）。Caudovirales目现在分为9个科，但其中3个科（Myoviridae、Podoviridae、Siphoviridae）非常广，被证明是多系群的，因此没有用于科级的分类分配，而其余6个科仅代表GenBank中可以找到的噬菌体的一小部分。我们的Caudovirales噬菌体编码的TerL系统发育树显示只有34个肠道噬菌体属于这6个ICTV接受科中的一个。其余这些未分类的噬菌体很可能发现目前完全由未培养的噬菌体组成的新科，或属于具有培养代表物的尚未在新的病毒多级分类下定义的科。

图1：在人类肠道宏基因组中发现的Caudovirales噬菌体的三个候选科。

（A）由GenBank（=3931）和肠道宏基因组（=1298）中的Caudovirales噬菌体基因组编码的大终止酶亚基的系统发育树。分支根据目前的ICTV科着色，除了Myoviridae，Podoviridae，和Siphoviridae，它们是橙色的。最外层环表示本研究提出的候选科的位置：1、“Quimbyviridae”噬菌体；2、“Flandersviridae”噬菌体；3、“Gratiaviridae”噬菌体（见正文）。

（B）Urovicota噬菌体的基因共享网络。将人类肠道宏基因组（蓝色节点）中鉴定的噬菌体基因组与GenBank数据库中的噬菌体进行比较（颜色见图1，添加了棕色的crAss-like噬菌体和黑色的本研究中提出的新Caudovirales科）。

（C）人类肠道病毒组中噬菌体丰度。x轴描述了给定噬菌体在所有病毒中的平均丰度（n=1258）；y轴是给定噬菌体至少有一次read的病毒组的分数。每个噬菌体基因组（n=7,888）按目（C）或科（仅Uroviricota）（D）的分类水平着色。

3 为比较基因组学研究的候选科的选择

        基于噬菌体标志蛋白的分类学分析表明，人类肠道中很少有噬菌体属于目前被接受的ICTV家族。为了优先考虑候选家族进行深入分析，接下来我们用全基因组比较和每个噬菌体的丰度计算来补充具有标记的基于基因的分类分析。

        利用从宏基因组中恢复的噬菌体和储存在GenBank的噬菌体建立了一个基因共享网络。根据具有显著序列相似性的ORF的数量，在两个节点（节点代表病毒基因组）之间绘制边缘。大多数宏基因组恢复的噬菌体在网络内与GenBank噬菌体有多个连接，证明这些通过人工优化的contigs是真正的噬菌体基因组（图1B）。然而，两个大噬菌体类群（暂时标记为“Flandersviridae”和“Gratiaviridae”）与较大的网络连接较弱，反映了它们不同的基因组含量。这些噬菌体的基因含量与先前已知的噬菌体的差异以及它们在系统发育树中的不同位置(图1A，见下文)表明它们代表了新的属，也可能是新的科。

        为了量化人类肠道病毒群落中每个噬菌体的分数丰度，将1265种人类肠道病毒组的reads与包含宏基因组恢复和GenBank噬菌体的数据库进行竞争性覆盖。大多数基因组没有从超过2%的病毒组（Q1-Q3，0%-2%的病毒组）中招募到任何reads（“检测到”）（图1C-D），这与先前报道的人类肠道病毒素组的个体性差异一致。一个值得注意的例外是crAss-like噬菌体，它们招募了来自三分之一的病毒的至少一个read（Q1-Q3，9-28%），与之前关于它们广泛分布的报道一致。在人类肠道病毒组的收集中经常观察到一个未被特征的Caudovirales基因组（54%，图1C），这表明这个噬菌体也是广泛性的。为了排除观察到的频率来自到一个或几个位点的非特异性read覆盖，而不是完整基因组的可能性，我们检查了整个基因组的测序读取范围（accession OMAC01000147.1）。该基因组在病毒组中的广泛覆盖证实了它的频繁检测不是一个假象，尽管在参考序列中几个基因座在病毒组中缺席了。在人类肠道病毒群落中对这种没有被特征的噬菌体（以下称为Quimbyvirus，以《辛普森一家》中的昆比市长这个角色命名）缺少检测，这种情况需要进行详细的检查。

          因此，根据它们在标记基因系统发育树中的不同位置（全部的三个类群），结合基因内容不同（“Flandersviridae”和“Gratiaviridae”）和在人体肠道病毒群落的高丰度（“Flandersviridae”和“Quimbyviridae”），我们选择三组噬菌体进行深入分析。下面是对每个候选科的比较基因组分析。

4 “Quimbyviridae”噬菌体丰度高，序列多变，且感染拟杆菌Bacteroides

         在TerL系统发育树中，Quimbyvirus属于一组噬菌体，其近亲包括Vequintavirinae科和Ounavirinae亚科，属于现在已不存在的Myoviridae。为了阐明昆比病毒的分类学关系，对邻近分支的基因组进行了检查（图2）。类Quimby噬菌体的基因组长度为75.2kb，接近类Quimby分支底部分支的基因组大小“4986组”（72kb），但小于邻近分支其他噬菌体的基因组，包括Ounavirinae（88kb）和Vequintavirinae（145kb）。尽管基因组大小相似，但门户蛋白和MCP的系统发育重建还是将Quimby样噬菌体从4986组中分离出来。此外，大多数类Quimby噬菌体编码一个DnaG家族primase酶和DnaB家族解旋酶，这些酶在4986组中都不存在。然而，在类Quimby噬菌体的一个分支中，primase酶从复制模块中丢失了。该分支的基因组编码一种邻近DnaB家族螺旋酶的蛋白，与RepA有翼螺旋转螺旋结构域显著相似关系（HHpred概率，96.5）（图2）。RepA家族蛋白通过与宿主DnaG primase酶相互作用来介导质粒的复制，这表明RepA类蛋白在复制过程中结合宿主primase酶，触发该谱系中噬菌体编码的dnaG的丢失。与RepA介导的染色体外复制策略一致，这个分支的基因组中不能检测到编码整合酶，但噬菌体编码许多抗抑制剂，参与温和噬菌体裂解-溶原命运的蛋白质。其余的类Quimby的噬菌体都有一个完整长度的整合酶，这表明这些噬菌体可以整合到它们的宿主细胞基因组中。基于TerL、门户、MCP和DnaG树的拓扑结构，我们提出类Quimby噬菌体代表了科水平上的一个新的分类群（此后，名为“Quimbyviridae”）。复制策略（染色体外和整合到染色体中）的潜在差异，结合标记蛋白的系统发育树的拓扑结构，表明“Quimbyviridae”分裂为两个不同的亚科。

图2：大终止酶亚基的系统发育树和类Quimby噬菌体的基因组图。

（A）每个分支的右侧显示了类Quimby噬菌体和ICTV分类噬菌体的个体基因组图。ORF根据功能赋予颜色：大终止酶亚基（红色）、结构蛋白（蓝色）、DNA复制修复（橙色）、溶原性（粉红色）、一般功能（绿色）和未知（灰色）。

（B）扩展的4个类Quimby噬菌体和来自邻近分支的一个肠道噬菌体基因组（“4986组”）。生成多样性的逆转录原件和超高多样性的ORF用虚线框和星号突出显示。两个图的个体基因组图的核苷酸尺度不同。

          Quimbyvirus病毒基因组与拟杆菌Bacteroides dorei的神秘原噬菌体（CP011531.1）一致，其92%的长度具有95%的核苷酸序列相似度，这表明 B. dorei这个人类肠道微生物组的常见组成部分，携带一种与Quimbyvirus密切相关的原噬菌体。对排列的检查表明，Quimbyvirus位点特异性地整合到的B. dorei的tRNA- Asp基因中，这是原噬菌体整合的典型地点。通过CRISPR-spacer分析确定的其它Quimbyviridae噬菌体的寄主包括Bacteroidetes 门下的Prevotella、 Bacteroides和Parabacteroides，Firmicutes门下的Lachnospiraceae。相对应地，4986组的宿主不包括任何Bacteroidetes。推断宿主范围的差异支持将4986组与“Quimbyviridae”噬菌体分离，提示4986组可能代表一个新科，但这些基因组未被进一步研究。

          一些“Quimbyviridae”噬菌体携带有产生多样性的逆转录元件（DGR），这是一盒选择性的突变基因位点的基因盒；被突变的位点是可变重复，是C型凝集素或免疫球蛋白样结构域的一部分。这些结构域的靶向突变会产生具有结合亲和力和特异性改变的蛋白质。Bordetella噬菌体BPP-1中的DGR基因盒是噬菌体中唯一经过实验研究的DGR系统，其中含C型凝集素结构域的尾纤维基因的多样化使得能噬菌体够吸附到不同的宿主细胞受体。在Quimbyvirus中，DGR的RT组件通过所有三个密码框中重叠的ORF编码（ORF52-54 ），这表明有活性的RT是由两个编程的密码框位移而产生的。尽管在许多紧凑的带尾噬菌体基因组中已经发现了重叠的ORF和编程的密码框位移，但到目前为止，DGR RT只被预测为由一个ORF编码。为了辨别这些密码框位移是否使RT不活跃，可变重复用腺嘌呤特异性取代来检验，这是DGR介导变异的特征。两个可变重复序列位于Quimbyvirus基因组的ORF47和80中，它们均编码包含C型凝集素结构域的蛋白，这是DGR的典型靶点（图2）。基于与几乎相同的Quimbyvirus基因组（>95%平均核苷酸相似度）的可变重复的同源模板排列在可变重复中检测到22个腺嘌呤位点，而所有其他碱基几乎完全保守。总的来说，这些结果表明密码框位移的RT具有特征DGR介导的超变异的选择性不亲和。

          第一个可变重复基因位于ORF51的C端，它位于尾纤维基因的下游，这表明该基因编码病毒子的结构成分，类似于噬菌体BPP-1的超变异尾纤维。第二个DGR靶位点位于ORF84中，它位于噬菌体结构基因模块的远端，并从相反的DNA链表达，提示它是一种非结构蛋白。ORF80的基因组邻域包括在7kb内编码一个核酸酶、4个甲基转移酶和一个tRNA连接酶的基因。该核酸酶与大肠杆菌mutY（HHpred概率97.3）有显著的序列和结构相似，这是一种参与碱基切除修复的DNA糖化酶。甲基转移酶最类似于腺苷酸和胞苷酸修饰酶（HHpred概率分别为100和99.9），可能防止宿主限制性内切酶剪切。同样地，tRNA连接酶可能会修复由宿主抗密码子内核酸酶剪切的tRNA。总的来说，ORF84与防御和反防御相关基因的邻接意味着这种超变异噬菌体蛋白在噬菌体宿主冲突中发挥作用；然而，DGR和超变异目标蛋白在“Quimbyviridae”噬菌体在生命周期中的确切功能仍有待研究。

5 “Flandersviridae”噬菌体在整个群落的宏基因组中很常见和很丰富

        对TerL的系统发育树的分析确定了29个肠道噬菌体的深层分支（从196个基因组中去重得到）和Ackermannviridae汇集到了一起（图3A）。对29个代表性contigs编码的ORF的注释表明基因组是共线性的，证实它们属于一个内聚群（图3B）。该基因组的内聚性被基因共享网络证实，这些基因组形成了一个相干的集群，与更大的网络几乎没有联系（图1B），反映了这些噬菌体编码的大多数蛋白质和GenBank噬菌体蛋白质之间遥远的（如果有的话）相似性。该组噬菌体的平均基因组大小为85.2kb，而Ackermannviridae噬菌体的平均基因组大小为157.7kb。有一个保守的结构基因模块，编码MCP、门户、鞘和基底板蛋白、TerL和病毒子成熟蛋白酶。收缩尾鞘的存在表明，这些病毒具有与Ackermannviridae相似的收缩尾巴，与TerL系统发育一致。结构块内的一些基因包含免疫球蛋白样或C型凝集素结构域（分别为BACON和GH5），预计在病毒子粘附细菌细胞或宿主相关的粘膜糖中发挥作用。结构块的下游是一个参与DNA复制的基因模块，包括DnaB家族解旋酶、DnaG家族primase酶和DNA聚合酶I（PolA）。polA基因广泛分布在dsDNA噬菌体中，因此可以作为描述噬菌体复制模块多样性的有用标记。由这些噬菌体编码的polA和dnaG的系统发育重建证实了它们的单系性。复制模块之后是一个长度约20kb长的位点，其中包含ORF，与功能特征蛋白没有可检测到的相似性。其中两个噬菌体与Bacteroides和Parabacteroides编码的CRISPR spacers相匹配，表明这些细菌是宿主。基于大终止酶和polA系统发育、其基因组的共线性以及在基因组大小和内容上与已知噬菌体的差异，我们认为这些感染Bacteroides的噬菌体代表了一个新的分类学群，可能以科水平命名为“Flandersviridae”（以一些宏基因组取样的地区命名）。

图3：大终止酶亚基系统发育树和“Flandersviridae”完整基因组图。

（A）“Flandersviridae”科和选定的ICTV分类噬菌体成员的基因组图，其颜色如图3所示。

（B）“Flandersviridae”科中三个属的基因组图。虚线框突出了一个“Flandersviridae”噬菌体复制模块中licD和ispD家族酶的插入。

        虽然“Flandersviridae”的所有成员都是共线性的，但其中一些包含了在DNA复制模块中编码核苷酰转移酶超家族酶的两个相邻基因。一种酶属于参与异戊二烯类生物合成的ispD家族，另一种是一种licD家族酶，负责向细菌细胞壁中的磷壁酸中添加磷酰胆碱（图3）。据我们所知，这两种酶以前都没有在噬菌体中被报道过。考虑到只有“Flandersviridae”的一些成员拥有这些基因，它们不太可能在噬菌体生殖中发挥基本功能，而可能涉及噬菌体与宿主的相互作用。licD家族酶可能会修饰磷壁酸，以防止其他噬菌体的超感染，因为这些多糖是某些噬菌体吸附到宿主的受体。ispD的作用不太清楚，因为ispD家族酶催化异戊烯基焦磷酸的生物合成一步，这是多种异戊二烯类物质的基石。噬菌体操纵宿主代谢网络，包括中枢碳代谢、核苷酸代谢和翻译；在“Flandersviridae”噬菌体基因组中发现ispD可能会增加异戊二烯类生物合成途径。

        从249个全群落的人类肠道宏基因组组合中恢复了完整的“Flandersviridae”噬菌体基因组。它们频繁组装成封闭的contigs表明，这些噬菌体可能作为染色体外环形DNA分子存在于宿主细胞中，类似于噬菌体P1的存在。然而，在“Flandersviridae”基因组中，参与DNA分裂的基因和裂解溶源开关的基因都不容易识别到。因此，这组噬菌体可能是专性裂解型，尽管仅从基因组序列中区分噬菌体的生活方式是充满挑战的。无论它们的生活方式如何，这些噬菌体都能频繁地从整个群落的宏噬菌体中恢复过来，这意味着它们是人类肠道病毒体中的常见成员。事实上，“Flandersviridae”噬菌体与crAss-like噬菌体检测频率相似（图1D），尽管数据库中类Flanders噬菌体的检测频率很低。与“Quimbyviridae”一样，“Flandersviridae”基因组中平均的read覆盖度（accession OLOC0100071.1）证实其检测不是人工的。病毒组中的类Flanders噬菌体的高分度丰度和高检测率通常与它们在整个群落宏基因组的频繁组装相一致，尽管它们不是最丰富的（见讨论）。总的来说，类Flanders噬菌体代表了一个以前未被发现的噬菌体群，广泛分布在人类肠道病毒组中。

6 “Gratiaviridae”，一个假定感染Bacteroides的噬菌体新科

        一个由18个基因组组成的深层分支群（从45基因组中去重得到）是TerL系统发育树上相对Autographiviridae、Drexlerviridae和Chaseviridae的底层分支（图4A）。虽然在肠道病毒组中并不常见（图1D），但这些contigs与已建立的噬菌体家族之间的深层关系促使人们对这些假定的噬菌体进行深入的基因组分析。所有18个基因组编码DnaG 家族primase酶和DnaE家族聚合酶，这些基因的系统发育重建表明这些噬菌体是单系的，唯一的例外是噬菌体phiST的dnaE基因属于多系，phiST是一种海洋Cellulophaga感染噬菌体，属于目前不存在的Siphoviridae。dnaG和dnaE基因嵌套在其他复制相关基因的模块中，包括超家族I和II螺旋酶、SbcCD外核酸酶和一个RecA家族ATP酶（图4B）。该结构模块由编码MCP、衣壳成熟蛋白酶、门户蛋白、底板蛋白和可收缩的尾鞘蛋白的基因组成。虽然这些基因组不像“Flandersviridae”噬菌体那样具有严格的共线性，但这些噬菌体编码的蛋白质的整体相似性在基因共享网络中很明显，它们形成一个相干簇，与crAss-like噬菌体共享一些边缘（图1B）。与crAss-like噬菌体相似，CRISPR-spacer的预测宿主是Bacteroides和Parabacteroides。综上所述，这些噬菌体的系统发育和基因组组织表明，它们代表了一个新的科，暂时被命名为“Gratiaviridae”（以先驱噬菌体生物学家Dr. Andre Gratia名字命名）。

图4：大终止酶亚基的系统发育树和“Gratiaviridae”噬菌体的基因组图。

（A）ICTV分类噬菌体的基因组图，着色如图3中所示。

（B）Gratiaviridae中四个属基因组图。虚线框突出显示了一个包含HipA家族激酶结构域的蛋白、AAA家族ATP酶和糖基转移酶（参见正文）。

        除了结构蛋白和复制蛋白外，“Gratiaviridae”噬菌体还编码类铁素二铁羧酸盐超家族的几种酶。这些噬菌体编码的类铁素酶属于两个家族，即DNA保护蛋白（DPS）和锰过氧化氢酶类。锰过氧化氢酶还没有在噬菌体基因组中被记录，类DPS酶只在7个感染Lactobacillus的噬菌体中被观察到。这两种酶都参与了厌氧菌对氧化应激的耐受性。过氧化氢酶将过氧化氢解毒为氧和水，促进厌氧Bacteroides在有氧压力下生存。DPS酶催化氧和游离铁之间的反应，产生不溶性氧化铁，降低细胞内氧和游离铁的浓度，否则它们将与过氧化氢反应并产生羟基自由基，一中最具毒性的活性氧基团。“Gratiaviridae”噬菌体可能会在感染过程中部署过氧化氢酶和类DPS酶，以增强其严格厌氧的Bacteroides寄主对氧化损伤的耐受性。值得注意的是，这些酶并不限于“Gratiaviridae”，可以在包括“Flandersviridae”的196（锰过氧化氢酶）和36（DPS）的其他噬菌体基因组中鉴定出来。在肠道噬菌体基因组中对这些酶的频繁鉴定凸显了细胞内铁和活性氧物种浓度对厌氧环境中生产性感染的重要性。

        五个“Gratiaviridae”噬菌体编码一种含有丝氨酸/苏氨酸蛋白激酶结构域的蛋白质，其序列与HipA家族激酶较远但相似（HHpred99）。虽然HipA家族激酶作为多样的II型毒素抗毒素系统[95,96]的毒素成分存在于很多系统发育上不同的细菌基因组中，但噬菌体编码的蛋白激酶只有两个有特征的例子。在感染早期，类T7噬菌体磷酸化RNA聚合酶和RNA聚合酶III作为接管宿主细胞转录和翻译的机制。相比之下，大肠杆菌噬菌体933W的蛋白激酶在溶源过程中表达， [100]并介导HK97噬菌体对宿主细胞超感染的感染流产。类HipA激酶在感染期间不太可能像类T7噬菌体的激酶那样在感染早期发挥作用，因为在所有五种“Gratiaviridae”噬菌体中，激酶编码在门户蛋白和MCP基因之间，它们在许多培养的噬菌体的感染后期表达。相反，激酶可能赋予对异型噬菌体感染的免疫力，类似于933W编码的激酶。为了支持免疫相关作用，在所有五种噬菌体基因组中，AAA家族ATP酶和糖基转移酶就在激酶得上游编码（图4）。糖基转移酶编码在荚膜多糖生物合成位点中，荚膜多糖的相位变化赋予了对依赖这些分子吸附的噬菌体的免疫力。HipA家族激酶、ATP酶和糖基转移酶的具体作用尚不清楚，但是总体来讲，这些酶可能会改变宿主细胞荚膜，在“Gratiaviridae”后代病毒子的形态发生完成的同时，对异型噬菌体感染给予暂时免疫。

        对人类肠道宏基因组的搜索发现了3738个假定的完整噬菌体基因组。为了恢复完整的噬菌体基因组，这种分析将搜索空间限制在具有直接末端重复的宏基因组contigs中，它们存在于一些噬菌体基因组的末端，从而形成环形的组装基因组。具有不同复制和DNA包装策略的噬菌体，如Preplasmiviricota门或大肠杆菌Mu噬菌体，缺乏直接重复而不会产生环形的组装基因组，因此在这里没有检测到。因此，通过这种策略恢复的噬菌体基因组集既是有偏见的，也低估了实际存在的噬菌体基因组。研究结果也偏向于一个较小的基因组，它们更有可能组装成一个单一的噬菌体，尽管在一个宏基因组中确定了一个294kb的噬菌体基因组。系统发育和比较基因组分析表明，这组contigs包括许多以前未被注意到的噬菌体谱系，有些，最有可能在科的水平。

        根据两组噬菌体在宏基因组和病毒组中的频繁恢复，我们选择了这些噬菌体进行深入分析。“Flandersviridae”和“Quimbyviridae”噬菌体的完整基因组分别聚集在249个和20个全群落宏基因组中。然而，Quimbyvirus比任何Flandersviridae噬菌体在病毒组中的检测频率更高（图1D）。这种差异可归因于几种因素，包括抽样偏差、参考数据库中更多的“Flandersviridae”基因组“稀释”了每个基因组的覆盖reads数，或者可变位点的存在断开了contigs的拼接（例如DGR的可变重复序列）。无论如何，这两组噬菌体都包含了大量的人类肠道病毒的成员。可以预测到，这些噬菌体的宿主包括Bacteroides，它是人类肠道中最主要的细菌类群，是其他常见的人类肠道噬菌体的宿主。这些噬菌体基因组中许多没有特征的“暗物质”很可能致力于防止这些噬菌体对Bacteroides寄主细胞的超级感染，并对抗寄主的防御。虽然在Bacteroidetes的一般防御系统研究不多，但大多数细菌含有具有活性的CRISPR-Cas系统，并检测到许多针对这里分析的噬菌体的CRISPR-spacer。这意味着许多如果不是绝大部分感染Bacteroidetes的噬菌体会编码Acrs。然而，目前在先前确定的Acrs序列上训练的预测方法只在少数这些噬菌体中检测到假定的Acrs。其余的噬菌体可能编码不同的Acrs或采用替代的抗CRISPR策略。

        一些噬菌体属拥有DGRs，包括Quimbyvirus。宏基因组调查表明，DGRs在居住在胃肠道环境的病毒中富集。结合将携带DGR的噬菌体从人类肠道菌中分离出来，这些观察结果反映了在胃肠道环境中基于超变异性的噬菌体-宿主互作的显著作用。值得注意的是，Quimbyvirus和另一种携带DGR的噬菌体（Hankyphage，BK010646.1）溶原性寄宿在同一Bacteroides中，这两种噬菌体经常在人类肠道病毒组中检测到。除了共性之外，Quimbyvirus的DGRRT由三个重叠的阅读框编码，并针对两个蛋白质，一个在结构模块中，另一个在与防御相关的基因岛中。DGRs在蓝藻和Gamma变形菌门基因组中与假定的防御和信号系统有关，但除了C型凝集素折叠的存在，超可变蛋白很少具有其他可识别的域来模糊它们的精确作用。

        本研究分析的第三组“Gratiaviridae”并不丰富，但在TerL树中相对于Autographiviridae、Chaseviridae和Drexlerviridae占据较深的位置。对“Gratiaviridae”基因组的分析将有助于这些科未来组织成更高的分类等级，可能是目级。此外，对“Gratiaviridae”基因组的分析表明存在过氧化氢酶和DPS家族酶来仲裁细胞对氧化应激的反应。氧浓度随胃肠道长度而变化，其中远端比近端肠道的浓度较低。氧也从组织径向扩散到流腔，并结合其他因素，这些梯度影响胃肠道微生物群的结构和组成。通过“Gratiaviridae”和其他肠道噬菌体获取氧解毒酶，表明需要补充其宿主细胞对感染期间氧化损伤的耐受性，这对位于氧暴露率较高的组织表面附近的细胞可能尤为重要。

        一些“Gratiaviridae”噬菌体的一个独特特征是一种HipA家族蛋白激酶。类T7噬菌体（来自Autographiviridae）和噬菌体933W（目前在科水平未分类）分别编码在宿主细胞接管和流产感染期间发挥作用的PKC家族蛋白激酶。第三种CotH家族蛋白激酶结构域偶尔在噬菌体基因组中观察到，它融合到一个超变异的C型凝集素结构域，但这些蛋白质目前尚未被研究。“Gratiaviridae”噬菌体家族招募自第四个蛋白质激酶家族，这些激酶和噬菌体编码的糖酵酸转移酶，可能改变宿主细胞荚膜，有助于防止超感染。

        综上所述，对本文所描述的噬菌体的比较基因组分析，以及对crAss-like噬菌体的补充分析，显著地增加了噬菌体的特征多样性，主要是感染Bacteroidetes的细菌，这些噬菌体是人类肠道微生物组的主要组成部分。这些发现还扩大了噬菌体基因功能的范围，特别增添了异戊二烯类代谢途径、类过氧化氢酶、HipA家族蛋白激酶和与防御相关的高变异基因，为实验研究开辟了多个方向。