The complete reference genome for grapevine (Vitis vinifera L.) genetics and breeding
葡萄(Vitis vinifera L.)遗传学和育种的完整参考基因组
摘要
葡萄是全球最具经济重要性的作物之一。然而,以往版本的葡萄参考基因组通常由成千上万个片段组成,缺失着丝粒和端粒,限制了重复序列、着丝粒和端粒区域的可及性,以及这些区域中重要农艺性状的遗传研究。在此,我们利用PacBio HiFi长读长序列为品种PN40024组装了一个从端粒到端粒(T2T)无间隙的参考基因组。该T2T参考基因组(PN_T2T)比12X.v0版本长69 Mb,且鉴定出更多的9018个基因。我们注释了67%的重复序列、19个着丝粒和36个端粒,并将之前版本的基因注释整合到PN_T2T组装中。我们检测到总共377个基因簇,这些基因簇与复杂性状如香气和抗病性相关。尽管PN40024源自九代自交,我们仍发现了九个与生物过程(如氧化还原过程和蛋白质磷酸化)相关的杂合位点基因组热点。因此,全面注释的完整参考基因组构成了葡萄遗传研究和育种项目的重要资源。
介绍
自2000年首个人类基因组发布以来,各种物种的数百个参考基因组相继被组装【1-3】。参考基因组对于生物学和遗传学研究至关重要,因此获得高质量的基因组一直是研究的目标。尽管如此,由于基因组中高度重复的序列,特别是端粒、着丝粒和核糖体DNA(rDNA)这三大代表性区域,许多片段仍然缺失【3-5】。
着丝粒承载CENPA/CENH3变异体核小体,是纺锤体微管形成和附着的位置,在细胞分裂过程中起着关键作用。它由α卫星DNA组成,α卫星由被称为高级重复序列(HORs)的单体重复序列组成,单体长度在100至200 bp之间【6-8】。尽管着丝粒的功能在不同物种中保持一致,但其结构和序列在物种内和物种间可以快速变化,并且在不同物种中可以观察到多样的组织形式。然而,着丝粒在基因组内表现出协同进化【7, 9-11】。目前,研究人员对着丝粒仍知之甚少。
端粒也大多未知。端粒由位于真核生物染色体末端的串联重复的相对保守的微卫星序列组成【12, 13】。端粒在细胞分裂过程中保护染色体末端序列【14-17】。核糖体DNA(rDNA)是基因组中最丰富的重复元素之一,在核糖体形成中起关键作用,推动细胞生长和增殖【18-20】。
由于之前组装的基因组中缺失的信息,过去二十年中对着丝粒、端粒和rDNA的研究极为有限。幸运的是,得益于测序技术和计算算法的改进,基因组组装迎来了一个新的时代:端粒到端粒(T2T)测序【21】。与碎片化的基因组相比,T2T基因组几乎没有间隙。它基于第三代测序平台,包括PacBio高保真长读(HiFi)、超长Oxford Nanopore Technologies(ONT)和Hi-C数据。此外,T2T基因组几乎包含了端粒、着丝粒和rDNA区域的完整信息【22, 23】。有望的是,T2T基因组使我们能够进入这些区域,了解这些区域的结构以及这些区域中基因的功能。自2020年发布首个人类X染色体的完整序列以来,T2T组装迅速成为研究热点【22, 23】。在植物中,首个T2T基因组于2021年在拟南芥中报道【7, 24】。目前,几种物种的T2T基因组组装已被获得,如水稻、香蕉和西瓜,使研究人员对基因组结构和功能及其与作物育种性状的关系产生了浓厚兴趣【25-28】。
葡萄(Vitis vinifera ssp. vinifera)是一种起源于近东的果树,是全球最广泛栽培和经济价值最高的作物之一【29】。驯化的葡萄通常具有高度杂合的基因组【30】,这极大地阻碍了高质量基因组的获取。例如,‘霞多丽’基因组中约有15%的基因是半合子的【31】。幸运的是,PN40024基因型,一种由‘Helfensteiner’自交产生的高度纯合的品种【31】,成为葡萄的参考基因组,首次于2007年获得(8X),是首个被测序的水果作物【32】。随后,几个更新版本相继发布:2017年的12X.v2版本及其升级注释VCost.v3,2021年的PN40024.v4.1版本【33】。葡萄基因参考目录现在包括所有注释版本的完整对应关系【34】。此外,近年来生产了各种葡萄品种的碎片化基因组组装,如‘黑科林斯’【35】、‘品丽珠’【36, 37】、‘赤霞珠’【37-39】、‘佳美娜’【40】、‘霞多丽’【30, 41】、‘梅洛’【35】和‘内比奥罗’【42】。作为果树中具有代表性的双子叶植物,葡萄的高质量基因组将极大地促进Vitis和双子叶物种的基因功能、遗传结构和进化研究。
尽管有大量的葡萄基因组序列可用,但这些基因组组装在重复区域、着丝粒和端粒方面仍然不完整。在此,我们生成了PN40024参考的T2T级别无间隙葡萄基因组,并进行了四项主要分析。第三代测序和组装技术在高保真长读的应用有助于无间隙基因组组装【43, 44】。因此,我们的第一个问题是,是否可以使用这些新的测序和组装方法完成葡萄参考基因组。其次,由于对着丝粒、端粒和rDNA的研究长期被忽视,我们基于组装的无间隙葡萄基因组分析了这些区域的特征、结构和分布。第三,基于T2T基因组,我们改进了高度重复区域中的转座子(TEs)和基因注释,这可以进一步提高我们对其生物功能,特别是基因簇功能的理解。最后,尽管PN40024基因组几乎完全纯合【32】,但在九代自交后,仍有一些位点保持杂合。研究这些杂合位点的基因组分布和遗传效应是值得的。
结果
无间隙的端粒到端粒葡萄参考基因组
我们使用源自‘Helfensteiner’的高度纯合自交系PN40024进行T2T基因组组装。总共由PacBio平台生成了21 Gb(21,024,461,524 bp,约42×覆盖)的HiFi读数。对于初步组装,使用hifiasm组装HiFi读数。然后,我们使用MUMmer和12X.v0基因组版本(V. vinifera基因组组装12X.v0;Vitis vinifera genome assembly 12X - NCBI - NLM)将38个片段排序到19条染色体中(图1)。在初始组装成片段后,仅剩一个间隙(补充数据图S1)。在用PN40024.v4的连续长读填补间隙后,最终生成了无间隙的PN_T2T基因组(494.87 Mb),比12X.v0(426.18 Mb,表1)长69 Mb,使用相同的统计方法。k-mer度量用于评估基因组的纯合度,估计为99.8%(补充数据图S2A-D)。BUSCO(基准通用单拷贝直系同源物)用于评估基因组的完整性;在基因组组装中发现了98.5%的核心保守植物基因(补充数据图S2E),比12X.v0(93.7%,表1)多出4.8%。
图1 葡萄参考基因组的T2T无间隙组装。(A) 基因组组装概述(12X.v0,右侧柱;PN_T2T,左侧柱)。染色体3和染色体5上的红色虚线框表示两个版本基因组组装之间的大型倒位差异。(B) (A)中染色体3上的红色虚线框区域的放大部分。(C) 显示染色体18上12X.v0和PN_T2T组装之间的1 Mb共线性区域的图。灰色条带连接对应的共线性区域,底部的红色框表示12X.v0中的缺口。(D) 在PN_T2T基因组中检测到的不同转座子家族的类型和百分比。
Table 1
Comparison of genomic features of 12X.v0, 12X.v2, PN40024.v4, and PN_T2T assemblies.
12X.v0 | 12X.v2 | PN40024.v4 | PN_T2T | |
---|---|---|---|---|
Total sequence length (bp) | 426 176 009 | 458 815 822 | 462 158 227 | 494 873 210 |
Number of chromosomes | 19 | 19 | 19 | 19 |
Contig N50 (bp) | 102 700 | 102 674 | 25 934 928 | |
Maximum length (bp) | 30 274 277 | 34 568 450 | 34 942 157 | 36 684 271 |
Number of gaps | 9429 | 5106 | 3391 | 0 |
Centromeres annotated | 19/ 19 | |||
Telomeres annotated | 36/38 | |||
Bases masked (bp) | 303 719 475 | 328 929 883 | ||
Retroelements (bp) | 217 819 122 | 241 027 616 | ||
LTR (bp) | 212 117 752 | 235 245 099 | ||
Number of genes | 28 516 | 41 182 | 35 256 | 37 534 |
Number of TEs | 942 096 | 935 783 | ||
BUSCO (%) | 93.70 | 97.70 | 98.20 | 98.50 |
12X.v0, Genome - NCBI - NLM. 12X.v2, https://urgi.versailles.inra.fr/files/Vini/Vitis%2012X.2%20annotations/12Xv2_grapevine_genome_assembly.fa.zip. PN40024.v4, https://grapedia.org/genomes/https://integrape.eu/resources/genes-genomes/genome-accessions.
图2 PN_T2T参考基因组中的重复序列注释。(A) 着丝粒和端粒预测的数据流。(B) 端粒、着丝粒和不同类型转座子的染色体分布。虚线表示预测的着丝粒中心位置。
与12X.v0基因组相比,我们的PN_T2T组装在多个指标上有显著改进。PN_T2T的contig N50长度比12X.v0高约250倍(25.93 Mb 对比 102 kb),12X.v0中的所有9429个缺口和PN40024.v4中的3391个缺口都在PN_T2T基因组中填补了(表1,补充数据表S1,图1A)。如图1C所示,12X.v0中的28个缺口在PN_T2T中得到了填补,最大的缺口位于染色体18的1 Mb共线性区域,长度为16,951 bp(图1C)。12X.v0中的方向错误也得到了纠正,例如与PN_T2T相比的倒位和易位(图1A,补充数据图S3)。例如,两个大型倒位分别位于染色体3的着丝粒周围和染色体5的末端,长度分别为4.9 Mb和1.9 Mb,在两个组装版本之间被观察到(图1A和B,图S8)。此外,PN_T2T基因组组装中检测到了19个着丝粒和38个端粒中的36个,除了染色体15上的一个端粒和染色体17上的一个端粒。总共注释了37,534个基因和41,064个转录本,其中24,526(86.01%)、27,696(78.83%)和27,717(78.75%)与较早版本PN40024.v2.1(Phytozome v13)和PN40024.v4.1(Grapevine Genomes - Grapedia)共享(补充数据表S2)。总共有5472个(14.58%)基因在任何三个版本中都没有对应。通过BUSCO分析评估了97.9%的完全组装基因,并在40,307个独特序列中的35,508个序列(88.1%)中检测到结构域,而PN40024.v4.1有38,364个独特序列,其中29,688个序列中检测到结构域(77.4%,补充数据表S2)。
基于RepeatModeler2构建的物种特异性pan-TE数据库,使用图2A所示的管道检测重复序列。最终,我们的无间隙葡萄基因组中有66.47%被标记为重复序列(图1D)。相比之下,使用相同管道在12X.v0基因组中识别出62.47%的重复序列(补充数据表S3)。在PN_T2T基因组中预测的重复序列中,最大部分为转座子(63.90%),总长度为316 Mb(在12X.v0中为59.96%和292 Mb)。转座子主要由长末端重复序列(LTR)类型组成,主要是Gypsy(20.22%)和Copia(19.67%)元件。总共检测到276个rDNA序列,占基因组的0.019%。
图3 PN_T2T基因组中着丝粒重复单元的示意图。(A) 整个基因组中不同重复单元长度的分布。图表上部显示不同重复单元拷贝数,图表下部显示不同重复单元在染色体中的百分比。(B) 每条染色体中107 bp重复单元拷贝的总长度。(C) 19条染色体中107 bp重复单元的比对。(D-H) 染色体16、17、3、14和18中不同重复单元的总长度。
端粒和着丝粒的识别
为了访问PN_T2T中的端粒和着丝粒区域,我们使用图2A所述的流程识别了端粒和着丝粒。对于端粒,我们检查了每条染色体两端的150 kb序列,端粒重复单元的长度设置为5至12 bp。最终检测到的端粒重复单元(TTTAGGG/CCCTAAA)在整个基因组中最为丰富,且存在于所有染色体中。Melters等人【11】和Castro等人【45】在葡萄中也报告了相同的端粒重复单元。我们进一步预测了PN_T2T基因组中38个端粒中的36个,除了染色体15和染色体17的短臂(图1A和2B,补充数据表S4)。其中,最长的端粒(31 kb)位于染色体8的短臂,有4479个重复,而最短的端粒(1260 bp)位于染色体7的长臂,仅有180个重复。
为了检测着丝粒区域,我们沿着基因组扫描了30至500 bp的候选重复序列。串联重复序列查找器(TRF)在PN_T2T基因组中发现了470个不同的重复单元。107 bp重复是整个基因组中最丰富的单元,有182,620.5(拷贝数≥2)次重复,占总基因组序列的约3.95%,其次是321 bp(2.45%)、214 bp(1.94%)和135 bp(1.05%)(图3A)。有趣的是,我们发现214 bp和321 bp重复单元的序列分别由两个和三个107 bp重复单元组成。转座子(TE)分析也支持了这个107 bp重复单元的着丝粒特征(图2)。因此,着丝粒主要基于107 bp重复单元识别,并定位于所有19条染色体上(图1A和2B,补充数据表S5)。如图3B所示,107 bp重复单元的总长度从1.4 kb到3.8 Mb不等,但107 bp重复单元的序列在染色体间高度保守(图3C)。除了染色体3、14和18,107 bp重复单元在所有染色体中最为丰富(图3D-H,补充数据表S6)。我们发现187 bp是染色体14中的主要重复单元,分布在整个染色体中,51 bp、56 bp、105 bp和107 bp重复单元在着丝粒中高度重叠和富集,IGV可视化显示染色体中有一个核心区域(补充数据图S4)。染色体3中的着丝粒重复单元是135 bp重复及其整数倍(270和405 bp)。对于染色体18,66 bp及其整数倍132 bp是主要重复单元(补充数据图S4)。
图4 着丝粒重复单元拷贝的特征和分布。(A) 整个基因组中基因、TEs和不同重复单元的分布。(B) IGV中可视化的染色体16上预测的着丝粒区域。(C) 着丝粒中捕获的基因的GO功能注释。MF,分子功能;CC,细胞成分;BP,生物过程。显著的富集P值:*P < 0.05,**P < 0.01。(D) 三角形显示每个单倍型内的序列相似性,按身份着色。
定位着丝粒重复序列
为了定位着丝粒重复序列,我们进一步检查了TEs与着丝粒之间的关系。LTR逆转座子或着丝粒逆转座子(CRs)通常与串联重复混合并富集在植物的着丝粒区域【46, 47】。我们发现(图4A),在综合基因组查看器(IGV)中查看染色体中特富集的巨大着丝粒串联重复序列时,基因和TE重复序列,如LTR(Gypsy和Copia)、DNA TE(MULE-MuDR)和RC(Helitron),在特殊区域的密度较低(图2,补充数据图S4)。然后,我们推断具有着丝粒重复序列和低TE密度的区域为着丝粒区域,通过逐一放大观察(补充数据图S4,补充数据表S5)。107 bp的模式是目标,与葡萄的着丝粒区域高度相关。然而,在染色体3、14和18上可能出现不同的重复单元和模式(图3F-H)。转座子的分布和着丝粒的分布显示特定的序列定义的重复超家族在不同程度上与着丝粒的接近度相关或不相关(图2B和4A),形成密度梯度,这些是主要的染色体级别的重复相关特征,推测反映了整体染色质结构(补充数据图S4)。
检测捕获的基因
为了检测捕获的基因,我们筛选了这些区域中所有与着丝粒区域高度相关的基因。有趣的是,我们发现了343个基因(补充数据表S7和S8)捕获在着丝粒中,其中179个基因通过BLASTP具有Uni-Prot ID。通过GO(基因本体论)功能注释,12个基因在蛋白质结合(分子功能,MF)中富集,例如参与乙烯、赤霉素和脱落酸信号通路的VviAMP1(Uni-Prot ID Q9M1S8)【48, 49】。此外,我们发现10个基因在细胞成分(CC)中富集于细胞质、线粒体和细胞浆,包括影响植物生长和发育的一般运输蛋白VviBIG(Uni-Prot ID Q9SRU2)【50】;催化线粒体克雷布斯循环相关酶活性的延胡索酸水解酶1 VviFUM1(Uni-Prot ID P93033)【51】;以及在雄配子体发育和花粉管与胚珠之间的相互作用中起关键作用的6-磷酸葡萄糖酸脱氢酶,脱羧2 VviPGD2(Uni-Prot ID Q9FWA3)【52】。此外,RNA修饰、蛋白质自磷酸化、DNA整合、DNA重组和光形态发生等生物过程(BP)相关术语也表现出富集(图4C)。
葡萄参考基因组中的基因簇
为了推断葡萄基因组中的基因簇,PN40024蛋白编码基因之间的蛋白质对蛋白质比对揭示了在基因组位置和功能方面丰富的重复结构。通过可视化全对全比对(补充数据图S5)可以看到突出且复杂的串联类似高相似性基因块。我们在葡萄参考基因组中发现了总共377个基因簇(补充数据表S9)。这些重复通常涉及局部重排,并可能扩展到包含几十到数百个基因的兆位(图5)。在染色体16(23-27 Mb)上,有599个富集域基因,主要包括WAKs(壁相关受体激酶聚糖结合)、PPR重复、富亮氨酸重复、ABC转运蛋白、整合酶域、肽酶家族、蛋白激酶和逆转录酶(图5A)。在染色体18(25-36 Mb)上,有1,237个基因富集在主要包括整合酶域、C-JID域、NB-ARC域、富亮氨酸重复、多铜氧化酶、逆转录酶、萜烯合酶和TIR的域中。我们的结果显示,许多强烈富集的结构域是植物抗病基因(R基因)的结构域的一部分,包括NB-ARC、TIR和由Colis数据库识别的结构。我们分析了41,064个PN_T2T PCG的结构域架构,并识别了3,381个可能的R基因。总的来说,这些R基因和葡萄中的基因簇为探索植物防御机制提供了巨大的机会。
自交九代后仍存在的杂合区域
基于PN_T2T基因组组装,我们从NCBI下载并分析了四个PN40024克隆的重测序数据【32, 53】。总共检测到244,215个SNPs,其中208,330个SNPs(85.3%)在所有四个样本中共享,而另外35,886个SNPs仅在一到三个样本中存在(图6A)。有趣的是,我们在染色体1、2、3、4、7、10、11和16上发现了九个杂合SNPs热点(图5A,补充数据图S6)。为了进一步研究高度杂合的区域,我们检查了前5%的杂合性窗口,并确定了总共九个大的连续片段(染色体1,1.1-1.3 Mb;染色体2,4.2-7.2 Mb;染色体3,9.4-9.9 Mb;染色体4,21.8-22.9 Mb;染色体7,15.3-26.2 Mb;染色体10,0.7-6.5 Mb,17.6-18.3 Mb;染色体11,7.1-7.8 Mb;染色体16,13.0-13.5 Mb)。对这些区域内基因的GO富集分析显示,最显著的富集术语是对缺水的反应、蛋白质磷酸化、细胞分裂、对氧化应激的反应和对盐胁迫的反应,这些与植物的关键生理活动密切相关(补充数据表S10和S11,图6C,补充数据图S7)。我们进一步在PN_T2T参考基因组上对这些九个杂合区域的热点进行了分相(补充数据文件2)。
图5 已识别基因簇的示意图。(A) 染色体16上22-27 Mb的基因簇。(B) 染色体18上25-36 Mb的基因簇。右侧的图表表示左侧橙色框内的区域。不同颜色表示不同的基因簇。None,Pfam数据库未识别出结构域的基因;Other,含有少量基因的其他基因簇。
图6 PN40024中杂合区域的特征。(A) 在所有四个PN40024样本中共享的杂合位点。灰色条表示着丝粒区域,橙色线表示存在于所有样本中的杂合位点。蓝色框突出显示了大的杂合片段。(B) 使用四个样本中不重叠的100 kb窗口计算PN40024基因组中的杂合性。(C) (A)中显示的杂合位点所含基因的GO富集分析。显著的富集P值:P < 0.05,**P < 0.001。
讨论
一个完整的参考基因组对于作物遗传学研究和育种目的至关重要。最新版本的PN40024.v4组装通过包括长读序列和收集标准注释来改进参考资源【31】。尽管如此,这些先前的版本仍然存在数千个缺口,并且缺乏重复区域、着丝粒和端粒,这些都限制了对这些区域内变异的访问。有时,这些无法到达的区域是重要农艺性状的数量性状位点(QTLs),例如染色体2上的果实颜色和性别决定【30, 54-56】以及染色体14上的抗病性【57, 58】。因此,一个完整的参考基因组具有揭示重要多基因农艺性状遗传力缺失的巨大潜力,从而提高葡萄育种的遗传增益。
越来越多的研究表明基因簇的重要功能,在PN_T2T中检测到总共377个基因簇。由于葡萄基因组在双子叶植物进化中的重要系统发育位置,葡萄基因组也广泛用于植物进化和比较基因组学研究【32】。T2T版本可以广泛用于植物进化基因组学,特别是重复序列、着丝粒和端粒。T2T无间隙参考基因组包含了先前版本的基因注释,并具有更准确的TE注释(占基因组的约67%),这将成为葡萄功能基因组学和育种的重要资源。
植物着丝粒的结构和背景
着丝粒区域长度从千碱基到十亿碱基不等,包括超过90%的串联重复序列【59】。着丝粒是基因组中最后一个未知的领域之一,因为先前的测序技术无法访问它。由于着丝粒区域高度重复的性质,组装通常会崩溃。我们为葡萄基因组的所有19条染色体组装并注释了着丝粒(图1)。大多数染色体有一个着丝粒,而其他染色体可能有多个着丝粒区域,即所谓的全着丝粒【60, 61】。在染色体16和18上,我们在多个区域发现了串联重复,而在其他染色体上仅检测到一个峰值(图2B),这表明着丝粒区域的结构可能更加复杂,需要进一步研究。
在PN40024葡萄参考基因组中,19条染色体上有三种主要的重复模式,表明不同的染色体进化历史(图3D-H)。在染色体3、14和18上,我们分别发现了135 bp、56 bp和66 bp的串联重复(补充数据图S4),而在其他染色体上,串联重复的主要单元是107 bp(图3D-H和4B和D)。每条葡萄染色体的着丝粒进化历史仍是一个悬而未决的问题,需要用所有Vitis基因组来解决。先前的比较基因组分析表明,着丝粒在染色体数量恒定的近缘物种中是保守的【9】。在染色体分裂和融合期间,着丝粒结构发生变化,当染色体数量在进化过程中变化时。麝香葡萄(Vitis rotundifolia)有20条染色体,其中染色体7和13与亚属Vitis染色体7共线,这与染色体融合事件有关【62】。在我们的葡萄参考基因组中,染色体7上只剩下一个着丝粒区域(图2B,补充数据图S4),这表明在Vitis属的进化过程中丢失了一个着丝粒。
着丝粒结构塑造了基因组内的内容、种群内的遗传多样性以及物种间的遗传分化。先前的种群遗传分析揭示了着丝粒区域中的遗传变异高度关联,其遗传多样性远低于染色体臂【63】。着丝粒捕获了数十到数千个高度与着丝粒串联重复序列相关联的基因。这些基因与着丝粒区域一起,作为超级基因发挥功能【64】。总的来说,我们在葡萄参考基因组的着丝粒区域中发现了343个捕获的基因(补充数据表S7)。有趣的是,这些基因主要参与乙烯、赤霉素和脱落酸信号通路【48, 49】。
几乎纯合基因型中的杂合热点
用于构建葡萄参考基因组的当前植物起源于“赫芬斯坦”品种,通过九代自交产生了99.8%的纯合基因组(补充数据图S2A-D)。剩余的杂合位点仍然引人关注,因为它们可能代表需要杂合性的热点,如果处于纯合状态会导致致命后果。因此,我们收集了在不同国际实验室维护的四个PN40024克隆的Illumina重测序读数。有趣的是,当映射到PN_T2T时,杂合SNP和结构变异(SVs)在特定区域中富集。总共发现208,330个在四个样本中共享的杂合SNPs,35,886个SNPs仅存在于一个至三个样本中。前者更可能是PN40024在九代自交后保留的原始变异,而后者可能是在不同实验室分布和组织培养期间产生的体细胞变异。有趣的是,我们发现常见变异的热点富集于中心生物过程,包括氧化还原过程和蛋白质磷酸化。染色体2上的热点还覆盖了性别决定QTL区域(图6),这使得性别决定基因的挖掘更加复杂【30, 56】,因为候选基因在旧版参考基因组中不存在。据报道,在果树的克隆繁殖过程中,这些杂合有害变异在基因组中积累【30, 65】。克隆过程隐藏了杂合的隐性有害变异,包括小的SNPs和插入缺失变异(indels)以及大的结构变异,在自交周期中暴露于致命或强隐性选择。葡萄育种中常见自交和外交衰退,因为在克隆繁殖过程中增加的隐藏杂合隐性有害变异在有性繁殖过程中暴露。
总体而言,并且仍然承认所有先前的测序工作,我们的工作代表了葡萄参考基因组的完整T2T序列的完成。这一组装,加上先前手动校对的注释,目前正在转移到PN_T2T中,应代表葡萄界的黄金标准。按照这一预测,T2T组装及其更新的注释可在葡萄基因组百科全书(GRAPEDIA;Home - Grapedia)下载,其中将与不同的应用程序接口一起使用,包括基因卡、转录组数据可视化和用于变异-基因表达-表型关联的软件。
材料和方法
样本采集和基因组测序
PN40024是通过连续自交步骤从‘Helfensteiner’品种中衍生出的近纯合系之一【31】,通过SSR标记测试估计其纯合度接近97%【32】。我们根据材料转移协议(MTA)从INRAE获得了这种自交材料,并将其移植到属于深圳中国农业科学院农业基因组研究所(AGIS)的温室中进行后续实验。从PN40024的幼叶和胚珠中快速冷冻至液氮中。使用DNeasy Plant Mini Kit(Qiagen)按照制造商说明分离基因组DNA和RNA。对于PacBio HiFi测序,在PacBio Sequel II平台上测序了两个单分子实时细胞,使用CCS(https://github.com/PacificBiosciences/ccs)生成了总计21 Gb的HiFi读长数据。对于RNA-seq,从总RNA中分离出的10 μg poly(A) mRNA用于制备每个样本的Illumina RNA-seq文库。这些文库然后按照制造商的说明使用Illumina HiSeq™ 2000系统进行测序。
端粒到端粒的基因组组装
最初,通过整合PacBio单分子实时长读序列组装PN40024基因组。由PacBio Sequel II平台生成的读数通过hifiasm(https://github.com/chhylp123/hifiasm)进行自我校正、修剪和组装,使用默认参数【43】。hifiasm(v.0.13)的初始输出生成了p_ctg草图组装。使用GenomeScope 2.0【69】的基于k-mer的方法估计基因组杂合度,估计其纯合度接近99.8%(补充数据图S2A-D)。然后,使用MUMmer(v.4.0.0)【70】的‘scaffold’工具生成基于同源性的框架,使用12X.v0参考基因组(补充数据图S3)。通过应用MUMmer工具,我们将contig级别的组装排序和定向到19条染色体,并将相邻的contig连接生成具有100 N的框架。最后,通过将以前版本的PN40024的基因组测序数据映射到新的组装基因组上,并在IGV(v.2.12.3)软件中可视化以观察缺口区域是否有读数支持来手动调整组装(补充数据图S1)。通过将缺口周围的50 bp序列映射到PN40024.v4的连续长读段上填补和关闭缺口,最终获得了所有19条葡萄染色体的无间隙T2T PN40024组装。基于BUSCO【71】的完整性和重复评分对组装进行了检查。对于高度杂合区域的分相,使用minimap2将所有读数对齐到PN_T2T组装。由hifiasm和ragtag组装的主contigs用于将这些contigs分相为两个单倍型。
基因和转座子的注释
我们使用了一种自开发的方法进行基因组注释。首先使用转录本和Uni-Prot作为证据搜索假定基因。然后为假定基因构建初步基因模型,并使用AUGUSTUS(v.3.4.0)【72】进行进一步搜索。然后过滤所有找到的假定基因片段,包括涉及重复区域的基因、编码序列长度小于90的基因以及没有任何证据支持的基因。我们尝试补充缺失的基因,并对完整基因进行可变剪接分析。最后,通过从Pfam数据库下载的隐马尔可夫模型检查所有结果,以获得最终的基因模型。Interproscan(v.5.56-89.0)【73】用于我们的组装功能注释,Pfam(v.34.0)【74】和Coils(v.2.2.1)【75】用于结构域识别(https://github.com/unavailable-2374/Genome-Wide-Annotation-Pipeline)。
主要的重复序列分析概述如图2A所示,首先通过RepeatModeler(open-2.0.3)【76】和一系列脚本构建了一个pan-Vitis重复家族数据库,然后使用RepeatMasker(open-4.1.2)进行应用。为了构建这个pan-Vitis重复数据库,我们从NCBI下载了17个Vitis基因组,然后使用RepeatModeler2识别TE家族。之后,我们获得了17个TE家族的共识fasta文件,并通过去除单拷贝和失败注释来聚合这些文件。我们使用NCBI-BLAST+2.9.0【77】去除一些冗余序列(-i 80%,-l 80%)。接下来,我们获得了重复序列的最终文件,然后使用deepTE【78】的植物模型分类未分类的重复元素。最后,通过RepeatMasker注释了完整参考基因组的重复序列。
葡萄参考基因组不同版本的比较
为了比较先前版本的葡萄基因组与PN_T2T,我们使用minimap2对基因组进行了比对,并使用SAMtools索引比对BAM文件(minimap2 -ax asm5 -t 4 --eqx A.fa B.fa | samtools sort -O BAM -> A_B.bam, samtools index A_B.bam)。接下来,为了检测基因组之间的结构变异,我们需要找到基因组之间的共线性和结构重排。为此,我们使用SyRI(syri -c A_B.bam -r A.fa -q B.fa -F B --prefix A_B)。最后,使用Plotsr生成图表(plotsr --sr A_Bsyri.out --sr B_Csyri.out --sr C_Dsyri.out --genomes genomes.txt -o output_plot.pdf, https://github.com/schneebergerlab/plotsr)。MUMmer(v.4.0.0)用于通过全基因组比对比较12X.v0基因组与参考基因组PN_T2T【70】。首先,我们使用nucmer对两个基因组序列进行比对(nucmer --mum),然后过滤最小比对长度为10,000 bp的单对单比对(delta-filter -i 95 -l 10,000)。
使用SAMtools(v.1.7)提取12X.v0中染色体18(25.0-26.0 Mb)的序列,并比对PN_T2T中的序列。使用python脚本(getgaps.py)检测缺口信息,最后使用LINKVIEW2(https://github.com/YangJianshun/LINKVIEW2)可视化比对结果。
端粒和着丝粒的识别
通过使用TIDK(v.0.2.0)(https://github.com/tolkit/telomeric-identifier)探索端粒重复单元,选项为tidk explore -f genome.fa --minimum 5 --maximum 12 -o tidk_explore -t 2 --log --dir telomere_find --extension TSV。然后使用以下参数搜索整个基因组:tidk search -f genome.fa -s TTTAGGG -o tidk_search --dir telomere_find。最后,我们基于TIDK图完成了端粒的快速统计,并使用R脚本可视化端粒峰。
对于着丝粒注释,使用TRF(v.4.09)【79】完成串联重复注释,参数为trf genome.fa 2 7 7 80 10 50 500 -f -d -m,然后使用TRF2GFF(https://github.com/Adamtaranto/TRF2GFF)合并注释结果。为了完成数据统计和可视化,我们使用Linux系统中的awk命令提取信息,并在IGV(v.2.12.3)【80】中分析结果。我们使用四个软件展示着丝粒区域的更多细节:Iqtree(v. 2.1.4-beta)【81】用于构建系统发育树(选项:-m GTR+I+G -bb 1000 -bnni -alrt 1000);itol(v.6)【82】用于可视化系统发育树;GeneDoc(v.2.7.0)(https://github.com/karlnicholas/genedoc)用于实现多序列比对;R脚本用于绘制数据统计和排版细节。
为了检测捕获在着丝粒区域中的基因功能,我们下载了Swiss-Prot的蛋白质序列库(2022/08/30, https://ftp.ncbi.nlm.nih.gov/blast/db/FASTA/)进行本地blast。之后,我们提取了所有由diamond(v.2.0.15)(参数:-k 1 -e 0.00001, https://github.com/python-diamond/Diamond)比对的PN_T2T蛋白质序列。我们进一步将Swiss-Prot ID上传到DAVID(https://david.ncifcrf.gov/tools.jsp)完成GO富集和注释。最后,通过我们的R脚本完成数据可视化。
基因簇的识别
为了定义参考基因组中的基因簇,使用gffread提取蛋白质序列,然后通过BLASTP进行全对全比对,筛选条件为e值<1e-5,相似性>30%。筛选后的比对结果结合功能注释过滤掉不共享相同结构域的比对结果。最终,通过识别三个连续相同Pfam登录号(https://www.ebi.ac.uk/interpro/entry/pfam/#table)确定基因簇的存在,使用这些Pfam登录号作为种子,上下30个基因找到具有相同Pfam登录号的基因。总共发现了377个基因簇(补充数据表S8)。
PN40024克隆中的杂合性
从NCBI数据库下载了四个重测序样本(SRR6156373, SRR8835144, SRR8835157, SRR8835168)并映射到新组装的PN_T2T基因组进行SNP调用。质量控制后的读数使用bwa(v.0.7.15)默认参数映射到基因组。SAMtools(v.1.4)和GATK(v.4.1.8)用于对BAM文件进行排序和索引,不去重。将gvcf文件在GATK中合并并用于联合调用所有样本的SNP。为了获得高质量的SNP,我们对SNP调用进行了严格筛选,基于以下标准:(i)在所有样本中具有两个以上等位基因的SNP被移除,参数为vcftools --min-alleles 2 --max-alleles 2;(ii)我们移除了质量分数(GQ)<30(--minGQ 30)和缺失率为0的SNP(--max-missing 1);(iii)SNP的小等位基因频率(MAFs)≥0.01,以去除不变位点。