閩楠轉(zhuǎn)錄組分析及基因功能注釋

劉 丹 曾欽朦 劉 斌 李 煜 陳世品

(福建農(nóng)林大學(xué)林學(xué)院，福州 350002)

閩楠(Phoebebournei(Hemsl.) Yang)俗稱楠木，為樟科(Lauraceae)楠屬(Phoebe)植物，是我國珍貴的用材樹種之一，具有較高的經(jīng)濟(jì)和生態(tài)價值[1]。天然閩楠現(xiàn)存野生資源日益減少，且受到其自身生物學(xué)特性影響，生長發(fā)育緩慢，在遺傳育種方面也大多采用傳統(tǒng)育種方法，周期長且效率低[2～4]。研究閩楠群體遺傳多樣性，不僅可以避免物種滅絕，還可以為遺傳群體的管理提供參考。目前為止，閩楠遺傳多樣性的研究相對較少，僅有李娟等[5]對閩楠5個野生群體進(jìn)行ISSR分析，利用10對引物對96個樣本的遺傳多樣性和遺傳結(jié)構(gòu)進(jìn)行研究，得出閩楠在基因水平上遺傳多樣性較高，種群間遺傳變異強(qiáng)于種群內(nèi)。江香梅等[6]利用RAPD分子標(biāo)記對福建和江西兩省的8個天然群體進(jìn)行遺傳多樣性分析，得出其遺傳變異主要存在于種群內(nèi)，并提出通過種群擴(kuò)繁和回歸自然的方式增強(qiáng)種群間的基因流，以維持群體遺傳多樣性水平。通過轉(zhuǎn)錄組測序開發(fā)分子標(biāo)記，可以研究閩楠的群體遺傳多樣性，對閩楠群體管理和種質(zhì)資源保護(hù)具有重要意義，也可以用于閩楠分子標(biāo)記輔助育種研究，獲得的轉(zhuǎn)錄組信息將為關(guān)鍵基因功能的研究提供基礎(chǔ)。

基于二代測序技術(shù)的轉(zhuǎn)錄組測序技術(shù)(RNA-seq)，具有高準(zhǔn)確性、快速和低成本等優(yōu)勢，能夠在缺乏基因資源的條件下獲得物種的生長和代謝規(guī)律，并揭示其生物學(xué)特性與基因內(nèi)在關(guān)聯(lián)，同時可獲得物種絕大多數(shù)的轉(zhuǎn)錄產(chǎn)物。通過轉(zhuǎn)錄組測序開發(fā)的分子標(biāo)記與傳統(tǒng)標(biāo)記相比，成本低廉，簡單迅速，多態(tài)性高且在植物物種間可轉(zhuǎn)移，不僅可以大大提高工作效率，還可以對基因結(jié)構(gòu)與功能進(jìn)行更深層次的分析與探究，是一種性價比高的基因序列研究手段[7～9]。目前，轉(zhuǎn)錄組測序技術(shù)已經(jīng)在林木遺傳育種中得到廣泛的應(yīng)用。例如：李清瑩[10]等通過對火力楠進(jìn)行轉(zhuǎn)錄組測序，獲得了97 503個Unigene基因；田曉明[11]等對半楓荷葉片進(jìn)行轉(zhuǎn)錄組測序，共獲得93 602條Unigenes，比對鑒定獲得19 279個基因，經(jīng)KEGG分析，有16 621個Unigenes參與了131條代謝通路；趙德剛[12]等對杜仲雌雄株進(jìn)行轉(zhuǎn)錄組測序，獲得平均長度為288 nt的雌株159 434個Unigene片段和平均長度為231 nt的雄株257 288個Unigene片段，共有48 761個表達(dá)序列標(biāo)簽(EST)。

