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

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

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

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

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

1 材料與方法

1.1 試驗(yàn)材料

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

1.2 愛文和金煌全基因組重測(cè)序

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

1.3 基因組變異檢測(cè)與注釋

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

1.4 變異基因分析

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

2 結(jié)果與分析

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

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

2.2 SNP檢測(cè)與注釋

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

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

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

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

表3 芒果愛文、金煌基因組中SNP位點(diǎn)注釋

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

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

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

表5 芒果愛文、金煌基因組中InDel位點(diǎn)注釋

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

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

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

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

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

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

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

3 討論與結(jié)論

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

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

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

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