云南金花茶轉錄組序列分析及功能注釋

， ， ， ，b，，培堯

（1.西南林業大學 a.西南山地森林資源保育與利用教育部重點實驗室；b.國家林業和草原局西南地區生物多樣性保育重點實驗室，云南 昆明 650224；2.畢節市林業科學研究所，貴州 畢節 551700；3.云南省大圍山國家級自然保護區河口管理分局，云南 河口 661399）

云南金花茶Camellia fascicularisH.T.Chang為山茶科Theaceae 山茶屬Camellia金花茶組的常綠木本植物[1]。其現有野生資源僅600 株左右，屬于我國二級保護植物，同時也被云南省認定為特有極小種群植物[2]。云南金花茶間斷分布于河口及其周邊個舊市和馬關縣的溝谷或山坡海拔400 ～1 000 m 濕潤肥沃的土壤中，耐貧瘠[3]。云南金花茶花形嬌艷多姿，呈金黃色，除具有園林觀賞價值外，還具有降血脂、抗氧化、抑制癌細胞等功效[4-6]。近年來，有關云南金花茶的研究報道主要集中在繁育技術、離體培養體系建立、保健功效、營養成分分析以及化學元素研究[7-11]等方面，而在分子水平的研究少有報道，僅見本課題組利用SSR 分子標記技術對云南金花茶遺傳多樣性的分析[12]。

隨著分子生物學相關知識和技術的快速升級，利用轉錄組測序（RNA-seq）技術獲得基因信息，對植物進行生物信息學分析的方法得到廣泛應用。使用RNA-seq 技術可在獲得大量物種注釋序列信息的同時，挖掘出生物重要的功能基因，因此該技術成為研究缺乏功能組學和基因組學信息的物種優良性狀和基因功能的重要手段[13-15]。伴隨高通量測序技術的不斷進步與完善，RNA-seq 成本不斷降低，核苷酸的檢測也更精準和快捷[16-19]。目前，該技術在前基因組學和后基因組學的研究中均有應用。魏開發等[20]對火龍果Hylocereus undatus進行RNA-seq，研究不同發育階段中花芽、果實和枝條基因的表達情況，并將組裝得到的Unigene 進行注釋，結果發現，在火龍果不同組織中分別有多個特異表達的Unigenes。宋猜等[21]對采于花芽萌動前后3 個發育時期的仁用杏Armeniaca vulgaris花芽進行轉錄組測序，試驗數據可提供大量仁用杏開花相關基因的信息，可為研究仁用杏成花分子機制及解決仁用杏花期凍害問題提供理論依據。張賢等[22]通過RNA-seq 對芒草Miscanthus sinensis的74 134 條Unigenes 進 行不同功能領域數據庫注釋，其結果可為進一步鑒定中國芒草的功能基因提供參考。另有學者對百香果Passiflora edulis[23]、云錦杜鵑Rhododendron fortunei[24]、桐花樹Aegiceras corniculatum[25]進行了相關研究。目前，僅見葉鵬等[12]對云南金花茶轉錄組微衛星序列的分布與特征進行了探討，而有關其基因功能注釋、代謝途徑和代謝通路方面的分析未見報道。本研究中對云南金花茶進行高通量轉錄組測序，獲得大量原始數據的同時，對其進行拼接與組裝處理，建立RNA-seq 數據庫，并將處理后的數據在7 大公共數據庫進行功能注釋、代謝途徑和相關通路分析，以期為云南金花茶乃至山茶屬植物功能基因的探索提供一定的理論參考。

1 材料與方法

1.1 試驗材料

在云南省河口縣（103°97′E，22°52′N）海拔1 036 m 的陽坡地帶采集云南金花茶植株，引種到西南林業大學溫室大棚中。選用1 份云南金花茶幼嫩葉片作為試驗材料，用錫箔紙包好放入液氮中，備用。

1.2 試驗方法

1.2.1 轉錄組測序

首先提取云南金花茶RNA，然后構建cDNA數據庫，繼而利用Illumina HiseqTM 2000 平臺對云南金花茶進行轉錄組測序。此部分工作由北京諾禾源科技股份有限公司完成。

1.2.2 序列組裝