目前國內(nèi)有關(guān)楠屬開展了浙江楠(P.chekiangensis)[13]；楨楠(P.zhennan)[14]；紫楠(P.sheareri)[15]；白楠(P.neurantha)[16]的轉(zhuǎn)錄組測序分析研究工作，但未見閩楠基因組和轉(zhuǎn)錄組的相關(guān)報道，且未發(fā)現(xiàn)其他楠屬樹種開發(fā)的分子標(biāo)記在閩楠上進(jìn)行過通用性驗證，遺傳信息匱乏。為此，結(jié)合閩楠優(yōu)質(zhì)木材的特性，選取3個重要部位的組織，包括木質(zhì)部、韌皮部和葉片，使用二代高通量測序技術(shù)對其進(jìn)行轉(zhuǎn)錄組測序，獲取大量閩楠轉(zhuǎn)錄本信息，利用生物信息分析從海量的轉(zhuǎn)錄本數(shù)據(jù)中篩選候選基因，以期初步了解閩楠基因表達(dá)與分布情況，為深入開展閩楠遺傳育種及分子生物學(xué)相關(guān)研究奠定基礎(chǔ)。

1 材料與方法

1.1 材料

2017年6月，在福建農(nóng)林大學(xué)校園內(nèi)選取1株約40年生的閩楠作為試驗植株，該植株生長健康、樹干通直、長勢良好且無病蟲害。采集試驗植株的木質(zhì)部、韌皮部、葉片，洗凈后經(jīng)液氮快速冷凍，于-80℃冰箱中保存?zhèn)溆谩?/p>

1.2 方法

1.2.1 總RNA提取

分別使用PureLinkTMPlant RNA試劑盒(Invitrogen)提取植物總RNA，再用RQ1 RNase-Free DNase(Promega)去除總RNA中的雜質(zhì)DNA。采用1.25%瓊脂糖凝膠電泳檢測RNA完整性，利用Agilent 2100核酸蛋白檢測儀(Agilent,Santa Clara，CA)檢測3個RNA樣品的濃度和完整度。

1.2.2 文庫構(gòu)建及測序

通過Oligo(dT)磁珠富集總RNA中帶有polyA結(jié)構(gòu)的mRNA，采用離子打斷的方式，將mRNA打斷到200～300 bp片段。以mRNA為模板，用6堿基隨機(jī)引物和逆轉(zhuǎn)錄酶合成cDNA第一鏈，并以第一鏈cDNA為模板進(jìn)行第二鏈cDNA的合成，第二鏈cDNA合成時，其中的堿基T被替換成U，從而達(dá)到鏈特異性文庫的目的。PCR擴(kuò)增富集文庫片段，擴(kuò)增完成后進(jìn)行文庫片段大小的選擇，文庫大小在300～400 bp。通過Agilent 2100 Bioanalyzer對文庫進(jìn)行質(zhì)檢，再進(jìn)行文庫總濃度檢測及有效文庫濃度檢測。采用第二代測序技術(shù)(Next-Generation Sequencing，NGS)，基于Illumina HiSeq測序平臺，對這些文庫進(jìn)行雙末端(Paired-end，PE)測序。

1.2.3 數(shù)據(jù)分析

