CHM：成年同卵双胞胎的病毒组多样性与肠道微生物组多样性相关

Virome Diversity Correlates with Intestinal Microbiome Diversity in Adult Monozygotic Twins

翻译：秋芒树 英国帝国理工学院

责编：刘永鑫 中科院遗传发育所

原文链接：https://www.sciencedirect.com/science/article/pii/S1931312819300599

Cell Host and Microbe    [IF:17.872]

DOI: 10.1016/j.chom.2019.01.019

Resource 2019-02-13

图形摘要：病毒组仍然是微生物组中相对未被探索的组分。Moreno-Gallego和Chou等人，在没有考虑宿主遗传变量的情况下，研究了同卵双胞胎的病毒组，探究微生物多样性如何与病毒组多样性相关。按微生物组一致性的高或低对同卵双胞胎进行分类，揭示了与病毒组的相关性。

热心肠日报

链接：https://www.mr-gut.cn/papers/read/1087944415

Cell子刊：肠道菌群可预测病毒组结构

创作：米见对 审核：小肠君 03月01日

原标题：病毒多样性与成年同卵双胞胎肠道菌群多样性的关联

在21对有相似或不同肠道菌群多样性的成人同卵双胞胎中，肠道病毒组具有高度特异的个体差异型；

噬菌体在肠道病毒组中占优势地位，主要为长尾噬菌体目、微小噬菌体科、短尾病毒科和长尾病毒科；

菌群的丰度和多样性与病毒的丰富度和多样性存在对应性，与菌群相似的双胞胎相比，菌群不同的双胞胎有更多不同的病毒，菌群的多样性越高，病毒组亦然；

病毒组的多样性模式由噬菌体驱动，而非真核病毒；

肠道菌群结构可用于预测病毒组结构。

主编评语：人类肠道病毒组具有高度变异性。本研究通过比较同卵成年双胞胎的肠道菌群与病毒组多样性之间的关联，发现尽管病毒组具有高度宿主特异性，其结构和多样性与肠道菌群存在对应关系，并强调了噬菌体在肠道病毒组中的核心作用，值得参考。

研究总结

病毒组是人类肠道微生物子中最可变的组分之一。在双胞胎中，病毒组是相似的，但成年人却不一样，这表明随着双胞胎年龄的增长，他们的环境和微生物组发生分歧（diverge），病毒组也一样发生分歧（diverge）。微生物组在多大程度上驱动了巨大的病毒组多样性尚不清楚。在这个研究中，我们研究了21对成年同卵双胞胎中微生物组和病毒组多样性之间的关系，这些双胞胎按微生物组一致性的高或低来分类。源自类病毒颗粒（virus-like particles）的病毒组对于每个个体都是独一无二的，Caudovirales 和 Microviridae 占优势，并显示出包含crAssphage的小核心（small core）。与微生物组具有一致性的双胞胎相比，在不具有一致性的双胞胎中显示出，微生物组越丰富，病毒组就越丰富。这些模式是由噬菌体而不是真核病毒驱动的。总的来说，这些观察支持微生物组在病毒组模式中的强大作用。

研究背景

