可可泛基因组-文献精读112
Genomic structural variants constrain and facilitate adaptation in natural populations of Theobroma cacao, the chocolate tree
基因组结构变异在可可树（Theobroma cacao）自然种群中的适应性限制与促进作用
意义
基因组结构变异（SVs）是适应和物种形成的重要因素，但我们对其整体适应性后果的理解仍然有限，目前的数据和分析主要集中在人体和寿命较短的驯化物种上。在本研究中，我们利用31个高质量的基因组组装，研究了结构变异在可可树（Theobroma cacao）自然种群中的进化影响。我们发现，大多数SVs是有害的，从而限制了适应性。这些不利效应可能是由于基因功能受损以及重组受到抑制所直接或间接造成的。然而，我们也发现了一些SVs，可能通过与病原抗性相关的性状促进局部适应。总体而言，我们为理解结构变异在自然种群中的适应性效应过程提供了重要的见解。
摘要
基因组结构变异（SVs）在适应和物种形成中可能起到重要作用。然而，SVs的整体适应性效应尚不清楚，部分原因是准确的种群级别SVs鉴定需要多个高质量的基因组组装。本文中，我们利用31个染色体尺度的、解决了单倍型的可可树（Theobroma cacao）基因组组装——一种异交的长寿命树种，是巧克力的来源——研究SVs在自然种群中的适应性后果。在这31个样本中，我们发现了超过160,000个SVs，它们覆盖的基因组区域是单核苷酸多态性和短插入缺失（SNPs和短indels）的8倍（分别为125Mb和15Mb）。我们的结果表明，这些SVs中绝大多数是有害的：它们以低频率分布，并且在基因组的功能区域中被抑制。我们展示了SVs影响基因表达，可能通过破坏基因功能来导致SVs的有害效应。我们还为理论预测提供了实证支持，特别是倒位SVs通过抑制重组积累有害核苷酸变异，从而增加基因负担。尽管总体上SVs表现为有害效应，我们还是鉴定出了若干可能促进局部适应的个体SVs，其中几个与种群间差异表达的基因相关。涉及病原抗性的基因在这些候选基因中表现出明显的富集，突显了SVs在这一重要局部适应性性状中的作用。除了揭示SVs在进化上的重要性外，这31个全新组装的基因组也为可可树的遗传和育种研究提供了宝贵资源。
引言
超过一个世纪以来，基因组结构变异（SVs）被认为是功能变异的重要来源（1-3）。这种变异影响核苷酸序列的存在、数量、位置和/或方向，通常影响比单核苷酸多态性（SNPs）更大比例的基因组区域（4-8）。SVs还可以对表型产生显著影响（9-12），并有助于适应和物种形成（13, 14）。然而，大规模基因组分析表明，SVs通常以低频率分布，并且在基因组的功能区域中被抑制，表明它们经历了强烈的净化选择（4, 7, 15, 16）。尽管有这种普遍的观察，导致这些适应性效应的过程仍然不太清楚。SVs的适应性后果可能是由于直接影响基因功能，如破坏编码区或调控元件（17-19），或者它们可能通过抑制重组间接产生效应（20, 21）。
SVs通过抑制重组可能在局部适应中发挥重要作用，因为理论上SVs可以通过阻止或减少在染色体杂合体内形成有效的交叉，从而使局部有利等位基因免受基因流的影响（22, 23）。与这一预期一致，SVs——特别是倒位——已被证明与多种物种中的局部有利表型相关，包括金鱼草（Mimulus guttatus）（24）、玉米（Zea mays）（25）、黄蓍（Boechera stricta）（26）和向日葵（Helianthus annuus）（27）。然而，抑制重组也有其弊端，因为它降低了安排的有效种群大小（Ne）。较低的Ne反过来可能削弱净化选择的效果，从而增加SVs中的有害突变积累（28, 29）。
为了研究SVs的适应性后果，我们为31个野生收集的可可树（Theobroma cacao，以下简称“可可”）样本构建了染色体尺度、解决了单倍型的全新基因组组装。大多数先前关于SVs进化作用的研究都集中在短命的、自交的和/或驯化的物种上。相比之下，可可树是一种主要通过异交繁殖、寿命较长的多年生植物，且样本来源于未驯化种群（30, 31）。结合野生收集的样本和基于二倍体组装的SV鉴定——这是SV鉴定的黄金标准（32）——使我们的数据集在研究自然种群中SVs的进化影响时具有独特优势。我们提出了以下问题：SVs的整体适应性效应是什么？SVs是否影响基因表达，提示基因功能受损？SVs是否由于抑制重组而积累有害核苷酸变异？哪些SVs可能是局部适应的重要贡献者？
结果
为了鉴定在31个野生收集的可可样本中存在的SVs（图1A），我们使用10x Genomics链接读技术构建了全新基因组组装（SI附录，图S1）。我们62个单倍型特异性的基因组组装的总大小在341到387 Mb之间，scaffold N50介于35和41 Mb之间，contig N50介于94和188 Kb之间（数据集S1和S2）。这些组装还捕获了超过96%的普遍保守的单拷贝基准（BUSCO）基因（数据集S3），并且与两个可可参考基因组，Matina 1-6版本1.1（33）和Criollo B97-61/B2版本2.0（34），表现出高度的共线性（SI附录，图S2和S3）。总体而言，我们的染色体尺度、单倍型解决的基因组组装在质量上优于之前发布的两个参考组装（图1B和SI附录，表S1）。
SVs通过62个高质量基因组组装鉴定。 (A) 研究种群的大致位置。 (B) 两个已发布参考基因组（黑线）与我们的31个组装（每个样本一个单倍型；彩色线）之间的序列连续性（基于5000个最长的contigs）。对于每一个我们的组装，累计序列长度的增加速度比参考基因组更快，表明具有更高的连续性。 (C) 为每个样本鉴定的SV数量。 (D) 不同基因组特征中重叠的SV数量，每个样本单独计数。 (E) 基因组不同区域的唯一定位SV数量，以2 Mb为窗口进行计数。
通过将每个62个组装与Criollo参考基因组比对，我们鉴定了五种类型的SVs：插入（INS）、缺失（DEL）、串联重复（DUP）、倒位（INV）和易位（TRA）。我们进一步使用二倍体组装来区分杂合子（SV在一个单倍型中发现）和纯合子（SV在两个单倍型中发现）。我们鉴定了36,303个定位于50 bp以上的唯一SV（总计163,423个SV），每个样本中的变异数为4,610到5,963个（图1C和SI附录，表S2）。每个样本中的杂合SV比例与SNP数据一致（皮尔逊相关系数r = 0.97），尽管SV的纯合性通常高于SNP（SI附录，表S3）。总的来说，这些SVs覆盖了比SNPs和短插入缺失（INDELs，<50 bp）更大比例的可可基因组（分别为125 Mb和15 Mb）。正如预期，较短的变异比较长的变异更常见：15,829个SV小于1 Kb，14,104个SV在1到10 Kb之间，6,186个SV在10到100 Kb之间，152个SV在100 Kb到1 Mb之间，32个SV大于1 Mb（SI附录，图S4和表S4）。虽然大多数SV位于基因间区（图1D和SI附录，表S5），但Criollo基因组中64%的基因有SV与编码区域或推定的调控元件重叠（≤5 Kb的基因上游或下游）。SVs几乎存在于基因组的所有部分，但某些区域的SVs数量显著高于平均水平（图1E），表明这些区域可能是SV热点。总体而言，我们在31个样本中检测到丰富的结构变异，表明SVs可能对可可的进化产生显著影响。
SVs总体上具有不利影响。
为了确定SVs的总体适应性效应，我们检查了它们的频率模式。我们对每种SV类型进行了分析，但由于易位数量非常少（图1 C和D），我们在这些及后续分析中将其排除。我们还将插入和缺失（INS和DEL）合并为一种INDEL类型，因为INDELs的突变不能根据参考基因组定义（例如，查询中的插入在参考中可能是缺失），这可能会影响适应性推断。次等位基因频率谱（AFS）显示，SVs的分布频率明显低于同义或非同义SNPs（图2A），这表明SVs在适应性上的负面效应通常大于SNPs（35）。然而，我们注意到，偏斜的AFS也可能是由于与有限的有效种群大小（Ne）、较高的点突变率以及SV的突变率显著较低相关的突变-选择平衡的变化（SI附录，图S5）。尽管如此，基于种群动态和突变率的模拟表明，85%以上的新SVs具有强烈的不利效应，而非同义SNPs的不利效应比例为56%（SI附录，图S6）。SVs的不利适应性效应也得到了大量频率差异的支持，这些SVs重叠在功能性和（假定）非功能性元素之间，这一差异约为非同义和同义SNP之间差异的6倍（图2B）。
*SVs的适应性效应。 (A) 每种SV类型的等位基因频率谱（AFS），与同义（sSNP）和非同义（nSNP）核苷酸变异进行比较。 (B) 重叠不同基因组特征的SVs的AFS。 (C) 受SVs影响的基因的选择性约束度量。展示的是SV内基因和重叠SV断点的基因的非同义/同义核苷酸多样性比率（πN/πS）和非同义/同义核苷酸分化比率（dN/dS）。灰色水平线表示对照基因的中位数。 (D) 在基因5 Kb范围内的变异与其表达相关的概率。SV结果与随机数据（RAND）和SNPs进行比较。常见变异：MAF > 0.05，稀有变异：MAF ≤ 0.05。误差条表示95%置信区间（CI）。
SVs重叠功能性元素的强烈去除表明基因功能受损，接下来我们检查了选择性约束的放松是否允许SVs在这些基因中被保留。为此，我们通过计算可可树内的非同义与同义核苷酸多样性比率（πN/πS）和可可与其近亲Herrania umbratica之间的非同义与同义核苷酸分化比率（dN/dS）来估算净化选择的强度。由于净化选择往往会降低功能性变异的频率，增加的πN/πS和dN/dS表明选择性约束放松。我们发现，位于SVs内的基因和重叠SV断点的基因的πN/πS和dN/dS比预期的要高（图2C；P ≤ 0.0001，Wilcoxon秩和检验）。这一效应在不同SV类型之间有所不同：重叠INV断点的基因具有比INVs内基因更高的πN/πS和dN/dS，而INDELs和DUPs则恰恰相反（图2C）。这些结果表明，平均而言，重叠SVs的基因受到的选择性约束较弱，可能是因为这些基因大多数在生理和发育中扮演非关键角色。
SVs影响基因表达。
为了进一步了解SVs的潜在适应性效应，我们研究了它们对基因表达的影响。我们首先检查了常见SVs（次等位基因频率[MAF] > 0.05）如何影响附近基因（≤5 Kb）的表达。通过一种简单的定量方法，我们检测到214个基因（占测试基因的3%）在携带主型和次型安排的样本之间表现出差异表达（SI附录，图S7）。然而，在控制由于SVs附近映射精度下降和样本间种群结构可能引起的偏差后，影响附近基因表达的SVs显著减少：在测试的8,004对中，只有12对（0.1%）通过了名义假发现率（Q值）阈值0.1。对于这12个基因，SV基因型解释了39%到84%样本间的表达差异（SI附录，图S8）。尽管只有少数基因显示出强烈的表达受SVs影响的信号，但通过随机分配基因型到基因上，我们发现常见SVs比预期的更频繁地与基因表达相关（图2D；P = 1 × 10−9，似然比检验[LRT]）。然而，使用SNPs重复分析时，结果表明常见SVs对基因表达的影响大约是常见SNPs的三分之一（图2D；P = 4 × 10−6，LRT）。
接下来，我们评估了稀有SVs（MAF ≤ 0.05）影响基因表达的潜力。与常见变异相比，稀有SVs更可能是最近出现的，并对基因表达产生有害影响。通过对比次型安排的表达与主型安排的表达分布，我们检测到405个SVs（占测试SVs的3%）对基因表达有潜在影响（SI附录，图S8）。与常见变异的模式不同，稀有SVs对基因表达的影响是稀有SNPs的四倍左右（图2D；P < 2 × 10−16，LRT）。尽管我们仅检测到少数倒位（INVs）与基因表达相关（SI附录，图S8），但等位基因特异性表达（ASE）数据揭示，在倒位杂合子中具有ASE的基因呈过度代表（SI附录，图S9），这表明倒位可能影响其中基因的表达。综上所述，这些结果表明，常见SVs可能是中性或有益的，它们对基因表达的影响较为微弱。相比之下，稀有SVs更倾向于影响基因表达，这可能会损害基因功能并加剧SVs的有害效应。
倒位增加遗传负担。
除了直接损害基因功能外，SVs的有害效应还可能通过有害核苷酸变异的积累间接产生。这一预测特别适用于倒位（INVs），因为它们通常导致重组单倍型的完全丧失（14, 20, 21）。与抑制重组的效应一致，我们发现倒位单倍型之间的核苷酸分化增加（图3A；P < 2 × 10−16，Wilcoxon秩和检验），并且在这些单倍型内的基因链接性也增加（图3B；P < 2 × 10−16，Wilcoxon秩和检验）。
倒位中的遗传负担 (A) 主型和次型纯合子（INV）之间的绝对核苷酸分化（dXY），与随机共线区域（CTRL）相比。 (B) 在共线区域和主型及次型倒位纯合子中的连锁不平衡衰退的平均情况，作为物理距离的函数（使用相同样本量估算）。 (C) 在共线区域和主型及次型倒位纯合子中的衍生核苷酸等位基因的百分比。结果分为基因间区域（>5 Kb与每个基因的距离）、同义位点（4倍）和非同义位点（0倍）。误差条表示95%置信区间（CI）。 (D) 预测为有害的非同义SNPs的百分比。误差条表示95%置信区间（CI）。
抑制重组预计会降低这些单倍型的有效种群大小（Ne），可能导致有害核苷酸变异的积累增多（28）。实际上，倒位（INVs）在功能性位点的衍生等位基因比例高于共线区域（图3C；P < 2 × 10−16，LRT），表明遗传负担增加（35）。在推定中性位点（图3C；P < 2 × 10−16，LRT）也发现了类似的衍生等位基因增加，可能反映了倒位区域背景选择较弱，与共线区域相比（见SI附录，图S10的模拟结果）。我们注意到，如果倒位主要捕获非必需基因，即使没有SVs，它们也可能包含更多的衍生等位基因。然而，次型倒位单倍型的衍生等位基因频率也高于主型倒位单倍型（图3C；P < 2 × 10−16，LRT），这表明这些区域的Ne因抑制重组而降低（28）。通过SIFT4G（36）预测的突变效应进一步支持了遗传负担增加的结果，SIFT4G鉴定出了倒位区域的有害SNPs比共线区域更多（图3D；P < 2 × 10−16，LRT），并且次型倒位单倍型的有害SNPs数量高于主型倒位单倍型（图3D；P = 1 × 10−5，LRT）。综上所述，我们的结果支持理论预测，即抑制重组会导致倒位区域的遗传负担增加，并且这种负担在次型倒位单倍型中大于主型倒位单倍型，可能是由于次型倒位单倍型的频率较低，并且由于其更多地处于杂合状态，导致重组水平较低（SI附录，图S10）。
SVs对可可树局部适应的贡献
尽管SVs作为一类变异通常表现为有害，因此限制适应性，但其中一些可能有助于局部适应。本研究中的31个样本采自四个环境不同的地点（SI附录，图S11和S12），提供了识别有助于适应性分化的SVs的机会。我们注意到，SVs与SNPs表现出相似的种群结构模式，尽管不那么明显（图4A和SI附录，图S13）。与SNPs（FST：INDEL = 0.22，DUP = 0.20，INV = 0.13）相比，SVs的种群结构信号较弱（FST：SNPs = 0.24），这似乎不是由基因流引起的，而是由于SVs的次等位基因频率（MAF）普遍低于SNPs。通过在MAF上对变异进行条件限制，我们发现，SVs的种群结构通常比SNPs更为显著（SI附录，图S14）。
SVs与局部适应 (A) 通过主成分分析（PCA）得到的SVs的前两个特征向量的变异。 (B) 三种SV类型的FST估计分布，与模拟的中性样本（SIM）进行比较。 (C) 由330 Kb倒位（INV）引起的一个非重组单倍型区块示例（阴影区域）。展示了Iquitos和Nanay种群之间基于SNP的FST估计。 (D) 受选择异常SVs影响的基因的表达水平。展示了五个基因，其中最大的表达方差比例由SVs解释（R2）。 (E) 126个重新测序的可可样本（30）与我们基于组装的SVs进行基因分型的关系。这些样本根据它们贡献最大的祖先种群进行了标注。展示了基于全基因组SNP数据的t分布随机邻域嵌入（t-SNE）投影（SI附录，图S19展示了基于SV数据的t-SNE）。 (F) 基于祖先的基因组扫描的P值，使用了为126个样本基因分型的3,011个SVs。红色虚线：Q < 0.1。
为了识别可能有助于局部适应的SVs，我们将SVs的FST估计与基于种群的动态模型模拟的中性变异进行比较，后者基于这些种群的有效种群大小（Ne）、分化时间和迁移率估算（图4B）。在每一对种群比较中，136到291个SVs的FST估计超过了它们模拟分布的第99百分位（总共846个独特SVs）。我们将这些SVs视为可能有助于局部适应的候选基因。绝大多数SVs为INDELs，但五个INV也被识别为异常值。这些INVs可能促进了Iquitos和Nanay（三个INVs）以及Marañon和Guiana（两个INVs）之间的适应性分化，提示抑制重组的INV单倍型区块（图4C）可能将局部有益的等位基因隔离于基因流之外。这些INVs涉及四个与压力反应、细胞分化和蛋白质合成相关的基因（数据集S4）。
为了理解这些异常SVs的潜在表型效应，我们检查了864个候选SVs附近5 Kb内的769个基因的基因本体（GO）注释。这些基因富集了22个GO术语（Q < 0.1，超几何检验），包括与对抗细菌和卤霉菌（SI附录，图S15）相关的术语。在这些基因中，SV基因型解释了样本间较高比例的表达方差（图4D和数据集S5）。使用所有测试的SVs（MAF > 0.05）数据，并位于基因5 Kb内，我们发现选择异常的SVs比其他SVs更可能影响基因表达，约为6倍（SI附录，图S16；P = 2 × 10−16，LRT）。在140个差异表达的异常基因中，五个GO术语显著富集，其中“防御反应（细菌）”（Q = 1 × 10−5，超几何检验）是最强的富集术语。为了确定通过SNP分析是否能够找到相同的基因，我们重复使用SNPs（为每个基因估算加权FST）进行异常值分析。尽管通过SV分析识别的基因具有高于平均水平的SNPs加权FST估计（SI附录，图S17），但只有47个基因在两种方法下都被识别为选择异常基因。我们还发现SV和SNP异常基因的GO术语之间没有重叠（SI附录，图S18）。
为了进一步研究SVs在局部适应中的作用，我们利用了126个可可样本（图4E）的公开全基因组数据，这些样本采自该物种的原生范围（30）。在这些样本的组装基因组中，3,311个SVs（MAF > 0.05）是分化的，我们的祖先基因组扫描方法识别了45个作为局部适应的候选基因（图4F）。与这些SVs≤5 Kb内的基因相关的“防御反应（细菌）”富集（Q = 8 × 10−6，超几何检验）（数据集S6）。综合来看，我们的结果表明，SVs对可可树的局部适应有所贡献，这些有益的适应效应可能通过与病原抵抗相关的性状体现出来。
材料与方法
我们使用10x Genomics linked-read技术构建了31个野生采集的可可样本的染色体尺度、单倍型分辨的基因组组装。每个基因组大约有65倍的原始读取覆盖度，我们将linked reads组装成相位化的contig和scaffold，使用Supernova（47），然后使用RaGOO（48）将它们锚定到可可栽培品种Matina 1-6（版本1.1）参考染色体（33）。我们使用MUMmer4（49）将每个62个单倍型特异性组装与Criollo B97-61/B2（版本2.0）参考基因组（34）对齐，并通过MUM&Co（50）从比对中识别SVs。我们通过估算覆盖不同基因组区域的SV的等位基因频率谱，并使用fit∂a∂i（51）推断适应效应的分布，来评估作用于SV的净化选择强度。
为了量化基因表达，我们从在哥斯达黎加的共同花园中维持的树木中收集了叶片样本。我们收集了与基因附近（≤5 Kb）的SV数据，并测试了常见（MAF > 0.05）和稀有（MAF ≤ 0.05）SVs是否与这些基因的表达相关。常见的SVs通过线性回归模型分析，使用DNA测序（DNA-seq）读取计数和基于基因型的10个主成分作为协变量，分别用于控制SV周围的映射准确度和样本间的种群结构。我们使用绝对中位Z标准化和正态累积分布函数来测试稀有SVs是否影响基因表达。对于这两种分析，我们使用基于假发现率（FDR）的Q值（52）来考虑多重检验问题。
我们利用标准Illumina测序的短读数据识别位于倒位中的SNPs。我们使用dXY（53）估算主型与次型倒位单倍型间的遗传分化，使用Mangin等人的方法（54）量化倒位内的连锁不平衡，并通过估算功能性位点上衍生的和有害的等位基因的比例来确定倒位单倍型中的遗传负担（35）。祖先与衍生等位基因的推断是使用与倒位相关的近缘物种H. umbratica的基因组序列作为外群进行的。我们使用SIFT4G（36）预测突变效应，以区分倒位区域中的有害和耐受的SNPs。
我们通过寻找种群间高度分化的变异来研究SVs如何促进局部适应。种群分化通过Hudson等人的FST估计量（55）进行量化，观察到的估计与使用msprime（56）模拟的中性样本进行比较，以检测异常值。我们通过线性模型测试局部适应候选SVs是否与基因表达相关，从而评估它们如何影响基因功能，同时使用DNA-seq读取计数控制技术性变异。我们还使用Paragraph（57）对我们基于组装的SVs进行了基因分型，使用来自126个样本（30）的公开重新测序数据，并在更大的数据集中寻找分化异常值。为此，我们使用多元线性回归和马哈拉诺比斯距离寻找最能解释样本间祖先比例分布的SVs。祖先比例通过ADMIXTURE（58）估算。然后，我们检查了在相同生物过程中的SVs是否在物种范围内受到选择。