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

(福建農林大學林學院，福州 350002)

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

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

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

1 材料與方法

1.1 材料

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

1.2 方法

1.2.1 總RNA提取

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

1.2.2 文庫構建及測序

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

1.2.3 數據分析

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

2 結果與分析

2.1 測序數據過濾

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

2.2 轉錄本拼接結果

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

2.3 Unigene功能注釋結果

對聚類得到的Unigene進行基因功能注釋(見表3)。基因功能注釋所用到的數據庫包括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數據庫注釋了45 036條Unigene，獲得的信息最多，KEGG數據庫獲得的信息最少，僅有5 623條Unigene得到注釋，在所有數據庫中均被注釋到的Unigene有4 646條，僅占3.06%。

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

注：Q30(bp)：堿基識別準確率在99.9%以上的堿基總數；N(%)：模糊堿基所占百分比；Q20(%)：堿基識別準確率在99%以上的堿基所占百分比；Q30(%)：堿基識別準確率在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 轉錄本拼接結果

注：N50(bp).將所有序列從長到短排列，將序列長度按照該順序依次相加，當相加的長度達到序列總長度的50%時，最后一條序列的長度；N90(bp).將所有序列按照長度從長到短排列，將序列長度按照該順序依次相加，當相加的長度達到序列總長度的90%時，最后一條序列的長度；N50 Sequence No..長度大于N50的序列總數；N90 Sequence No..長度大于N90的序列總數；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功能注釋結果

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

2.3.1 閩楠轉錄組Unigene的NR功能分類

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

表4 閩楠轉錄組GO功能分類的Unigene數量分布

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

2.3.2 閩楠轉錄組Unigene的GO分類注釋

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

2.3.3 閩楠轉錄組Unigene的eggNOG功能注釋

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

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

2.3.4 閩楠轉錄組Unigene的KEGG代謝途徑分類分析

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

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

圖3 閩楠轉錄組Unigene的eggNOG功能分類 A.RNA加工與修飾；B.染色質結構與動力；C.能量生成與轉化；D.細胞周期調控，細胞分裂，染色體分配；E.氨基酸運輸及代謝；F.核苷酸運輸及代謝；G.糖類運輸及代謝；H.輔酶運輸及代謝；I.脂類運輸及代謝；J.翻譯，核糖體結構和生物合成；K.轉錄；L.復制，重組與修復；M.細胞壁生源/細胞膜生源/囊膜生源；N.細胞動力；O.翻譯后修飾，蛋白質轉換與分子伴侶；P.無機離子的運輸及代謝；Q.次生代謝產物生物合成，運輸及代謝；R.一般功能(預測)；S.未知功能；T.信號轉導機理；U.胞內運輸，分泌及小泡運輸；V.防御機制；W.胞外結構；Y.細胞核結構；Z.細胞骨架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 閩楠轉錄組SSR搜索結果

2.4 閩楠轉錄組Unigene的SSR分布特征

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

3 討論

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

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

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