大流行尺度系统基因组学揭示SARS-CoV-2重组格局

重组是许多流行病原体(包括β -冠状病毒)新遗传变异的主要原因⁶包括SARS-CoV-2在内的分支。通过混合来自不同基因组的遗传物质，重组可以产生具有潜在重要表型效应的新的突变组合⁷．例如，重组被认为在中东呼吸综合征的近期进化史中发挥了重要作用⁸严重急性呼吸综合征冠状病毒(SARS-CoV)^{9，10，11，12}．重组也有可能在未来产生具有人畜共患病潜力的病毒¹³．因此，准确和及时地鉴定重组是理解进化生物学和已建立和新出现的病原体在人类、农业和自然种群中的感染潜力的基础。

现在，在SARS-CoV-2人群中存在着大量的遗传多样性¹⁴已知有时会发生与不同的SARS-CoV-2变体合并感染的情况¹⁵重组有望成为大流行期间新基因变异的重要来源。自大流行早期以来，SARS-CoV-2基因组中是否存在可检测到的重组事件信号一直存在激烈的争论¹³．尽管如此，一些明显真正的重组世系已经确定使用特别的方法¹⁶以及通过减少对可能的重组祖先对的搜索空间来处理大量SARS-CoV-2数据集的半自动方法^16，17．由于在持续的SARS-CoV-2大流行期间及时和准确监测病毒遗传变异的重要性，需要检测和表征重组单倍型的新方法，以便在新的变异基因组序列获得后尽快评估它们。这种快速转变对于推动在知情和协调的情况下对新型SARS-CoV-2变种的公共卫生应对至关重要。

我们开发了一种检测大流行尺度系统发生中的重组的新方法，利用系统发生位置进行重组推断(涟漪，图。1)。因为重组违背了许多系统发育方法的中心假设，即单一的进化史在整个基因组中共享，从不同基因组中产生的重组谱系通常会在“长分支”上发现，这是由于适应两个亲本单倍型的不同进化史(图2)。1)。请注意，只要重组相对不常见，即使分支长度人为扩大，系统发育推断也有望保持准确¹⁸．waves通过首先识别综合SARS-CoV-2突变注释树上的长分支来利用这一信号^19，20.．然后，waves将潜在的重组序列详尽地分解为不同的片段，并使用最大限度的简约性将每个片段替换到一个全局系统发育中。ripple报告了两个亲本节点(下文称为供体和受体)，相对于全球系统发育的原始位置，这两个亲本节点导致了最高的简约分数改善1)。因此，我们的方法利用了每个亲本谱系的系统发育信号和基因组上标记的空间相关性。我们使用基于推断的特定位点突变率的零模型来建立显著性(补充文本)2而且3.)。

通过模拟的大量测试表明，涟漪是有效的，敏感的，可以自信地识别重组世系(补充文本4- - - - - -6)。正如预期的²¹，当重组发生在基因组边缘或基因相似序列之间时，使用waves(扩展数据图)更难检测。1而且2)。尽管如此，waves检测模拟重组的灵敏度为75.8%。在作为重组体检测的模拟样本中，waves准确地识别了90%的模拟断点(扩展数据表1及补充文本6)。此外，waves能够检测到在以前的分析中确定的所有高置信重组¹⁶(补充文本6)。利用ripple对全球约160万个SARS-CoV-2基因组进行了重组分析，结果显示，部分测序的SARS-CoV-2基因组属于可检测的重组谱系。为了减轻测序和组装错误的影响，我们排除了只有一个后代的所有节点，我们应用保守滤波器从由waves标记的重组集中去除潜在的虚假样本，并且我们使用原始序列读取数据手动确认假定重组样本子集中的突变(补充文本)7而且8，扩展数据表2和扩展数据图。3.)。在此之后，我们保留了589个独特的重组事件，这些事件总共有43104个后代样本(补充表1)。这意味着总样本SARS-CoV-2基因组中约2.7%被推断属于可检测的重组谱系。事后统计分析得出我们的统计阈值的经验错误发现率估计为11%9及扩展数据表3.)。此外，供体节点和受体节点的后代之间的地理位置和日期元数据的过度相似性支持了重组基因组的许多祖先在人类群体中共同循环的观点(补充文本)10而且11和扩展数据图。4而且5)。因为在基因相似的病毒谱系之间发生的重组事件难以检测(扩展数据图。2)，我们的估计可能大大低估了重组的总体频率。因此，相对于SARS-CoV-2人群的全球重组频率，涟漪效应估计可能是保守的。

waves发现了SARS-CoV-2基因组中重组断点位置的显著不均匀分布，这与之前对β -冠状病毒的分析一致^11，22．特别是，在假定的重组事件中，相对于基于随机断点位置的预期，朝向SARS-CoV-2基因组3'端的重组断点过多(P< 1 × 10⁻⁷；排列测试;补充的文本12)。值得注意的是，当我们模拟遵循均匀分布的重组断点时，这种偏差是不明显的13和扩展数据图。1)。突变点分析确定了穗蛋白区域立即5'的重组断点频率的增加(20875个碱基对;补充的文本14)，当我们将自己限制在具有最大数量后代的假定节点和不同数据源之间时，这种模式是一致的，进一步表明它不是人为的(补充文本15及扩展数据表4)。假定的重组中断点的速率在变更点的3'处比5'区间高约3倍(图2)。2)，这与其他人类冠状病毒基因组的相对重组率相似¹¹．

