大豆重测序二（同一领域互竞）-文献精读58

High-quality genome of a modern soybean cultivar and resequencing of 547 accessions provide insights into the role of structural variation

现代大豆品种的高质量基因组及对547个种质资源的重测序揭示结构变异的作用

大豆重测序-文献精读53

摘要

大豆提供蛋白质、油脂以及多种与健康相关的化合物。理解结构变异（SV）对现代育种中经济性状的影响对大豆改良至关重要。在此，我们组装了现代品种农大豆2号（NDD2）的高质量基因组，并在与已报道的29个基因组的比较中鉴定出25,814个SV-基因对，其中在547个深度重测序（平均覆盖度=18.05倍）的种质资源中验证了13个NDD2特有的SV，这促进了我们对基因组变异生物学的理解。我们发现了一些与种子蛋白质和重量形成相关的插入/缺失变异，一个与干旱适应相关的倒位，以及一个与大豆关键分化事件相关的大型染色体互易易位。在547个种质资源中鉴定出的749,714个SV中，有6,013个与在十个不同地点和年份环境下确定的22个与产量和种子品质相关的性状显著关联。我们还发现了1,761个关联SV落在基因或调控区域，其中12个位于GmMQT，影响油脂和异黄酮含量。我们的研究为大豆改良中的SV作用提供了资源和见解。

正文

大豆（Glycine max）是一种重要的豆科作物，提供了全球超过一半的油料作物产量和四分之一的蛋白质，以及许多与健康和药物相关的化合物1,2,3。理解大豆产量、种子品质及其他重要特性的遗传基础，对于科学界来说，高质量的参考基因组具有重要的生物学意义。尽管越来越多的大豆基因组已经被组装，这对于揭示植物驯化和育种的见解至关重要，但几乎所有已组装的种质资源都在2010年之前获得批准（参考文献1,4,5,6,7,8,9,10,11）。自2011年起（中国大豆育种和种子工业的第三阶段），大豆品种的产量和种子质量改良取得了巨大进展；然而，现代品种基因组的表征仍然十分有限。

结构变异（SV），包括插入、缺失、倒位和易位，已被证明在植物进化、驯化和育种中扮演着越来越重要的角色12,13。泛基因组的构建揭示了野生和栽培大豆中大量非冗余的SV11,14，其中两个SV被发现决定了种皮颜色和光泽11。在POWR1的CCT结构域中的插入/缺失被证明可以增加大豆种子的重量和油含量，同时降低蛋白质含量15。许多SV被证明影响了番茄果实的特性及其对各种病害的抗性16。在棉花中，通过参考现代品种基因组，挖掘出了数百个SV，这些SV与纤维质量和产量相关性状在数十年的育种过程中紧密相关12。与其他作物相比，SV对大豆重要性状如产量、种子质量和营养成分的遗传影响仍需进一步研究。已经验证，增加种质资源和提高测序深度是关联研究中进行变异探索的有效策略17,18。因此，有必要基于现代品种基因组并通过深度重测序更多种质资源，全面揭示大豆SV的作用。

在本研究中，我们报告了现代品种农大豆2号（NDD2）的高质量参考基因组，该品种于2014年发布并自那时以来被广泛种植（补充图1）。NDD2展现了许多优良特性，如大粒种子（百粒重（HSW）=26.9克）、高产量（在河北省大豆区域试验中较对照显著提高4.12%）、高种子蛋白质含量（SPC；43.56%）、对大豆花叶病毒的高抗性（SC3和SC7株系的病害指数分别为12%和6%），并成功克服了产量、种子质量和病害抗性之间的权衡19（补充表1）。此外，我们对来自中国九省和美国大豆种植区的547个种质资源进行了迄今为止最深入的重测序（~20倍覆盖度），大约70%的资源不同于之前的研究11。基于三年七个地点的十个环境中31个植物学和经济学重要性状的表征，我们鉴定了大量未在先前基因组中报道的SV和基因。现代大豆品种的高质量基因组以及群体的深度重测序为大豆改良中SV的作用提供了见解。

结果

高质量的现代品种NDD2基因组增加了基因组资源

我们通过整合五种测序技术对NDD2基因组进行了测序，生成了总计138.51 Gb的PacBio单分子实时序列，N50为26.74 kb；121.85 Gb的Nanopore序列，N50为26.64 kb；473.35 Gb的Bionano光学图谱数据，平均长度为267 kb；102.57 Gb的Illumina双端测序碱基以及362.86百万对染色质构象捕获（Hi-C）相互作用对（补充表2-5）。我们应用了一种改进的策略20，组装了NDD2的染色体级基因组，总大小为1,013.66 Mb，contig N50值为27.16 Mb，与29个泛基因组种质资源中表现最好的基因组相当1,7,9,11（表1，图1a,b及补充图2）。为了评估组装基因组的质量，我们将高质量的双端测序读段重新比对到组装上，达到了99.75%的比对率和99.56%的基因组覆盖率（补充表6）。进一步的评估显示，组装的质量值（QV）为41.83，超过了脊椎动物基因组计划的QV40标准（参考文献21）。此外，基因组中约99.70%的胚基植物基准通用单拷贝直系同源基因（BUSCOs）和98.79%的超保守真核基因存在（补充表6）。我们还与其他29个已发表的染色体级大豆基因组进行了共线性分析，得到了平均共线性率96.34%的结果（补充表7）。我们在NDD2基因组的所有染色体中恢复了Cent91/92大豆特有的着丝粒重复序列22。值得注意的是，我们成功生成了端粒序列（即CCCTAAA/TTTAGGG重复序列23），这在除四条染色体外的单端组装中较为困难（图1c）。与29个泛基因组资源相比1,7,9,11，我们生成的端粒组装（40条中的36条）为最佳（补充表8）。综上所述，我们组装了一个具有良好准确性、完整性和连续性的现代大豆基因组。

