滑皮金柑全長轉(zhuǎn)錄組測序分析*

李智龍，熊志偉，尹 輝，高玉霞,?

(贛南師范大學(xué) a.生命科學(xué)學(xué)院； b.國家臍橙工程技術(shù)研究中心，江西 贛州 341000)

金柑(Kumquat)屬于蕓香科(Rutacease)柑橘亞科(Aurantioideae)金柑屬(Fortunella)植物，常見品種為4種：圓金柑(Fortunellajaponica)、金桔(Fortunellamargarita)、金彈(Fortunellacrassifolia)和山金柑(Fortunellahindsii)[1].金柑為常綠灌木或小喬木，樹形美觀，枝條密生有短刺，葉多為卵狀橢圓形，基部寬楔形.果實(shí)顏色黃亮，且含有豐富的維生素、果糖、果膠等營養(yǎng)成分，比其他柑橘類水果小，其整個(gè)水果包括果皮均可食用，也可以用糖漿腌制或保存，被用作傳統(tǒng)的民間藥物來治療呼吸道炎癥，為我國重要的經(jīng)濟(jì)水果[2-3].據(jù)文學(xué)作品記載，金柑在12世紀(jì)起源于中國，后傳入日本，在中國和日本都有相當(dāng)長的栽培歷史，1684年由倫敦園藝協(xié)會(huì)收藏家Robert Fortune帶入歐洲，不久后傳入北美[4].融安金柑在中國廣西已有300年種植歷史，滑皮金柑為融安金柑的實(shí)生變異株.對比普通融安金柑，滑皮金柑味更甜少辣而更受消費(fèi)者喜愛，是廣泛的栽培品種[5].滑皮金柑已成為我國重要的經(jīng)濟(jì)水果，僅在融安縣的種植面積達(dá)到6 272公頃，年產(chǎn)值達(dá)3.48億，是當(dāng)?shù)刂涡援a(chǎn)業(yè)[6].目前滑皮金柑的分子生物學(xué)研究幾乎空白，對滑皮金柑進(jìn)行組學(xué)研究有助于該品種的抗逆抗病相關(guān)研究，確保該產(chǎn)業(yè)的健康發(fā)展.

1 材料及方法

1.1 實(shí)驗(yàn)材料

2年齡滑皮金柑的葉、莖和根部樣品采自贛南師范大學(xué)國家臍橙工程技術(shù)研究中心苗圃，經(jīng)液氮冷凍后送至北京諾禾致源科技股份有限公司進(jìn)行總RNA的提取及文庫建立.

1.2 方法

1.2.1 總RNA樣品檢測及文庫的構(gòu)建

利用Trizol法對分別對金柑葉、莖和根樣品進(jìn)行總RNA提取[7].分別利用瓊脂凝膠電泳、Nanodrop儀、Qubit熒光定量儀和Agilent 2100分析儀對總RNA的降解程度、污染情況、純度、數(shù)量和完整性進(jìn)行檢測，經(jīng)檢測結(jié)果達(dá)到要求后混合均勻，用于文庫的建立.

利用Oligo(dT)富集總RNA中含有poly(A)的mRNA，后根據(jù)SMARTer PCR cDNA Synthesis Kit使用說明將mRNA反轉(zhuǎn)錄為cDNA，PCR擴(kuò)增富集cDNA.取其中部分cDNA利用BluePippin進(jìn)行片段的篩選，富集大于4 Kb片段，并將篩選片段進(jìn)行大規(guī)模PCR，以獲得足夠的cDNA總量.對全長cDNA進(jìn)行損傷修復(fù)、末端修復(fù)、連接SMRT啞鈴型接頭、消化去除端末接頭，構(gòu)建未篩選片段和大于4 Kb片段等摩爾混合構(gòu)建混合文庫，結(jié)合引物、綁定DNA聚合酶，形成完整的SMRT bell文庫.

1.2.2 全長轉(zhuǎn)錄組測序及序列分析

文庫檢測合格后運(yùn)用PacBio Sequel平臺(tái)進(jìn)行測序.原始下機(jī)數(shù)據(jù)用軟件包SMRT link處理獲得Subreads序列，校正后獲得環(huán)形一致性序列(circular consensus sequence, CCS)，再根據(jù)序列是否包含5′端引物、3′端引物以及poly(A)將序列分為全長序列(Full-Length Read，F(xiàn)L)與非全長序列；對全長序列進(jìn)行聚類，獲得聚類一致序列，利用二代高質(zhì)量數(shù)據(jù)對三代序列進(jìn)行測序錯(cuò)誤校正后獲得高質(zhì)量的一致性序列進(jìn)行后續(xù)分析.

