芒果炭疽病抗感品种全基因组重测序分析

刁兴旺， 吴莉君， 何 红， 姚全胜， 柳 凤

(1.广东海洋大学滨海农业学院，广东湛江 524088； 2.中国热带农业科学院南亚热带作物研究所/海南省热带园艺产品采后生理与保鲜重点实验室/农业部热带果树生物学重点实验室，广东湛江 524091)

芒果(MangiferaindicaL.)属漆树科(Anacradiaceae)芒果属(Mangifera)常绿乔本，原产于印度和东南亚，目前主要分布在热带和亚热带地区。芒果树具有速生易长、抗逆性强、结果早、产量高、易栽培管理、经济寿命长等优点，在我国热带农业经济中占有重要地位[1]，但芒果炭疽病对产业危害严重。芒果炭疽病是芒果产业中普遍且危害性最大的一种真菌病害，其病原菌主要是胶孢炭疽菌[Colletotrichumgloeosporioides(Penz.) Sacc.]。该病害主要危害芒果的新梢、嫩叶、花和果实，也是芒果最严重的采后病害之一，在贮运期会引起果实的腐烂，占果实总病害量的70%以上，严重影响芒果的品质和经济价值[2]。选育抗病品种是防控芒果炭疽病经济有效的方法。与传统选育技术相比，分子标记辅助选择育种以及基因组辅助选择育种等技术是基于全基因组范围内发掘变异位点，在果树遗传育种中不仅选择准确性高而且不受环境影响，可以极大提高选择效率并缩短育种年限[3]。

全基因重测序可为分子标记等育种技术提供大量突变位点信息，随着测序技术的发展及成本降低，已有至少20个果树全基因组测序并公布基因草图[4]，包括苹果(Malusdomestica)[5]、梨(Pyrusspp)[6]、枣(ZiziphusjujubaMill.)[7]、龙眼(DimocarpuslonganLour.)[8]、柑橘(CitrusreticulataBlanco)[9-10]、芒果[11]等。以此为基础开展的全基因组重测序分析也在果树重要经济性状研究领域取得重要成果。2020年高立志教授团队选用芒果栽培种红象牙(杂合度约为1.8%)，结合单分子实时测序(SMRT)、Illumina测序技术和遗传图谱，生成了一个包含 34 529 个预测蛋白质编码基因 371.6 MB 大小的芒果基因组，并公开发布了染色体级别的芒果参考基因组序列，为基于全基因组重测序技术发掘芒果全基因组范围内重要性状相关的单核苷酸多态性(SNP)、插入/缺失(InDel)位点和进一步开发相关分子标记进行辅助育种奠定了可行的基础。

目前，关于芒果炭疽病抗感品种遗传群体构建以及全基因组重测序相关分析尚未见报道。本研究选取芒果炭疽病高感品种爱文和炭疽病高抗品种金煌，通过全基因组重测序，并以品种“红象牙”的基因组为参考比对后分析检测SNP、InDel位点，对多态性进行差异分析及功能注释。一方面可用于后期开发多态性分子标记，芒果抗炭疽病品种快速鉴定、种质资源评价及遗传作图选择利用；另一方面可为全基因组关联分析、功能基因定位及分子标记辅助育种提供研究基础。

1 材料与方法

1.1 试验材料

2019年秋季，在中国热带农业科学院南亚热带作物研究所(110.28°E，21.17°N)将采自我国芒果主产区广东省、广西壮族自治区、云南省、四川省的芒果炭疽病高感品种爱文和高抗品种金煌的叶片进行初步处理。选取健康嫩叶送至广州基迪奥生物科技有限公司进行全基因组重测序。

1.2 爱文和金煌全基因组重测序

经质检合格的样品用于构建DNA文库，用FASTP (版本0.18.0)[12]对Illumina平台的原始数据进行过滤，最终得到质量较好的有效数据用于后续分析。

1.3 基因组变异检测与注释