RNA-seq 測序完成后，統計raw reads 的數量和長度。原始數據中含有低質量序列、重復冗余序列、接頭和無法確定堿基信息的序列，必須將上述序列去除，以獲得clean reads，繼而統計clean reads 的數量、總長度、處理后不確定序列所占比例、GC 堿基所占比例、N50（拼接轉錄本不小于總長50%的長度）以及Q20（處理后質量高于20的堿基）所占比例等。對clean reads 通過Trinity Software（http:///trinityrnaseq.Github.io/） 進 行de novo 組裝。首先利用clean reads 之間的overlap 將其向兩邊伸展形成序列重疊克隆群（contig），再依據序列雙末端的信息對contig 進行再次連接，得到該樣品的Transcript，去除Transcript 冗余reads 獲得Unigene 后，進行Transcript 和Unigene的分布和長度分析[26]。

1.2.3 基因功能注釋、分類及生物學通路分析

將處理得到的Unigene 在7 個不同功能領域的公共數據庫中進行基因功能注釋和分類分析，從而獲得較全面的云南金花茶基因功能信息。數據庫包括：Nr（Non-Redundant Protein Database，非冗余蛋白數據庫）、Nt（Nucleotide Sequence Database，核酸序列數據庫）、Pfam（Pfam Protein Sequence Database，Pfam 蛋白序列數據庫）、KOG（euKaryotic Ortholog Groups，真核生物蛋白直系同源數據庫）、Swiss-Prot（Swiss-ProtProtein Database，Swiss-Prot 蛋白質序列數據庫）、KEGG（Kyoto Encyclopedia of Genes and Genomes，基因組百科全書）、GO（Gene Ontology，基因本體數據庫）[27]。使用BLAST 軟件將Unigene 在Nr、Nt、Swiss-Prot 等數據庫中進行比對（e-value ＜1×10-5），獲取相關基因注釋。比對到Nr 數據庫中[28]，從而獲取云南金花茶基因序列相似性和物種分布信息。依據Nr 中注釋的結果，在Blast2GO 數據庫比對，得到GO 功能注釋信息[29]。GO 數據庫包括3 大類別，分別為生物過程、分子功能與細胞組分，以此可以宏觀解讀云南金花茶基因功能的分布及特征[30]。將Unigene比對到KOG 數據庫中，并按可能的功能對獲得結果的Unigene 進行分類與統計；另外，對Unigene進行KEGG 數據庫相關通路（包括細胞過程、遺傳信息處理、新陳代謝、環境信息處理、有機系統5 大類別）分析，了解云南金花茶的代謝通路以及各通路之間的關系[31]。

1.2.4 云南金花茶轉錄組Unigene 的CDS 預測

將Unigene 序列依次比對到Nr（https://www.ncbi.nlm.nih.gov/）、Swiss-Prot（http://www.ebi.ac.uk/uniprot/）、KEGG（http://www.genome.jp/kegg/）、KOG（http://www.ncbi.nlm.nih.gov/COG/）等蛋白數據庫中，對于未比對上或未預測到結果的序列，使用ESTScan（3.0.3）軟件進行預測。

2 結果與分析

2.1 云南金花茶RNA-seq 及de novo 組裝結果

通過RNA-seq，共得到云南金花茶57 051 836條原始序列。將原始序列中的接頭（dadpter）、低質量reads、重復冗余以及不確定堿基含量超過10%的讀序經處理后，獲得54 817 600條有效序列，總長為8.22 Gb，Q20、Q30（處理后質量高于30的堿基）高質量序列分別占96.39%和91.28%，GC 含量占總堿基數的44.54%，堿基錯誤率為0.02%，說明由高通量測序平臺獲得了較高數量和質量的云南金花茶序列，有利于后續數據的組裝，滿足后期生物信息學的研究。得到的clean reads經de novo 組裝后，共獲得155 011 條Transcript，這些Transcript 經進一步組裝之后，得到95 979 條Unigenes，序列信息達107 907 727 nt。對Transcript的序列長度分析結果表明，其平均長度是807 nt，N50是1 411 nt。其中，以200～500 nt的短序列居多，有85 904 條，占總數的55.42%；500 ～1 000 nt長度的序列為30 871 條，占總數的19.92%；1 000 ～2 000 nt 長度的序列為23 853 條，占總數的15.39%；大于等于2 000 nt 長度的序列占總數的9.28%（圖1A）。Unigene 分析統計結果表明，其平均長度為1 124 nt，N50 為1 660 nt，其中1 000 ～2 000 nt 的序列占總序列的24.84%，超過2 000 nt 的序列占14.99%。通過對高通量RNA-seq得到的大量序列進行處理，經組裝后Unigenes 數據的完整性明顯提高，可進行下一步的分析統計（圖1B）。

