苦茶全長轉錄組測序及基因結構分析

龐丹丹，劉玉飛，孫云南，田易萍，宋維希，陳林波

(1.云南省農(nóng)業(yè)科學院 茶葉研究所，云南勐海 666201；2.云南省茶學重點實驗室，云南勐海 666201)

苦茶(Camelliasinensis)是中國的一種特異茶樹種質資源[1]，由其制成的茶葉感官品質表現(xiàn)極苦，且與其他茶類的苦味不同，因此，被命名為苦茶。苦茶多分布于云南、四川，以及廣東、湖南、江西毗鄰區(qū)，尤以云南以及南嶺山脈兩側最多[2-3]。在苦茶的原生地，苦茶成為當?shù)厝嗣裆畹谋仨毱?，將它作為一種藥物飲用，認為長期飲用苦茶具有“退火發(fā)汗、解毒、治病”的功效[2]。研究已證實，苦茶嘌呤生物堿的組分與常規(guī)茶樹差異較大，其以苦茶堿(1,3,7,9-四甲基尿酸)為主，其次是咖啡堿、可可堿[3-4]，苦茶堿具有鎮(zhèn)靜催眠[5]、抗抑郁[6]等藥理活性，并且毒理試驗鑒定苦茶堿為無毒[6]，其還可以削弱咖啡堿的興奮作用[7]，這使得苦茶的價值不斷提升，消費需求連年增加?？嗖璩懈吆康目嗖鑹A外，還有非兒茶素類茶多酚[8]。同時，苦茶中還具有高含量的表沒食子兒茶素(EGC)，EGC含量與紅茶中關鍵品質成分茶黃素的含量呈極顯著正相關[2,9]，所以苦茶可以作為選育和培育特異高級紅茶茶樹品種的育種材料；以往的研究表明，苦茶茶葉中存在一種具有丁香香氣的特殊物質丁子香酚甙[2]。這說明苦茶還具有多種不同于常規(guī)栽培茶樹的特異性狀，這些性狀形成的機理、性狀相關基因標記的開發(fā)以及遺傳機制的研究，都需要借助于苦茶資源。

當前，大多數(shù)轉錄組測序是基于Illumina平臺的第二代測序(SGS)技術生成的[10]。但是，SGS技術不能產(chǎn)生較長的轉錄本，并且SGS無法獲得可變剪接等信息，從而限制了該技術在轉錄組測序中的利用[11]。而全長轉錄組測序是基于第三代測序(TGS)平臺進行的，該方法可以產(chǎn)生長讀段，因此可以直接測序全長轉錄本[12]，在這方面全長轉錄測序優(yōu)于短讀測序；另外，它還可用于選擇性剪接事件及初級-前體-成熟RNA結構的分析，以幫助更好地理解RNA的加工過程；此外，其可以在轉錄水平上全面分析由交替剪接和基因融合產(chǎn)生的結構變異，比SGS更準確地檢測結構變異(SV)，使其適合于多倍體物種或具有高重復序列和高雜合性的物種的轉錄組分析[13]。由于茶樹基因組具有高的重復序列(至少包含64%的重復序列)和高的雜合度[14]，以及苦茶特異性狀形成機制需進一步解析。因此，本研究以苦茶為材料，利用PacBio平臺進行Iso-Seq，對所獲得序列與參考基因組比對、同時對其進行功能注釋、基因結構分析等，為苦茶特征形狀遺傳分析提供基礎數(shù)據(jù)，為更好地了解和利用這一重要的特異茶樹資源提供幫助。

1 材料與方法

1.1 試驗材料

用于測序分析的苦茶資源‘老曼娥苦茶’來源于云南省農(nóng)業(yè)科學院茶葉研究所試驗基地，其6個不同組織(芽、葉、花芽、花蕾、莖和幼果)的樣品采摘后用液氮冷凍，并置于-80 ℃冰箱，之后送至諾禾致源科技股份有限公司(北京)進行測序。

1.2 RNA提取與構建文庫

參照張亞真等[15]的方法，分別提取苦茶芽、葉、花芽、花蕾、莖、幼果不同部位的總RNA，檢測合格后進行等量混勻。利用帶有Oligo(dT)的磁珠從混勻后的RNA樣品中分離出mRNA，并利用SMARTer PCR cDNA Synthesis Kit將其反轉錄為cDNA，接著用BluePippin對cDNA進行片段篩選，對全長cDNA進行損傷修復、末端修復，并連接SMRT啞鈴型接頭，最終使用核酸外切酶消化獲得文庫。

