芒果炭疽病抗感品種全基因組重測序分析

刁興旺， 吳莉君， 何 紅， 姚全勝， 柳 鳳

(1.廣東海洋大學濱海農(nóng)業(yè)學院，廣東湛江 524088； 2.中國熱帶農(nóng)業(yè)科學院南亞熱帶作物研究所/海南省熱帶園藝產(chǎn)品采后生理與保鮮重點實驗室/農(nóng)業(yè)部熱帶果樹生物學重點實驗室，廣東湛江 524091)

芒果(MangiferaindicaL.)屬漆樹科(Anacradiaceae)芒果屬(Mangifera)常綠喬本，原產(chǎn)于印度和東南亞，目前主要分布在熱帶和亞熱帶地區(qū)。芒果樹具有速生易長、抗逆性強、結(jié)果早、產(chǎn)量高、易栽培管理、經(jīng)濟壽命長等優(yōu)點，在我國熱帶農(nóng)業(yè)經(jīng)濟中占有重要地位[1]，但芒果炭疽病對產(chǎn)業(yè)危害嚴重。芒果炭疽病是芒果產(chǎn)業(yè)中普遍且危害性最大的一種真菌病害，其病原菌主要是膠孢炭疽菌[Colletotrichumgloeosporioides(Penz.) Sacc.]。該病害主要危害芒果的新梢、嫩葉、花和果實，也是芒果最嚴重的采后病害之一，在貯運期會引起果實的腐爛，占果實總病害量的70%以上，嚴重影響芒果的品質(zhì)和經(jīng)濟價值[2]。選育抗病品種是防控芒果炭疽病經(jīng)濟有效的方法。與傳統(tǒng)選育技術(shù)相比，分子標記輔助選擇育種以及基因組輔助選擇育種等技術(shù)是基于全基因組范圍內(nèi)發(fā)掘變異位點，在果樹遺傳育種中不僅選擇準確性高而且不受環(huán)境影響，可以極大提高選擇效率并縮短育種年限[3]。

全基因重測序可為分子標記等育種技術(shù)提供大量突變位點信息，隨著測序技術(shù)的發(fā)展及成本降低，已有至少20個果樹全基因組測序并公布基因草圖[4]，包括蘋果(Malusdomestica)[5]、梨(Pyrusspp)[6]、棗(ZiziphusjujubaMill.)[7]、龍眼(DimocarpuslonganLour.)[8]、柑橘(CitrusreticulataBlanco)[9-10]、芒果[11]等。以此為基礎開展的全基因組重測序分析也在果樹重要經(jīng)濟性狀研究領域取得重要成果。2020年高立志教授團隊選用芒果栽培種紅象牙(雜合度約為1.8%)，結(jié)合單分子實時測序(SMRT)、Illumina測序技術(shù)和遺傳圖譜，生成了一個包含 34 529 個預測蛋白質(zhì)編碼基因 371.6 MB 大小的芒果基因組，并公開發(fā)布了染色體級別的芒果參考基因組序列，為基于全基因組重測序技術(shù)發(fā)掘芒果全基因組范圍內(nèi)重要性狀相關(guān)的單核苷酸多態(tài)性(SNP)、插入/缺失(InDel)位點和進一步開發(fā)相關(guān)分子標記進行輔助育種奠定了可行的基礎。

目前，關(guān)于芒果炭疽病抗感品種遺傳群體構(gòu)建以及全基因組重測序相關(guān)分析尚未見報道。本研究選取芒果炭疽病高感品種愛文和炭疽病高抗品種金煌，通過全基因組重測序，并以品種“紅象牙”的基因組為參考比對后分析檢測SNP、InDel位點，對多態(tài)性進行差異分析及功能注釋。一方面可用于后期開發(fā)多態(tài)性分子標記，芒果抗炭疽病品種快速鑒定、種質(zhì)資源評價及遺傳作圖選擇利用；另一方面可為全基因組關(guān)聯(lián)分析、功能基因定位及分子標記輔助育種提供研究基礎。

1 材料與方法

1.1 試驗材料