Genomic features	NDD2
Assembled genome size (Gb)	1.01
Percentage of anchoring (%)	97.65
Number of contigs	256
Number of scaffolds	187
Contig N50 (Mb)	27.16
Scaffold N50 (Mb)	50.87
Gap number	69
GC content (%)	35.03
Repeat ratio (%)	51.46
Predicted PCG model number	58,899
BUSCOs (%)	99.70

a. 与其他27个contig N50大于10 Mb的大豆基因组相比，NDD2基因组的contig N50值增加。我们将contig按从长到短的顺序累积，当contig累积长度达到总基因组长度的x%时，相应contig的长度为n(x)。图中的黑色虚线对应每个基因组组装的contig N50和N90值。 b. NDD2组装中20条染色体的染色质相互作用。每个热图的分辨率为500 kb。深红色点表示高概率的相互作用，浅色点表示低概率的相互作用。 c. NDD2基因组中着丝粒和端粒的密度分布与另外三种代表性大豆（W05、Wm82和ZH13）相比。 d. 根据将每条染色体划分为等宽的1,000个窗口，Gypsy、蛋白编码基因（PCG）和Copia的平均密度分布。

随后，我们预测了NDD2中560.92 Mb（占55.34%）的转座元件，其中长末端重复序列（LTRs）占比最高（44.39%），主要是Gypsy（42.85%）和Copia（20.20%；补充表9）。我们鉴定了58,899个蛋白编码基因（PCGs），其中96.70%具有功能注释，8,503个基因显示了Gypsy和Copia的插入（补充表10-12）。此外，我们估计96.51%的同源PCGs在与已发表的大豆基因组中的蛋白质序列相似度超过80%（补充表13），并且根据166个转录组数据，91.93%的基因表达得到了验证（补充表14），表明NDD2的98.73%的PCGs得到了支持。此外，2.38%（1,404个）的PCGs显示出小于20%的相似度，定义为新预测的基因（补充表15）。值得注意的是，我们发现48个新预测的基因位于至少一个之前发表的基因组中有缺口的序列中（补充数据1），其中17个基因得到了转录组或功能数据库的支持（补充表16），并通过PCR扩增和Sanger测序验证了8个随机选择的基因（补充图3）。

接下来，我们通过将每条染色体划分为1,000个等宽窗口，分析了每条染色体上PCGs的密度分布，发现在靠近染色体端粒的20%窗口中PCGs分布最为集中，与其他区域相比增加了0.57倍（P < 2.20 × 10−10，Wilcox检验），而与Gypsy和Copia元件的模式相反，分别减少了0.24倍和2.37倍（P < 2.20 × 10−10，Wilcox检验）（图1d）。该现代品种基因组为大豆重要性状的遗传解析和改良提供了基因组资源。

在NDD2中探索SV提供了对现代育种的见解

高质量的NDD2组装为探索现代大豆育种中的基因组SV提供了潜力。我们使用NDD2基因组和29个已发表的大豆基因组1,7,9,11构建了一个基于图的泛基因组，鉴定出47,058个非冗余SV，其中包括37,304个插入/缺失（INS/DEL；≥50 bp）、3,071个倒位（INV；1.01–29.14 kb）和6,683个易位（TRANS；1.01–19,040.53 kb）（补充数据2）。在INS/DEL中，平均93.94%和93.46%分别通过Illumina短读段和PacBio长读段成功进行基因分型（补充表17和18）。值得注意的是，我们根据NDD2参考基因组鉴定了25,814个SV-基因对，包括23,119个INS/DEL-基因对、719个INV-基因对和1,976个TRANS-基因对（补充表19-21），这些SV可能通过调控基因表达影响相关性状12,16。例如，在NDD2和W03之间鉴定出的8,609个SV-基因对中，71.79%（6,180个）表现出表达变化，NDD2和W03在性状上存在显著差异（扩展数据图1和补充表22）。总体而言，这些SV（包括新探索或未描述的SV以及在我们对547个种质资源的深度重测序群体中验证的SV）为研究SV在作物改良中的作用提供了资源平台。

首先，我们鉴定了13个NDD2特有的INS/DEL（补充表23），这些SV通过对NDD2自身基因组的长读段比对以及26个已发表的大豆基因组进行确认，并通过PCR实验和Sanger测序验证（补充数据3）。这13个SV（238-6,231 bp）也得到了547个重测序种质资源中的31-69个一致性资源的支持（补充表23），其中7个特有的INS/DEL能够以极显著的水平将547个资源分为两种纯合类型，分别对应产量相关性状（百粒重（HSW）、株高和每株豆荚数（PN））和种子品质性状（蛋白质、油脂、异黄酮和生育酚），表明这些SV可能与大豆育种中的性状分化相关（补充表24）。我们分析了其中一些反映了SV对相关基因及其对应性状影响的7个SV-基因对12,16。

a. NDD2基因组中DEL238与29个泛基因组样本的对比。灰色矩形表示29个样本中的del238。蓝色矩形表示Glyma.NDD2.06G308200（B细胞受体蛋白）的外显子。黄色矩形和紫色箭头分别表示5′ UTR和3′ UTR。DEL238位于Glyma.NDD2.06G308200的内含子中。 b. 根据DEL238和del238的547个重测序样本的SPC（种子蛋白质含量）和HSW（百粒重）小提琴图。在小提琴图中，红色中心线表示中位数，两侧虚线表示上四分位数和下四分位数；n表示具有DEL238或del238的样本数量。差异显著性通过双尾t检验分析。 c. 对不同大豆样本的DEL238及相应的SPC和HSW的PCR验证。样本1–11表示具有DEL238且SPC和HSW较高的样本；样本12–22表示具有del238且SPC和HSW较低的样本。SPC和HSW为相应11个样本的平均值。具体来说，样本1–11为NDD2、Ji13BB9、Jidou12、Cang0749、Cangdou0734、Ji06B9、Wuxing4、Lang2013-6、Langlunxuan2、Ji12B9和CangM0707；样本12–22为Wm82、ZH13、Zhonghuang58、Zhongpin036025、Yudou25、Jindou40、PI547875、PI547872、Zhongzuo036052、PI547843和Zhonghuang15。M表示DNA标记DL2000。比例尺为1厘米。 d. DEL238-高SPC-HSW的NDD2和del238-低SPC-HSW的LD14（Liaodou14）大豆种子中Glyma.NDD2.06G308200的差异表达。数据来自大豆种子的转录组。差异显著性通过双尾t检验分析。数据表示为平均值±标准误差（n = 2）。 e. DEL238和del238在大豆育种中对SPC和HSW的不同影响。

