朱興正，夏麗飛，陳林波*，孫云南，田易萍，宋維希，蔣會兵

?

保護品種云茶1號茶樹全長轉(zhuǎn)錄組測序分析

朱興正1,2，夏麗飛1,2，陳林波1,2*，孫云南1,2，田易萍1,2，宋維希1,2，蔣會兵1,2

1. 云南省農(nóng)業(yè)科學院茶葉研究所/云南省茶樹種質(zhì)資源創(chuàng)新與配套栽培技術工程研究中心，云南 勐海 666201； 2. 云南省茶學重點實驗室，云南 勐海 666201

為探討云茶1號茶樹品種優(yōu)異性狀的遺傳基礎，采用PacBio平臺進行全長轉(zhuǎn)錄組測序，最終獲得Polished consensus序列213?389個，預測到CDS有223?120個，檢測到195?062個SSR位點。在NR數(shù)據(jù)庫有170?264個同源序列比對到980個物種；有103?124個在KOG數(shù)據(jù)庫得到注釋，根據(jù)其功能各分為26類；有65?524個得到GO注釋分為細胞組分、分子功能及生物學過程等三大類的55個功能組；根據(jù)KEGG數(shù)據(jù)庫，105?972個得到了注釋，涉及到216個代謝途徑分支，包括茶葉品質(zhì)、活性物質(zhì)代謝以及抗逆等相關基因等；還預測到隸屬于60個轉(zhuǎn)錄因子家族的轉(zhuǎn)錄因子有5?785個。這些結(jié)果為進一步開展云茶1號茶樹特異性狀基因的標記性引物開發(fā)、遺傳研究以及品質(zhì)形成和抗逆機制研究奠定了基礎。

云茶1號；全長轉(zhuǎn)錄組；基因分析；功能注釋

云茶1號茶樹品種是云南省農(nóng)業(yè)科學院茶葉研究所于1993—2005年從云南元江細葉糯茶群體品種中采用單株育種法育成，其芽葉生育力強、持嫩性好，發(fā)芽密，茸毛特多，芽葉肥壯，一芽二葉百芽重達144.0?g。春茶一芽二葉蒸青樣水浸出物44.6%，茶多酚31.3%，氨基酸3.4%，咖啡堿4.3%，兒茶素總量158.7?mg·g-1，適制紅茶、綠茶、普洱茶、白茶。抗逆性鑒定發(fā)現(xiàn)其抗寒、抗旱、抗小綠葉蟬和抗茶餅病能力強，在云南的西雙版納、普洱、保山等州（市）有種植，湖南、廣西等省有引種，是一個性狀優(yōu)異的優(yōu)良品種[1-3]。2005年11月國家林業(yè)局植物新品種保護辦公室授予植物新品種權，品種權號20050030[4]。

基于PacBio平臺的第三代測序技術是一種集高通量、快速度以及超長的讀長優(yōu)勢和低成本等多種優(yōu)點于一身的新型測序技術。它最大特點是無需進行PCR擴增，可直接讀取目標序列，因此假陽性率大大減少，同時避免了堿基替換及偏置等常見PCR錯誤的發(fā)生，精準度可達到99.9%[5-6]。因此，本試驗采用PacBio平臺的第三代測序技術對云茶1號進行全長轉(zhuǎn)錄組測序，通過生物信息學方法對所測的序列進行序列分析、功能注釋、功能分類及代謝途徑分析等，旨在為進一步挖掘云茶1號茶樹品種次生代謝相關功能基因、抗逆基因及開發(fā)分子標記奠定基礎。

1 材料與方法

1.1 試驗材料

采摘云南省農(nóng)業(yè)科學院茶葉研究所試驗基地樹齡為13年的云茶1號茶樹的芽、葉、莖、花芽、花蕾、幼果，利用液氮迅速固樣后送至北京諾禾致源科技股份有限公司進行RNA提取及測序等分析。

1.2 RNA的提取與文庫構建