利用FastQC軟件(http://www.bioinformatics.babraham.ac.uk/projects/fastqc)進(jìn)行測序數(shù)據(jù)的質(zhì)量評估，采用Cutadapt[17](Version 1.2.1)去除3′端的接頭序列，使用Trinity[18](版本r20140717，K-mer 25 bp)軟件進(jìn)行轉(zhuǎn)錄組De Novo組裝，利用轉(zhuǎn)錄組表達(dá)定量軟件RSEM[19]計算每個基因的FPKM值，采用DEGSeq[20]分析表達(dá)差異Unigene，用Varscan[21](Version 2.3.7)程序獲取cSNP位點(diǎn)。使用MISA(http://pgrc.ipk-gatersleben.de/misa/misa.html)程序在序列中搜索SSR位點(diǎn)。

2 結(jié)果與分析

2.1 測序數(shù)據(jù)過濾

將RNA片段化，建立PCR富集文庫，文庫插入長度為380 bp，采用Illumina MiSeq測序平臺測序，測序模式為雙末端的2×150 bp(見表1)。測序過濾后的數(shù)據(jù)平均質(zhì)量較好，可以進(jìn)行后續(xù)分析。

2.2 轉(zhuǎn)錄本拼接結(jié)果

提取每個基因下最長的轉(zhuǎn)錄本作為該基因的代表序列，稱為Unigene。對拼接得到的Contig、Transcript和Unigene序列進(jìn)行統(tǒng)計(見表2)。

2.3 Unigene功能注釋結(jié)果

對聚類得到的Unigene進(jìn)行基因功能注釋(見表3)?；蚬δ茏⑨屗玫降臄?shù)據(jù)庫包括NR(NCBI non-redundant protein sequences)、GO(Gene Ontology)[22]、KEGG(Kyoto Encyclopedia of Genes and Genome)[23]、eggNOG(evolutionary genealogy of genes:Non-supervised Orthologous Groups)[24]、Swiss-Prot。由表3可知，NR數(shù)據(jù)庫注釋了45 036條Unigene，獲得的信息最多，KEGG數(shù)據(jù)庫獲得的信息最少，僅有5 623條Unigene得到注釋，在所有數(shù)據(jù)庫中均被注釋到的Unigene有4 646條，僅占3.06%。

表1 測序數(shù)據(jù)統(tǒng)計

注：Q30(bp)：堿基識別準(zhǔn)確率在99.9%以上的堿基總數(shù)；N(%)：模糊堿基所占百分比；Q20(%)：堿基識別準(zhǔn)確率在99%以上的堿基所占百分比；Q30(%)：堿基識別準(zhǔn)確率在99.9%以上的堿基所占百分比；GC(%)：GC含量

Q30(bp). The total number of bases whose recognition accuracy is above 99.9%; N(%). Percentage of fuzzy bases; Q20(%). The percentage of bases whose recognition accuracy is above 99%; Q30(%). The percentage of base recognition accuracy above 99.9%;GC(%). GC content

表2 轉(zhuǎn)錄本拼接結(jié)果

注：N50(bp).將所有序列從長到短排列，將序列長度按照該順序依次相加，當(dāng)相加的長度達(dá)到序列總長度的50%時，最后一條序列的長度；N90(bp).將所有序列按照長度從長到短排列，將序列長度按照該順序依次相加，當(dāng)相加的長度達(dá)到序列總長度的90%時，最后一條序列的長度；N50 Sequence No..長度大于N50的序列總數(shù)；N90 Sequence No..長度大于N90的序列總數(shù)；GC(%).序列的GC含量

Note:N50(bp).All sequences are arranged from long to short,and the sequence lengths are added in this order. When the added length reaches 50% of the total sequence length,the length of the last sequence is; N90(bp).All sequences are arranged from long to short in length,and the sequence length is added in this order. When the added length reaches 90% of the total sequence length,the length of the last sequence is; N50 Sequence No..Total number of sequences longer than N50; N90 Sequence No..Total number of sequences longer than N90; GC(%).The GC content of the sequence

表3 Unigene功能注釋結(jié)果

圖1 NR注釋結(jié)果統(tǒng)計圖Fig.1 NR annotated species distribution of Unigenes of transcriptome for P.bournei

2.3.1 閩楠轉(zhuǎn)錄組Unigene的NR功能分類

通過Blastx將閩楠轉(zhuǎn)錄組中的Unigene與NR數(shù)據(jù)庫進(jìn)行比對，可以獲取本物種基因序列與近緣物種基因序列的相似性以及本物種基因的功能信息，其中有45 036條Unigene有注釋結(jié)果。通過NR庫比對注釋的結(jié)果，統(tǒng)計并繪制比對上的物種分布圖(見圖1)。閩楠與黃瓜、野草莓、大豆、歐洲大葉楊、碧桃、蓖麻、葡萄都有同源序列分布，其中與葡萄有34%的同源序列，與歐洲大葉楊、碧桃、蓖麻分別有8%、7%、6%的相似序列，在黃瓜、野草莓、大豆中的同源序列較少，均為3%。此外，還有35%的Unigene屬于其他序列，可能包含了閩楠自身特有的與大多數(shù)物種不同的序列。