接下来，我们检查了与大豆育种分化相关的倒位事件。尽管已有研究关注种皮颜色驯化中的倒位事件9,11，但育种过程中倒位对大豆重要性状的影响较少被描述。我们首先在第5号染色体上鉴定了一个新的倒位（倒位05 (INV05)，3.06 kb），该倒位可能与大豆适应不同地理区域有关（图3）。我们发现，在我们的重测序样本中，INV05的频率从地方品种的37.70%（61个样本中的23个）增加到改良品种的69.75%（486个样本中的339个），表明INV05分化对大豆育种的贡献。我们发现一个与脱水反应元件结合蛋白（DREB）转录因子相关的基因Glyma.NDD2.05G253650位于INV05中，其断点位于启动子区域和3′ UTR。这导致了NDD2与其对应的W05（带有inv05）的基因在上游3 kb处有11个特定的调控元件（干旱、ABA和光反应等）和下游273 bp处（可能导致调控变化31,32）的巨大差异（75.83%），C02（Tiefeng18，带有INV05）和C01（Xudou1，带有inv05）的结果也证实了这一点。DREB在植物中响应各种非生物胁迫，例如低温或高温、盐度、低氮，尤其是干旱33,34,35,36，例如通过过表达大豆GmDREB1提高小麦抗旱性37。我们推测INV05与大豆在种植区域对干旱和温度的适应有关。由于95%的重测序样本来自中国黄淮海生态区，这是中国最重要的大豆驯化中心38，我们分析了该地区改良品种中的INV05分布。INV05在南北自然屏障——秦岭-淮河线附近呈现出有规律的变化，导致了气候条件的差异，南方气温和降水较高（湿润，年降水量超过800毫米），而北方为半湿润，降水量低于800毫米39。北方样本中INV05的显著更高频率（76.90%）与南方的较低频率（50.00%，P = 5.5 × 10−4）一致，表明大豆适应性与育种过程中的人工改良相符。同样，北方地区的省份中，山西为82.86%，河北为82.08%，山东为75.00%，河南为56.25%，这一趋势与从山西到河南干旱程度的下降相符（图3）。超过80%的频率与山西位于干旱易发的黄土高原地区，河北位于干旱易发的华北平原漏斗形区域，地下水短缺且降水较少相符。

最后，通过对NDD2与泛基因组样本的染色体共线性分析【11】并将这些大豆的Hi-C数据与NDD2作为参考进行比对，我们将两个染色体易位事件从约0.12 Mb和0.14 Mb【9】扩大到野生大豆W05和NDD2在第11和第13号染色体上的19.1 Mb和4.5 Mb（图4a）。这两个小的易位事件包含在我们识别的大型易位事件中，且具有相同的前置断点。除了W02之外，泛基因组的其余样本都表现出与NDD2相同的模式，W02被新识别为与W05具有内部易位（扩展数据图3）。基于这些信息，我们还通过序列分析发现了W05与野生大豆PI483463之间的染色体易位【7】。由于这些易位，野生大豆被划分为以下两种类型：类型I（W01、W03和PI483463）具有较长的13号染色体（45.7-47.6 Mb）和较短的11号染色体（39.2-41.9 Mb），而类型II（W02和W05）则表现出较长的11号染色体（53.5-54.7 Mb）和较短的13号染色体（31.5-31.9 Mb），这表明易位可能发生在野生大豆中，而NDD2及泛基因组中的其他代表性改良品种在育种过程中继承了类型I的特征。基于上述发现，我们提出了一个野生大豆的两种类型分化模型，这一模型从大型结构变异的角度阐明了栽培大豆的祖先，揭示了大豆中的一个关键分化事件（图4b）。

a. 交互信号图和共线性图。每个面板展示了其他大豆Hi-C数据与NDD2基因组比对后的交互信号热图，并显示了其他大豆基因组与NDD2基因组之间的共线性图。 b. NDD2与五个其他大豆基因组（W01、W02、W03、W05和ZH13）在Gm11和Gm13上的易位痕迹分析。NDD2与PI483463基因组之间的共线性在扩展数据图3中展示。

此外，通过分析NDD2与其I型祖先野生大豆之间的易位区域的基因组成变化，我们发现NDD2在13号和11号染色体上分别有261个和76个特有基因（同源性低于20%，见补充表25），这清楚地反映了育种过程中基因的分化。在这些基因中，我们发现一些可能与栽培大豆的特性相关，例如与花青素积累相关的F-box/kelch重复蛋白基因Glyma.NDD2.11G275400【40】，与确定型和不确定型植物生长习性相关的自我修剪样基因Glyma.NDD2.13G154100【41】，与病原识别和植物抗病性相关的富含亮氨酸重复序列受体样蛋白激酶Glyma.NDD2.13G343800和Glyma.NDD2.13G343900【42,43】，以及与种子萌发和多种非生物胁迫耐受性相关的类胚蛋白基因Glyma.NDD2.13G154300【44,45】，这些都符合栽培大豆的进化和改良方向。此外，我们还发现NDD2与三种野生大豆相比，在13号和11号染色体上分别有35个和9个基因拷贝数变异（CNV，见补充表26），其中43个基因的拷贝数增加，1个基因的拷贝数减少，表明这些变异可能是性状变化的来源，因为CNV通常会导致基因表达及相应性状的改变【46,47】。

重测序大豆中的结构变异（SVs）与重要性状相关