采用Trizol法[7]提取分別提取云茶1號芽、葉、莖、花芽、花蕾、幼果總RNA，利用瓊脂糖凝膠電泳、Nanodrop、Agilent 2100、Qubit等檢測合格后，進行等量混勻。混勻后的RNA利用帶有Oligo(dT)的磁珠富集mRNA，使用SMARTer PCR cDNA Synthesis Kit將mRNA反轉(zhuǎn)錄為cDNA，再用BluePippin進行片段篩選，對全長cDNA進行損傷修復、末端修復、連接SMRT啞鈴型接頭和核酸外切酶消化后獲得文庫。

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

對庫檢合格文庫運用Pacbio Sequel平臺進行測序，原始下機數(shù)據(jù)采用PacBio官方軟件SMRTlink處理獲得Subreads序列，由Subreads之間的校正獲得環(huán)形一致性序列（CCS），再根據(jù)序列是否包含5'端引物、3'端引物以及polyA尾將序列分為全長序列與非全長序列，然后利用同種類型聚類（ICE）對全長序列進行聚類，獲得Cluster consensus序列，最后使用非全長序列對得到的一致序列進行校正（Polishing），獲得高質(zhì)量的序列進行后續(xù)分析。再利用ANGEL軟件對獲得的序列進行編碼蛋白產(chǎn)物的序列（Coding sequence, CDS）預測分析[8]。采用MISA軟件對單核苷酸、二核苷酸、三核苷酸、四核苷酸、五核苷酸和六核苷酸按照最少重復次數(shù)分別為9、6、5、5、5和5對Gene進行SSR檢測。

1.4全長轉(zhuǎn)錄組序列的功能注釋

對獲得的高質(zhì)量序列，利用公共數(shù)據(jù)庫包括非冗余蛋白數(shù)據(jù)庫（Non-Redundant Protein Database, NR）、蛋白質(zhì)真核同源數(shù)據(jù)庫（Eukaryotic Orthologous Groups, KOG）、蛋白質(zhì)序列數(shù)據(jù)庫（Swiss Prot Protein Database, Swiss Prot）、蛋白質(zhì)原核同源數(shù)據(jù)庫（Cluster of Orthologous Groups of Proteins, COG）、基因本體論數(shù)據(jù)庫（Gene Ontology, GO）、蛋白質(zhì)家族域數(shù)據(jù)庫（Protein Families Database, Pfam）、東京基因與基金組百科全書（Kyoto Encyclopedia of Genes and Genomes, KEGG）進行基因功能注釋。使用iTAK軟件[9]以及采用database中分類定義好的轉(zhuǎn)錄因子（Transcription Factor, TF）family及規(guī)則，通過hmmscan鑒定TF。

2 結(jié)果與分析

2.1 全長轉(zhuǎn)錄組測序與組裝分析

2.1.1 測序結(jié)果與數(shù)據(jù)組裝

采用PacBio Sequel測序平臺對云茶1號芽、葉、莖、花芽、花蕾、果等混合樣進行全長轉(zhuǎn)錄組測序，共獲得6?282?894個Subreads（10.62?Gb），平均subreads長度為1?690?bp，N50為2?646。通過檢測其中環(huán)形一致性序列（Circular Consensus Sequence, CCS）Read有499?432個，帶有5'端引物的reads數(shù)為454?699個，帶有3'端引物的reads數(shù)有464?035個，帶有PolyA尾的reads數(shù)有451?234個；全長（Full-Length, FL）reads數(shù)有400?959個；全長非嵌合（Full-Length non-chimericRead, FLNC）reads數(shù)有388?370，F(xiàn)LNC平均長度為2?469?bp；使用ICE算法將同一轉(zhuǎn)錄本的FLNC序列進行聚類，得到consensus序列，并采用非全長的序列對得到的consensus序列進行校正，最終獲得Polished consensus序列213?389個，用于后續(xù)分析。

2.1.2 CDS預測