表4 閩楠轉(zhuǎn)錄組GO功能分類的Unigene數(shù)量分布

圖2 閩楠轉(zhuǎn)錄組Unigene的GO功能分類Fig.2 GO analysis of P.bournei transcriptome Unigene

2.3.2 閩楠轉(zhuǎn)錄組Unigene的GO分類注釋

對閩楠轉(zhuǎn)錄組Unigene進(jìn)行GO功能分類，有26 463條Unigene獲得了基因注釋，占總Unigene數(shù)的17.44%。將注釋的到的Unigene劃分為生物過程、細(xì)胞成分、分子功能3大類共計52個分支(見圖2)，對每一類的基因數(shù)量進(jìn)行統(tǒng)計發(fā)現(xiàn)(見表4)，在生物過程這一類中，代謝過程、細(xì)胞過程、單一有機(jī)體過程占比較大；在細(xì)胞成分一類中，細(xì)胞、細(xì)胞器、膜、細(xì)胞組分占比較大；在分子功能一類中，催化活性和連接占比較大。

2.3.3 閩楠轉(zhuǎn)錄組Unigene的eggNOG功能注釋

將閩楠轉(zhuǎn)錄組的Unigene與eggNOG數(shù)據(jù)庫進(jìn)行比對，可得到43 388個Unigene具有功能信息，根據(jù)功能可將其劃分為25類(見圖3:A～W，Y，Z)。由圖3可知，一般功能預(yù)測(R)的Unigene數(shù)量最多，其次為未知功能(S)和信號傳導(dǎo)機(jī)理(T)，細(xì)胞動力(N)最少，僅有7條Unigene注釋。

通過對每一類eggNOG注釋的基因功能進(jìn)行統(tǒng)計，可獲得閩楠轉(zhuǎn)錄組中各類功能基因的數(shù)量分布信息(見表5)。其中，此生代謝產(chǎn)物生物合成、運(yùn)輸及代謝(Q)有1 617條Unigene。

2.3.4 閩楠轉(zhuǎn)錄組Unigene的KEGG代謝途徑分類分析

根據(jù)KEGG數(shù)據(jù)庫的注釋信息進(jìn)一步將閩楠Unigene進(jìn)行pathway注釋，共有5 623條Unigene獲得注釋，參與的代謝通路可歸為5個大類、35個子類，又可分為176個信號代謝分支(見圖4)。其中，新陳代謝(Metabolism)一類中獲得注釋最多的是碳代謝(Carbon metabolism)，有164條Unigene；遺傳信息處理(Genetic Information Processing)一類中核糖體(Ribosome)獲得注釋最多，有221條Unigene；環(huán)境信息處理(Environmental Information Processing)一類中獲得注釋最多的是植物激素信號轉(zhuǎn)導(dǎo)(Plant hormone signal transduction)，有118條Unigene，細(xì)胞進(jìn)程(Cellular Processes)和有機(jī)系統(tǒng)(Organismal Systems)兩類中獲得注釋最多的分別為內(nèi)吞作用(Endocytosis)96條、生成信號通路(Neurotrophin signaling pathway)90條(見表6)。

表6 閩楠轉(zhuǎn)錄組的KEGG代謝途徑(>1.5%注釋基因)