香港金柑(Fortunellahindsii)基因組已公布，是目前與滑皮金柑親緣關(guān)系最近的物種，因此被定為全長轉(zhuǎn)錄組分析的參考基因組.利用GMAP軟件將兩者序列進(jìn)行比對分析[8].根據(jù)比對結(jié)果，使用Tapis軟件對一致性序列進(jìn)行進(jìn)一步的校正、聚類去冗余得到最終的高質(zhì)量的轉(zhuǎn)錄本(isoform)[9].

1.2.3 全長轉(zhuǎn)錄組的基因結(jié)構(gòu)分析

利用SUPPAV2.3軟件進(jìn)行可變剪接事件分析[10].使用iTAK軟件進(jìn)行植物轉(zhuǎn)錄因子預(yù)測[11].利用Tapis軟件進(jìn)行多聚腺苷酸化分析[9].長鏈非編碼RNA的預(yù)測：利用PLEK和CNCI兩種軟件同時(shí)對轉(zhuǎn)錄本進(jìn)行編碼潛能預(yù)測，后利用CPC軟件將上一步預(yù)測得到的轉(zhuǎn)錄本與NCBI蛋白庫進(jìn)行BLAST對比，進(jìn)一步預(yù)測轉(zhuǎn)錄本[12-14].使用GMAP軟件，以參數(shù)——max-intron length-ends 50000；-f 4；-z sense_force；-n 0進(jìn)行融合基因的挖掘.

1.2.4 抗病基因的挖掘

結(jié)合比較基因組學(xué)和系統(tǒng)發(fā)生學(xué)方法系統(tǒng)挖掘滑皮金柑轉(zhuǎn)錄組中中含NB-ARC結(jié)構(gòu)域抗病蛋白.具體方法是將全長基因組進(jìn)行tblastn比對，篩選e-value值小于1*10-5的轉(zhuǎn)錄組序列進(jìn)行翻譯，獲得含有NBS結(jié)構(gòu)域蛋白.挖掘出的所有的含NBS結(jié)構(gòu)域蛋白和已報(bào)道的1 601條真核及原核生物中的含NBS結(jié)構(gòu)域蛋白序列一起利用hmmalign程序進(jìn)行多重序列比對，得到比對序列[15].利用PhyML3.0[16]對比對文件建立系統(tǒng)發(fā)生關(guān)系，利用外群法確定樹根.系統(tǒng)發(fā)生樹上，與植物的含NB-ARC結(jié)構(gòu)域序列形成單系群的柑橘含NB-ARC結(jié)構(gòu)域序列為柑橘的NBS-LRR類抗病基因.利用HMMER(https://www.ebi.ac.uk/Tools/hmmer/)網(wǎng)站的phmmer功能對抗病蛋白進(jìn)行結(jié)構(gòu)域注釋.

2 結(jié)果與分析

2.1 全長轉(zhuǎn)錄組測序及序列分析

通過測序，獲得原始數(shù)據(jù)13.81 G.過濾原始下機(jī)數(shù)據(jù)，去除接頭和長度小于50 bp的原始下機(jī)數(shù)據(jù)，得到序列，稱為subreads.subreads共有10,429,790條(13.02 Gb)，經(jīng)測定其平均長度為1 249 bp，N50為1 556 bp.通過檢測發(fā)現(xiàn)環(huán)形一致序列CCS共238 054個(gè)，全長非嵌合序列(Full-Length non-chimeric Read, FLNC)197 790個(gè)，占CCS的83.09%，平均長度是1 403 bp，N50為1 653 bp.使用hierarchical n*log(n)算法將同一轉(zhuǎn)錄本的FLNC序列進(jìn)行聚類去冗余再進(jìn)行校正，最終獲得22 759個(gè)高質(zhì)量一致序列用于后續(xù)分析.

利用GMAP軟件將一致序列比對到參考基因組上，比對到基因組未注釋區(qū)域的序列認(rèn)為是新基因位點(diǎn)，比對到已知基因區(qū)域不同外顯子的序列認(rèn)為是新轉(zhuǎn)錄本.經(jīng)比對分析發(fā)現(xiàn)總共12 040個(gè)轉(zhuǎn)錄本，對應(yīng)6 984個(gè)基因，其中已知基因的已知轉(zhuǎn)錄本有5 033條，占比41.80%；已知基因的新轉(zhuǎn)錄本有1 418條，占比11.78%；新基因的新轉(zhuǎn)錄本有5 589條，占比46.42%，對應(yīng)878個(gè)新基因(圖1A：轉(zhuǎn)錄本分類圖)，由此可見，全長轉(zhuǎn)錄組發(fā)現(xiàn)了大量新轉(zhuǎn)錄本，將金柑轉(zhuǎn)錄本數(shù)量提升到原有的2.39倍.分析全長轉(zhuǎn)錄組和參考基因組轉(zhuǎn)錄本的長度數(shù)據(jù)，繪制長度密度圖(圖1B：轉(zhuǎn)錄本密度分布圖).結(jié)果發(fā)現(xiàn)全長轉(zhuǎn)錄組利用PacBio 平臺(tái)超長讀長優(yōu)勢，與參考基因組已注釋轉(zhuǎn)錄本相比得到長度更長的轉(zhuǎn)錄本序列.將參考基因已知基因的轉(zhuǎn)錄本和全長轉(zhuǎn)錄組基因的轉(zhuǎn)錄本數(shù)目進(jìn)行比較分析，繪制成柱狀圖(圖1C：每個(gè)基因?qū)?yīng)轉(zhuǎn)錄本數(shù)目分布圖).結(jié)果發(fā)現(xiàn)全長轉(zhuǎn)錄組發(fā)現(xiàn)更多含有2個(gè)以上轉(zhuǎn)錄本基因，說明全長轉(zhuǎn)錄組結(jié)構(gòu)更完整，與參考基因組比對發(fā)現(xiàn)更多新轉(zhuǎn)錄本，補(bǔ)充基因注釋結(jié)果.將參考基因組已知轉(zhuǎn)錄本和全長轉(zhuǎn)錄組轉(zhuǎn)錄本含有的外顯子數(shù)量進(jìn)行比較(圖1D：每個(gè)轉(zhuǎn)錄本對應(yīng)的外顯子數(shù)目分布圖)，結(jié)果發(fā)現(xiàn)與參考基因組比對發(fā)現(xiàn)更全面的轉(zhuǎn)錄本外顯子數(shù)量，說明全長轉(zhuǎn)錄組可以得到轉(zhuǎn)錄本更真實(shí)結(jié)構(gòu)特征.