使用比对软件BWA(版本 0.7.12)[13]将2个样本全基因组重测序数据对比到参考基因组上(https：//bigd.big.ac.cn/search/？dbId=gwh&q=PRJCA002248)，结果使用软件Picard (版本1.129)(Http：//sourceforge.net/projects/picard/)标记重复读长(reads)，统计标记后的reads深度及覆盖度。使用BEDTools(版本 2.25.0)[14]软件进行覆盖度统计。使用变异检测软件GATK(版本 3.4-46)[15]进行群体SNP检测，并使用ANNOVAR(版本 2)[16]对检测出的变异(variant)进行功能注释。对检测、过滤得到最终的SNP位点集进行统计分析。使用BreakDancer(版本 1.1.2)[17]进行大片段结构变异的检测。使用CNVnator(版本0.3.2)[18]来进行拷贝数变异的检测。

1.4 变异基因分析

在获得的基因组变异检测分析基础上，查找爱文和金煌与参考基因组之间可能存在功能差异的基因序列，进而通过Blastx将基因序列比对到蛋白数据库NR、SwissProt、KEGG和KOG(E值<0.000 01)，得到给定基因具有最高序列相似性的蛋白，从而得到该基因的蛋白功能注释信息，进行统计分析。

2 结果与分析

2.1 全基因组重测序数据分析

提取芒果爱文、金煌样本的DNA经质检合格后进行全基因组重测序建库，分别获得4.73、5.58 Gb原始数据通过测序得到的原始图像数据转化序列数据后进行测序评估，芒果爱文、金煌重测序数据经质量控制后分别得到37 207 272、31 501 486个过滤后的数据，正确识别率大于Q30 (测序正确率为99.90%)的碱基分别占92.27%、92.76%，基因组GC含量爱文略低于金煌，分别为36.18%、37.80% (表1)。测序深度为15.01×和13.73×，基因组覆盖度分别为92.23%、92.27%。测序质量良好，所得数据量和序列长度可准确反映基因组信息。

2.2 SNP检测与注释

通过与参考基因组比对(https：//ngdc.cncb.ac.cn/search/？dbId=gwh&q=PRJCA002248)，分别检测到爱文SNP位点 3 316 161个，金煌SNP位点 3 075 003 个，其中爱文的转换数是2 288 222个，颠换1 027 939个；金煌的转换数是2 129 307个，颠换945 696个，2个芒果品种的SNP位点类型多为转换。金煌的杂合率为57.12%，略高于爱文的杂合率(表2)。

表1 芒果爱文、金煌的重测序数据质量对比

表2 芒果爱文、金煌SNP位点数量统计

对芒果爱文、金煌2个样本基因组SNP进行注释统计(表3)，经分析发现，芒果爱文、金煌的SNP变异主要分布在基因间区，分别检测到1 928 078、1 775 226 个；2个品种在编码区的非同义突变分别为126 159、123 930个；分别因SNP突变获得终止子 1 927、1 866个，丢失终止子470、460个；分别在可变剪切位点2 bp以内的基因组区域内检测到SNP位点 1 608、1 573个。

表3 芒果爱文、金煌基因组中SNP位点注释

2.3 芒果爱文、金煌基因组 InDel检测与注释

结合芒果爱文、金煌2个样本全基因组重测序数据与参考基因组比对结果，检测发现2个样本分别获得244 686、201 520个InDel位点，其中2个样本纯合基因型的InDel数量分别为154 530、117 125个，杂合率分别为36.85%、41.88%。分别统计芒果爱文、金煌2个样本基因组InDel位点进行注释并对基因组位置信息和编码信息进行统计发现，芒果爱文、金煌2个样本InDel位点主要分布在基因间区分别有118 467、94 428个。2个样本InDel位点主要为插入类型分别有4 048、3 861个；2个芒果品种在外显子区域InDel位点存在显著差异分别有 8 312、7 605个(表4、表5)。

表4 芒果爱文、金煌基因组中InDel编码信息统计

表5 芒果爱文、金煌基因组中InDel位点注释

2.4 芒果爱文、金煌在DNA水平变异基因分析

将芒果爱文、金煌间DNA水平变异基因进行KEGG、KOG、NR、SwissProt数据库注释，发现爱文、金煌基因组间存在多种类型的变异基因，与功能数据库对比得出2个样本间对应变异基因数量分别为28 981、17 151、30 013、22 910个(图1)。

为进一步研究含有突变位点的基因功能，向GO数据库各术语标签(term)进行映射(图2)。对芒果爱文、金煌基因组重测序结果进行GO 注释统计，二者差异基因主要存在于参与生物过程基因中，其中参与水杨酸(SA)介导的信号通路(GO：0009863)的有354个(1.93%)，参与水杨酸代谢过程(GO：0009696)的有269个(1.47%)，富集到对细菌的反应(GO：0009617)相关基因为656个(3.58%)，对真菌的反应(GO：0009620)相关基因有523个(2.85%)，对线虫的反应(GO：0009624)相关基因9个(0.05%)。

蛋白相邻类的聚簇(KOG) 变异基因分析结果显示，芒果爱文、金煌2个样本间重测序数据在以下几类功能基因存在变异，其中信号转导途径类基因2 529(12%)个，能量生产和转化类基因778(4%)个，氨基酸转运与代谢类基因798(4%)个，防御机制类基因181(1%)个(图3)。

根据芒果爱文、金煌基因组重测序结果进行KEGG注释(图4)可知，含有差异位点的基因可注释到139条代谢途径，其中与植物抗病性紧密相关的代谢途径：植物激素信号转导代谢途径(ko04075)有差异基因数395 (5.93%)，植物-病原互作代谢途径(ko04626)差异基因数为303(4.55%)，淀粉和蔗糖代谢途径(ko00500)基因差异数为200个(3.00%)。

2.5 芒果爱文、金煌抗病相关基因差异性分析

通过对芒果爱文、金煌2个样本全基因组重测序数据进行样本间差异基因统计共获得33 466个差异基因，对其中功能性变异基因进行初步分类和筛选，得到2个样本间与E3泛素蛋白连接酶相关差异基因数为349个，与WRKY转录因子相关差异基因数为100个，与植物抗性相关差异性基因数为383个，其中包括抗病蛋白、耐药性、抗线虫蛋白等基因。在所获得的抗病蛋白差异基因中包括抗病蛋白RGA2差异基因48个、叶锈病10号抗病基因座受体差异基因31个、抗锈激酶基因24个，抗病蛋白RPM1差异基因24个、抗病蛋白DSC1差异基因19个、增强抗病性蛋白基因19个、NBS-LRR抗病蛋白基因18个、含NB-ARC结构域的蛋白质变异基因16个、抗病蛋白At4g27190差异基因15个，抗霜霉病蛋白基因7个、抗稻瘟病蛋白基因9个。含BTB/POZ结构域和锚蛋白重复序列的NPR1蛋白差异基因5个。这些抗病蛋白与植物抗病性密切相关，2个样本中抗性相关基因的变异可能是导致二者间抗病性差异的原因。

3 讨论与结论

植物对病原菌的抗性由多方面原因组成，其中最重要的原因是遗传因素中的抗病基因，决定抗病基因的产物的结构、功能、信号传递途径的因素是多方面的[19]。Li等公开发表了全球首个栽培芒果染色体级别的参考基因组序列，使开展芒果种质资源的全基因组重测序分析成为可能[20]。芒果品种爱文肉质腻滑，纤维少，味甜，品质较好，但是易感炭疽病；芒果品种金煌树势强，果实耐贮藏，中熟，对炭疽病表现为高抗。对芒果炭疽病感病品种爱文和芒果炭疽病抗病品种金煌的全基因组重测序，通过全基因组变异分析技术可以深度挖掘二者存在抗感差异的原因，在基因水平阐明芒果抗炭疽病机制，并且开展芒果优异资源挖掘和基因功能分析。

近年来，基于全基因重测序的SNP检测分析已经成为植物学科尤其是果树学遗传研究的高效便捷的手段[21]。陈璇等对我国野生型大麻和栽培型大麻进行测序深度为10×全基因组重测序分析，共检测到2 264 150个单核苷酸多态性位点，2个样本的杂合度分别为0.12%、0.11%，在一定程度上反映我国野生大麻和栽培大麻在基因组水平上的差异性[12]。邵勤等在盘菜品种温盘2号测序深度为46×全基因组重测序中共检测到1 424 852个SNP位点，298 244 个Small InDel位点，3 068个结构变异，共导致了42 334个基因变异，并发现了大量参与肉质根形成发育的关键基因发生突变[22]。非同义SNP位点变异会使氨基酸序列发生改变进而影响蛋白质序列。洪森荣等对上饶梨2个品种全基因组重测序分析发现，在六月雪、黄皮消2个品种中分别检测到6 171 357、6 140 603个SNP位点，共有 2 282 个SNP 关联了2 067 个基因，2个样本InDel位点分别为800 388、799 603个，共5 115个差异InDel关联到了3 682个基因，并检测到耐药蛋白和抗病蛋白等重要基因发生变异[23]。通过单核苷酸多态性以及结构变异分析可以发掘同一物种间存在表型差异的变异本质。杨赟等在红叶杜仲、小叶杜仲的10×深度全基因组中发现严重影响蛋白质功能的SNP位点有1 516、1 640个，中度影响功能的SNP位点有 41 328、47 192个，筛选了12个特异性的SNP位点，并从SNP位点的差异出发探讨了导致杜仲颜色差异的原因[24]。刘冰浩等通过对桂柚一号、沙田柚进行全基因组重测序，分别检测到了 2 985 308、2 997 229 个SNP差异位点，并且发现结构变异(SV)在两者变异中起到了很大作用[25]。

本研究通过对芒果易感炭疽病品种爱文、炭疽病高抗品种金煌进行全基因组重测序，以“红象牙”的基因组为参考，分别检测到3 316 161、3 075 003个单核苷酸多态性位点，通过对比到数据库检测到二者间SNP突变共导致了33 466个基因变异，这些变异基因导致了2个品种间植物-病原互作代谢途径、激素介导信号通路、抗病蛋白表达等多种生物过程的差异。植物激素不仅可以作为内源信号分子参与植物逆境胁迫应答，也可以作为外援信号分子诱导植物抗病反应[26]。在2个样本间检测到大量与E3泛素蛋白连接酶相关基因变异，E3泛素蛋白连接酶不仅可以将泛素与底物蛋白链接参与细胞功能调节[27]，还可以作为正负调节因子参与植物对病原菌防御中的效应器诱导免疫反应[28]。蛋白质泛素化修饰广泛地参与植物免疫应答的调控，尤其广泛地参与了茉莉酸、水杨酸和乙烯等植物激素防御途径中的信号传导和生物合成等环节。泛素介导的蛋白酶水解在水杨酸的信号传导和生物合成等多个过程中发挥作用，水杨酸是系统获得性免疫中的关键激素，在植物抵御真菌病害侵染中发挥重要作用[29]，在爱文、金煌2个芒果品种间SNP突变而引起的差异基因中，发现有354个参与水杨酸介导的信号通路，有269个参与水杨酸代谢过程，水杨酸以及其类似物可激活抗病相关基因表达并对病原菌的侵害产生明显抗性[30]。张庆雨等研究发现，水杨酸能够提高草莓对炭疽病侵染的抗性，水杨酸可能参与了表达水平和炭疽病抗性水平相关的FaNBS20基因的重要响应途径[31]。水杨酸还可以作为诱导信号激活NPR1表达，NPR1的过量表达转基因植物株系增强了对细菌和真菌等病原菌的抗性[32]，Silva 等通过草莓过表达NPR1发现其对炭疽病的抗性明显增加[33]。NPR1包含的BTB结构域增加了其自身参与泛素化的可能性，另外NPR1包含的锚定蛋白重复序列可以介导蛋白之间互作，SA很可能参与调节此过程[34-35]。同时NPR1也是水杨酸途径中关键的辅基活因子，研究表明，NPR1位点突变会导致植物不能激活系统性获得免疫而易感病[36]。本研究中2个样本间检测到的5个NPR1基因变异位点分布在芒果5条不同染色体上，可能是导致2个样本表现对炭疽病抗性差异的关键。

炭疽病在不同植物中的抗性遗传表现不同，一般是由单基因控制的，大多数受显性基因控制，少数情况下是由隐性抗病基因或者不完全显性抗病基因控制，不同抗病基因间相互独立或存在互作关系[37]。通过对芒果炭疽病抗感品种金煌、高感品种爱文全基因组重测序数据分析，深度挖掘与炭疽病抗性相关的差异基因并解析其所参与的代谢途径，可为进一步阐明其遗传特点，进行抗病基因定位、克隆、转化及特异分子标记开发并且选育芒果抗病品种提供理论研究基础和良好的数据支撑。