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

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

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

1 材料及方法

1.1 實驗材料

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

1.2 方法

1.2.1 總RNA樣品檢測及文庫的構建

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

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

1.2.2 全長轉錄組測序及序列分析

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

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

1.2.3 全長轉錄組的基因結構分析

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

1.2.4 抗病基因的挖掘

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

2 結果與分析

2.1 全長轉錄組測序及序列分析

通過測序，獲得原始數(shù)據(jù)13.81 G.過濾原始下機數(shù)據(jù)，去除接頭和長度小于50 bp的原始下機數(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個，全長非嵌合序列(Full-Length non-chimeric Read, FLNC)197 790個，占CCS的83.09%，平均長度是1 403 bp，N50為1 653 bp.使用hierarchical n*log(n)算法將同一轉錄本的FLNC序列進行聚類去冗余再進行校正，最終獲得22 759個高質量一致序列用于后續(xù)分析.

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

圖1 全長轉錄組序列分析

2.2 全長轉錄組的基因結構分析

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

圖2 全長轉錄組的基因結構分析圖

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

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

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

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

2.3 柑橘抗病基因挖掘

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

3 討論

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

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

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