圖1 全長轉(zhuǎn)錄組序列分析

2.2 全長轉(zhuǎn)錄組的基因結(jié)構(gòu)分析

可變剪接，也稱為選擇性剪接，是某些基因的一個(gè)mRNA前體通過不同的剪接方式(選擇不同的剪接位點(diǎn))產(chǎn)生不同的mRNA剪接異構(gòu)體的過程.可變剪接是調(diào)節(jié)基因表達(dá)和產(chǎn)生蛋白質(zhì)組多樣性的重要機(jī)制，是導(dǎo)致真核生物基因和蛋白質(zhì)數(shù)量較大差異的重要原因.經(jīng)鑒定，6 984個(gè)基因中共發(fā)現(xiàn)3 031個(gè)基因發(fā)生AS事件，包括1 099個(gè)基因發(fā)生A3(Alternative 3′ splice-site)、601個(gè)基因發(fā)生A5(Alternative 5′ splice site)、569個(gè)基因發(fā)生SE(Skipped exon)、339個(gè)基因發(fā)生RI(Retained intron)、273個(gè)基因發(fā)生AF(Alternative first exon)、123個(gè)基因發(fā)生AL(Alternative last exon)和18個(gè)基因發(fā)生MX(Mutually exclusive exon)(圖2A：基因結(jié)構(gòu)分析圈圖，由外及內(nèi)依次為：染色體序列；可變剪接位點(diǎn)(堆積柱狀圖，不同的可變剪接類型不同的顏色來表示,淺藍(lán)色表示RI，綠色表示A3，黃色表示A5，紫色表示SE，紅色表示MX，棕色表示AF，深藍(lán)色表示AL)；APA位點(diǎn)；新轉(zhuǎn)錄本分布圖 ,顏色越接近于紅色密度越大；新基因分布圖 ,顏色越接近于紅色密度越大；lncRNA密度分布；融合基因，紫色線(相同)，黃色線(不同)染色體上的基因發(fā)生了融合.)，說明可變剪接事件在滑皮金柑基因組中頻繁發(fā)生.

圖2 全長轉(zhuǎn)錄組的基因結(jié)構(gòu)分析圖