人类肠道微生物组由大量的细菌，连同少数古菌和真核细胞，共同形成一个密度非常高的微生态系统(每克粪便10¹¹ - 10¹²个细胞) 。微生物组的细胞和病毒组的成分(每克粪便中10⁹ - 10¹²个类病毒颗粒（virus-like particles）的比例大致相等。病毒组主要由噬菌体和溶源噬菌体组成，也包括更罕见的真核病毒和内源性逆转录病毒。目前，大多数噬菌体在数据库中没有匹配，它们的宿主仍有待阐明。将噬菌体与宿主相匹配是一项挑战：例如，最常见的人类肠道噬菌体——crAssphage——的宿主最近才被鉴定为Bacteroides的物种。除了对宿主的鉴定，还存在其他问题，比如塑造病毒组最重要的因素，以及微生物组的细胞部分对病毒组的预测程度。

普遍认为噬菌体及其宿主在时间序列上的种群动态可能是相关联的。事实上，研究者描述过水生系统中的病毒及其细菌宿主的群体振荡，这些研究表明病毒在调节细菌种群（population）方面起关键性作用。但是这种捕食者-被捕食者模式（predator-prey）对于人类肠道病毒组和微生物组中并不典型。为了清楚起见，从这里开始，我们使用“微生物组（microbiome）”来指代微生物组的细胞部分（cellular fraction of the microbiome），例如，主要是细菌细胞。尽管如此，病毒组和微生物组确实显示出相似的跨宿主的多样性模式，例如人与人之间的（interpersonal）高度的差异和相对长期的稳定性。与无关（unrelated）个体相比，有联系（related）的个体的微生物组更相似，这可能是由于共同的饮食习惯推动了微生物组之间的相似性。相应地，饮食与病毒组多样性相关，很可能这是通过饮食对微生物组的影响。在婴儿中，同卵双胞胎之间的病毒组比无关（unrelated）个体之间的病毒组更相似。在成年双胞胎中没有观察到这种模式，可能是因为它们的微生物组已经发生了分歧（divergence）。由于宿主相关性等混杂因素，微生物组本身在多大程度上驱动宿主间病毒组的多样性模式难以评估。

在这个研究中，我们关注成年同卵双胞胎的肠道微生物组，以进一步探索微生物组与病毒组多样性之间的关系。通过研究同卵双胞胎的病毒组，我们控制了宿主的遗传相关性。尽管与双卵双胞胎或不相关的个体相比，同卵双胞胎通常具有更相似的微生物组，但是同卵双胞胎仍然可以显示出大量双胞胎内的微生物组多样性。我们以前从TwinsUK队列（cohort）中收集了双胞胎的粪便微生物，基于这一信息，我们选择了具有高度一致性或高度不一致性的微生物组的双胞胎。我们从类病毒颗粒（virus-like particles）中产生病毒，这些类病毒颗粒（virus-like particles）是从与微生物组来源相同的样品中获得的。结果表明，微生物组多样性和病毒组多样性的指标呈正相关关系。

研究结果

选择微生物组一致或者不一致的同卵双胞胎

Selection of Microbiome-Concordant and -Discordant Monozygotic Twin Pairs

我们选择了具有相似体重指数( BMI )的双胞胎，基于之前获得的16S rRNA基因测序数据，它们的微生物组样本间多样性( β-多样性)要么一致，要么不一致。这项研究中的成年同卵双胞胎不共享一个家庭（household），我们假设双胞胎之间的其他环境变量相似。我们基于三个β-多样性距离确定同卵双胞胎微生物组之间的一致或不一致程度，这三个距离是Bray-Curtis、加权（weighted）UniFrac和未加权（unweighted）UniFrac。与预期一致，β-多样性是相关的（Pearson相关系数>0.4）。基于成对距离(pairwise distance)的分布，我们从所有三种分布的边界中选择了21对同卵双胞胎（图1A），同时保持年龄和体重指数在整个组中的平衡。在所选的21对双胞胎中，微生物组具有一致性的同卵双胞胎之间比不具有微生物不一致的更相似(p = 6.31 X 10 ^-12，两类相似度分组统计)。不具有微生物组一致性的同卵双胞胎的微生物组在所有的分类学水平上都有不同的组成，特别是在门的水平上，Firmicutes和Bacteroidetes是两个主要的门，对同卵双胞胎之间的差异贡献最大(图1B/C)。

图1. 双胞胎的微生物组差异

( A )来自之前研究中的354对单卵双胞胎微生物组的β-多样性。每个点代表一对双胞胎的β-多样性，使用加权（weighted）UniFrac (x轴)、未加权（unweighted）UniFrac (z轴)和Bray-Curtis (y轴)方法测量的。平面是最小二乘法（the least squared）拟合的Bray-Curtis ~ Weighted UniFrac + Unweighted UniFrac平面。从两个边缘（edges）中选择具有微生物组具有一致性(蓝色)和微生物组不具有一致性(橙色)的同卵双胞胎子集。黑点表示用于病毒组和整个粪便宏基因组比较的样本。

( B )比较21对同卵双胞胎在门水平上的分类学概况(相对丰度)，1 – 9组的微生物组具有一致性，10 – 21组的微生物不具有不一致性。

( C )微生物组具有一致性(蓝色点，样本量n = 9 )和不具有一致性(橙色点，样本量n = 12 )同卵双胞胎的主要门水平的相对丰度的差异。使用了Mann-Whitney’s U test. **表示p < 0.0005, 表示p = 0.055。

类病毒颗粒（VLPs）的鸟枪法宏基因组测序

Shotgun Metagenomes of VLPs

我们从用于16S rRNA基因多样性分析的相同粪便样品中分离出类病毒颗粒（virus-like particles）。从类病毒颗粒（virus-like particles）中提取的DNA用于全基因组扩增，然后进行鸟枪法宏基因组测序。第一个大片段文库（“large-insert-size library”）选择平均插入大小为500 bp (总共34,325,116对序列reads；质量控制后每个样品817,265 ± 249,550对序列)，用于病毒重叠群(contigs)的重新组装（de novo assembly）。第二个小片段文库平均插入大小为300 bp的较小片段在(“small-insert-size library”)中纯化并测序。最终的pair-end 序列合并为2,5324,163条较长的过滤过质量的序列，以提高比对（mapping）精确度。

识别推定的细菌污染

Identification of Putative Bacterial Contaminants

从类病毒颗粒（virus-like particles）制备和测序的病毒可能被细菌DNA污染。然而，鉴于噬菌体是水平基因转移的主要媒介，在溶源噬菌体状态下，温和病毒通常占细菌基因组的10 %，去除潜在的细菌污染也可能去除病毒序列。为了评估细菌DNA的污染，我们把病毒序列比对（mapped）到一组8163个完全组装的细菌基因组上。我们的策略是评估每个基因组全部长度上每个100kb滑窗的覆盖率，那些覆盖率中值大于100的基因组被认为是污染。比对到短区域的序列被认为是溶源噬菌体或水平转移的基因而被保留下来(图2A)。比对（map）到基因组的序列且被确定为潜在污染物，去除它们之后再进行进一步的分析。

我们识别了65个细菌基因组是潜在污染，每个样品的序列中1 % ± 1.125 %比对（mapping）到细菌基因组。大多数( 37 / 68 )属于Firmicutes门；在物种水平上，在至少50 %的样品中检测到Bacteroides dorei, B. vulgatus, Ruminococcus bromii, Faecalibacterium prausnitzii, B. xylanisolvens, Odoribacter splanchnicus和B. caecimuris。如果微生物组中最丰富的细菌种类是最可能的污染源，那么它们作为污染物的相对丰度应该与它们在微生物组中的相对丰度相对应。然而，我们观察到污染的DNA和微生物组中代表分类群的相对丰度之间没有显著的相关性(图2B)。

图2. 类病毒颗粒（VLP）在准备过程中的细菌污染

( A ) 单个样品（4A）在去除被认定为污染的read之前(上图)与之后(下图)的类病毒颗粒（VLP）比对到细菌基因组的热图。细菌基因组用竖条表示，按长度排序，分成了100,000bp大小的bins。中值覆盖率超过了100的序列被认为是污染。左边的色标显示了bin覆盖率，散点图显示了每个基因组bin覆盖率的中位数。

( B )从所有类病毒颗粒（VLP）提取物中识别出来的65个污染基于NCBI的进化树图。右图: 类病毒（VLP）提取物中细菌基因组的丰度与微生物组中16S rRNA基因图谱之间的Spearman相关系数（rank correlation co-efficient (rho)）。左图:所有个体中细菌基因组的总丰度。

功能组成支持病毒在类病毒颗粒纯化中的富集

Functional Profiles Support Viral Enrichment in VLP Purification

为了评估病毒组的功能性内容，我们使用整合基因集Integrated Gene Catalog，IGC)的KEGG对“short-insert-size”文库的原始序列进行了注释。根据以前的研究，大部分宏基因组类病毒颗粒（virus-like particles）的序列( 85.43 % ± 5.74 % )被比对(mapped)到具有未知功能的基因(图3A)。