几行证据表明，重组断点位置的倾斜分布不是宿主间传播动力学水平上的正向选择的结果。首先，这些重组进化支中的许多只存在了相对较短的时间，并且可能已经灭绝了。检测到的重组淋巴结的观察后代最早和最晚的日期之间的平均时间跨度仅为37天。第二，我们推断的重组事件子集发生在关注的变量之间(VOC;谱系B.1.1.7, B.1.351, B.1.617.2和P.1(参考。²³)和其他谱系，VOCs对刺突蛋白突变的贡献平均略低于非VOC谱系(125个VOC/非VOC重组中有60个，P= 0.48，符号检验)。第三，重组分支的大小与其他分支的大小差别不大，这是预期的，如果重组谱系经历了强选择(P= 0.8470，排列检验)。因此，尽管重组世系宿主间传播动态的自然选择也会影响观察到的重组断点位置分布¹¹，我们的数据表明，其他偏见影响了SARS-CoV-2基因组重组事件的分布。这可能包括影响重组断点分布的中性机械偏差。

尽管重组尚未在流行的SARS-CoV-2基因组中广泛传播，但重组已显著促进了SARS-CoV-2谱系的遗传多样性。由重组所贡献的可变位置比率(R)与由从头突变产生的基因(米)，R/米，通常用来概括这两种变化源的相对影响²²．利用我们的假定重组事件数据集，我们进行了估计R/米SARS-CoV-2 = 0.00264(补充文本16)。这对于冠状病毒人群来说是低的(例如，对于中东呼吸综合征，R/米估计为0.25-0.31(参考。²²))，这大概反映了在大流行的最初阶段，可能的重组祖先之间的遗传多样性极低，以及我们的方法的保守性质。随着新冠病毒群体不断积累遗传多样性，并与其他病毒共同感染宿主，重组在产生功能性遗传多样性方面将发挥越来越大的作用，这一比例可能会增加²⁴．因此，随着大流行的发展，waves将在检测新的重组谱系和量化其对病毒基因组多样性的影响方面发挥主要作用。

我们广泛优化的waves实现使其能够搜索整个系统进化树，并检测SARS-CoV-2谱系内部和之间的重组，而无需先验地定义一组谱系或分支定义突变。这是我们的方法相对于其他方法的一个关键优势，这些方法通过减少可能的重组事件的搜索空间来处理SARS-CoV-2数据集的规模(例如，参考文献。^16，17，25)。waves在同一Pango谱系的分支中发现了223个重组事件。我们的结果还包括366个谱系间重组事件(补充表1)。此外，我们发现重组影响了穿山甲SARS-CoV-2命名系统的证据²³．具体来说，我们发现B.1.355谱系的根可能是由属于B.1.595和B.1.371谱系的节点之间的重组事件引起的(图2)。3.及补充表1)。这些不同的重组事件突出了在《涟漪》中所采取的方法的多功能性和优势。

在含有刺突蛋白的SARS-CoV-2基因组3'部分检测到重组率增加，突出了持续监测的实用性。刺突蛋白是病毒谱系功能新颖的主要位置，因为它们适应人类宿主内部和之间的传播。我们发现，特别是在刺突蛋白周围的重组事件过多，目前循环中的重组蛋白水平相对较高，这强调了通过实时分析病毒基因组，监测通过突变或重组产生的新病毒谱系进化的重要性。我们的工作还强调了明确考虑系统发育网络将对准确解释SARS-CoV-2序列产生的影响¹¹．

除了SARS-CoV-2，重组是推动病毒和微生物适应的主要进化力量。它可以推动抗生素耐药性的传播⁷，耐药¹，免疫和疫苗逃逸²．重组的鉴定是病原体进化分析管道的重要组成部分，因为重组会影响系统发育、传播和系统动力学推理的质量^3.．基于这些原因，近年来检测微生物重组的计算工具已经变得非常流行和重要⁴．SARS-CoV-2大流行推动了病原体基因组测序和数据共享的空前激增，这反过来凸显了当前软件在研究大型基因组数据集方面的一些局限性⁵．涟漪是为大流行规模的数据集构建的，并经过充分优化，可以在40分钟内彻底搜索有史以来推断的最大的系统发育之一的重组(补充文本)17)。我们预计waves在密集采样的基因组数据集上表现最好，这可能会成为许多全球分布的病原体的标准，但我们警告说，它还没有在其他物种上得到验证。为了方便实时分析全球不同研究小组每天产生的数万个新SARS-CoV-2序列的重组^26，27，28， ripple提供了一个选项，可以在几分钟内评估任何用户提供的样本中重组血统的证据17)。因此，waves为在大量采样和快速进化的病原体群体中快速分析重组打开了大门，并为在大流行期间实时调查重组提供了工具。

