全球微生物组超突变的生态学和分子靶点

原名：Ecology and molecular targets of hypermutation in the global microbiome

译名：全球微生物组超突变的生态学和分子靶点

期刊：Nature Communications

IF：14.919

发表时间：2021.5.24

通讯作者：Simon Roux, Emiley A. Eloe-Fadrosh

通讯作者单位：美国能源部联合基因组研究所劳伦斯伯克利国家实验室

DOI号：10.1038/s41467-021-23402-7

实验设计

结果

1.大规模宏基因组数据挖掘揭示DGRs的广泛多样性

为了识别候选DGRs，研究者搜索了异常重复序列1kb范围内的逆转录酶（RT）基因，并根据系统发育分析中所处的位置和错配模型来过滤假阳性。利用这种方法，他们从81,404个公共基因组数据和来源于163个环境类型的9467个公共宏基因组数据获得了全球范围的DGRs多样性数据（补充数据1）。在分析的基因组中，他们一共检测到1314个DGRs，数目和多样性跟以往数据挖掘所发现的相差无几。在分析的宏基因组数据中，他们一共检测到31,007个DGRs，比此前发现的总DGRs多了近15倍。总体而言，DGRs分布在1500多个细菌和真菌属和90多种环境类型中（补充数据2，补充说明1）。值得注意的是，由于宏基因组数据通常都是片段化的（编码DGRs的contig的平均长度为9584bp），他们的方法无法检测出DGRs和靶标基因距离较远的情况（例如可变重复区VR不在RT附近的DGRs），因此这一数据依然少于实际DGRs情况。

然后，研究者根据平均氨基酸相似度（AAI）对发现的RT序列进行聚类，先基于95% AAI可以分为13,415个OTUs，再基于50% AAI可以分为1318个clusters（附图1，补充数据3）。每个OTUs和clusters的成员都与某一基因组类型（例如，病毒基因组或细胞型基因组），物种类型和生物群类型构成相关关系，表明这一分组代表了连续且主要是垂直传播的DGRs进化单元（补充说明2，附图2）。为了获得DGRs多样性的全面概览，研究者构建了包含每一RT cluster的代表序列以及二组内含子（Group Ⅱ introns）、反转座子和未被表征的RT等其他RT基因的系统发生树（图1b）。结果表明，DGRs构成了一个跟其他类型RT区分开来的单系进化枝（monophyletic clade），这表明这些元件具有单独的进化起源，而这一发现跟以前的猜想相吻合。总体而言，75%的clusters只含有宏基因组来源的DGRs序列。此次研究中DGRs的系统发生学多样性是此前已知DGRs的近6倍（573%），这一发现也强调了宏基因组和宏转录组数据对于发现DGRs序列的突出作用。

图1.DGRs多样性在生物体、生物群和生物物种中的分布情况。a.DGRs突变逆转录回巢错过程模式图。不同颜色代表DGRs的主要组成部分；图注指出该过程的三个主要阶段。b.利用DGRs和非DGRs的逆转录酶（RT）构建的系统发育树。先对所有RT蛋白序列进行聚类，每一簇中选取代表序列进行系统发育分析。分支颜色代表每一簇中RT的类型。所有支持度<50%的节点均收缩起来。外层的圆圈由内到外分别代表基因型类型、物种分类和生物群落类型。缩写解释：CPR: Candidate Phyla Radiation. DPANN Diapherotrites, Parvarchaeota, Aenigmarchaeota, Nanoarchaeota, Nanohaloarchaeota, FCB Flavobacteria, Fibrobacteres, Chlorobi, Bacteroides, PVC Planctomycetes, Verrucomicrobia, Chlamydiae, Aq aquatic, Te Terrestrial, En Engineered, H-a Host-associated。NA代表特征无法估计的情况。c.不同DGRs clade在RT OUT水平的一些特征分布。颜色代表内容与图1b相同。不含NA值。d.不同物种编码DGRs的基因组的丰度。根据单拷贝标记基因（见“方法”）计算出整个宏基因组组装结果中观察到的基因组总数，并从整个数据集中得出DGR的平均频率。然后为每个组计算每个基因组的DGR检测频率，并与总体频率相比较，得出对数富集率。除了蓝细菌组（P值=0.21），图中所有的对数比都具有统计学意义（卡方检验校正后的P值<1E-10）。e.不同生物群落编码DGRs的基因组的丰度。对于每个生物群落，在估计的基因组总数和每个宏基因组中检测到的DGRs的数量之间计算出线性回归（见附图5）。然后，回归斜率被认为是对每个基因组中DGR是、平均数量的估计并此显示出来，误差条代表斜率估计的标准误差。每基因组0.05和0.25个DGRs的截止点用垂直虚线突出显示。在这些计算中，病毒和低复杂性的宏基因组被排除在外（见方法和补充数据1）。圆点根据生物群落类型着色（蓝色：水生，棕色：工程，粉色：与宿主相关，绿色：陆生）。

  2.DGRs的传播受到严重限制并体现在整体的谱系划分