为了进一步验证序列来源于类病毒颗粒（virus-like particles）而不是微生物组，我们生成并且比较了四个补充个体的类病毒颗粒（virus-like particles）和大量粪便样品的宏基因组数据。与预期一致，来自相同样品的病毒组和微生物组的功能分布是不同的。比对（map）到注释基因的病毒组序列在两个类别中富集：遗传信息过程（Genetic Information Process）( 48.87 % ± 12.12 % )和核苷酸代谢（Nucleotide Metabolism）( 17.59 % ± 8.81 % )，相比之下微生物宏基因组中分别为24.31 % ± 1.28 %和5.47 % ± 0.4 %（图3B）。基本上，细菌宏基因组中存在的大多数其他功能类别在病毒组中都不存在。此外，相对于微生物组，病毒组的功能注释显示样本间的更高的变异，更低的组内相关系数（intraclass correlation coefficient，图3B）。

图3.  粪便宏基因组和病毒组的功能基因（gene content）比较

( A ) 粪便宏基因组和病毒组中KEGG类别的相对丰度，包括整合基因集(Integrated Gene Catalog，IGC)的所有命中（hits），不考虑注释。

( B )粪便宏基因组和病毒组在KEGG注释第二级类别中相对丰度的热图，不包括带有未知注释的整合基因集基因。

色标显示相对丰度的平方根。A.V.表示附加的病毒组；A.M.表示附加的微生物组(全基因组提取)。A.M.的Intra-class coefficient(ICC) 值为0.99，A.V.的ICC值为0.85。具有微生物一致性的同卵双胞胎的ICC值0.69，不具有一致性的ICC值为0.68。

病毒组是个体独特的

Viromes Are Unique to Individuals

我们从“large-insert-size library”中组装了序列，总共产生了大于500 nt 的107307条重叠群 (最大值：79,863 nt；平均值，1,118±1741 nt)。为了评估病毒组的结构和组成，建立了一个比对到去重复的重叠群的序列补充矩阵（a matrix of the recruitment of 序列 against dereplicated 重叠群 was built）。这个补充的矩阵包括14584条，既长(＞1300 nt )，覆盖面又广( ＞ 5X )的重叠群，他们被称为“病毒类型（virotypes）”。对补充矩阵的分析表明，每个个体都有一组独特的病毒类型（virotypes）： 3415种病毒类型只存在于一个个体中(占总数的23.41 % )， 413种病毒类型存在于至少50 %的个体中 ( 2.83 % )，只有18种病毒类型存在于所有个体中 ( 0.1 % )。

微生物组具有一致性的双胞胎共享病毒类型

Twins with Concordant Microbiomes Share Virotypes

我们检查了双胞胎之间共享的病毒类型（virotypes），并观察到同卵双胞胎共享的病毒类型（virotypes）并不比无关的个体多。然后，我们分别评估了微生物组具有一致和不具有一致性的双胞胎：微生物组不具有一致性的双胞胎与无关个体相比不共享更多的病毒类型( p = 0.254 )，而微生物组具有一致性的双胞胎确实比无关个体共享更多的病毒类型( p = 0.048 )。此外，我们还发现微生物组具有一致性的双胞胎比微生物组不具有一致性的双胞胎共享更多的病毒类型( p = 0.015 )。

在肠道病毒组中噬菌体占优势

Bacteriophage Dominance of the Gut Virome

为了描述病毒组的分类组成，我们试图使用投票系统方法（a voting system approach that）对所有66,446个去重复并且有良好覆盖率的重叠群进行注释，该方法利用组装后的重叠群及其编码蛋白中的信息。此外，我们对两个高度丰富的肠道相关噬菌体家族（families）进行了定制注释: ( 1 )crAssphage；( 2 ) Microviridae科。为此，我们使用隐马尔可夫模型( Hidden Markov Models ，HMMs )来搜索crAssphage (双链DNA，dsDNA )病毒和Microviridae科(单链DNA，ssDNA )病毒重叠群。

HMMs使我们能够识别远源的同源物（distant homologs），然后我们将它们与已知的参考序列整合到系统发育树中，以确认注释并更好地解析分类。我们注释了108个重叠群( 19个crAssphage，90个Microviridae科)，验证了68个重叠群在科水平的分配（assignments），并为97个重叠群分配了一个亚科。对于Microviridae科，只有11个重叠群以前有分配的物种分类（taxonomic assignment），都属于Gokushovirinae，我们确认了这些分配，另外有23个重叠群为Gokushovirinae，54个重叠群为Alpavirinae，1个重叠群为Pichovirinae。对于crAssphage，11个重叠群与原始的crAssphage聚集在一起，3个重叠群与Chlamydia噬菌体参考基因组聚集在一起，5个重叠群与IAS病毒聚集在一起。

在校对了投票系统注释和HMM注释之后，总共有12751个重叠群( 29.62 % )被分配物种分类。病毒组以噬菌体为主，只有6.42 %的重叠群被注释为真核病毒。与预期一致，大多数重叠群( 96.98 % )是dsDNA病毒，而只有2.43 %的重叠群被标注为ssDNA病毒。Caudovirales是最丰富的目，其三个主要科分别为: Myoviridae ( 20.22 % ± 4.83 % )、Podoviridae ( 10.54 % ± 3.27 % )和Siphoviridae ( 35.25 % ± 7.19 % )。crAssphage家族平均构成重叠群的13.26 % ( ± 12.24 % )，在一个病毒组中最大的分布可以达到55.80%，Microviridae占病毒群的3.87 % ± 2.57 %。有趣的是，我们观察到Phycodnaviridae超过平均丰度的1 % ( 1.77 % ± 1.12 % 图4A)，以及与核质大DNA病毒(nucleocytoplasmic large DNA viruses, NCLDV)相关的重叠群的平均相对分布为3.99 % ± 2.22 %。在所有样品中都存在的18个重叠群包括10个标记为crAssphage的重叠群，2个标记为“未分类的Myoviridae”，2个标记为“未分类的Caudovirales”，1个标记为Microviridae，3个未分类。在每个样本定义的分类概况中，我们在所有分类水平上寻找具有微生物组一致和不具有微生物组一致性的双胞胎的病毒组组成的差异。任何分类群在目和科水平上没有显著差异，包括crAssphage和Microviridae科(图4B)。