可變多聚腺苷酸化是真核生物中一種廣泛存在的基因轉(zhuǎn)錄后調(diào)控過程.APA在mRNA 成熟過程中對前體 mRNA 3′ 端進(jìn)行加工修飾，通過調(diào)控 3′ 非翻譯區(qū) (3′ UTR) 長度而影響 mRNA 穩(wěn)定性、翻譯效率和定位等重要過程[17].利用Tapis軟件對全長轉(zhuǎn)錄組進(jìn)行APA檢測，結(jié)果表明6 106個(gè)基因中至少含有一個(gè)poly(A)位點(diǎn)，其中42.25%的基因含有2個(gè)或2個(gè)以上的poly(A)位點(diǎn)(圖2B：多聚腺苷酸化分布圖)，說明APA在滑皮金柑轉(zhuǎn)錄組中是廣泛分布的.

轉(zhuǎn)錄因子是一群能與基因5′端上游特定序列專一性結(jié)合，從而保證目的基因以特定的強(qiáng)度在特定的時(shí)間和空間表達(dá)的蛋白質(zhì)分子.使用iTAK軟件對轉(zhuǎn)錄組的轉(zhuǎn)錄因子預(yù)測分析，結(jié)果表明：預(yù)測共得596個(gè)轉(zhuǎn)錄因子分布于82個(gè)轉(zhuǎn)錄因子家族，其中轉(zhuǎn)錄因子較多的家族可見圖2C(圖2C：轉(zhuǎn)錄因子分析結(jié)果圖).與植物生長發(fā)育相關(guān)的常見轉(zhuǎn)錄因子中NAC家族有24個(gè)、bHLH家族23個(gè)、AUX/IAA家族18個(gè)、WRKY家族21個(gè)、MYB-related家族35個(gè)，這些轉(zhuǎn)錄因子家族的預(yù)測為之后進(jìn)行金柑生長發(fā)育及次生代謝產(chǎn)物的合成，機(jī)體的調(diào)控提供數(shù)據(jù)基礎(chǔ).

長鏈非編碼RNA是一類轉(zhuǎn)錄本長度超過200 nt，不編碼蛋白質(zhì)的RNA分子.LncRNA因具有非常重要的調(diào)控功能，且?guī)缀鯀⑴c到了各種生物學(xué)過程和通路，與各種疾病的發(fā)生發(fā)展緊密關(guān)聯(lián)，從而成為過去幾年和將來的研究熱點(diǎn)和重點(diǎn).通過PLEK、CNCI和CPC軟件和pfam數(shù)據(jù)庫的嚴(yán)格篩查，我們在全長轉(zhuǎn)錄組中共發(fā)現(xiàn)360個(gè)LncRNA(圖2D：LncRNA預(yù)測結(jié)果維恩圖).統(tǒng)計(jì)出mRNA的數(shù)量為11 680個(gè)，說明滑皮金柑轉(zhuǎn)錄組的RNA中主要是mRNA.根據(jù)LncRNA在基因組上的位置可被分為4類：基因間區(qū)的lncRNA(Intergenic lncRNA，LincRNA)，主要產(chǎn)生于兩個(gè)編碼基因的中間區(qū)域；反義lncRNA (Antisense_lncRNA)，主要產(chǎn)生于編碼基因的反義鏈；內(nèi)含子lncRNA (Sense_intronic lncRNA) 主要產(chǎn)生于編碼基因的內(nèi)含子區(qū)域；Sense_overlaping lncRNA主要位于基因內(nèi)與內(nèi)含子有重疊[18].對分析全長轉(zhuǎn)錄組中的360個(gè)lncRNA類型鑒定，發(fā)現(xiàn)LincRNA為122個(gè)，Antisense_lncRNA有76個(gè)，Sense_intronic lncRNA為130個(gè)和Sense_overlaping lncRNA是32個(gè).

兩個(gè)或多個(gè)基因的編碼區(qū)首尾相連，置于同一套調(diào)控序列(包括啟動(dòng)子、增 強(qiáng)子、核糖體結(jié)合序列、終止子等)控制之下構(gòu)成的嵌合基因稱為融合基因.過去的諸多研究不斷表明，基因融合與各種疾病，特別是癌癥的發(fā)生發(fā)展緊密相關(guān)，甚至是一些癌癥的直接誘因，所以融合基因也成為了當(dāng)前組學(xué)大數(shù)據(jù)分析中的一項(xiàng)重要研究內(nèi)容[10].依據(jù)融合基因鑒定原理，共鑒定出79個(gè)融合基因(圖2A).

2.3 柑橘抗病基因挖掘