2.2 云南金花茶轉錄組Unigene 的功能注釋及分類

將獲得的95 979 條Unigenes 通過BLAST軟件在7 大數據庫進行比對， 共有63 888（66.56%）條Unigenes 獲得注釋。其中，在Nr（e-value ≤1×10-5）注釋成功 的Unigenes 有58 830 條，占Unigenes 總數量的61.29%；在Nt（e-value ≤1×10-5）注釋成功43 623 條，占總數的45.45%；在KEGG（e-value ≤1×10-10）注釋成功的有23 214 條，占總數的24.18%；在Swiss-Prot（e-value ≤1×10-5）注釋成功的有44 315 條，占總數的46.17%；在Pfam（e-value ≤0.01）注釋成功的有41 096 條，占總數的42.81%；在GO（e-value ≤1×10-6）注釋成功的有41 905 條，占總數的43.66%；在KOG（e-value ≤1×10-3）注釋成功的有23 499 條，占總數的24.48%。在7 大數據庫中均能得到成功注釋的序列數目為11 933條，占總數的12.43%，其中63 888 條序列至少在1 個數據庫中注釋成功，占總數的66.56%。

圖1 云南金花茶轉錄組組裝序列長度分布Fig.1 Length distribution of assembly transcript andunigenes for C.fascicularis

2.2.1 云南金花茶轉錄組Unigene 的Nr 功能注釋

通過Nr 庫比對，云南金花茶有58 830 條Unigenes 在Nr 數據庫中找到相似序列，注釋匹配的物種主要有葡萄Vitis vinifera、中粒咖啡Coffea canephora、可可樹Theobroma cacao、荷花Nelumbo nucifera、芝麻Sesamum indicum這5類，其中獲得注釋基因最多的是葡萄，有29.9%，中粒咖啡、可可樹、荷花、芝麻分別僅占5.6%、5.3%、4.8%、4.7%，其余49.7%的注釋基因分布于其他物種。從這些注釋的信息中可以得出，云南金花茶的大部分序列均可以在被子植物中得到相應的匹配。從e-value 分布（圖2A）可以看到，有44.2%的e-value 分布于l×10-100～l×10-45，有30%的e-value 分布于l×10-45～l×10-5，當e-value為0 時占25.8%。此外，有49.5%的序列相似度可達80%～95%，甚至有9.3%的序列相似度達到95%～100%，僅7.7%的序列相似度在60%以下，可以看出物種的序列相似度較高（圖2B）。總體而言，從e-value 和序列相似度分布情況可看出，云南金花茶在Nr 數據庫中比對的匹配度較高，但是由于缺乏云南金花茶一些基因組及轉錄組信息，導致部分Unigene 在數據庫中未得到匹配。

2.2.2 云南金花茶轉錄組Unigene 的GO 功能注釋

根據Nr 注釋成功的基因進行GO 功能分類注釋，其結果如圖3所示。分析結果表明，共有41 905 條Unigene 注釋了224 129 個GO功能，占Unigenes 總數量的43.66%。按3 大功能類別劃分，生物過程功能類別基因序列為107 044 條，占總數的47.76%；細胞組分功能類別65 990 條，占總數的29.44%；分子功能類別51 095 條，占總數的22.80%。由此可知，在生物過程功能類別中所注釋的基因比例最大。3 個功能大類進一步可劃分為56 個GO 功能亞類，分別包括25、21 和10 個亞類。在生物過程包含的25 個功能亞類中，獲得注釋偏多的分別是代謝過程、細胞過程、單一有機體過程，分別占該類型的20.970%、22.570%和16.750%，細胞聚合過程所得到注釋的比例最少，僅有0.008%。在細胞組分類別中，細胞和細胞部分所得到的注釋居多，分別為12 867 和12 865 條，均約占細胞組分這一大類的19.500 0%，而細胞外基質組分、擬核、共質體得到注釋較少，分別占0.006 1%、0.007 6%、0.006 1%。分子功能類別中，結合、催化活性得到的注釋較多，各占所屬分類總數的46.830 0%和38.100 0%，而金屬伴侶蛋白活性所得注釋最少，僅占0.005 9%。上述GO 功能注釋結果顯示了云南金花茶葉片中基因表達的基本情況，可以看出，在56 個功能亞類中，生物過程中代謝活動相關的基因量較多，說明云南金花茶有著較強的代謝能力。