将DGRs来源的物种类型（病毒型或细胞型）、物种分类和生物群信息比对到系统发生树上，发现所有DGRs主要分为6个clade（clade1-6，图1b，c，补充说明3）。其中，三个clade基本都是病毒编码的DGRs，而且大多是被预测可以感染Firmicutes、Bacteroidetes、Proteobacteria等肠道细菌的噬菌体（clade 1，4，6，图1c，附图3）。Clade 2和clade 5 几乎全由细胞型生物编码的DGRs组成，主要是水生生物，其中clade2主要分布在Patescibacteria（也被称为CPR），而clade 5主要分布在Patescibacteria和Proteobacteria等多种生物中。Clade 3的DGRs在病毒型和细胞型生物中均匀分布，主要包括Bacteroidetes和Proteobacteria。

虽然在这里腺嘌呤位置的突变偏移没有作为检测DGRs的指标，模板重复区（TR，template repeat）和可变重复区（VR, variable repeat）之间的配对在TR区腺嘌呤位点处的错配率大于75%（图1c）。在人工检查离群值之后（补充数据4），发现实际上仅7个cluster的DGRs的腺嘌呤位点错配率低于75%（附图4，补充数据4），这与之前的比较基因组研究和生化研究结果一致，证实了腺嘌呤的突变倾向是DGRs RT的固有特征。结合系统发生树中DGRs的单系群特征和腺嘌呤突变的普遍性可知，目前已知的DGRs都起源于同一个祖先，并在之后在其他生物之间散发转移，并导致了现在看到了6个clade。这一进化过程与DGRs最开始由可移动基因元件所编码并被传播至各个物种和生态系统的猜想相吻合。

在分析的9467个宏基因组中，有几种物种和生物群含有明显多的DGRs。首先，DGRs在CPR, Firmicutes和Flavobacteria-Bacteroidetes-Chlorobi（FCB）的成员中明显更常见（图1d）。然后，研究发现人肠道、咸水湖、垃圾堆、地下水库等环境样本中的基因组也携带更多DGRs（图1e，附图5）。系统发育逻辑回归分析进一步证实种系发生和生态因素共同驱动DGRs的分布（补充说明5）。DGRs分布和特定病毒型生物和细胞型生物的单系进化枝具有相关关系。在解释完这一系统发育信息后，病毒编码的DGRs依然跟特定的生物群具有明显正相关，而细胞型生物编码的DGRs仅与一种生物群成负相关（附表1）。综上所述，这些结果表明DGRs具有漫长复杂的进化历史，并可以在不具有系统发育相关性的物种中进行传播，但只能在靶点突变是选择优势的特定位点或支持DGRs突变逆转录归巢 (mutagenic retrohoming) 的物种才能保留下来。而DGRs在基因组和生物群中分布广泛且不均匀的特点支持了这一猜想。

3.DGRs的靶标基因具有相同的组织结构特征

为了解识别的DGRs可能的功能，研究者接下来分析了这些DGRs的36,611个靶点基因的多样性。跟之前报道的一样，大部分（68%）的靶点基因单独比对到参考数据库时无法进行功能注释。但是，从头聚类分析（de novo clustering）结果表明，大部分（92%）DGRs靶标基因编码蛋白可以分为24个蛋白质簇（PCs，protein clusters）并可见明显的基因组类型和DGRs 进化枝划分（图2a，附图6，补充数据5，补充说明6）。