基于深度重测序和更多样本的全基因组关联研究（GWAS）能够更精确地解析SVs对重要性状的遗传影响。因此，我们对547个代表性样本进行了高基因组深度重测序（平均覆盖度为18.05倍，见补充表27），并在2019年至2021年中国的10个不同地点×年份环境中评估了31个性状，包括6个与产量相关的性状（株高、分枝数、主茎节数、单株荚数（PN）、百粒重（HSW）和底部豆荚高度）以及16个种子质量性状（油脂、蛋白质、苷元、染料木苷、大豆苷、总异黄酮、α-生育酚、δ-生育酚、γ-生育酚、总生育酚、棕榈酸、硬脂酸、油酸、亚油酸、亚麻酸和不饱和脂肪酸含量），并在单一环境中评估了9个植物学性状（花色、茸毛、种皮颜色、种脐颜色、种子光泽、生长习性、株型、豆荚开裂和植物倒伏，见补充图4）。

通过NDD2基因组作为参考，我们从这些样本中探测到749,714个插入/缺失SVs（见补充表28和29）。群体结构和遗传关系分析表明没有群体偏差，使该样本集适合用于GWAS（见补充图5）。基于这些SVs、各环境中的31个性状均值、最佳线性无偏预测（BLUP）值以及十个环境中的总均值，进行了GWAS分析。总共发现14,237个非冗余SVs与性状显著关联，其中4,458个与产量相关性状、6,552个与种子质量性状、3,333个与植物学性状相关（见补充数据4和补充图6）。值得注意的是，我们鉴定了6,013个同时与BLUP和总均值匹配的非冗余SVs，并发现其中1,043个与产量相关性状，4,970个与种子质量性状相关（见扩展数据图4-10、补充数据5-6和补充图7-21）。此外，我们发现这些6,013个非冗余SVs的69.73%的关联信号与SNP-GWAS的关联信号一致（见补充表30），表明SV-GWAS可能提供额外的关联信号，类似于InDel-GWAS【48】。我们发现989个SVs位于调控区域，772个SVs位于基因体内，这些SVs可能直接改变调控元件和基因的功能【12】。此外，我们发现上述三类性状中共享的GWAS基因频繁出现（见补充表31和32）。

这些SVs为大豆育种中基因型与表型关系的全面理解提供了依据。我们首先发现位于第6号染色体15.25–15.35 Mb区域的21个SV与重要的产量相关性状PN（单株荚数）相关（图5）。每个SV都能将547个样本分为两种纯合类型，且PN均值显著不同（见补充数据5）。特别是，当我们整体分析这21个SV时，能够将样本群体分为两种单倍型，其中Hap-Alt的PN较大（55.7），而Hap-Ref较小（47.3，图5）。此外，我们还比较了这21个相关SV连锁区域中36个基因在影响豆荚形成的相关器官（如花、未成熟豆荚和茎）中的表达差异，比较了NDD2（Hap-Ref，约50个PN）和Yudou22（Hap-Alt，约71个PN）的差异，发现候选基因Glyma.NDD2.06G174200（钾/氢离子逆向转运蛋白）在Yudou22的花中表现出特异性高表达，而在Yudou22和NDD2的其他相关器官中未表达（图5）。这一现象类似于拟南芥中钾/氢离子逆向转运蛋白基因的突变，影响花粉管的伸长，进而导致卵细胞无法受精，种子无法形成【49,50】。由于正常的授粉和卵细胞受精是有效豆荚形成的必要条件，我们推测Glyma.NDD2.06G174200在决定PN方面具有重要功能。因此，我们推断这些21个SV可能是不同大豆样本中影响PN的遗传变异。

a. 与PN（单株荚数）关联的第6号染色体峰值区域的局部曼哈顿图和连锁不平衡（LD）热图。统计分析使用了双尾Wald检验。红色虚线表示候选基因组区域。 b. 基于重测序群体中21个SV的PN箱线图。在箱线图中，中心线表示中位数，箱线表示上四分位数和下四分位数，须线表示数据范围；n表示具有相同基因型的样本数量。差异显著性通过双尾t检验分析。 c. NDD2和Yudou22之间21个SV区域附近180 kb范围内基因的差异表达情况。红色标记的是候选基因Glyma.NDD2.06G174200。Chr表示染色体。

此外，我们首次在第5号染色体上鉴定了与六个种子质量性状相关的基因组区域，该区域有一个强信号，包含181 kb内的137个SV。在与这些SV相关的连锁区域中有25个候选基因，我们关注了编码FT和TFL1母本蛋白（MFT）的基因Glyma.NDD2.05G269100，该基因受到位于基因体内或调控区域的12个SV的影响（图6）。MFT通常被认为是FT和TFL1的祖先【51】。与FT和TFL1在分别正向和负向调控植物开花方面的明确功能不同，MFT的功能尚不清楚【52,53】，特别是本研究中多个与种子质量性状相关的MFT功能。因此，我们将此基因命名为GmMQT（multiple quality traits），并发现其在大豆种子和豆荚中的表达远高于叶片、根瘤和根毛中（见补充图22），表明它在种子发育和成分形成中的作用。GmMQT中的12个相关SV将547个样本分为两种单倍型（Hap-Ref和Hap-Alt）。Hap-Ref的油脂含量（19.05%）显著高于Hap-Alt（17.11%），γ-生育酚和生育酚含量也表现出一致的趋势。相反，Hap-Alt的异黄酮含量（2313.50 μg/g）显著高于Hap-Ref（1741.91 μg/g），其两个成分（染料木苷和大豆苷）也呈现出类似趋势。我们观察到通过毛状根转化的GmMQT在低异黄酮大豆品种Ji11B9中过表达，导致染料木苷、大豆苷和总异黄酮含量显著高于野生型，分别增加了1.87和2.93倍的染料木苷、1.93和0.76倍的大豆苷以及1.91和2.64倍的总异黄酮。进一步实验显示，十个Hap-Alt样本中的GmMQT表达和异黄酮含量显著高于十个Hap-Ref样本（见补充图23）。在相同条件下，带有pGmMQTAlt::GUS（β-葡糖醛酸酶）的转基因毛状根比带有pGmMQTRef::GUS的根显示出更强的GUS染色和显著更高的GUS表达，表明上游的两个SV对基因表达的影响（见补充图23）。此外，两株过表达GmMQT的拟南芥转基因株系的油脂和生育酚含量显著高于野生型，种子的油脂含量分别增加了6.71%和4.63%，生育酚含量增加了51.89%和47.74%（图6），这与之前的研究发现一致，即通过SNP-GWAS鉴定的大豆GmST05（属于MFT，我们的研究将其命名为GmMQT）影响了油脂含量【53】。此外，利用两个同源CRISPR-Cas9基因敲除的大豆株系（CR#1和CR#2）的GmST05【53】，我们发现这两株CR株系的种子中异黄酮含量极显著下降（分别减少94%和83%），相比于野生型（见补充图24）。因此，我们推测通过SV-GWAS鉴定的GmMQT可以影响更多种子质量性状，包括首次在本研究中确认的异黄酮含量以及之前研究中的油脂含量【53】。