圖2 云南金花茶轉錄組Unigenes 的Nr 注釋分類Fig.2 Annotation category on unigenes Nr in transcriptome of C.fascicularis

圖3 云南金花茶轉錄組Unigene 的GO 功能分類Fig.3 GO classification of unigenes in transcriptome for C.fascicularis

2.2.3 云南金花茶轉錄組Unigene 的KOG 功能注釋與分類

將獲得的Unigene 進行KOG 蛋白數據庫分類注釋，其結果如圖4所示。分析結果表明，在KOG 中有23 499 條Unigenes 能夠匹配，占總數的24.480%，共獲得KOG 功能注釋信息26 430 個，根據比對結果可分為26 個功能大類，包括能量產生、轉化，次生代謝物的生物合成、加工、運輸等不同類別的基因表達。其中，主要為一般功能預測基因，占總數的16.620%；其次是翻譯后修飾、蛋白質轉化和分子伴侶基因，占總數的11.630%；信號傳導機制（8.310%），轉錄（5.420%），RNA 加工和修飾（5.360%），胞內運輸、分泌和膜泡運輸（5.310%）的功能基因也占有較高的比例。而胞外結構所獲得的功能注釋信息最少，僅有54 個，占總數的0.020%。可見，云南金花茶在轉錄、翻譯和蛋白質運輸等方面的基因表達量較多。此外，還有1 個未知蛋白，不能獲知其具體生物學功能，占總數的0.004%。

2.2.4 云南金花茶轉錄組Unigene 的KEGG 代謝通路分析

將獲得的Unigene 比對到KEGG 數據庫中，共有23 214 條Unigenes 獲得注釋，占總Unigenes數量的24.18%。根據其涉及的代謝通路可將云南金花茶Unigenes 歸為5 大類別和19 個亞類，其結果如圖5所示。

圖4 云南金花茶轉錄組Uingenes 的KOG 功能注釋分布Fig.4 KOG functional annotation distribution of unigenes in transcriptome for C.fascicularis

圖5 云南金花茶轉錄組Unigene 的KEGG 分類Fig.5 KEGG classification of C.fascicularis Unigene

通過對圖5中相關通路分類下Unigenes 的具體分析統計發現：5 大類別中代謝通路所占比例最多，有9 584 條Unigenes，占總數的55.94%；其次是遺傳信息處理相關的通路，占27.46%；而環境信息處理和細胞過程中相關的通路，皆占6.60%；有機系統相關的通路最少，僅占4.93%。將5 大類別進一步細分為亞類，其中代謝相關的通路可分為11 個亞類，以糖類代謝居多，占注釋Unigenes 總數的19.99%。其次為整體映射相關的通路，占總數的14.80%，萜類化合物和聚酮化合物的代謝相關通路最少，所占比例僅為4.22%。另外，遺傳信息處理相關的通路分為4 個亞類，翻譯相關的通路最多，占總數的38.41%；其次為折疊、分選和降解通路，占32.87%；復制和修復通路所占比例最少，僅為11.99%。在環境信息處理通路中，僅包括2 個亞類，以信號傳導通路居多，占87.92%。細胞過程和有機系統相關的通路均僅1 個亞類，分別占6.60%和4.93%。在KEGG 代謝通路分析結果中，代謝通路類別所獲得的注釋基因最多，表明云南金花茶在這一時期有較強的代謝活動。

2.3 云南金花茶轉錄組Unigene 的CDS 預測

蛋白數據庫的比對結果顯示，有60 939 條Unigene 比對到蛋白庫中，另預測到有26 428 條CDS，其長度分布如圖6所示。從blast 比對得到的CDS 長度分布（圖6A）可看出，達到1 000 nt以上長度的CDS 序列，占CDS 總數的30.59%，其中1 000 ～2 000 nt 長度的CDS 序列占總數的22.940%，2 000 ～9 000 nt 長度的CDS 序列占總數的7.630%，9 000 nt 以上長度的CDS 序列占0.260%。通過ESTScan 預測的CDS 長度集中分布于100 ～500 nt，占82.100%，有極少數序列的長度在9 000 nt 以上，占0.004%（圖6B）。