通過tblastn軟件獲得6條含有NB-ARC結(jié)構(gòu)域蛋白，將此6序列和代表性1 601條R基因序列進(jìn)行系統(tǒng)發(fā)生學(xué)分析，發(fā)現(xiàn)此6條序列均落在植物R蛋白的多樣性里面，說明此6條序列均為滑皮金柑R基因，值得說明得是其中2條為新基因.利用phmmer程序?qū)基因進(jìn)行結(jié)構(gòu)域注釋，發(fā)現(xiàn)6條R基因蛋白的結(jié)構(gòu)域結(jié)構(gòu)分別為RPW8-NBS、Rx_N-NBS、NBS、NBS、NBS、NBS.全長轉(zhuǎn)錄組中的R基因不能代表滑皮金柑中所有的R基因，而是正在發(fā)生轉(zhuǎn)錄的基因，所以滑皮金柑的R基因遠(yuǎn)不止6個(gè).

3 討論

本項(xiàng)目利用三代結(jié)合二代測序策略對滑皮金柑(金彈種)進(jìn)行全長轉(zhuǎn)錄組測序，將2019年報(bào)道的香港金柑(山金柑種)基因組[19]作為參考基因組進(jìn)行轉(zhuǎn)錄組分析，共得到12,040個(gè)轉(zhuǎn)錄本序列，其中41.8%為已知基因轉(zhuǎn)錄本，46.42%為已知基因的新轉(zhuǎn)錄本，11.78%新基因轉(zhuǎn)錄本，表明我們共挖掘出新轉(zhuǎn)錄本總占比為58.2%，超過香港金柑總轉(zhuǎn)錄本數(shù)，將金柑的總轉(zhuǎn)錄本數(shù)提升到原先的2.15倍.通過轉(zhuǎn)錄本序列分析，全長轉(zhuǎn)錄組利用PacBio Sequel平臺(tái)超長讀長優(yōu)勢，與參考基因組已注釋轉(zhuǎn)錄本相比，本研究獲得長度更長的轉(zhuǎn)錄本序列，具有更高的轉(zhuǎn)錄本密度；更真實(shí)轉(zhuǎn)錄本結(jié)構(gòu)特征，以及更全面的轉(zhuǎn)錄本外顯子數(shù)量.

本研究對全長轉(zhuǎn)錄組的可變剪切、多聚腺苷酸化、 轉(zhuǎn)錄因子、長鏈非編碼RNA和融合基因等5個(gè)方面進(jìn)行結(jié)構(gòu)注釋.我們發(fā)現(xiàn)6 984個(gè)基因中有3 031個(gè)基因發(fā)生了可變剪切事件，占比為43.4%；6 106個(gè)基因至少含有一個(gè)poly(A)位點(diǎn)，且其中含有2個(gè)和2個(gè)以上的poly(A)位點(diǎn)的基因占比為42.3%；596個(gè)轉(zhuǎn)錄因子分布于82個(gè)轉(zhuǎn)錄因子家族；360個(gè)LncRNA和75個(gè)融合基因.結(jié)果說明在滑皮金柑全長轉(zhuǎn)錄組中，可變剪切和多聚腺苷酸化發(fā)生頻繁.我們發(fā)現(xiàn)轉(zhuǎn)錄因子預(yù)測結(jié)果多分布在金柑生長發(fā)育及次生代謝產(chǎn)物合成的表達(dá)部分，說明轉(zhuǎn)錄因子在調(diào)節(jié)金柑生長發(fā)育中發(fā)揮作用.

柑橘作為世界第一大水果，一直以來遭受嚴(yán)重的病害，例如柑橘黃龍病，培育柑橘抗病品種一直科學(xué)家們的研究目標(biāo).R基因在植物抵御病原入侵過程中發(fā)揮重要作用，本項(xiàng)目利用比較基因組學(xué)結(jié)合系統(tǒng)發(fā)生學(xué)方法系統(tǒng)挖掘了全長轉(zhuǎn)錄組中的R基因并對其進(jìn)行結(jié)構(gòu)域注釋.我們共發(fā)現(xiàn)6個(gè)R基因，其中2條序列除含有NBS結(jié)構(gòu)域外還含有rpw8和Rx_N結(jié)構(gòu)域，其他四條序列只含有NBS結(jié)構(gòu)域.全長轉(zhuǎn)錄組中的R基因不能代表滑皮金柑中所有的R基因，而是正在發(fā)生轉(zhuǎn)錄的基因，根據(jù)此思路可初步篩選和某一病害相關(guān)抗病基因，為滑皮金抗病機(jī)制和抗病品種的培育打下基礎(chǔ).