a. 第5号染色体上与多种种子质量性状（油脂、异黄酮和生育酚）相关的峰值区域的局部曼哈顿图和LD热图。统计分析使用了双尾Wald检验。箭头指向多效基因GmMQT（Glyma.NDD2.05G269100）。红色虚线表示候选基因组区域。 b. GmMQT的基因结构。蓝色矩形表示外显子，灰色矩形和黄色箭头分别表示5′ UTR和3′ UTR。垂直线表示相关的SV。 c, e. 基于12个相关SV的种子油脂含量(c)和总生育酚含量(e)的箱线图。在箱线图中，中心线表示中位数，箱线表示上四分位数和下四分位数，须线表示数据范围；n表示具有相同基因型的样本数量。差异显著性通过双尾t检验分析。 d, g. GmMQT的过表达（OE1和OE2）显著提高了转基因拟南芥的种子油脂含量(d)和总生育酚含量(g)。WT表示野生型。差异显著性通过双尾t检验分析。数据以平均值±标准误差（n=3）表示。 f. 转基因拟南芥中GmMQT的PCR扩增。M-1表示DNA标记DL5000。 h. 基于12个相关SV的种子总异黄酮含量的箱线图。在箱线图中，中心线表示中位数，箱线表示上四分位数和下四分位数，须线表示数据范围；n表示具有相同基因型的样本数量。差异显著性通过双尾t检验分析。 i. 转基因毛状根中观察到的绿色荧光，显示GmMQT的表达。比例尺为500 μm。 j. 转基因毛状根中GmMQT的PCR扩增。M-2表示DNA标记DL2000。 k. qRT-PCR分析转基因毛状根中GmMQT的表达。差异显著性通过双尾t检验分析。GmActin11（Glyma.18G290800）作为内参基因。数据以三次重复的平均值±标准误差表示。 l. GmMQT的过表达显著提高了毛状根中的种子总异黄酮含量。差异显著性通过双尾t检验分析。数据以平均值±标准误差（n=3）表示。

此外，我们在第11号染色体上鉴定了种子大豆苷含量的最高关联峰，该区域包含1.04 Mb内的363个SV。在之前的研究中，几个QTL被定位到此区域【54】，但因果基因尚不清楚。我们重点关注了候选基因体内的相关SV，发现未知基因Glyma.NDD2.11G106800（命名为GmSGI，seed glycitin increasing）中存在一个内含子插入。在参考基因型的384个样本中，糖苷含量（217.62 μg/g）显著高于161个替代基因型样本（178.11 μg/g），表明GmSGI-SV与糖苷含量之间存在密切关系。此外，通过比较泛基因组中SoyC04、C06、C07、C08和C11的基因组序列【11】，验证了该内含子SV，C04、C06、C08和C11与NDD2（参考基因型）相同，而C07（替代基因型）不同于NDD2。参考品种种子中GmSGI的表达水平与其糖苷含量大致呈现一致的增加趋势，与替代基因型相比（见补充图25）。我们从高糖苷样本Wuxing4中克隆了开放阅读框（ORF）。在低糖苷样本Ji11B9中，两个通过毛状根转化过表达GmSGI的样本糖苷含量显著高于野生型，分别增加了45.4%和21.4%。此外，Hap-Ref样本中GmSGI表达和糖苷含量显著高于Hap-Alt样本（见补充图26），表明SV影响了基因表达进而影响糖苷含量。因此，我们得出结论，GmSGI是大豆种子糖苷含量的一个新鉴定基因。

a. 第11号染色体上与大豆苷含量相关的峰值区域的局部曼哈顿图和LD热图。统计分析使用了双尾Wald检验。箭头指向基因GmSGI（Glyma.NDD2.11G106800）。红色虚线表示候选基因组区域。 b. GmSGI的基因结构。蓝色矩形和箭头分别表示外显子和3′ UTR。红色垂直线表示GmSGI中的插入位点。 c. 基于相关SV的大豆苷含量的箱线图。在箱线图中，中心线表示中位数，箱线表示上四分位数和下四分位数，须线表示数据范围；n表示具有相同基因型的样本数量。差异显著性通过双尾t检验分析。 d. 在转基因大豆毛状根中观察到的绿色荧光，显示GmSGI的表达。比例尺为500 μm。 e. 转基因毛状根中GmSGI的PCR扩增。M表示DNA标记DL2000。 f. 通过qRT-PCR分析转基因毛状根中GmSGI的表达。差异显著性通过双尾t检验分析。GmActin11（Glyma.18G290800）作为内参基因。数据以平均值±标准误差（n=3）表示。 g. GmSGI的过表达显著提高了大豆毛状根中的大豆苷含量。差异显著性通过双尾t检验分析。数据以平均值±标准误差（n=3）表示。

讨论

在本研究中，我们组装了现代品种NDD2的高质量参考基因组，表现出更高的准确性、完整性和连续性，具有更高的contig N50值、最少的缺口以及清晰定义的着丝粒和端粒（补充表33和图1），为进一步生物学和遗传学研究挖掘更多变异提供了独特的基因组资源。例如，使用NDD2作为参考基因组，对547个重测序样本鉴定了64,529个SV和1,332个相关基因，再加上由于不同发布年代和来源国导致的不同遗传成分，补充了以Williams82为参考的资源（补充表34）。