CDS（Coding sequence）是編碼一段蛋白產(chǎn)物的序列，與蛋白質(zhì)的密碼子完全對應。利用ANGEL軟件[8]進行CDS預測分析，得到其編碼區(qū)的氨基酸序列和核酸序列。將獲得的Polished consensus在蛋白質(zhì)數(shù)據(jù)庫中進行Blast比對以后，發(fā)現(xiàn)有223?120個Coding Sequence（CDS），序列長度范圍分布在200~ 2?500?bp之間，主要集中在300~1?100?bp（圖1），說明Unigenes的序列質(zhì)量較好。

2.1.3 SSR分析

采用MISA對云茶1號全長轉(zhuǎn)錄組Gene進行SSR檢測，搜索標準為：1~6個核苷酸基序（motif）重復次數(shù)分別為大于等于9、6、5、5、5和5，結(jié)果見表1。有195?062個序列檢測到SSR位點，單核苷酸和二核苷酸重復類型占優(yōu)勢，分別占總SSR的57.58%和31.14%；其他４種重復類型所占比例相對較少：三核苷酸重復占9.28%，四核苷酸重復占0.92%，五核苷酸重復占0.46%，六核苷酸重復占0.62%。在檢測到的SSR位點中，單核苷酸中出現(xiàn)頻率最高的核苷酸基序為T/A（94?710）和A/T（45?750個）；二核苷酸以CT/AG最多有13?993個，其次是TC/GA和AG/CT，分別有13 492個和8?246個；三核苷酸以GAA/CTT（956個）、TTA/AAT（919個）、CCA/GGT（875個）占優(yōu)勢；四核苷酸出現(xiàn)頻率以TTTA/AAAT（296個）和TTAT/AATA（133個）為多；五核苷酸和六核苷酸分別以TGTGT/ACACA（159個）和TCCACC/ AGGTGG（25個）最多。云茶1號中兒茶素、氨基酸、咖啡堿等茶葉品質(zhì)成分含量高、抗逆性強以及芽葉大等優(yōu)良性狀，研究價值很高。因此，上述對SSR特征分析的結(jié)果，為進一步開展云茶1號及茶樹通用性引物設計及其遺傳圖譜構建等打下了良好基礎。

圖1 CDS長度分布

Fig.1CDS length distribution

表1 SSR不同重復基序分布及優(yōu)勢堿基組成

2.2 Unigene的功能注釋

2.2.1注釋結(jié)果統(tǒng)計

通過BLAST程序，對轉(zhuǎn)錄組序列進行NR數(shù)據(jù)庫比對。有170?264個Unigene在NR數(shù)據(jù)庫比對到980個物種上，比對較多的20個物種包括葡萄（28?588個）、茶樹（8?753個）、核桃（7?601個）、中?？Х龋??000個）、芝麻（5?892個）、可可（5?767個）、蓮（5?471個）、棗（4?170個）、木薯（4?043個）、胡蘿卜（3?962個）、麻瘋樹J（3?499個）、土瓶草（3?335個）、柑橘（3?156個）、牽?；ǎ??144個）、蓖麻（3?123個）、桃（2?899個）、野生煙草（2?855個）、楊樹（2?663個）、煙草（2?470個）、刺菜薊（2?330個）、青梅（2?243個）。從NR庫中的注釋情況看，注釋到葡萄的基因最多，可能是GenBank中登錄葡萄的基因信息多，且茶樹與葡萄的相對較為近緣。同時隨著茶樹生物學研究的不斷深入，登錄到GenBank的基因越來越多，獲得較多茶樹基因注釋。

2.2.2 Unigene的KOG分類統(tǒng)計

將云茶1號Unigene與KOG數(shù)據(jù)庫進行比對并根據(jù)其功能進行分類統(tǒng)計見表2。有103?124個Unigene得到注釋，按照功能一共分為26類，其中一般功能預測最多（19?620個），其次是翻譯后修飾、蛋白折疊和分子伴侶（13?209個），最少的為未知蛋白（2個）。在KOG注釋分類中，注釋到次生代謝物合成、運輸和代謝有5?183個，氨基酸運輸和代謝的有4?363個，說明了云茶1號茶樹中次生代謝和氨基酸代謝較為活躍，這為云茶1號茶樹茶葉品質(zhì)的形成奠定了基礎。