PCs的功能结构域和位点保守性分析表明DGRs的靶标基因具有近乎普遍的模块化组织（modular organization）。靶点蛋白都是明显的多结构域蛋白，VR区通常位于C末端（图2a），这通常是由于逆转录归巢需要顺式作用DNA元件，只有作用于C端的VR区时这些顺式作用元件才会位于基因间并不受氨基酸编码的限制。虽然也含有一些DGRs可变的位点，但这些C末端区域整体来说比PCs的其他区域更保守，这可能是由于跟DGRs诱导超突变有关的结构限制（附图7）。相应地，虽然在靶点蛋白N末端预测出一些折叠结构，但是含VR区跟C型凝集素折叠结构相关。少数VR区以往被认为跟Ig样折叠有关，但是通过扩展序列的范围对这些序列重新进行分析，发现这些靶标对应于含Ig样结构域的噬菌体尾部纤维，旁边是一个未被表征的非Ig样VR结构区（附图8和9，补充说明7）。由于新的C型凝集素折叠结构域突变体的发现仍然具有一定的误差，因此其他一些跟VR重叠的未知保守结构域也可能是属于C型凝集素折叠结构域（图2a）。基于DGRs靶标PCs的分布，这些C型凝集素折叠结构域最有可能跟新型病毒、未培养细菌（CPR）和真菌相关（附图9）。DGRs靶标蛋白中的模块性组织也表明DGRs靶标区域可能发生基因内的重组，这种重组可能将大范围独立折叠的结构域跟C端或C型凝集素编码区融合从而产生易于突变的嵌合靶标。

  图2. DGR靶标基因的多样性和主要类型。 a.丰度最高的DGRs靶点PCs的流行程度和序列特征。列出的24个PCs可以代表92%的DGRs靶点。这些靶点蛋白分成四种类型（左边图；用红色星星高亮出目标蛋白的预测未知），右侧图列出了它们相关的DGRs RT OTUs的数目、基因组类型的相对百分比、DGR clade分布、在C端可检测到VR区凝集素折叠结构域的PCs占比，可能具有尾部结构蛋白或膜蛋白的PCs占比以及靶点序列中VR区的相对位置等特征。箱线图下线和上线分别对应于第一和第三四分位数，箱须延伸不超过四分位数范围的±1.5倍。对于C-Lec折叠VR和定位预测数据，仅考虑高质量目标（见“方法”）。对于混合PC（PC_00009和PC_00021），来源不明的目标被排除在外。b.匹配至少50个CPRISPR间隔序列的基因组编码或没有DGRs的估计病毒宿主多样性。纵轴数据代表每50间隔序列相关联的宿主物种数目。箱线图下线和上线分别对应于第一和第三四分位数，箱须延伸不超过四分位数范围的±1.5倍。c.含有一个跨膜结构域（上图）或VR区外含有一个以上细胞粘附/糖类结合相关功能域（下图）的细胞型生物靶点的百分比。目标序列根据其PC内容分为“膜结合”PC或“其他”PC。     4. DGRs的靶标基因主要参与病毒和细胞以及细胞和颗粒的相互作用

由于DGRs靶标蛋白具有近乎普遍的模块化组织结构，因此可以基于目标蛋白C末端VR区外的保守结构域或序列特征推测其可能的功能。在24个主要靶标PCs中，即使使用高灵敏度的注释工具（如HHblits）进行分析，大部分（71%）C末端VR区外的序列跟已知的蛋白质结构域也不存在有意义的序列相似性。因此，我们利用不基于相似度的蛋白质注释方法将这些PCs按照功能大类进行划分，也就是把病毒DGRs的PCs分为“结构蛋白”类和“未知”类，把细胞型生物DGRs的PCs分为“膜结合”类和“未知”类（补充说明8）。

病毒DGRs（clade1，4和6）的靶标PCs大部分被注释为尾部结构蛋白，参与宿主识别和粘附过程（图2a），这其中包括涉及先前报道的DGR介导鲍特氏菌属（Bordetella）噬菌体尾纤维突变的目标蛋白，而这种突变可使其改变宿主。我们推断宿主粘附蛋白发生的超突变可能广泛地使编码DGRs的病毒能够接触到更多种类的宿主细胞。于是，研究者通过将编码DGRs和不编码DGRs的病毒跟包含来自576,561个原核生物参考基因组的670万条CRISPR 间隔序列（spacer）的数据库进行比对来证实此猜想。结果表明，控制某些混杂因子后，编码DGR的病毒比不编码DGRs的病毒与更多种宿主存在相关关系（图2b，附图10，补充说明9）。尽管这些结果表明DGRs介导的超突变可以使噬菌体粘附并感染更多种宿主细胞，但这并不一定意味着这些噬菌体就能成功有效地进行复制周期。宿主体内的某些特征和防御机制可以使这些感染无效，而这些病毒要成功感染新的宿主（即使跟原宿主亲缘关系较近），可能需要一些额外的适应条件。