图4.病毒组组成

21对同卵双胞胎的微生物组在科水平的分类特征比较，1 - 9 组是具有一致性的微生物组，10 – 21组是不具有一致性的微生物组。

( A )同卵双胞胎的病毒组在科水平的组成。

( B )微生物组具有一致性(蓝色点，样本量n = 9 )和不具有一致性的(橙色点，样本量n = 12 )同卵双胞胎在每个科水平的相对丰度差异。

我们使用CRISPR比对微生物与噬菌体数据库(microbe-versue-phage，MVP)预测病毒类型（virotypes）和分类特征重叠群（taxonomically characterized）的细菌宿主。由于注释噬菌体的宿主，我们没有获得任何注释为真核病毒的重叠群信息。这些方法使我们能够识别910个重叠群推定的（putative）宿主。在这910个重叠群中，只有一个先前被注释为crAssphage，与预期一致，其宿主被推断为Bacteroidetes的物种。我们总共鉴定了1280个推定的（putative）细菌宿主菌株，包括来自多个门87个属的187个物种：其中大多数来自Firmicutes门（92），其次是Bacteroidetes门( 41 )和Proteobacteria门 ( 38 )。每个重重叠群的宿主数量中值为1（IQR = 1–2 )，而在菌株水平上，每个宿主的噬菌体数量中值为2 (IQR= 1–3 )。

病毒多样性与微生物多样性相关

Virome Diversity Correlates with Microbiome Diversity

为了评估病毒组和微生物组多样性之间的关系，我们使用从测序数据中获取的三个不同层次的信息检查了病毒组的样本内多样性(α-多样性)和β-多样性: ( 1 )病毒类型（virotypes），( 2 )分类注释的重叠群，( 3 )从短序列中注释的基因。

α-多样性

微生物组和病毒组的α-多样性在用于测试相关性的三层信息中的两层中呈正相关(病毒类型（virotypes）和分类注释的重叠群，图5A)。我们使用带注释的重叠群来查询病毒亚组（subgroups）(真核ssDNA、真核dsDNA、细菌ssDNA和细菌dsDNA)中的α-多样性。结果表明真核病毒的多样性与微生物组α -多样性无关。相反，ssDNA或dsDNA的噬菌体与微生物组α-多样性呈正相关关系。

图5. 噬菌体多样性与微生物组多样性相关，但与真核病毒多样性不相关

( A )病毒组香农α-多样性（Shannon α-diversity）与微生物组香农α-多样性（（Shannon α-diversity））的相关性( 样本量n = 42 )。绘制了线性回归95 %置信区间的最佳拟合线。病毒类型（Virotypes）：皮尔逊相关系数（Pearson correlation coefficient）= 0.406，m= 0.3，p = 0.007，R2 = 0.165。分类学：皮尔逊相关系数 = 0.389，m = 0.25，p = 0.010，R2 = 0.151。基因：皮尔逊相关系数 = 0.105，m = 0.11，p = 0.506，R2 = 0.011。

( B )根据注释为真核ssDNA病毒、ssDNA噬菌体、真核dsDNA病毒和dsDNA噬菌体的重叠群计算的病毒组香农α-多样性（Shannon α-diversity）与微生物组香农多样性（Shannon α-diversity）的相关性( 样本量n = 42 )。绘制了线性回归95%置信区间的最佳拟合线。真核病毒: 皮尔逊相关系数= 0.027，m= 0.034，p = 0.863，R2= 0.000751。sDNA噬菌体: 皮尔逊相关系数= 0.394，m= 0.35，p = 0.009，R2 = 0.155。dsDNA真核病毒: 皮尔逊相关系数= 0.143，m= 0.15，p = 0.368，R2 = 0.020。dsDNA噬菌体: 皮尔逊相关系数= 0.400，m = 0.25，p = 0.008，R2 = 0.16。

β-多样性

我们观察到，当使用Hellinger距离时，与不具有一致性的微生物组的双胞胎相比，具有微生物组一致性的双胞胎病毒组β-多样性较低；平均无权重Jaccard 距离和Bray-Curtis距离也显示出相同的趋势。与我们观察到的α-多样性相似，无论使用哪一层信息，具有一致性的微生物组的同卵双胞胎中病毒组的平均Hellinger距离明显低于不具有一致性的微生物组的同卵双胞胎（图6）。把具有微生物组一致性的双胞胎或不具有一致性的微生物组的双胞胎按性别分开时，我们没有观察到β-多样性的显著差异。尽管如此，任何关于性别影响的推断都是有局限的，因为每个群体的个体数量都减半了。此外，当使用注释重叠群时，在微生物组和病毒组β-多样性之间观察到显著正相关。这种关系是由噬菌体而不是真核病毒驱动的。

图6. 病毒组β-多样性模式反映了微生物组的β-多样性

根据分析的三个不同的信息层(病毒类型（virotypes）、基因（genes）和分类学（taxonomy）)，箱线图显示了微生物组和病毒组的Hellinger距离分布，对于微生物组具有一致性的同卵双胞胎(蓝色，样本量n = 9 )，微生物组不具有一致性的同卵双胞胎(橙色，样本量n = 12 )，微生物组具有一致的同卵双胞胎中的无关样本(蓝色，样本量n = 144 )，以及微生物组不具有一致性的同卵双胞胎中的无关样本(橙色，样本量n = 264 )。均值之间的显著差异(Mann-Whitney’s U test, p < 0.020) 用不同的字母表示。