2.2.3 Unigene的GO功能注釋

有65?524個Unigene注釋到GO數(shù)據(jù)庫，可分為細胞組分、分子功能、及生物學過程等3個大類的55個功能組（表3）。其中生物學過程包括25個功能組，以代謝過程（33?845個）、細胞過程（33?284個）、單生物過程（24?567個）、定位（9?572個）的最多，細胞聚集功能最少（1個），其次是節(jié)律過程（5個）；細胞組分包括20個功能組，以細胞部分（15?349個）、細胞（15?349個）、細胞器（9?900個）、膜結(jié)構（8?973個）、膜部分（8?487個）較多，以細胞連接、突觸以及突觸部分最少均為14個；分子功能包括10個功能組，以結(jié)合活性（39?544個）、催化活性（32?858個）最多，金屬伴侶蛋白活性為最少（3個）等。由于同個基因可能有多個注釋，因此，注釋到GO數(shù)據(jù)庫的基因總數(shù)大于實際注釋到的基因數(shù)。

表2 Unigene的KOG系統(tǒng)分類

2.2.4 Unigene的KEGG功能注釋

KEGG數(shù)據(jù)庫已建立了一套完整KO注釋的系統(tǒng)，整合了基因組、化學分子和生化系統(tǒng)等方面的數(shù)據(jù)，包括代謝通路（KEGG PATHWAY）、功能模型（KEGG MODULE）、基因序列（KEGG GENES）及基因組（KEGG GENOME）等，可完成新測序物種的基因組或轉(zhuǎn)錄組的功能注釋。云茶1號茶樹轉(zhuǎn)錄組有105?972個得到了注釋，根據(jù)參與的KEGG代謝通路分為三個層次，第一層次分為四大類，第二層為23小類（表4），第三層次為216個代謝途徑分支。涉及茶葉滋味相關代謝途徑的有氨基酸代謝（8?479個）、黃酮類化合物的生物合成（464個）、黃酮和黃酮醇的生物合成（86個）、咖啡因的代謝（41個）。涉及與香氣相關的單萜生物合成（80個）、倍半萜生物合成（168個）、萜類化合物骨架生物合成（378個）。涉及植物激素信號轉(zhuǎn)導（1?772個）、MAPK信號通路（1?194個）、AMPK信號通路（1?149個）等。這些Unigene為今后深入開展云茶1號茶滋味、香氣以及抗逆等相關的功能基因研究奠定了基礎。

表3 Unigene的GO分類

表 4 Unigene的KEGG分類

2.2.5 轉(zhuǎn)錄因子分析

轉(zhuǎn)錄因子（Transcription factor, TF）也稱為反式作用因子，是一群能與基因5'端上游特定序列專一性結(jié)合，從而保證目的基因以特定的強度在特定的時間與空間表達的DNA結(jié)合蛋白，通過它們之間以及與其他相關蛋白之間的相互作用，激活或抑制轉(zhuǎn)錄，在植物生長發(fā)育過程、組織分化、營養(yǎng)運輸、代謝合成及環(huán)境應答等生命過程中起著至關重要的轉(zhuǎn)錄調(diào)控作用[10]。使用iTAK軟件預測到轉(zhuǎn)錄因子有5?785個，隸屬于60個轉(zhuǎn)錄因子家族，其中比較多的轉(zhuǎn)錄因子家族見（圖2）。在獲得的轉(zhuǎn)錄因子家族中NAC家族有202個、AP2/ERF家族有185個、AUX/IAA家族有216個、MYB家族有92個、bHLH家族有235個、WRKY家族有172個等。這些轉(zhuǎn)錄因子家族成員的獲得，為后期研究參與云茶1號茶樹中次生代謝物的生物合成、生長發(fā)育以及抗逆調(diào)控奠定基礎。

3 討論