对细胞型生物目标蛋白来说，大多数PCs都包含至少一个N端跨膜结构域或信号肽以及参与蛋白质结合、糖类结合和细胞粘附的功能结构域（图2a，c，补充数据5），表明这些靶点蛋白大多是结合微生物胞体和病毒颗粒等细胞外成分的膜锚定蛋白。而且，最常见靶标分子 (PC_00001) 相关的宏基因组组装的基因组（MAGs）基因和功能注释呈现出富营养化和颗粒相关的生活方式（补充说明10），而作用于这一目标分子的DGRs主要分布在水环境，但在各水环境之间分布并不均匀（图1e，附图5，附图9）。因此，颗粒结合、细胞与细胞粘附或细胞与表面粘附所提供的选择优势在不同环境可能不同。例如，在大海中，生物的随机结合可能并不有利，因为它可能会导致粒子沉降而使细胞丢失增加，这可以解释为什么在这些样品中很少检测到DGRs（图1）。然而，重要的是，细胞膜蛋白也可以成为病毒进入细胞的受体，因此，作为病毒宿主军备竞赛（arms race）的一部分，其中一些目标蛋白可能受到DGRs驱动的多样化影响。Clade 2中CPR和clade 5真菌DGRs作用的目标分子跟其他有所不同，它们不含有典型可供识别的跨膜或其他保守结构域（图2c），这与之前报道的结果相同。而这是由于功能结构域在这些差异大的基因组中难以被发现，还是这些蛋白的确是非膜结合DGRs的，还需要进一步研究。综上所述，这项研究发现的大量DGRs靶点证明，DGRs主要跟细胞-颗粒、细胞-细胞和细胞-病毒之间的反应相关，在少数情况下也可能参与微生物的细胞调节过程。

  5.DGRs在生物中广泛活跃

接下来研究者评估了不同分类单元和生态环境DGRs位点的种群多样性。为此，他们分析了长度大于10 kb，覆盖度大于20×的基因组的6901个DGRs的单核苷酸突变和氨基酸突变。总体而言，70.1%的VR区含有单核苷酸突变（SNVs）（图3a，附图11，补充说明11），且主要是非同义SNVs。研究者首先通过非同义多态性和同义多态性的比率的比值（pN/pS）来评估这一点：几乎所有非靶点基因pN/pS都小于1（与纯化选择相吻合），含有多于1个SNVs的VR位点中，超80%位点的pN/pS大于1，意味着非同义突变的累积（图3b）。由于在没有同义SNV的情况下无法计算pN/pS以及36%的VR位点只表现出非同义SNV，研究者选择直接估计非同义SNVs中VR位点的富集程度并以此作为新近DGRs活性的标志（见方法和补充说明11）。他们推测当VR位点具有明显高于邻近基因的非同义SNV密度时，这是由近期或正在进行的DGRs驱动的突变事件造成的。 

对于所有DGRs分类，50-75%的DGRs呈现出近期活性的标志（图3c）。病毒编码的DGRs具有最高的活性，而clade 2和clade 3的DGRs 的活性明显低于平均值（补充说明12）。然而，在这些单样本变异分析中，强纯化选择会会降低种群多样性并可能会掩盖DGRs活性，即，活性DGRs可能产生新的变异，这些变异几乎会在瞬间从种群中清除，从而无法被检测出来。

图3. DGR靶基因座的多样性模式。a.DGRs RT 10k范围内非靶点基因（左）和VR位点（右）中，SNV不小于1的基因的百分比。“low coverage”类包括VR区的覆盖度显著低于附近基因覆盖度的情况，此时只利用到部分reads，无法可靠推测VR区的种群多样性。b.非靶点基因和VR位点的同义SNV不小于1的基因pN/pS值分布。虚线代表pN/pS=1的位置。箱线图下线和上线分别对应于第一和第三四分位数，箱须延伸不超过四分位数范围的±1.5倍。使用Kruskal–Wallis检验比较非靶基因和VR位点的pN/pS分布。c. 在不同 DGR 分类中，基于非同义 SNV 中 VR 基因座相对于周围基因的富集，“有活性”与“无活性” DGRs 的比例。活性序列占比明显低的组用星形进行标注。图3.d.VR基因座受到的两种竞争力（纯化选择和DGRs多样化）的示意图。DGRs可能活性水平的3种示例用不同的颜色表示，随时间发展的结局在右侧。e. 在两个时间序列数据集中观察到的一个 DGRs 靶点的多样性和变化示例。热图上标注了每一个位置相关的氨基酸，其在该群落中的频率用不同颜色标注。右边面板表示位置的类别，颜色设置与图d相同，基于样本内熵、样本间余弦距离和时间序列中氨基酸变化的数量（见补充说明12）。上图表示每个样本中所有位置的覆盖率中位数。为参考起见，在预测的VR之外的同一蛋白质的十个随机位置也包括在内。f.不同生物群落和clade中有活性的VR位点的分布。变异位置<50%和氨基酸变化<5%的病例被认为是 "低DGR活性"，用白色标示。只包括有≥10个DGRs的组别。