1.3 Iso-Seq及組裝分析

采用軟件SMRT Link v5.0對輸出進行過濾和處理，參數(shù)：--minLength=200，--minReadScore=0.75，最終得到的數(shù)據(jù)即為有效數(shù)據(jù)，校正后獲得CCS序列(Circular Consensus Sequence，參數(shù)：--minPasses=1，minPredictedAccuracy=0.8)；通過檢測CCS是否包含5′-primer、3′-primer、poly-A對CCS進行分類找出FLNC(full-length non chimera)序列，之后對FLNC進行聚類和校正，最終獲得polished consensus序列用于后續(xù)分析。使用MISA(http://pgrc.ipk-gatersleben.de/misa/misa.html)軟件對單、二、三、四、五和六核苷酸的最低重復次數(shù)分別設置為10、6、5、5、5和5，進行SSR搜索。

1.4 參考基因組比對及功能注釋

使用GMAP(http://research-pub.gene.com/gmap/)將三代測序reads比對到物種的參考基因組[16]，進行飽和度曲線及轉錄本密度分析，并對未比對到參考基因組上轉錄本進行七大數(shù)據(jù)庫(NR、NT、Pfam、KOG/COG、Swiss-Prot、KEGG、GO)的功能注釋。

1.5 基因結構分析

根據(jù)轉錄本與參考基因組比對結果，對轉錄本進行再次校正，然后進行利用SUPPA(https://bitbucket.org/regulatorygenomicsupf/suppa)、Tapis(https://bitbucket.org/comp_bio/tapis)等軟件進行可變剪切分析、APA預測、新基因和新轉錄本鑒定、新基因數(shù)據(jù)庫注釋、轉錄因子分析。并利用CNCI(https://github.com/www-bioinfo-org/CNCI)、PLEK(https://sourceforge.net/projects/plek/)、CPC(http://cpc.cbi.pku.edu.cn/)軟件以及Pfam數(shù)據(jù)庫對PacBio測序數(shù)據(jù)進行編碼潛能預測，并將最終得到的LncRNA進行后續(xù)分析。

2 結果與分析

2.1 Iso-Seq分析

2.1.1 測序結果與組裝全長 轉錄組測序共獲得8 730 808個Subreads(10G)，平均Subreads的長度為1 145 bp，N50為1 741。通過校正，CCS序列(環(huán)形一致性序列)有418 424個，含有5′端引物的reads數(shù)有338 306個，含有3′端引物的reads數(shù)目為368 458個，含有PolyA尾的reads數(shù)目為351 952個；全長序列數(shù)目為 295 803個，F(xiàn)LNC序列的數(shù)目有217 701個，且FLNC的平均長度為1 836 bp；最終獲得Polished consensus序列132 334個，用于后續(xù)分析。

2.1.2 SSR分析 利用MISA對苦茶全長轉錄組測序獲得的25 735條(大于500 bp)轉錄本序列進行搜索，共在12 926條轉錄本序列中發(fā)現(xiàn)符合標準的21 106個SSR位點，其中5 227個轉錄本含有大于1個SSR位點，SSR的發(fā)生頻率為50.22%，SSR出現(xiàn)率為82%，平均2.4 kb出現(xiàn)1個SSR。在檢測到的苦茶轉錄組SSR中，單、二和三核苷酸重復類型的比例較大，分別占總SSR的32.49%、47.91%和16.18%；其他3種重復類型所占比例相對較少，只占總SSR的3.42%，結果見表1?？嗖柁D錄組SSR重復單元的重復次數(shù)分布在5～46，其中5～10次重復的SSR位點有13 625個，占位點總數(shù)的64.56%；11～20次重復的SSR位點有7 299個，占總數(shù)的34.58%；20次重復以上的SSR位點有0.86%(表1)。

表1 苦茶SSR的類型、數(shù)量及分布頻率Table 1 Type, number and frequency of SSRs in Kucha

檢測到的苦茶全長轉錄組SSR核苷酸基序類型中(表2)，二核苷酸重復基序中AG/CT (7 720個)，AT/AT(1 653個)的出現(xiàn)頻率較高；三核苷酸中出現(xiàn)最多的重復基序是AAG/CTT(926個)，其次是ACC/GGT(474個)、ATC/ATG(446個)。四核苷酸中以AAAT/ATTT(126個)、AAAG/CTTT(47個)、AAAC/GTTT(41個)占優(yōu)勢；五核苷酸中分別以AAAAC/GTTTT(35個)、AAAAG/CTTTT(34個)和AAACC/GGTTT(27個)出現(xiàn)頻率較高；AAAAAC/GTTTTT(16個)是六核苷酸中出現(xiàn)頻率最高的，占其總數(shù)的12.5%，其次是AAAAAG/CTTTTT(14個)和AACCCT/AGGGTT(13個)。上述分析結果為下一步苦茶種質資源鑒定、遺傳多樣性的分析和遺傳圖譜的構建等研究工作提供參考。

2.2 參考基因組比對與Unmapped轉錄本的功能注釋

分析與茶樹的參考基因組的比對結果，發(fā)現(xiàn)能比對到基因組上的一致性序列的reads共 105 898個，依據(jù)比對情況將序列分為五種類型：Unmapped、Multiple mapped、Uniquely mapped、Reads map to ‘+’、Reads map to ‘-’，同參考基因組比對結果見圖1。其中不能比對到基因組上(Unmapped)的reads有26 436個，占 19.98%；在參考序列上有多個比對位置(Multiple mapped)的reads有1 359個，占1.03%；在參考序列上有唯一比對位置(Uniquely mapped)的reads有104 539個，其中比對到基因組上正鏈(Reads map to ‘+’)的reads有101 781個，占76.91%；比對到基因組上負鏈(Reads map to ‘-’)的reads有2 758個，占2.08%。

為了獲得轉錄本全面的功能信息，將未比對到參考基因組上的轉錄本進行七大數(shù)據(jù)庫(NR，NT，Pfam，KOG/COG，Swiss-Prot，KEGG，GO)的功能注釋，結果有11 152個轉錄本在NR數(shù)據(jù)庫中比對到311個物種上，比對上較多的5個物種包括茶樹Camelliasinensis(1 519個)、葡萄Vitisvinifera(1 413個)、胡桃Juglansregia(441個)、中粒咖啡Coffeacanephora(376個)、水芙蓉Nelumbonucifera(354個)。從NR數(shù)據(jù)庫的注釋結果來看，注釋到茶樹的基因最多，其次是葡萄；對基因進行KO注釋的結果顯示，苦茶轉錄組中有10 976個得到注釋，與茶葉品質相關的KEGG通路主要有氨基酸代謝途徑(409個)、苯丙烷物質的生物合成(70個)，與茶葉香氣有關的KEGG通路主要有萜類物質骨架生物合成(38個)、倍半萜生物合成(114個)、黃酮及黃酮醇的生物合成(2個)，涉及苦茶中苦茶堿合成相關的嘌呤代謝(74個)，與植物激素信號轉導相關的有64個，MAPK信號通路相關的則有31個等。同時在KOG數(shù)據(jù)庫中6 296個轉錄本得到注釋，其中注釋到次生代謝物合成、運輸和代謝有425個，氨基酸運輸和代謝的有357個，另外，還有6 221個轉錄本注釋到GO數(shù)據(jù)庫，由于同一個轉錄本對應到多個GO條目下，使得在GO數(shù)據(jù)庫中得到注釋的基因數(shù)要比注釋到的轉錄本數(shù)目多。

表2 苦茶全長轉錄組中SSR基序的分布Table 2 Distribution of SSR motif in full-length transcriptome of Kucha

2.3 基因結構分析

結構分析是三代全長轉錄組中的一個重點內容，不同的樣本轉錄物不盡相同，進行轉錄本結構分析，可以統(tǒng)計所有轉錄本的結構差異，能夠準確識別二代測序中無法區(qū)分的同源基因或者同源異構體等。

2.3.1 可變剪切分析 利用SUPPA軟件對數(shù)據(jù)進行分析，結果(圖2)共獲得4 892個基因發(fā)生了可變剪切事件，其中外顯子跳躍事件(Skipped exon，SE)的基因有640個，占比為5.21%；外顯子互斥事件(Mutually exclusive exon，MX)的基因有43個，占比為0.35%；內含子滯留事件(Retained intron，RI)的基因有1 635個，占比為 13.32%；5′UTR區(qū)可變事件(Alternative 5′ splice site，A5)的基因有947個，占比為7.71%；3′UTR區(qū)可變事件(Alternative 3′ splice site，A3)的基因有1 355個，占比11.04%；起始外顯子可變事件(Alternative first exon，AF)的基因數(shù)目有202個，占比1.65%；終止外顯子可變事件(Alternative last exon，AL)的基因數(shù)目有70個，占比 0.57%，其中內含子滯留事件占最大比例，其次是3′UTR區(qū)可變事件。

2.3.2 新基因和新轉錄本分析 從與茶樹參考基因組比對結果來看，全長轉錄本可分為三類，其中比對到已知基因的已知轉錄本有10 842個，占比24.56%；已知基因的新轉錄本有26 184個，占比59.32%；比對到參考基因組未注釋區(qū)域的新基因有7 115個，占比為16.12%(圖3)。為了探究新基因的功能信息，同樣將新基因的轉錄本進行七大數(shù)據(jù)庫(NR、NT、Pfam、KOG/COG、Swiss-Prot、KEGG、GO)的功能注釋。NR數(shù)據(jù)庫比對結果顯示，3 627個轉錄本被比對到142個物種上，比對上的數(shù)目較大的物種分別為葡萄Vitisvinifera(866個)、胡桃Juglansregia(176個)、可可樹Theobromacacao(173個)、中粒咖啡Coffeacanephora(136個)、芝麻Sesamumindicum(135個)、茶樹Camelliasinensis(121個)、水芙蓉Nelumbonucifera(120個)。以上的分析結果顯示，注釋到葡萄的基因數(shù)量最多，其原因有可能是茶樹與葡萄的親緣關系較近。

與GO數(shù)據(jù)庫進行比對，結果發(fā)現(xiàn)有575個轉錄本得到了注釋，在生物學過程類別中，代謝過程(259個)、細胞過程(243個)、單生物過程(176個)最多；細胞組分類別中，細胞部分(120個)、細胞(120個)、細胞器(96個)較多；分子功能類別中，結合活性(382個)、催化活性(237個)較多。根據(jù)它們參與的KEGG進行分類，結果發(fā)現(xiàn)新基因中有5 626個成功得到注釋，苯丙烷生物合成相關的有18個、其中與苦茶中苦茶堿合成相關的嘌呤代謝有17個，與茶葉品質相關的氨基酸代謝有69個，與制成的茶葉相關的萜類化合物代謝有31個，信號轉導相關的途徑有150個，其中植物激素信號轉導有43個。在KOG數(shù)據(jù)庫中1 896個轉錄本被注釋，按照功能一共分成24類，其中一般功能預測最多，其次是翻譯后修飾、蛋白折疊和分子伴侶，之后為信號轉導機制。

2.3.3 轉錄因子分析 本研究預測到的轉錄因子(TF)有1 918個，隸屬于84個TF家族，將注釋到轉錄本數(shù)目最多的前30個TF家族進行柱形圖展示，分析結果如圖4。在獲得的TF家族中，bHLH家族有96個、C3H家族有80個、bZIP家族有70個、MYB-related家族有66個、C2H2家族有66個、MYB家族有57個。

2.3.4 LncRNA分析 經(jīng)過篩選后，依據(jù)在基因組的位置對最終得到LncRNA進行分類，并對占比情況進行展示，其中基因間區(qū)(LincRNA)最多，含有569個，占比73.14；完全位于蛋白編碼基因的intron區(qū)(Sense-intronic LncRNA)次之，有141個，占比為18.12%；與蛋白編碼基因的exon區(qū)有overlapping(Sense-overlapping)有42個，占比5.4%；反義鏈(Anti-sense_LncRNA)有26個，占比3.34%(圖 5-A)；各軟件預測為noncoding的轉錄本條數(shù)畫成維恩圖，結果展示(圖5-B)。

3 討論與結論

目前，PacBio-Iso-Seq測序手段已在多種植物中得到了廣泛的應用，該方法相比于二代測序具有讀長長、無需組裝轉錄組就可以直接獲得全長轉錄本、更準確地檢測結構變異和適合于具有高重復序列和高雜合性的物種的轉錄組分析等優(yōu)勢[11-13]。本試驗利用Iso-Seq技術，結合生物信息學分析，獲得了苦茶全長轉錄組數(shù)據(jù)。

本研究共得到8 730 808個Subreads，平均Subreads的長度為1 145 bp，N50為1 741，可以看出全長轉錄組測序讀長長且連續(xù)性較高。全長非嵌合(FLNC)reads有217 701個，通過對FLNC-reads的聚類及校正分析，最終獲得polished consensus序列132 334個。龐丹丹等[17]對不同葉色的6份茶樹材料進行二代測序，共獲得112 233個Unigene，平均長度為759 bp，N50為 1 081 bp。Li等[18]采用RNA-seq對龍井43的根、莖、4個不同嫩度的葉片及花蕾和果實等13個組織進行測序分析，共獲得347 827條Unigene(>200 bp)，平均長度為791.2 bp，N50為1 342 bp。上述結果表明在獲得的序列質量和基因數(shù)等方面，三代測序結果均優(yōu)于二代測序。

苦茶的全長轉錄組測序在很大程度上補充了現(xiàn)有茶樹基因資源，并為發(fā)現(xiàn)新的或以前未被識別的蛋白質的編碼基因和轉錄亞型提供了優(yōu)勢。與茶樹參考基因組比對，分析發(fā)現(xiàn)大多數(shù)PacBio轉錄本是已知基因的新亞型(59.32%)，有 16.12%的轉錄本為7 115個新基因提供了證據(jù)(圖3)。分析表明，新的轉錄組數(shù)據(jù)對探究苦茶的注釋具有巨大的潛力。Unmapped轉錄本在NR數(shù)據(jù)庫中的比對結果顯示，共比對上311個物種，比對到最多的物種仍是茶樹，其次是葡萄，這表明比對所用到的參考基因組可能還需進一步完善。比較KEGG(10 976+5 626)、KOG (6 296+1 896)、GO(6 221+575)3大數(shù)據(jù)庫功能統(tǒng)計情況，發(fā)現(xiàn)注釋到KEGG的轉錄本數(shù)最多(括號內數(shù)值為未比對到基因組的轉錄本和新轉錄本數(shù)量，下同)，這些轉錄本中包含多個轉錄本與茶葉品質和苦茶特征性狀形成相關的代謝途徑，如苦茶堿等嘌呤生物堿代謝(74+17)、氨基酸代謝(409+69)、苯丙烷生物合成(70+18)和萜類物質的生物合成(152+31)等。Wang等[19]對3個苦茶品種進行轉錄組測序，通過與嘌呤生物堿代謝的KEGG通路富集分析，篩選出一批與苦茶堿和咖啡堿生物合成相關的關鍵基因。同樣，本研究新注釋到的一些轉錄本亦能夠為后期對苦茶性狀形成的研究提供基因資源。

可變剪切(Alternative splicing，AS)，即某些基因的一個mRNA前體以不同的剪接方式，形成不同的mRNA異構體[20]；其能夠調控基因的表達水平和促進蛋白質組的多樣性，在植物生長、發(fā)育和脅迫響應中發(fā)揮著關鍵作用[21-23]。選擇性剪接的亞型具有基于組織或時間的優(yōu)先表達特征，并且它們也受環(huán)境條件的影響[24]。在茶樹中，研究發(fā)現(xiàn)可變剪切亞型在高溫、干旱、低溫等不同環(huán)境脅迫條件下的表達模式存在明顯差異[25-26]，同時可變剪切事件還參與了黃酮、類黃酮和花青素等重要物質的次生代謝過程中[24, 27-28]。因此，對基因結構進行分析，發(fā)現(xiàn)4 892個基因發(fā)生了選擇性剪接，其中內含子滯留事件占比最大，這為進一步研究可變剪切在苦茶特異成分代謝中的作用奠定了基礎。

以往的研究表明，轉錄因子(TF)參與植物多種生長代謝過程中。在茶樹中，發(fā)現(xiàn)它們不僅參與葉片花等器官的發(fā)育[29-30]，以及脅迫反應[31]，還參與兒茶素[17]、茶氨酸[32]和花青素[33-34]等多種次生代謝過程。本研究還預測得到1 918個TF，在獲得的TF家族中，bHLH家族有96個、C3H家族有80個、bZIP家族有70個、MYB-related家族有66個、C2H2家族有66個、MYB家族有57個。這些轉錄因子的預測與分析，為之后研究苦茶特異性狀(苦茶堿[3-4]、丁子香酚甙[2]和高EGC[9]等)相關基因表達以及詳細的基因家族分析提供數(shù)據(jù)。此外，本研究還預測得到778個LncRNA，它們可能在茶樹次生代謝的調控中起著特殊的作用。

本研究為苦茶品質和特異性狀形成機制的研究提供了基礎數(shù)據(jù)，上述結果將有助于進一步開展苦茶特異性狀相關基因標記的開發(fā)、并為其形成的機理及遺傳機制的研究奠定基礎。