最后，我们比较了相关(同卵双胞胎)和无关个体之间的病毒组和微生物组的成对距离（pairwise distance）。成对距离矩阵显示病毒组和微生物组β-多样性测量值之间的正相关，不仅在双胞胎之间，在所有个体之间都是如此。这些结果表明，不管宿主之间的遗传相关性如何，个体的微生物组越相似，病毒组也会越相似。

研究讨论

同卵双胞胎，像其他兄弟姐妹一样，通常比无关个体拥有更相似的肠道微生物组。此外，同卵双胞胎总体上比异卵双胞胎具有更相似的微生物组，尽管在整个微生物组水平上，这种影响很小，这主要是由一小组可遗传的微生物驱动的。然而，在同卵双胞胎群体中，微生物组中双胞胎内部差异可能和异卵双胞胎的一样大。我们利用同卵双胞胎β-多样性的巨大差异，选择了肠道微生物组高度一致或不一致的同卵双胞胎。我们对它们的病毒环境的分析表明，尽管个体间肠道病毒环境差异很大，而且不管宿主的相关性如何，他们的微生物环境越不相似，他们的病毒环境就越不相似。这种模式是由病毒组的噬菌体驱动的。

通过从微生物组的β-多样性分布中选择同卵双胞胎，我们去除了宿主遗传相关性作为可能影响病毒组的变量。以前对婴儿双胞胎病毒组和微生物组的研究表明，同卵双胞胎的微生物组和病毒组比无关个体的更相似，这表明共享的宿主基因型和/或环境是关键。相反，一项对成年双胞胎病毒组的研究表明，成年同卵双胞胎没有比不相关的个体拥有更相似的病毒组；然而，根据目前的研究结果，这可能是一个统计功效（power）问题。事实上，在我们的数据集中，我们观察到无论双胞胎的微生物组一致还是不一致，同卵双胞胎比无关个体都具有更多相似的病毒组病毒类型（virotypes）和物种分类（taxonomy）。

之前的研究表明，与成年双胞胎相比，年轻双胞胎的病毒组相似性更高，这与婴儿双胞胎共享更大的环境有关，特别是在饮食方面。Minot等人的研究也表明，相同饮食的个体比不同饮食的个体具有更相似的肠道病毒组。众所周知，饮食是日常微生物组波动的强大驱动力，所以饮食对病毒组的影响可能是由微生物组介导的。然而，我们没有控制饮食，所以我们观察到的微生物组不一致可能是由于双胞胎在取样时饮食不同造成的。不管影响微生物组一致性的差异的基础是什么，它都与病毒组的一致性密切相关。

以前没有直接讨论过成年人的病毒组丰富度（richness）和微生物组丰富度之间的关系。我们观察到，使用描述病毒多样性的三层信息中的两层，微生物组和病毒组的α-多样性正相关。具体来说，这种模式在病毒类型（virotypes）和物种分类学（taxonomy）上被观察到，但在基因（gene）上没有。然而，由于观察到病毒基因仅在两类中得到富集，即遗传信息处理（Genetic Information Processing）和核苷酸代谢（Nucleotide Metabolism），我们预计受试者之间病毒组基因的多样性不会有差异。分类注释层面（taxonomic annotation layer）的信息表明，是病毒的噬菌体组分而不是真核病毒推动了这种α-多样性相关的模式。

病毒组和微生物组多样性之间的正相关关系表明，宿主的可获取性（availability）越强，病毒的多样性就越高。这些观察与“piggy back the winner”模型一致，该模型假设在密度大的环境中，噬菌体选择进入溶源循环，并与其宿主一起复制。事实上，对人类肠道病毒组的纵向研究已经报道了与溶源性相关的基因，温带型的重叠群随时间（over time）的低突变率，以及病毒组的长期稳定性，表明了对溶源性周期的偏好。然而，噬菌体捕食行为（phage predation）被认为是维持高度多样性和高效生态系统的一个重要因素，并可能在快速变化的生态系统中扮演维持多样性的角色，如人类肠道中。对病毒组-微生物组相互作用的短期时间序列（time series）分析，以及对噬菌体繁殖中溶源-裂解转变（lysogenic-lytic switch）的更好理解，将有助于解释在人类肠道病毒组中观察到的模式。

这里描述的病毒组的组成与以前报道的成人粪便病毒组的组成有相似性。从经过注释的部分来看，Caudovirales目及其科Siphoviridae、Myoviridae和Podoviridae以及crAssphage是所有样本中的优势噬菌体。Manrique等人将婴儿肠道的噬菌体定殖过程总结如下：真核病毒首先占据新生儿肠道，随后是Caudovirales，在2.5岁前Microviridae开始成为优势噬菌体。的确，在我们的样本集中观察到了丰富的Microviridae，但是Caudovirales是优势群体。年龄与这里研究的成人受试者的多样性模式无关。

尽管这里描述的每个病毒组具有高度的多样性和独特性，但我们在受试者中发现了一组病毒：所有样本中都存在18个重叠群。这些重叠群中有一半以上被注释为crAssphage，这与最近这种噬菌体广泛存在的报道一致。我们的数据集中其他共有的病毒类型（virotypes）被分类为Myoviridae和Microviridae。我们还发现了比对到以核质大DNA病毒、Phycodnaviridae和Mimiviridae为代表的科水平的重叠群。这些类型的病毒越来越多地被报道为人类肠道病毒组的成员。一组核心噬菌体由九个代表组成，包括crAssphage。广泛共享的病毒类型（virotypes）可能表示了，个体之间特定宿主的广泛共享，或者这些病毒在人类微生物组中具有广泛的宿主范围。