越来越多的研究揭示了SV在作物改良中的关键作用【12,13,55,56】，然而SV在大豆改良中的作用，特别是对产量和种子质量性状的影响，尚不清楚。我们发现了大规模的SV和SV-基因对，它们可能影响大豆的产量和种子质量性状，以及对多样化环境的适应性，并在NDD2中识别了许多有利的等位基因，如DEL1815、DEL238和INV05（图2和图3，扩展数据图2和补充数据7）。此外，我们发现了许多SV-GWAS中检测到但SNP-GWAS中未检测到的位点，其中253个SV与产量相关性状相关，233个影响种子质量相关性状（补充数据8）。这些发现促进了我们对SV在大豆育种中重要性状基因型-表型关系中的作用的理解。

本研究中的约95%的重测序样本来自黄淮海大豆种植区，由于多样的地理和气候因素，这里栽培大豆具有丰富的遗传多样性【38】。由于约70%的样本不同于之前的研究，且具有最高的重测序深度，我们鉴定了更多准确的SVs用于大豆改良，包括使用BLUP值和总均值同时与22个重要性状相关的许多SV。除了与单株荚数（PN）、种子油脂、异黄酮和生育酚含量相关的SVs外，我们还发现了与种子蛋白质含量（SPC）相关的61个SV和与百粒重（HSW）相关的154个SV（补充数据4）。此外，使用额外的429个样本对三个随机选择的SV在GmMQT相关基因中的测定进一步证明了SV与重要性状之间的密切关系，NDD2特有的DEL238也显示出在该群体中SPC和HSW的显著增加（补充图27–30和补充数据9）。我们新鉴定的SV不仅更新了相关性状遗传基础解析的分子信息，还为大豆育种中的基因组辅助选择提供了有用的分子靶标。

方法

基因组组装和重测序材料

选择了现代大豆品种NDD2进行基因组组装，因为它在农艺、产量相关性状和种子质量性状方面具有多项优良特性（补充表1）。来自中国和美国的共547份大豆种质资源（61个地方品种和486个改良品种）进行了重测序（补充表35）。

NDD2和547份大豆种质资源的采样

对于NDD2基因组组装，在幼苗期采集NDD2的新鲜幼叶，而基因注释则采集根、茎、叶、花和豆荚组织。对于547份大豆种质资源的重测序，在幼苗期从每个样本的单株上采集新鲜叶片。

两种长读段测序

从NDD2叶片中提取基因组DNA，采用CTAB（十六烷基三乙基溴化铵）方法按照制造商提供的标准操作程序进行。测试DNA质量，确保DNA适用于所有类型的基因组文库。随后，我们构建了两种DNA文库，在PacBio Sequel II平台（Pacific Biosciences）和PromethION平台（Oxford Nanopore Technologies）上进行长读段测序。对于PacBio数据，使用默认参数过滤子读段；对于Nanopore数据，移除质量低于Q7的低质量读段。

短双端测序和质量控制

首先，将约1.5 μg的基因组DNA通过超声波处理成350 bp的短插入片段。其次，这些DNA片段经过末端抛光和全长接头的加尾连接后进行PCR扩增。第三，使用AMPure XP纯化PCR产物，并通过Agilent 2100生物分析仪检查DNA文库的大小分布，并使用实时PCR进行定量。最后，在Illumina NovaSeq平台上对DNA文库进行测序，生成原始测序数据。通过过滤接头污染和低Phred质量的低质量读段来清理数据。

Hi-C文库构建和测序

Hi-C文库由NDD2的叶片细胞构建。首先，用甲醛固定细胞并裂解。其次，用DpnII酶切交联的DNA，黏性末端进行生物素化，并进行近端连接以形成嵌合连接。第三，物理剪切后富集300–500 bp的DNA片段。表示原始交联长距离物理相互作用的嵌合片段被处理为双端测序文库。最后，双端读段在Illumina NovaSeq平台上进行测序。使用fastp软件【57】过滤原始测序数据中的低质量读段。

Bionano光学图谱数据的生成

从NDD2的新鲜幼叶中提取高分子量DNA。按照标准Bionano协议进行标记后，用Bionano Irys系统对标记的DNA分子进行成像。将原始图像数据转换为BNX文件，并使用Bionano Solve软件包（v.3.5.1）中的AutoDetect工具将基础标记和DNA长度信息进行转换（https://bionanogenomics.com/support/software-downloads/）。过滤分子长度和标记密度后生成NDD2的光学图谱数据。

新的组装流程

我们使用了之前改进的组装流程【20】来生成完整的NDD2基因组。该流程包括以下四个主要步骤，详述如下。

组装并校正的contig

我们使用高质量的Nanopore序列，通过NextDenovo软件包（v.2.4.0；https://github.com/Nextomics/NextDenovo）中的“correct-then-assemble”策略，使用参数“read_cutoff=1k, seed_cutoff=32k, blocksize=3g”对NDD2基因组进行组装。随后，使用NextPolish软件包（v.1.3.1）【58】与Illumina双端读段和Nanopore长读段结合使用“best”算法模块对初始contig进行校正。基于上述分析策略，我们获得了ContigV1基因组组装。

通过Hi-C交互对contig进行聚类

Hi-C读段对使用Bowtie2软件【59】的单端模型映射到ContigV1。我们获得了唯一映射的对，并使用HiCUP流程（v.0.8.0）【60】丢弃了无效的自连接和未连接片段。最后，我们使用有效的交互对计算所有contig之间的连接频率，并使用聚合层次聚类算法对contig进行聚类。基于Hi-C信号密度建议的链接，我们将链接的contig聚类，假设它们位于同一条同源染色体上，并设置了分区数。

PacBio读段分类以进行局部组装