圖6 云南金花茶轉錄組Unigene 的CDS 序列長度分布Fig.6 CDS length distribution of transcriptome for C.fascicularis

3 結論與討論

采用Illumina HiSeq 2000 高通量測序可以同時完成前基因組學研究（測序和注釋）以及后基因組學（基因表達及調控，基因功能，蛋白和核酸相互作用）研究。目前，此項技術已被應用于鐵皮石斛Dendrobium officinale[32]、馬尾松Pinus massoniana[33]、云南松Pinus yunnanensis[34]、紅 豆杉Taxus chinensis[35]等多個物種基因組的分析。鑒于此，本研究中采用RNA-seq 技術對云南金花茶進行測序，獲得95 979 條Unigene，平均長度為807 bp，結果與其他茶科植物如油茶Camellia oleifera[36]、紫芽茶樹Camellia sinensis(Linn.) O.Kuntze[37]、“紫鵑”茶樹Camellia sinensisvar.Zijuan[38]等相比，拼接完整性較好。從整體上看，通過對云南金花茶Unigene 總長、GC 含量、堿基正確率、序列Q20 的分析，測序獲得的序列質量較高，數量較多，可為后續云南金花茶基因功能分析、分子標記開發、代謝通路等方面的研究提供參考。

通過Nr 數據庫比對，有58 830（占總數的61.29%）條序列在不同物種均有相應的注釋，其中注釋匹配的物種主要為葡萄，有29.9%，其次為中粒咖啡、可可樹等。從這些成功匹配的物種可以看出，云南金花茶的大部分序列均可以在被子植物中得到匹配。通過KOG 數據庫比對，注釋成功的KOG 功能信息為26 430 個，將其分為26 個功能大類，一般功能參與的基因最多，其次是翻譯后修飾、蛋白質轉化和伴侶基因。這與蔣會兵等[37]經研究得出的紫芽茶樹轉錄組KOG 功能注釋（26 個KOG 功能類別）分布大體一致。在KOG 注釋中有0.004%的未知蛋白，難以確定其具體生物功能，可能是注釋信息不完善所導致，這種情況在其他物種轉錄組分析中也有出現，如金錢松Pseudolarix amabilis[39]、云南松[34]、文冠果Xanthoceras sorbifolia[40]等。通過GO 數據庫比對，共獲得224 129 個GO 功能信息，按其具體的序列信息又可分為3 個大類和56 個亞類，其基因主要分布于細胞過程、代謝過程和單一有機體過程，且在56 個功能亞類中，生物過程中代謝活動相關的基因量較多，由此可看出，金花茶的代謝能力較強。根據KEGG 代謝通路分析結果，共有23 214（24.18%）條序列注釋成功，按照注釋結果可將其劃分為5 大類別和19 個亞類，其中代謝通路相關基因所占比例最高，表明云南金花茶在整個時期均有較強的代謝活動。其基因注釋分布特征與李明璽等[41]對靖安白茶轉錄組KEGG 注釋結果基本一致。將云南金花茶Unigene 通過ESTScan進行預測，共獲得26 428 條CDS，其長度集中分布于100 ～500 nt，而蔡年輝等[34]對云南松轉錄組Unigene 的CDS 預測結果顯示其長度集中分布于200 ～1 000 nt，說明不同物種間存在較大差異。

通過7 個不同功能領域的基因蛋白數據庫注釋，可以看出云南金花茶所含基因信息豐富，通過分析所有注釋信息，可以更深層次地探索其基因組信息和基因分布情況。盡管這些Unigenes 序列并未覆蓋整個云南金花茶蛋白編碼區，但所注釋成功的基因仍有助于云南金花茶功能基因的挖掘和利用，以及為山茶屬植物遺傳育種等方面的研究提供理論參考。在注釋中還有33.44%的Unigenes 未注釋成功，這些Unigenes 序列可能為其他未編碼RNA 序列和未含有蛋白質功能信息的序列，也可能是因為基因數據庫中信息不足所導致。本研究結果也可為后續云南金花茶基因組水平的研究及遺傳改良等方面的研究提供一定的參考，且可為云南金花茶的分子標記開發以及抗逆機理研究提供數據。