我们使用HMMs来注释病毒重叠群，这使得我们能够深入探讨病毒组的分类学的内容。除了与公共数据库的比较所揭示的内容之外，我们并确认了这些注释。因为每种类型的病毒(例如，科水平)都需要自己的HMM，所以我们将这种方法应用于几个关键组（key groups）。当应用于crAssphage时，HMM检索的重叠群仅与来自粪便病毒的序列聚在一组，而不与来自其他环境(例如陆地或海洋)的序列。这表明，尽管crAssphage是一个多样化的噬菌体群体，但其在人类肠道中的多样性仅限于crAssphage参考基因组，IAS病毒参考基因组或者Chlamydia噬菌体相关的序列。我们还将HHM应用于Microviridae科，这是一种ssDNA噬菌体。我们能够确认Gokushovirinae和Alpavirinae亚科的不同成员的存在。尽管有证据表明，所描述的Alpavirinae基因组构成了Microviridae科的第三组，它们对应于溶源噬菌体，这使得很难将它们整合到国际病毒分类委员会(International Committee on Taxonomy of Viruses， ICTV )的分类中；因此，在应用HMM图谱（profiles）前，没有任何重叠群被标注为Alpavirinae。

对于每一个病毒分类组（taxonomic group），都有一组相应的细菌宿主。从我们用来选择双胞胎的16S rRNA基因多样性数据中，很明显是哪一个门水平的细菌对具有一致性的微生物组双胞胎和不具有一致性的微生物组双胞胎的微生物组的差异贡献最大。但是与细菌不同的是，我们无法根据病毒中的目或者科来辨别出这样清晰的模式。事实上，大部分噬菌体多样性只归入一个目Caudovirales，及其三个科：Myoviridae, Podoviridae和Siphoviridae。这些科的病毒可以感染不相关的宿主。因此，我们没有必要期望特定的病毒目或病毒科显示在细菌门水平上观察到的这种模式（很明显是哪一个门水平的细菌对具有一致性的微生物组双胞胎和不具有一致性的微生物组双胞胎的微生物组的差异贡献最大）。

最后，我们注意到了一个有趣的模式，即在选定的细菌种类中，病毒组中完全覆盖细菌基因组。由于这些推定的污染物不是微生物组中最丰富的成员，它们不太可能代表大量（bulk）DNA的随机污染。目前还不清楚为什么某些细菌基因组显示如此高的覆盖率。一种可能性是，我们正在观察具有转座功能的噬菌体宿主的物种范围。噬菌体，比如Mu噬菌体可以随机整合到宿主基因组中，通过连续几轮复制转座扩增，然后可以包装宿主基因组的任何部分。有趣的是，在此检测到的几种污染也被报道为其他人类肠道病毒研究中的污染，这可能表明具有转座功能的噬菌体的宿主特异性。另一种解释包括囊泡产生、基因转移剂和/或广义转导过程。对不同病毒数据库中发现的完整细菌基因组的进一步比较，可能有助于揭示它们的来源，特别是在多项研究中回收的相同细菌物种。

研究方法

EXPERIMENTAL MODEL AND SUBJECT DETAILS

粪便样品

Fecal Samples

本研究中使用的粪便样本是作为之前的研究的一部分而获得的。从16S rRNA基因多样性中，先前共测量了354对同卵双胞胎，他们的粪便样本是在2013年1月28日至2014年7月14日期间收集的。我们基于双胞胎中的三个微生物组β-多样性距离，选择了9对微生物组具有一致性和12对微生物组不具有一致性的同卵双胞胎。这三种距离是，unweighted UniFrac, weighted UniFrac和Bray-Curtis。在微生物组具有一致性和不具有一致性组内都选择了能够在性别、年龄、BMI和BMI差异方面（BMI difference）保持平衡的双胞胎。一致组中的双胞胎年龄在23至77岁之间，包括5对男性和4对女性，微生物组不具有一致性的双胞胎年龄在29至81岁之间，包括5对男性和7对女性。所有涉及使用这些先前收集的样本的工作都得到康奈尔大学IRB的批准。

从人类粪便样品中分离类病毒样颗粒

Isolation of Virus-like Particles (VLPs) from Human Fecal Samples

类病毒颗粒（Virus-like Particles, VLPs）的分离基于之前描述过的方法。对于类病毒样颗粒分离，之前使用0.02 μm滤膜( Whatman )过滤，0.5 g粪便样品在15 mL PBS中涡旋5 - 10分钟，使之重新悬浮。将匀浆在4,500 xg下离心30分钟，上清液通过0.22 μm polyethersulfone（PES ) Express Plus Millipore Stericup ( 150 ml )过滤，以去除细胞碎片和细菌大小的颗粒。然后滤液在Millipore Amicon Ultra-15离心过滤单元上浓缩100K至1毫升。将浓缩物转移到5 Prime Phase Lock Gel中，并在室温下用200 ml氯仿（chloroform）温育10分钟。在15000 xg下离心1分钟后，将水层转移到新的微量离心管中，并在37℃用Invitrogen TURBO DNase (14 U )、Promega RNase One (20 U )和1 ml Benzonase Nuclease (E1014 Sigma Benzonase Nuclease) 处理3小时。温育后，向每个样品中加入0.04 体积（volumes）0.5 M EDTA。在进一步处理之前样品被储存在- 80℃环境下。

病毒DNA鸟枪法测序

Viral DNA Shotgun Sequencing

病毒DNA是用从Invitrogen购买的PureLink Viral RNA/DNA Mini Kit提取的。然后使用从Sigma-Aldrich购买的GenomePlex Complete Whole Genome Amplification (WGA2) Kit扩增每个病毒DNA样品。该步骤包括两个空白对照，但是产量非常低，无法构建文库。然后用Covaris S2 Adaptive Focused Acoustic Disruptor 对扩增产物进行片段化，参数设置如下：duty cycle为10 %，cycle per burst为200，强度为4, 时长60 s。每个病毒测序文库按照Illumina TruSeq DNA Preparation Protocol制备，每个样品有一个独特的标签。所有标签化文库汇集在一起。其中一半由BluePippin根据大小选择用来富集更长的插入物( 425 - 875 bp，包括adapters)。在Illumina HiSeq 2500仪器上，两个文库（pools），即“大片段文库large-insert-size library”和“小片段文库short-insert-size library”，在康奈尔生物技术资源 测序中心的Illumina HiSeq 2500上，以250 bp双端的快速模式运行，在独立lane上进行测序。

全部粪便样品的宏基因组鸟枪法测序

Whole Fecal Metagenome Shotgun Sequencing