我们使用Minimap2软件【61】将PacBio子读段重新对齐到ContigV1。去除次优对齐的读段后，提取每个contig组的映射读段。随后，我们对每个分类映射读段进行局部组装。与全局组装相比，该策略避免了组装过程中构建串图时基因组重复序列引起的错误重叠关系【62】。最后，我们再次对所有contig进行校正。根据“seq_stat”命令的结果，为每个局部组装设置合适的“seed_cutoff”参数。此步骤使用的组装和校正方法与步骤1中提到的相似。

contig锚定到染色体上

使用Bionano Solve软件包（v.3.5.1），使用参数“-B 1 -N 1”，我们将Bionano数据映射到contig级组装上，生成混合框架级组装。然后，我们将Hi-C数据映射到混合框架级组装上，从而生成连接信息。接下来，使用ALLHiC算法【63】对染色体级基因组进行锚定。最后，在由Juicebox（v.1.22；https://github.com/aidenlab/Juicebox）生成的可视化界面中手动调整具有明显离散染色质交互模式的定位和方向错误。

基因组评估

采用多种方法评估组装基因组的质量。具体来说，使用Merqury软件，通过无参考和k-mer方法评估基因组组装的质量和完整性【21】。我们进行了BUSCO分析【64】，通过搜索1,614个胚基植物BUSCO（v.5.2.1）来评估基因组的完整性，还使用BWA软件【65】将Illumina双端序列比对到组装基因组，以计算重比对率和覆盖深度。

基因组重复序列注释

通过同源性搜索和从头预测相结合的方法鉴定重复元件。对于基于同源性的预测，使用RepeatMasker【66】和RepeatProteinMask，在Repbase TE库中进行搜索，采用默认参数。对于从头预测，首先基于LTR FINDER【67】、PILER【68】和RepeatScout【69】的结果构建参考重复库，采用默认参数，然后使用RepeatMasker进行搜索。此外，使用Tandem Repeats Finder【70】鉴定基因组中的串联重复序列，参数设置为“2 7 7 80 10 50 2000 -d –h”。

基因结构和功能注释

基因预测通过同源性和从头预测相结合的方法进行，并辅以转录证据。

同源性方法

利用拟南芥、Williams82（G. max）、W05（Glycine soja）、SoyZH13（G. max）、NDM8（Gossypium hirsutum）、Medicago truncatula、水稻、菜豆和番茄的蛋白质组，通过TBLASTN【71】搜索目标基因组。然后使用SOLAR程序连接BLAST命中的局部比对结果，并使用GeneWise【72】定义每个BLAST命中中对应基因组区域的精确基因结构。

RNA测序数据

将RNA测序数据对齐到基因组，使用Tophat【73】鉴定假定的外显子区域和剪接位点。然后，使用Cufflinks【74】从比对的读段中组装基因模型（Cufflinks-set）。

从头预测方法

将NDD2基因组上的五种NDD2组织RNA测序数据的组装转录本使用PASA【75】程序进行对齐。通过PASA将转录本比对组装为基因结构模型，作为Augustus【76】、SNAP【77】和GlimmerHMM【78】的训练集。基于这些训练集，我们使用Augustus、GlimmerHMM和SNAP对重复掩蔽的基因组进行从头编码区预测。此外，GeneID【79】和GeneScan【80】用于直接预测重复掩蔽基因组中的基因模型。

在上述方法完成后，使用EVidenceModeler整合生成的基因模型。每种证据的权重设置如下：Homology-set = Cufflinks-set > Augustus > GeneID = SNAP = GlimmerHMM = GeneScan。最后，通过PASA2更新基因模型，以生成未翻译区和可变剪接变异信息。

随后，通过搜索已知数据库中的功能基序、结构域和可能的生物过程，对蛋白编码基因（PCG）进行功能注释，包括SwissProt【81】、Pfam【82】、NCBI的非冗余蛋白序列数据库（NR）、基因本体（Gene Ontology）【83】和京都基因与基因组百科全书（KEGG）【84】。

大豆基于图的基因组构建

我们使用NDD2和29个已发布的大豆基因组【1,7,9,11】，构建了包含大豆基因组SV图谱的基于图的基因组。首先，使用MUMmer软件包（v4.0.0）中的nucmer程序，将29个发布的大豆基因组比对到NDD2参考基因组，参数为“-c 1000 -maxgap = 500”。然后，应用delta过滤器程序，通过一对一比对过滤比对块。长度超过1,000 bp的比对块用于进一步检测。最后，使用SVMU流程鉴定插入、缺失、倒位和CNV，并应用SyRI软件【86】鉴定易位区域。

群体基因组测序和SNP识别

对547份大豆样本中的每个样本提取约1.5 µg基因组DNA，制备DNA文库，并在Illumina NovaSeq平台上生成测序数据。使用fastp软件【57】过滤低质量读段后，将高质量数据使用BWA（0.7.17-r1188）软件映射到NDD2基因组上，命令为“mem -t 4 -k 32 -M”【87】。然后，使用Picard软件包中的“MarkDuplicates”命令标记重复的BAM文件【v2.18.15】。

随后，我们根据最佳实践与HaplotypeCaller方法结合使用Genome Analysis Toolkit（v4.1.2.0）【88】，对个体gVCF进行调用，然后通过命令“GenomicsDBImport”和“GenotypeGVCFs”合并所有gVCFs进行群体SNP调用。为获得可信的群体SNP集，我们执行以下筛选过程：(1) 使用硬过滤命令“VariantFiltration”排除潜在的假阳性变异调用，参数为“--filterExpression 'QD < 2.0 || MQ < 40.0 || FS > 60.0 || SOR > 3.0 || MQRankSum < −12.5 || ReadPosRankSum < −8.0'”；(2) 使用Hardy–Weinberg平衡P值≥0.01筛选双等位基因变异；(3) 当群体内缺少变异的样本比例超过20%且次要等位基因频率（MAF）< 0.01时，过滤并移除变异。最终，鉴定出高质量SNP用于后续分析。根据NDD2参考基因组，使用ANNOVAR【89】对SNP进行注释。

群体SV识别