圖3 閩楠轉(zhuǎn)錄組Unigene的eggNOG功能分類 A.RNA加工與修飾；B.染色質(zhì)結(jié)構(gòu)與動力；C.能量生成與轉(zhuǎn)化；D.細(xì)胞周期調(diào)控，細(xì)胞分裂，染色體分配；E.氨基酸運(yùn)輸及代謝；F.核苷酸運(yùn)輸及代謝；G.糖類運(yùn)輸及代謝；H.輔酶運(yùn)輸及代謝；I.脂類運(yùn)輸及代謝；J.翻譯，核糖體結(jié)構(gòu)和生物合成；K.轉(zhuǎn)錄；L.復(fù)制，重組與修復(fù)；M.細(xì)胞壁生源/細(xì)胞膜生源/囊膜生源；N.細(xì)胞動力；O.翻譯后修飾，蛋白質(zhì)轉(zhuǎn)換與分子伴侶；P.無機(jī)離子的運(yùn)輸及代謝；Q.次生代謝產(chǎn)物生物合成，運(yùn)輸及代謝；R.一般功能(預(yù)測)；S.未知功能；T.信號轉(zhuǎn)導(dǎo)機(jī)理；U.胞內(nèi)運(yùn)輸，分泌及小泡運(yùn)輸；V.防御機(jī)制；W.胞外結(jié)構(gòu)；Y.細(xì)胞核結(jié)構(gòu)；Z.細(xì)胞骨架Fig.3 eggNOG analysis of P.bournei transcriptome Unigene A. RNA processing and modification; B. Chromatin structure and dynamics; C.Energy production and conversion; D. Cell cycle control,cell division,chromosome partitioning; E. Amino acid transport and metabolism; F. Nucleotide transport and metabolism; G.Carbohydrate transport and metabolism; H. Coenzyme transport and metabolism; I. Lipid transport and metabolism; J.Translation,ribosomal structure and biogenesis; K.Transcription; L. Replication,recombination and repair; M.Cell wall/membrane/envelope biogenesis; N. Cell motility; O.Posttranslational modification,protein turnover,chaperones; P. Inorganic ion transport and metabolism; Q.Secondary metabolites biosynthesis,transport and catabolism; R.General function prediction only;S.Function unknown;T.Signal transduction mechanisms;U.Intracellular trafficking,secretion,and vesicular transport; V.Defense mechanisms; W.Extracellular structures;Y.Nuclear structure; Z.Cytoskeleton

圖4 KEGG注釋統(tǒng)計圖Fig.4 KEGG annotation statistical chart

表7 閩楠轉(zhuǎn)錄組SSR搜索結(jié)果

2.4 閩楠轉(zhuǎn)錄組Unigene的SSR分布特征

SSR(Simple Sequence Repeats)，稱為短串聯(lián)重復(fù)序列或微衛(wèi)星標(biāo)記[25]，通過對閩楠的151 729條Unigene進(jìn)行SSR位點(diǎn)搜尋，共檢測到35 972個SSR位點(diǎn)，SSR重復(fù)類型從單核苷酸到六核苷酸均有不同數(shù)量的分布，所占比例有較大差異，其中單核苷酸數(shù)量最多，占60.5%，六核苷酸僅占0.04%(見表7)。

3 討論