PacBio第三代測序技術是基于單分子實時（Single molecule real time, SMRT）的單分子測序儀，由美國Pacific Biosciences（PacBio）[5]公司設計制造，最大特點是無需進行PCR擴增，可直接讀取目標序列，因此假陽性率大大減少，是轉(zhuǎn)錄組從頭測序的首選[11-12]。本研究采用PacBio平臺對云茶1號茶樹全長轉(zhuǎn)錄組測序及分析，獲得6?282?894個Subreads，平均subreads長度為1?690?bp，N50為2?646，利用公共數(shù)據(jù)庫NR、KOG、Swiss Prot、COG、GO、Pfam、KEGG對獲得的213?389個Polished consensus序列進行功能注釋和分類。有170?264個Unigene比對到980個物種上；有103?124個Unigene得到KOG注釋，按照功能一共分為26類；有105?972個Unigene注釋216個代謝途徑分支。Shi等[13]采用二代測序技術對龍井43茶樹的嫩葉、成熟葉、莖、幼根、花蕾以及成熟種子的總RNA等量混合進行轉(zhuǎn)錄組測序，共獲得127?094條Unigene，平均長度為355?bp，N50為506?bp，注釋到Nr、COG以及KEGG數(shù)據(jù)庫的分別為53?937個、15?701個和16?939個。陳林波等[14]以紫娟茶樹的芽、第二葉、開面葉、成熟葉為材料，利用二代測序技術對其轉(zhuǎn)錄組進行測序分析，獲得206?210條Unigene，注釋到NR的為71?799個、注釋到GO的為61?141個、注釋到KEGG的為30?995個。這些研究結(jié)果表明，不論是從獲得的序列質(zhì)量、基因數(shù)以及注釋到的基因信息來看，第三代測序結(jié)果均優(yōu)于第二代。

圖2轉(zhuǎn)錄因子家族分析

茶樹為葉用經(jīng)濟作物，云茶1號茶樹其芽葉大，產(chǎn)量高，抗茶餅病和抗小綠葉蟬強，芽葉中兒茶素、黃酮以及氨基酸等生物活性物質(zhì)含量高，制作的茶葉品質(zhì)優(yōu)，是一個集多種生物性狀優(yōu)異的茶樹品種，是茶樹遺傳改良的重要親本。本研究對云茶1號茶樹全長轉(zhuǎn)錄組進行測序，通過7個數(shù)據(jù)庫進行功能注釋，在KEGG代謝通路注釋到的105?972個unigenes，歸屬于216條通路。其中一些代謝途徑與茶葉品質(zhì)形成相關，包括氨基酸代謝（8?479個），咖啡因的代謝（41個），黃酮和黃酮醇的生物合成（86個），黃酮類化合物的生物合成（464個），單萜的合成（80個），倍半萜的生物合成（168個），二萜類化合物的生物合成（349），這些數(shù)據(jù)為云茶1號茶葉品質(zhì)形成機制的研究提供了基礎數(shù)據(jù)。還有一些與植物體內(nèi)重要的抗逆途徑相關，包括植物激素信號轉(zhuǎn)導（1?772個），AMPK信號途徑（1?149個），鈣信號途徑（431個），MAPK信號途徑（1?194個），苯丙素的生物合成（826個）以及預測到的轉(zhuǎn)錄因子5?785個和所獲得SSR信息，這將有助于進一步開展云茶1號茶樹的抗病機理以及特異性狀基因的標記性引物開發(fā)和遺傳研究奠定基礎。

[1] 張俊, 田易萍, 徐丕忠, 等. 優(yōu)質(zhì)、抗病大葉茶新品種“云茶1號”選育[J]. 茶葉, 2008, 34(1): 39-41.

[2] 王深. 云茶1號主要特點與栽培管理[J]. 農(nóng)村實用技術, 2013(2):19.

[3] 冉隆珣, 玉香甩, 田易萍, 等. 不同茶樹品種對茶餅病抗性鑒定初探[J]. 遼寧農(nóng)業(yè)科學, 2017(1): 69-70.

[4] 田易萍, 張俊, 徐丕忠, 等. 茶新品種‘云茶1號’[J]. 園藝學報, 2009, 36(1): 153.

[5] Liao Y C, Lin S H, Lin H H. Completing bacterial genome assemblies: strategy and performance comparisons [J]. Scientific Reports, 2015, 5: 8747.