基因组DNA是使用PowerSoil - htp DNA提取试剂盒从每份样品中的大约100毫克样品中提取的。然后，每个测序文库都是按照Illumina TruSeq DNA标准建库方法，用500 ng DNA，使用gel-free方法，经过14次PCR循环制备的。在康奈尔生物技术资源 测序中心的Illumina HiSeq 2500上，以2x150 bp模式运行，在独立lane上进行测序。

细菌污染的评估

Assessment of Bacterial Contamination

2017年2月21日，从NCBI FTP上检索了8163个完整的细菌基因组。使用Bowtie 2 v.2.2.8比对（mapped）每个样本的序列到细菌参考基因组上，参数如下: -local-maxins 800 -k = 3。使用view和depth Samtools commands v.1.5，计算每个碱基的基因组覆盖率时，只考虑了比对质量超过20的序列。接下来，平均基因组覆盖率为100K bp bins。我们观察到均匀覆盖的基因组的median bin coverage至少是100；那些median bin coverage大于100的基因组被认为是污染，删除比对到这些基因组的序列。细菌基因组中可能有一个或多个溶源噬菌体；这些溶源噬菌体的爆发（bursting）事件可能会发生，产生几个类病毒样颗粒（Virus-like Particles, VLPs）。作为一项保守措施，为了避免源自溶源噬菌体而非细菌基因组本身的序列丢失，bins覆盖度超过细菌平均覆盖度三个标准差的bins也被识别并归类为类溶源噬菌体（prophages-like）区域。比对）到潜在污染基因组的read被标记为“污染”，并从进一步的分析中移除，而比对（mapping）到高覆盖度 bins的序列比对被标记为“可能的溶源噬菌体”。

使用内部Python脚本构建每个样品中每个潜在污染的丰度的矩阵，并采用PKM标准化方法。与此同时，来自古德里奇等人的数据，每个OTU的相对丰度被回收，并在物种层面使用summarize_taxa.py qiime脚本进行汇总。针对两组物种，计算了污染物相对丰度与其对应的16S rRNAs数据之间的Spearman相关性。

功能组成

Functional Profiles

经过joined和trimmed的来自“小片段文库”的序列被比对到整合基因集上，这是人类肠道微生物组中参考基因的总集。使用的方法是BLASTX DIAMONAL v.0.7.5，最大e-value值的截止值为0.001，要报告的最大目标序列数目设置为25。

比对到整合基因集后，使用内部Python脚本生成了丰度矩阵。然后根据整合基因集)提供的每一个基因，使用KEGG对矩阵进行注释。注释的丰度矩阵被稀释(rarefied，subsampling without replacement)为每个样本2000000 read命中（hits）。然后使用QIIME 1.9的命令summarize_taxa_through_plots.py生成KEGG功能图谱。使用R包Psych计算每组(附加的微生物组、附加的病毒组、具有一致性的微生物组样品病毒组和不具有一致性的微生物组样品病毒组)功能图谱的组内相关系数。

从头组装

De novo Assembly

经过质量剪切步骤的序列还是双端序列（正向和反向），使用针对于短序列的整合宏基因组装 (InteMAP)流程组装，插入（insert）片段大小设置为325 bp ± 100 bp。每个样品单独组装。组装第一次运行后，所有干净的序列都使用Bowtie 2 v.2.2.8比对到组装的重叠群，参数如下: -local-maxins 800。然后至少一次比对具有一致性的序列都提交到InteMAP以供二次组装。使用自写Perl脚本，将来自所有样本的大于500 bp的重叠群汇集在一起，并用两两比对方法进行比较。从这一分析中，有可能识别出潜在的环状基因组，并在另一个重叠群中将超过90%长度的重叠群去冗余。

根据Roux等人对覆盖度（coverage）和长度（length）的建议，对去重复后的宏基因组集合（assemblies）进行序列的补充，用于构建丰度矩阵。使用Rsubread v.1.28.0将序列比对到非冗余的重叠群。使用自写Python脚本将比对（mapping）的输出解析成丰度矩阵，该矩阵通过每个样品中相对丰度RPKM方法进行标准化，并转换为Log10(x+1)，x为标准化丰度。标准化覆盖率低于5x的重叠群被除去。最后，应用重叠群长度过滤获得病毒类型（virotypes）。绘制重叠群数量作为长度的函数的衰减曲线，长度阈值就是衰减曲线的转折点1,300 bp。

HMM注释

HMM Annotation

建立了独立的HMM组成(profiles，图谱/打分矩阵)，以识别类crAss（crAss-like）重叠群和Microviridae 重叠群。为了构建HMM 类- crAss组成特征，我们从NCBI网站上下载了类-crAss家族的主要衣壳蛋白( Major Capsid Protein，MCP )序列。多序列比对使用MUSCLE v.3.8.31进行的，并使用UGENE v.1.31.0进行检查。间隙（gaps）超过30 %的位置（positions）被移除。最后，HMM – crAsslike组成是使用软件包HMMER v.3.1b2 的hmmbuild构建的。对于Microviridae科的例子，还使用了Alves等人开发的病毒蛋白1 (viral protein 1，VP1 )的所有HMM组成。

使用hmmsearch查询组装好的重叠群的预测蛋白质，以匹配HMM -图谱。e-value值低于1x10^-5的蛋白质被认为是真正的同源物，其他的标准还包括，参考蛋白质的大小在参考蛋白质的大小范围内( crAsslike MCP : 450 - 510残基（residues）；Microviridae: 450 - 800个残基（residues）)，至少50 %的覆盖率（coverage）和至少40 %的同一性（identity）。覆盖率和同一性通过BLASTp，用真实同源物查询参考序列来确定。

通过上述步骤的真正同源物被用于系统发育分析。参考序列和同源序列使用MUSCLE v.3.8.31进行比对，使用UGENE v.1.31.0去除至少30 %缺口的位点。使用RAxML v.8.2.4进行最大似然( ML )系统发育分析，使用prottest v.3.4.2获得了最佳进化模型，ML树中节点的支持是通过具有100次伪复制（pseudoreplicates）的自举（bootstrap）获得的。