6.DGRs活性驱动靶点的频繁改变

鉴于DGRs的广泛活跃性，多数VR位点的进化可能处于DGRs多样化和纯化选择这两种相反的压力。只要不同突变体的适应值（the adaptive value）随时间波动，就可以从时间序列数据中分析这些进化压力的关系（图3d，附图13，补充说明12）。具体来说，研究者推测DGRs多样化造成了每个样本的高群体多样性和时间节点间优势等位基因的改变，但纯化选择会降低每个样本的群体多样性。同时，在多样化和净化选择的交替阶段，研究者期望观察到每个样本中的种群多样性低，但样本间的优势等位基因会发生变化（附图14）。

为了验证这一猜想并阐明自然条件下DGRs多样化作用和纯化选择的平衡状态，研究者对某一时间序列的宏基因组中发现的一部分DGRs进行了分析。他们找到适合分析的包含563个DGRs的130个纵向数据集（基因长度不少于10kb，覆盖度不低于10x，并在2个及2个以上时间节点可以被检测到。补充数据6）。总体而言，在这些DGRs中，大多数预测的VR区多样性程度高，氨基酸替换频繁，表明DGR的高活性强于纯化选择（图3e,f）。自然生长的生物群都可以观察到这种模式，而在培养的微生物（如实验室孵育）中无法观察到，这种情况下各类DGRs多样性和氨基酸改变频率都要明显降低，可能跟实验过程中的种群瓶颈效应（population bottleneck）、自然选择压力减少或短时间跨度有关（图3f）。总的来说，主要以结构蛋白为靶点的病毒编码DGRs（即clade1、4和6），比细胞型生物编码的DGRs（clade2和5）和“混合”DGRs（clade3）具有更高氨基酸突变率（图3f，附图15）。从这里观察到的平均突变率推断，即使DGRs仅针对约0.1%的氨基酸残基产生作用，保守估计DGRs驱动的突变平均造成了病毒基因组中6-16%的氨基酸突变（补充说明13）。除此之外，研究者还观察到，与人类微生物组项目中“未受干扰”的肠道微生物组样品相比，12个月减肥计划后的个人肠道样品中DGRs活性水平更高，其中细胞编码的clade3 DGRs活性增加更明显（图3f，附图16）。综上所述，DGRs驱动病毒结构蛋白比驱动病毒结构蛋白和细胞型靶点发生随时间更稳定的突变，后者可能跟选择压力下的适应更相关，而这是由DGRs的RT区活性控制更强还是选择压力更强导致的，仍有待确定。

讨论

此宏基因组DGRs比较分析表明，DGRs是微生物和病毒基因组进化过程的重要成分。靶点蛋白腺嘌呤突变偏倚和C端凝集素折叠的普遍保守性意味着DGRs的RT区在产生突变的类型和位置上受到机制上的限制，而C端凝集素折叠似乎能够适应大规模的序列变异的。在选定的生物群落和分类群中观察到的明显DGRs富集可能反映了特定的生态条件和生活方式，在这些条件下超突变是有利的。例如，在人类肠道微生物群中，丰富的资源和多种噬菌体的频繁感染相结合可能通过细胞壁修饰产生耐药性。这将反过来选择编码DGRs的病毒，它将利用宿主识别蛋白（如尾部纤维）的超突变来躲避宿主抗性，并在此过程中扩大其潜在宿主的范围。最后，DGRs广泛持久的活性意味着它们可以用于维持目标的高群体多样性而不是只在极端压力条件下被触发。综上所述，此次对DGRs多样性和活性的全球分析结果表明，DGRs可能在多个分类群和生物群落中塑造了长期的微生物-微生物和病毒-宿主相互作用；它们驱动了病毒颗粒和微生物包膜中关键成分的多样性和进化；它们可能是病毒和细胞微生物适应和应对不断变化的环境的重要机制。

打赏