[6] Shin S C, Ahn D H, Kim S J, et al. Advantages of single-molecule real-time sequencing in high-GC content genomes [J]. PLoS One, 2013, 8(7): e68824.

[7] Gao J P, Wang D, Cao L Y, Sun H F. Transcriptome sequencing of codonopsis pilosula and identification of candidate genes involved in polysaccharide biosynthesis [J]. PLoS One, 2015, 10(2): 117-134.

[8] Shimizu K, Adachi J, Muraoka Y. ANGLE: a sequencing errors resistant program for predicting protein coding regions in unfinished cDNA [J]. Journal of Bioinformatics and Computational Biology, 2006, 4(3): 649-664.

[9] Zheng Y, Jiao C, Sun H, et al. iTAK: a program for genome-wide prediction and classification of plant transcription factors, transcriptional regulators, and protein kinases [J]. Molecular Plant, 2016, 9(12): 1667-1670.

[10] Lutoval A, Doduevai E, Lebedeva M A, et al. Transcription factors in developmental genetics and the evolution of higher plants [J]. Russian Journal of Genetics, 2015, 51(5): 449-466.

[11] 曹晨霞, 韓琬, 張和平. 第三代測序技術在微生物研究中的應用[J]. 微生物學通報, 2016, 43(10): 2269-2276.

[12] 任毅鵬, 張佳慶, 孫瑜, 等. 基于PacBio平臺的全長轉(zhuǎn)錄組測序[J]. 科學通報, 2016, 61(11): 1250-1254.

[13] Shi C Y, Yang H, Wei C L, et al. Deep sequencing of thetranscriptome revealed candidate genes for major metabolic pathways of tea-specific compounds [J]. BMC Genomics, 2011, 12(1): 131.

[14] 陳林波, 夏麗飛, 周萌, 等. 基于RNA-Seq技術的“紫娟”茶樹轉(zhuǎn)錄組分析[J]. 分子植物育種, 2015, 13(10): 2250-2255.

Full-length Transcriptome Analysis of Protected Cultivation ‘Yuncha 1’ (Var)

ZHU Xingzheng1,2, XIA Lifei1,2, CHEN Linbo1,2*, SUN Yunnan1,2, TIAN Yiping1,2, SONG Weixi1,2, JIANG Huibin1,2

1. Tea Research Institute, Yunnan Academy of Agricultural Sciences/Yunnan Technology Engineering Research Center of Tea Germplasm Innovation and Supporting Cultivation, Menghai, 666201, China; 2.Yunnan Provincial Key Laboratory of Tea Science, Menghai, 666201, China

To explore the genetic basis for important traits, the full-length transcriptome of the ‘Yuncha 1’ () was sequenced by using PacBio Platform. A total of 213?389 polished consensus were generated, 223?120coding sequences were predicted and annotated, and 195?062 SSR loci were found. According to NR databases, 170?264 homologous sequences were mapped to 980 species, 103?124 unigenes were further annotated and grouped into 26 functional categories in KOG databases, 65?524 unigenes were annotated against GO database and divided into cellular component, molecular function and biological process categories with a total of 55functional groups. KEGG pathway analysis showed that 105?972 unigenes could be broadly classified into 216 metabolism pathways according to their function, and some of them were involved in quality, bioactive substances, and resistance gene, etc. It is also predicted that there were 5?785 transcription factors belonging to 60 transcription factor families. The experimental results will give important data for development of SSRs of specifictraits,genetic analysis and studies involved in quality formation and resistence mechanism in tea cultivar ‘Yuncha 1’.

tea cultivar ‘Yuncha 1’, full-length transcriptome, gene analysis, function annotation

S571.1；S324

A

1000-369X(2018)02-193-09

2017-10-31

2017-11-27

國家自然科學基金項目（31660224、31560220、31460216)、云南省人才培養(yǎng)計劃項目（2015HB105）

朱興正，男，副研究員，主要從事茶樹育種與生物技術方面的工作。

chenlinbo2002@sina.com