物种组成

Taxonomic Profile

为了推断组装的类病毒颗粒（Virus-like Particles, VLPs）的分类学关系（taxonomic affiliation），使用GeneMarkS v.4.32对所有大于500 bp的重叠群预测基因。预测基因的氨基酸序列随后被用于BLASTp搜索中，使用DIAMOND v.0.7.5对NR NCBI病毒数据库进行搜索最大e-value阈值（cutoff）为0.001，报告中目标序列输出最大数目设置为25。使用BLASTp结果，每个基因的分类由MEtaGenome Analyzer（MEGAN5 ）v.5.11.3中的最近共同祖先算法(lowest-common-ancestor algorithm)进行分配。参数如下: Min Support: 1, Min Score: 40.0, Max Expected: 0.01, Top Percent: 10.0, Min-Complexity filter: 0.44。独立地，使用CENTRIFUGE v.1.0.4查询NT NCBI病毒基因组数据库获得每个重叠群的分类结果。然后使用投票系分配每个重叠群的最终分类注释，其中每个蛋白质的分类注释和重叠群的CENTRIFUGE注释被视为投票。通过对重叠群的所有可能投票，构建了一个N元树（N-ary tree），每个节点（node）的权重（weight）是包括该节点的投票数。重叠群的分类注释将是遍历（traversing）通过权重最高节点（heaviest nodes）的树的结果，但有一个考虑：如果节点的所有子节点具有相同的权重，遍历（traversing）必须停止。分类概况被认为是补充矩阵的一个子集，包含所有由投票系统（voting system）注释或通过HMM概况标注的重叠群 (见上文)。

噬菌体-宿主相互作用的预测

Prediction of Phage-Host Interaction

使用PilerCR v.1.06，从用于评估细菌污染的8,163种细菌基因组中识别出有规则间隔的短回文重复序列（Clustered Regularly Interspaced Short Palindromic Repeats, CRISPRs)。间隔物（spacers）的预期大小为20 bp和72 bp，被用作对病毒类型和分类标注重叠群的查询，使用BLASTn (v.2.6.0+ )短查询参数: 覆盖至少90%间隔区且e值<0.001的匹配被认为是CRISPR间隔区-病毒关联。此外，使用LAST-959病毒类型和有分类注释的重叠群被比对到MVP数据库中病毒群的代表基因组上。由于MVP中的病毒簇包含至少95%同一性的序列，其长度的至少80%，因此只保留满足这些限制的匹配。重叠群的宿主是由其匹配的病毒簇确定的。

多样性指数

Diversity Indexes

样本内的香农多样性指数( a -多样性)和同卵双胞胎内的Hellinger距离( B-多样性)是使用生成的所有三个丰度矩阵(函数、分类和读取补充矩阵)的Vegan R包的多样性和距离vegdist函数计算的。使用皮尔逊相关系数测量病毒组α-多样性和微生物组α -多样性之间的相关性。病毒组β-多样性和微生物组β-多样性之间的相关性是通过使用皮尔逊相关系数的曼特尔检验（Mantel test）来计算的。此外，对微生物组具有一致性的同卵双胞胎之间的β-多样性与微生物组不具有一致性的同卵双胞胎之间的β-多样性进行了比较；p值使用Mann-Whitney U test来计算。

统计分析

QUANTIFICATION AND STATISTICAL ANALYSIS

p值是使用使用Python “scipy”包的Mann-Whitney U测试或Mantel测试获得的；使用Python “scipy” 包将相关系数测量皮尔逊相关系数；α和β多样性是用R的 “vegan” 包计算的；使用R包“psych”计算类别内系数（Intra-class coefficien）；使用RAxML进行最大似然系统发育分析。

拓展阅读

1. 数据分析过程中使用的脚本 https://github.com/leylabmpi/TwinsUK_virome
2. Human Genetics Shape the Gut Microbiome https://www.cell.com/fulltext/S0092-8674(14)01241-0
3. Environment dominates over host genetics in shaping human gut microbiota https://www.nature.com/articles/nature25973?platform=hootsuite
4. Early life dynamics of the human gut virome and bacterial microbiome in infants https://www.nature.com/articles/nm.3950
5. The Human Gut Phage Community and Its Implications for Health and Disease https://www.mdpi.com/1999-4915/9/6/141
6. Healthy human gut phageome https://www.pnas.org/content/113/37/10400

猜你喜欢

• 10000+: 菌群分析
  宝宝与猫狗 提DNA发Nature 实验分析谁对结果影响大  Cell微生物专刊 肠道指挥大脑
• 系列教程：微生物组入门 Biostar 微生物组  宏基因组
• 专业技能：生信宝典 学术图表 高分文章 不可或缺的人
• 一文读懂：宏基因组 寄生虫益处 进化树
• 必备技能：提问 搜索  Endnote
• 文献阅读 热心肠 SemanticScholar Geenmedical
• 扩增子分析：图表解读 分析流程 统计绘图
• 16S功能预测   PICRUSt  FAPROTAX  Bugbase Tax4Fun
• 在线工具：16S预测培养基 生信绘图
• 科研经验：云笔记  云协作 公众号
• 编程模板: Shell  R Perl
• 生物科普: 肠道细菌 人体上的生命 生命大跃进  细胞暗战 人体奥秘  

写在后面

为鼓励读者交流、快速解决科研困难，我们建立了“宏基因组”专业讨论群，目前己有国内外5000+ 一线科研人员加入。参与讨论，获得专业解答，欢迎分享此文至朋友圈，并扫码加主编好友带你入群，务必备注“姓名-单位-研究方向-职称/年级”。技术问题寻求帮助，首先阅读《如何优雅的提问》学习解决问题思路，仍末解决群内讨论，问题不私聊，帮助同行。

学习扩增子、宏基因组科研思路和分析实战，关注“宏基因组”

点击阅读原文，跳转最新文章目录阅读
https://mp.weixin.qq.com/s/5jQspEvH5_4Xmart22gjMA