2019年秋季，在中國熱帶農(nóng)業(yè)科學院南亞熱帶作物研究所(110.28°E，21.17°N)將采自我國芒果主產(chǎn)區(qū)廣東省、廣西壯族自治區(qū)、云南省、四川省的芒果炭疽病高感品種愛文和高抗品種金煌的葉片進行初步處理。選取健康嫩葉送至廣州基迪奧生物科技有限公司進行全基因組重測序。

1.2 愛文和金煌全基因組重測序

經(jīng)質(zhì)檢合格的樣品用于構(gòu)建DNA文庫，用FASTP (版本0.18.0)[12]對Illumina平臺的原始數(shù)據(jù)進行過濾，最終得到質(zhì)量較好的有效數(shù)據(jù)用于后續(xù)分析。

1.3 基因組變異檢測與注釋

使用比對軟件BWA(版本 0.7.12)[13]將2個樣本全基因組重測序數(shù)據(jù)對比到參考基因組上(https：//bigd.big.ac.cn/search/？dbId=gwh&q=PRJCA002248)，結(jié)果使用軟件Picard (版本1.129)(Http：//sourceforge.net/projects/picard/)標記重復讀長(reads)，統(tǒng)計標記后的reads深度及覆蓋度。使用BEDTools(版本 2.25.0)[14]軟件進行覆蓋度統(tǒng)計。使用變異檢測軟件GATK(版本 3.4-46)[15]進行群體SNP檢測，并使用ANNOVAR(版本 2)[16]對檢測出的變異(variant)進行功能注釋。對檢測、過濾得到最終的SNP位點集進行統(tǒng)計分析。使用BreakDancer(版本 1.1.2)[17]進行大片段結(jié)構(gòu)變異的檢測。使用CNVnator(版本0.3.2)[18]來進行拷貝數(shù)變異的檢測。

1.4 變異基因分析

在獲得的基因組變異檢測分析基礎上，查找愛文和金煌與參考基因組之間可能存在功能差異的基因序列，進而通過Blastx將基因序列比對到蛋白數(shù)據(jù)庫NR、SwissProt、KEGG和KOG(E值<0.000 01)，得到給定基因具有最高序列相似性的蛋白，從而得到該基因的蛋白功能注釋信息，進行統(tǒng)計分析。

2 結(jié)果與分析

2.1 全基因組重測序數(shù)據(jù)分析

提取芒果愛文、金煌樣本的DNA經(jīng)質(zhì)檢合格后進行全基因組重測序建庫，分別獲得4.73、5.58 Gb原始數(shù)據(jù)通過測序得到的原始圖像數(shù)據(jù)轉(zhuǎn)化序列數(shù)據(jù)后進行測序評估，芒果愛文、金煌重測序數(shù)據(jù)經(jīng)質(zhì)量控制后分別得到37 207 272、31 501 486個過濾后的數(shù)據(jù)，正確識別率大于Q30 (測序正確率為99.90%)的堿基分別占92.27%、92.76%，基因組GC含量愛文略低于金煌，分別為36.18%、37.80% (表1)。測序深度為15.01×和13.73×，基因組覆蓋度分別為92.23%、92.27%。測序質(zhì)量良好，所得數(shù)據(jù)量和序列長度可準確反映基因組信息。

2.2 SNP檢測與注釋

通過與參考基因組比對(https：//ngdc.cncb.ac.cn/search/？dbId=gwh&q=PRJCA002248)，分別檢測到愛文SNP位點 3 316 161個，金煌SNP位點 3 075 003 個，其中愛文的轉(zhuǎn)換數(shù)是2 288 222個，顛換1 027 939個；金煌的轉(zhuǎn)換數(shù)是2 129 307個，顛換945 696個，2個芒果品種的SNP位點類型多為轉(zhuǎn)換。金煌的雜合率為57.12%，略高于愛文的雜合率(表2)。

表1 芒果愛文、金煌的重測序數(shù)據(jù)質(zhì)量對比

表2 芒果愛文、金煌SNP位點數(shù)量統(tǒng)計

對芒果愛文、金煌2個樣本基因組SNP進行注釋統(tǒng)計(表3)，經(jīng)分析發(fā)現(xiàn)，芒果愛文、金煌的SNP變異主要分布在基因間區(qū)，分別檢測到1 928 078、1 775 226 個；2個品種在編碼區(qū)的非同義突變分別為126 159、123 930個；分別因SNP突變獲得終止子 1 927、1 866個，丟失終止子470、460個；分別在可變剪切位點2 bp以內(nèi)的基因組區(qū)域內(nèi)檢測到SNP位點 1 608、1 573個。

表3 芒果愛文、金煌基因組中SNP位點注釋

2.3 芒果愛文、金煌基因組 InDel檢測與注釋

結(jié)合芒果愛文、金煌2個樣本全基因組重測序數(shù)據(jù)與參考基因組比對結(jié)果，檢測發(fā)現(xiàn)2個樣本分別獲得244 686、201 520個InDel位點，其中2個樣本純合基因型的InDel數(shù)量分別為154 530、117 125個，雜合率分別為36.85%、41.88%。分別統(tǒng)計芒果愛文、金煌2個樣本基因組InDel位點進行注釋并對基因組位置信息和編碼信息進行統(tǒng)計發(fā)現(xiàn)，芒果愛文、金煌2個樣本InDel位點主要分布在基因間區(qū)分別有118 467、94 428個。2個樣本InDel位點主要為插入類型分別有4 048、3 861個；2個芒果品種在外顯子區(qū)域InDel位點存在顯著差異分別有 8 312、7 605個(表4、表5)。

表4 芒果愛文、金煌基因組中InDel編碼信息統(tǒng)計

表5 芒果愛文、金煌基因組中InDel位點注釋

2.4 芒果愛文、金煌在DNA水平變異基因分析

將芒果愛文、金煌間DNA水平變異基因進行KEGG、KOG、NR、SwissProt數(shù)據(jù)庫注釋，發(fā)現(xiàn)愛文、金煌基因組間存在多種類型的變異基因，與功能數(shù)據(jù)庫對比得出2個樣本間對應變異基因數(shù)量分別為28 981、17 151、30 013、22 910個(圖1)。

為進一步研究含有突變位點的基因功能，向GO數(shù)據(jù)庫各術(shù)語標簽(term)進行映射(圖2)。對芒果愛文、金煌基因組重測序結(jié)果進行GO 注釋統(tǒng)計，二者差異基因主要存在于參與生物過程基因中，其中參與水楊酸(SA)介導的信號通路(GO：0009863)的有354個(1.93%)，參與水楊酸代謝過程(GO：0009696)的有269個(1.47%)，富集到對細菌的反應(GO：0009617)相關(guān)基因為656個(3.58%)，對真菌的反應(GO：0009620)相關(guān)基因有523個(2.85%)，對線蟲的反應(GO：0009624)相關(guān)基因9個(0.05%)。

蛋白相鄰類的聚簇(KOG) 變異基因分析結(jié)果顯示，芒果愛文、金煌2個樣本間重測序數(shù)據(jù)在以下幾類功能基因存在變異，其中信號轉(zhuǎn)導途徑類基因2 529(12%)個，能量生產(chǎn)和轉(zhuǎn)化類基因778(4%)個，氨基酸轉(zhuǎn)運與代謝類基因798(4%)個，防御機制類基因181(1%)個(圖3)。

根據(jù)芒果愛文、金煌基因組重測序結(jié)果進行KEGG注釋(圖4)可知，含有差異位點的基因可注釋到139條代謝途徑，其中與植物抗病性緊密相關(guān)的代謝途徑：植物激素信號轉(zhuǎn)導代謝途徑(ko04075)有差異基因數(shù)395 (5.93%)，植物-病原互作代謝途徑(ko04626)差異基因數(shù)為303(4.55%)，淀粉和蔗糖代謝途徑(ko00500)基因差異數(shù)為200個(3.00%)。

2.5 芒果愛文、金煌抗病相關(guān)基因差異性分析

通過對芒果愛文、金煌2個樣本全基因組重測序數(shù)據(jù)進行樣本間差異基因統(tǒng)計共獲得33 466個差異基因，對其中功能性變異基因進行初步分類和篩選，得到2個樣本間與E3泛素蛋白連接酶相關(guān)差異基因數(shù)為349個，與WRKY轉(zhuǎn)錄因子相關(guān)差異基因數(shù)為100個，與植物抗性相關(guān)差異性基因數(shù)為383個，其中包括抗病蛋白、耐藥性、抗線蟲蛋白等基因。在所獲得的抗病蛋白差異基因中包括抗病蛋白RGA2差異基因48個、葉銹病10號抗病基因座受體差異基因31個、抗銹激酶基因24個，抗病蛋白RPM1差異基因24個、抗病蛋白DSC1差異基因19個、增強抗病性蛋白基因19個、NBS-LRR抗病蛋白基因18個、含NB-ARC結(jié)構(gòu)域的蛋白質(zhì)變異基因16個、抗病蛋白At4g27190差異基因15個，抗霜霉病蛋白基因7個、抗稻瘟病蛋白基因9個。含BTB/POZ結(jié)構(gòu)域和錨蛋白重復序列的NPR1蛋白差異基因5個。這些抗病蛋白與植物抗病性密切相關(guān)，2個樣本中抗性相關(guān)基因的變異可能是導致二者間抗病性差異的原因。

3 討論與結(jié)論

植物對病原菌的抗性由多方面原因組成，其中最重要的原因是遺傳因素中的抗病基因，決定抗病基因的產(chǎn)物的結(jié)構(gòu)、功能、信號傳遞途徑的因素是多方面的[19]。Li等公開發(fā)表了全球首個栽培芒果染色體級別的參考基因組序列，使開展芒果種質(zhì)資源的全基因組重測序分析成為可能[20]。芒果品種愛文肉質(zhì)膩滑，纖維少，味甜，品質(zhì)較好，但是易感炭疽??；芒果品種金煌樹勢強，果實耐貯藏，中熟，對炭疽病表現(xiàn)為高抗。對芒果炭疽病感病品種愛文和芒果炭疽病抗病品種金煌的全基因組重測序，通過全基因組變異分析技術(shù)可以深度挖掘二者存在抗感差異的原因，在基因水平闡明芒果抗炭疽病機制，并且開展芒果優(yōu)異資源挖掘和基因功能分析。

近年來，基于全基因重測序的SNP檢測分析已經(jīng)成為植物學科尤其是果樹學遺傳研究的高效便捷的手段[21]。陳璇等對我國野生型大麻和栽培型大麻進行測序深度為10×全基因組重測序分析，共檢測到2 264 150個單核苷酸多態(tài)性位點，2個樣本的雜合度分別為0.12%、0.11%，在一定程度上反映我國野生大麻和栽培大麻在基因組水平上的差異性[12]。邵勤等在盤菜品種溫盤2號測序深度為46×全基因組重測序中共檢測到1 424 852個SNP位點，298 244 個Small InDel位點，3 068個結(jié)構(gòu)變異，共導致了42 334個基因變異，并發(fā)現(xiàn)了大量參與肉質(zhì)根形成發(fā)育的關(guān)鍵基因發(fā)生突變[22]。非同義SNP位點變異會使氨基酸序列發(fā)生改變進而影響蛋白質(zhì)序列。洪森榮等對上饒梨2個品種全基因組重測序分析發(fā)現(xiàn)，在六月雪、黃皮消2個品種中分別檢測到6 171 357、6 140 603個SNP位點，共有 2 282 個SNP 關(guān)聯(lián)了2 067 個基因，2個樣本InDel位點分別為800 388、799 603個，共5 115個差異InDel關(guān)聯(lián)到了3 682個基因，并檢測到耐藥蛋白和抗病蛋白等重要基因發(fā)生變異[23]。通過單核苷酸多態(tài)性以及結(jié)構(gòu)變異分析可以發(fā)掘同一物種間存在表型差異的變異本質(zhì)。楊赟等在紅葉杜仲、小葉杜仲的10×深度全基因組中發(fā)現(xiàn)嚴重影響蛋白質(zhì)功能的SNP位點有1 516、1 640個，中度影響功能的SNP位點有 41 328、47 192個，篩選了12個特異性的SNP位點，并從SNP位點的差異出發(fā)探討了導致杜仲顏色差異的原因[24]。劉冰浩等通過對桂柚一號、沙田柚進行全基因組重測序，分別檢測到了 2 985 308、2 997 229 個SNP差異位點，并且發(fā)現(xiàn)結(jié)構(gòu)變異(SV)在兩者變異中起到了很大作用[25]。

本研究通過對芒果易感炭疽病品種愛文、炭疽病高抗品種金煌進行全基因組重測序，以“紅象牙”的基因組為參考，分別檢測到3 316 161、3 075 003個單核苷酸多態(tài)性位點，通過對比到數(shù)據(jù)庫檢測到二者間SNP突變共導致了33 466個基因變異，這些變異基因?qū)е铝?個品種間植物-病原互作代謝途徑、激素介導信號通路、抗病蛋白表達等多種生物過程的差異。植物激素不僅可以作為內(nèi)源信號分子參與植物逆境脅迫應答，也可以作為外援信號分子誘導植物抗病反應[26]。在2個樣本間檢測到大量與E3泛素蛋白連接酶相關(guān)基因變異，E3泛素蛋白連接酶不僅可以將泛素與底物蛋白鏈接參與細胞功能調(diào)節(jié)[27]，還可以作為正負調(diào)節(jié)因子參與植物對病原菌防御中的效應器誘導免疫反應[28]。蛋白質(zhì)泛素化修飾廣泛地參與植物免疫應答的調(diào)控，尤其廣泛地參與了茉莉酸、水楊酸和乙烯等植物激素防御途徑中的信號傳導和生物合成等環(huán)節(jié)。泛素介導的蛋白酶水解在水楊酸的信號傳導和生物合成等多個過程中發(fā)揮作用，水楊酸是系統(tǒng)獲得性免疫中的關(guān)鍵激素，在植物抵御真菌病害侵染中發(fā)揮重要作用[29]，在愛文、金煌2個芒果品種間SNP突變而引起的差異基因中，發(fā)現(xiàn)有354個參與水楊酸介導的信號通路，有269個參與水楊酸代謝過程，水楊酸以及其類似物可激活抗病相關(guān)基因表達并對病原菌的侵害產(chǎn)生明顯抗性[30]。張慶雨等研究發(fā)現(xiàn)，水楊酸能夠提高草莓對炭疽病侵染的抗性，水楊酸可能參與了表達水平和炭疽病抗性水平相關(guān)的FaNBS20基因的重要響應途徑[31]。水楊酸還可以作為誘導信號激活NPR1表達，NPR1的過量表達轉(zhuǎn)基因植物株系增強了對細菌和真菌等病原菌的抗性[32]，Silva 等通過草莓過表達NPR1發(fā)現(xiàn)其對炭疽病的抗性明顯增加[33]。NPR1包含的BTB結(jié)構(gòu)域增加了其自身參與泛素化的可能性，另外NPR1包含的錨定蛋白重復序列可以介導蛋白之間互作，SA很可能參與調(diào)節(jié)此過程[34-35]。同時NPR1也是水楊酸途徑中關(guān)鍵的輔基活因子，研究表明，NPR1位點突變會導致植物不能激活系統(tǒng)性獲得免疫而易感病[36]。本研究中2個樣本間檢測到的5個NPR1基因變異位點分布在芒果5條不同染色體上，可能是導致2個樣本表現(xiàn)對炭疽病抗性差異的關(guān)鍵。

炭疽病在不同植物中的抗性遺傳表現(xiàn)不同，一般是由單基因控制的，大多數(shù)受顯性基因控制，少數(shù)情況下是由隱性抗病基因或者不完全顯性抗病基因控制，不同抗病基因間相互獨立或存在互作關(guān)系[37]。通過對芒果炭疽病抗感品種金煌、高感品種愛文全基因組重測序數(shù)據(jù)分析，深度挖掘與炭疽病抗性相關(guān)的差異基因并解析其所參與的代謝途徑，可為進一步闡明其遺傳特點，進行抗病基因定位、克隆、轉(zhuǎn)化及特異分子標記開發(fā)并且選育芒果抗病品種提供理論研究基礎和良好的數(shù)據(jù)支撐。