waves使用突变注释树(MATs)的高效空间数据结构^20.，其中系统发育树的分支被注释为已经推断发生在它们身上的突变，以识别重组事件。数字1说明底层算法。waves识别假定的重组节点，包含用户指定的至少数量的突变，并通过计算从根开始的路径上分支上注释的所有突变，推断出在其相应序列上发生的突变集。然后，waves在突变位点上添加一个或两个断点，并使用与初始简约相比的部分放置来评估简约得分的提高。详情请参见补充文本1．为了确定假定的重组是否有意义，我们通过随机选择节点并添加，建立了一个零模型k从我们全局树的实际突变谱中提取的其他突变。然后，我们将这些样本放在树上，并使用waves来确定他们的节俭分数改进(补充文本)2)。对于我们全局树中的每个假定重组，我们将其简约分数改进与相同初始简约分数的零简约分数改进的分布进行比较3.)。我们首先采用2021年5月28日的公共树来开发我们的起始树^19，20.，屏蔽所有有问题的网站²⁹，以及修剪含有少于28000个非n核苷酸和含有两个或两个以上非[ACGTN-]核苷酸的样本(补充文本)5)。在此之后，我们通过运行matOptimize(补充文本)来优化该树4)两次，第一次子树修剪和再嫁接(SPR)半径为10，随后为40，并以蒙面变量调用格式(VCF)文件作为输入。使用waves的说明可在https://usher-wiki.readthedocs.io/en/latest/tutorials.html．我们在包含224个虚拟中央处理器(vcpu)的n2d-highcpu-224谷歌云平台实例上运行了waves(补充文本)18)。

为了测试waves的敏感性，我们从至少有10个后代的系统发育中随机选择两个内部节点，并在整个基因组中随机选择断点来模拟重组样本。我们分别为一个和两个断点重组生成了1000个模拟，在重组事件后，没有、一个、两个和三个额外的突变添加到序列中https://github.com/bpt26/recombination/．这些组合总共产生了2000个模拟重组谱系。然后，我们测量了作为断点位置函数的waves检测断点的能力以及重组节点到任一父节点的最小遗传距离(补充文本)6;遗传距离是根据推断出的分离病灶样本、谱系或节点的突变数量来估计的)。我们还通过确保ripple检测到Jackson等人的每一个高置信重组SARS-CoV-2聚类来评估其敏感性。¹⁶．

我们应用了几个事后过滤器来去除假定的重组节点，这些节点可能是由几个可能的错误源引起的假阳性。对于每个组成重组事件的三个内部节点(假定的重组节点、供体节点和受体节点)的每个内部节点，我们从COG-UK、GenBank、GISAID和中国生物信息学国家中心下载了每个节点最近的后代的共识基因组序列。然后，我们使用MAFFT对每个三人组的所有后代序列进行对齐^30.该方法特别关注重组信息位点，也就是说，重组节点的等位基因与一个父节点匹配，而与另一个不匹配。如果重组信息突变靠近索引位点或缺失碱基，或者重组的整个碱基是20个核苷酸范围内的单簇突变(补充文本)7)。我们还通过手动检查10个样本的原始读取数据来确认序列质量，我们可以自信地将原始序列读取数据与给定的共识基因组联系起来(补充文本)8)。为了估计与我们选择的特定方法和统计阈值相关的错误发现率(FDR)，我们计算了一个事后经验FDR。我们得到了我们测试的内部节点的数量，并与给定的简约分数相关联。然后，对于每个初始简约分数和简约分数的改进，我们获得了在null模型下显示简约分数改进的预期内部节点数。我们的扩展数据表3.)为给定初始和最终简约得分的期望节点数与检测到的具有相同初始和最终简约得分的重组节点数之比(补充文本)9)。

我们还使用样本元数据进行了事后分析，以确定重组节点的祖先是否具有比预期的更高的空间或时间重叠。我们将地理重叠计算为从供体节点和受体节点的后代中选择来自同一国家的样本的联合概率。对于时间重叠，我们分别记录来自供体和受体的最早到最近的样本间隔，并计算分隔这两个间隔的最小天数(重叠间隔为0)。对于每一个检测到的三元组，我们从树中选择两个随机的内部节点，它们的后代分别等于真正的供体和受体，从而为这两个类别生成一个空分布。然后，我们以同样的方式计算该随机集的地理和时间重叠(扩展数据图。4及补充文本10)。

为了确定已识别的重组断点是否显著地向基因组的3'端转移，我们进行了置换检验，比较均匀模拟断点分布的平均值与真实集中检测到的断点位置分布的平均值之间的差异(补充文本)12)。我们还使用变更点R包进行了变更点分析³¹并对重组预测区间中点计数进行泊松模型拟合。然后，我们计算了在已识别的变化点两侧的区间内重组断点的平均速率，以估计基因组3'部分重组率的倍增率(补充文本)13)。估计R/米，我们发现吝啬分数的下降与每个检测到的重组事件相关，作为估计R．然后我们计算米通过取这个值并从我们整个系统发育中观察到的突变总数中减去它(补充文本)16)。R/米是这些值的比值。

报告总结

有关研究设计的进一步资料，请参阅自然研究报告摘要链接到这篇文章。