使用两种SV识别方法来鉴定547份大豆样本的基因组SV。首先，我们使用基于图的基因组中的基因组SV图谱，并使用vg工具【90】对群体SV进行基因分型。然后将Illumina短读段映射到基于图的基因组上，vg工具以默认参数运行。vg工具能够准确地对基于图的基因组中的每个样本的变异进行基因分型，包括“1/1”（纯合变异基因型）、“0/1”（杂合变异基因型）、“0/0”（纯合参考基因型）和“-/ -”（缺失）。基于此，我们对所有547份样本中的所有变异进行基因分型，并生成547份样本的SV图谱矩阵。第二，每个样本使用BWA软件【88】对NDD2基因组进行比对，生成的所有547个样本的BAM文件用于通过LUMPY软件包【91】进行群体SV检测。

为进一步评估每个样本中SV基因分型结果的可靠性，我们首先拆分出547个样本的SV。其次，提取SV边界区域周围的已映射短读段，并使用Velvet【92】进行局部组装，以获得每个样本的SV相关contig序列。最后，我们将此SV相关的contig序列映射到NDD2序列中，以获得比对结果。我们计算解析该比对，以确定是否存在支持SV的证据。最后，我们将所有个体的调用重新合并为非冗余集，并确保每个调用至少存在于十个样本中。

连锁不平衡分析

使用Haploview软件计算高质量变异对之间的连锁不平衡系数（r²），参数设置如下：maxdistance 500, MAF > 0.05, miss ≤ 0.2。然后，使用这些结果估计LD衰减。所有大豆样本的LD衰减在最大值的一半时为189 kb。

农艺、产量和种子质量性状的表型测定

547个重测序样本分别在2019年、2020年和2021年10个环境中种植，环境标记为E1–E10（补充表36）。每个实验按照随机完全区组设计种植，三次重复。种植行长度设置为2.5米，行距和株距分别设置为0.5米和0.1米。

在单一环境中研究了9个植物学性状（花色、茸毛颜色、种皮颜色、种脐颜色、种子光泽、株型、生长习性、豆荚开裂和植株倒伏），在10个环境中

测定了6个产量相关性状（株高、分枝数、主茎节数、单株荚数、百粒重和底部豆荚高度）和16个种子质量性状（油脂、蛋白质、染料木苷、大豆苷、甘草苷、总异黄酮、α-生育酚、δ-生育酚、γ-生育酚、总生育酚、棕榈酸、硬脂酸、油酸、亚油酸、亚麻酸和不饱和脂肪酸含量）。

对于百粒重和16个种子质量性状，种子样本来自同一行中五株植物的合并样本；而对于5个产量相关性状，数据独立来源于每次重复中每个样本的5株植物。9个植物学性状和6个产量相关性状的评估方法见补充表37【93】。SPC（种子蛋白质含量）通过凯氏定氮法测定（GB 5009.5-2016），种子油脂含量通过核磁共振法测定。种子异黄酮及其成分含量通过高效液相色谱法（GB/T 26625-2011）测定，脂肪酸通过气相色谱法检测，生育酚及其成分含量通过高效液相色谱法测定。

通过Microsoft Excel 2017软件计算单一环境和10个环境中的产量和种子质量性状的均值，并通过R的lme4包计算这些性状的BLUP值。这些性状在单一环境和10个环境中的均值以及BLUP值用于GWAS分析。

全基因组关联分析 (GWAS)

我们使用高质量的SV集合，通过应用全基因组高效混合模型关联分析 (GEMMA) 软件【12】对31个性状进行了GWAS分析。GEMMA分析中使用了以下方程：y = Xα + Sβ + Kμ + e，其中y表示表型，α和β表示固定效应（标记和非标记效应），μ表示未知的随机效应，X、S和K分别为α、β和μ的发生矩阵，e为随机残差效应向量。在本研究中，使用前三个主成分构建S矩阵以进行群体结构校正。使用简单匹配系数构建K矩阵。位于GWAS信号附近约180 kb的基因（补充图31）被视为候选基因。

在拟南芥中过表达GmMQT

使用在样本JiNF58中设计的引物（补充表38），通过PCR扩增GmMQT的ORF。扩增片段插入二元载体pCHAC的SalI和KpnI限制性内切酶位点。通过花浸法将该构建体转化为拟南芥Colombia-0生态型【94】。转基因植株在含50 μg/ml潮霉素（H370; PhytoTechnology Laboratories）的1/2 Murashige和Skoog (MS) 培养基中筛选（M519; PhytoTechnology Laboratories）。通过基因组PCR进一步确认转基因植株中的GmMQT存在（补充表38）。T3转基因植株用于分析种子油脂、γ-生育酚和生育酚含量（分析方法见补充表37），实验进行了三次生物重复。

转基因大豆复合植株的生成

为大豆毛状根转化制作GmMQT和GmSGI过表达构建体，分别从样本Zhongdou32和Wuxing4中通过PCR扩增基因的全长ORF（补充表38）。扩增的ORF克隆到带有报告基因GFP的二元载体pB7WG2D中，所得质粒通过冻融法引入根癌农杆菌K599菌株【95】。通过A. rhizogenes介导的转化方法获得转基因大豆毛状根【96】。大豆样本Ji11B9的种子在28°C下萌发3天，胚轴被K599菌株感染。毛状根长出并生长12天后，通过荧光显微镜（Olympus BX51）检测转基因毛状根的GFP荧光。

分别提取带有GmMQT和GmSGI的转基因大豆毛状根的总RNA，合成cDNA。在CFX96TM实时系统（Bio-Rad）上使用SYBR Premix Ex Taq（Takara Bio）通过qRT-PCR测定GmMQT和GmSGI的表达水平（补充表38）。使用管家基因GmActin11（Glyma.18G290800）作为内参，并使用2−ΔΔCt方法【96】标准化每个转化体中GmMQT或GmSGI的表达倍数变化。使用GFP阳性毛状根进行染料木苷、大豆苷、甘草苷和总异黄酮含量检测【97,98,99】，实验进行了三次重复。