本研究對閩楠的木質(zhì)部、韌皮部、葉片3個不同部位的組織進(jìn)行高通量測序，獲得木質(zhì)部41 383 707條、韌皮部43 343 922條、葉片44 191 586條Clean Reads片段，經(jīng)轉(zhuǎn)錄本拼接后得到序列總長度達(dá)120 535 288bp的383 331條Conting片段，進(jìn)一步組裝得到平均長度為542 bp的151 729條Unigenes。相比于其他利用相同測序手段的物種來說長度較長，如茶(Camelliasinensis)(平均長度355 bp)[26]、杉木(Cunninghamialanceolata)(平均長度449 bp)[27]、棗(Ziziphuscelata)(平均長度408 bp)[28]等。在Swiss-Prot數(shù)據(jù)庫中共檢測到34 885條蛋白同源序列，與NR數(shù)據(jù)庫進(jìn)行比對發(fā)現(xiàn)，與葡萄科葡萄屬的葡萄(Vitisvinifera)比對的同源信息最多，占34%，出現(xiàn)此種情況，可能是由于葡萄與閩楠的進(jìn)化史和生活史較為相似，也可能是因為葡萄具有參考基因組；與黃瓜、野草莓、大豆等物種相比同源性較低，僅有3%，另有35%的Unigene未在數(shù)據(jù)庫中找到同源序列，可能是因為閩楠本身特有的基因存在，也可能是因為這些序列是一些非編碼的RNA序列，或是因為長度較短不包含蛋白質(zhì)功能域信息，也有可能是目前公共數(shù)據(jù)庫中的基因信息還不夠完善[29]。從進(jìn)化關(guān)系上來看，不同科屬的植物，其具有同源序列的數(shù)量較低，閩楠符合這一特點(diǎn)。將閩楠轉(zhuǎn)錄組Unigene注釋到GO數(shù)據(jù)庫，可劃分為生物過程、細(xì)胞成分、分子功能3大類共計52個分支，與eggNOG數(shù)據(jù)庫比對可分為25類，通過KEGG功能注釋可知轉(zhuǎn)錄組中涉及的基因共參與了176條代謝通路，其中涉及代謝通路和次生物質(zhì)的生物合成基因較多。此外，轉(zhuǎn)錄組序列中共檢測到35 972個SSR位點(diǎn)，其中單核苷酸占比最大，達(dá)到60%以上。

丁亞軍[13]等在2014年對浙江楠葉、木質(zhì)部、韌皮部、花、花柄5個組織的轉(zhuǎn)錄組進(jìn)行測序，共獲得了52 527 954條Reads片段，經(jīng)組裝拼接后得到平均長度為711.25 bp的111 250個Unigenes，在NR、Swiss-Prot、GO、KEGG等數(shù)據(jù)庫中均獲得較好的注釋。有3大類60個分支的36 370條Unigenes成功注釋到GO數(shù)據(jù)庫，有25類15 605條Unigenes具有具體的蛋白功能定義，16 135條Unigenes在KEGG數(shù)據(jù)庫映射到298條代謝通路中。值得注意的是，在與NR數(shù)據(jù)庫比對中發(fā)現(xiàn)，與浙江楠的Unigene同源信息最高的是葡萄科葡萄樹的葡萄，占比62%，其次為楊柳科的毛果楊和鐵莧菜亞科的蓖麻，各占12%，與豆科的青仁烏豆和苜蓿共有8%相似，這與閩楠在NR數(shù)據(jù)庫中的同源序列比對結(jié)果有一定的相似性，但也存在差異。由此可見，閩楠和浙江楠雖同為樟科楠屬植物，但是在基因的同源性分布上仍有不同，一方面可能是因為轉(zhuǎn)錄組測序結(jié)果的時空異質(zhì)性較高，不僅不同物種差異較大，同物種，乃至同植株在不同階段的測序差異結(jié)果也會較大，我們在轉(zhuǎn)錄組測序時用于提取RNA選取的組織來自于不同部位，導(dǎo)致了不同組織器官在不同時期具有不同的基因表達(dá)，另一方面可能是因為與其他物種進(jìn)行比對具有一定的局限性。

本研究通過生物信息學(xué)方法對閩楠進(jìn)行轉(zhuǎn)錄本拼接、基因功能注釋、SNP檢測以及SSR位點(diǎn)預(yù)測等，獲得的轉(zhuǎn)錄組數(shù)據(jù)不僅可以作為樟科楠屬樹種遺傳信息的重要組成部分，進(jìn)一步豐富該科樹種的基因數(shù)據(jù)庫，同時也為閩楠基因的分子克隆和功能鑒定提供數(shù)據(jù)支持，除了可以發(fā)掘候選基因外，還可以開發(fā)EST-SSR分子標(biāo)記，為后續(xù)閩楠的遺傳多樣性分析奠定基礎(chǔ)。