基于第三代納米孔測序技術的茶樹全長轉錄組分析

王 芬，木 仁，張 雯，馬 媛，郭治友，崔寶祿，李 靜

（黔南民族師范學院生物科學與農學院，貴州 都勻 558000）

0 引言

【研究意義】茶樹Camellia sinensis（L.）O.Kuntze屬山茶科、山茶屬，異花授粉植物，起源于中國西南一帶，為葉用經濟作物，作為受大眾歡迎的世界性飲品，具有非常強的保健作用[1-2]。茶產業作為中國的優勢特色產業，在助力貴州農業農村鄉村振興方面具有重要的經濟價值[3]。伴隨著茶葉貿易全球化，消費者對茶葉質量的追求日益提高，提高茶葉品質是發展茶產業的重中之重。都勻毛尖茶在貴州省黔南州大面積推廣，種植面積高達160萬畝，40萬群眾以茶為業。都勻毛尖原料大部分為都勻本地種與引進的福鼎大白茶，福鼎大白茶發芽率高、抗旱性強，是十分優良的茶樹遺傳資源[4]。茶中的茶多酚、氨基酸和咖啡堿是茶葉滋味和品質的主要組成成分，其中茶多酚和氨基酸具有多種保健功能[5-7]。目前基于二代轉錄組測序技術對茶樹的研究已經得到廣泛應用[8-10]，但由于其讀長短，拼裝困難，導致轉錄本的拼接組裝不完整，而三代牛津納米孔技術（Oxford Nanopore Technology, ONT）憑借其超長讀長的優勢可以彌補二代測序技術的不足，在轉錄本等信息的識別上更全面。基于ONT全長轉錄組測序可以從理論水平上更全面地挖掘與茶葉品質、獨特芳香以及產量相關的代謝通路和生物過程?！厩叭搜芯窟M展】目前，單分子ONT技術已經在動物、植物、微生物等領域得到應用。Jansen等[11]利用ONT技術對歐洲鰻鱺的基因組進行了測序，與之前的草圖相比準確性大大增加。Fellers等[12]通過ONT技術對感染小麥組織進行測序，證實了小麥條紋花葉病毒的存在，該結果表明ONT技術可以更準確地識別病原體。Giordano等[13]證明通過單分子ONT平臺測序的數據足以完整地組裝釀酒酵母S288C菌株。同時，一些茶樹的基因組相繼被測序，如云抗10號[14]、舒茶早[15-16]、小葉茶碧云[17]、古茶樹[18]等，為改良茶葉品質、提高茶葉產量提供大量數據支持。RNA-Seq分析是解鎖生命密碼非常重要的工具，但由于第二代高通量測序平臺[19]的轉錄組技術測得數據的拼接準確性不高，使得深入理解細胞生命活動困難重重。與二代技術相比，第三代測序技術具有通量高、讀長長、成本低等優點，讀長可達10 kb[20-21]。基于ONT[22]的全長轉錄組測序技術是一種單分子實時測序技術，并且在測序時無須打斷RNA片段，所測即所得[23]，大大縮減生物基因組重構和組裝的時間和成本，為轉錄組學的研究減少阻礙。第三代測序技術的錯誤率雖然比二代高，大約為15%[24]，但是利用更正軟件加大測序深度可以大大降低錯誤率，可使準確率達到99.9%[25]。此外，龐丹丹等[26]利用PacBio三代測序技術對苦茶全長轉錄組進行分析，研究結果為探索苦茶特異性狀相關基因標記的開發奠定基礎。夏麗飛等[1]通過PacBio平臺對紫娟茶樹全長轉錄組進行分析，為開展紫娟茶樹葉片呈色機理提供數據支持。【本研究切入點】目前，利用三代納米孔測序技術研究茶葉的品質和滋味鮮有報道。【擬解決的關鍵問題】本研究利用三代測序平臺ONT技術對都勻福鼎大白茶葉、根和莖進行全長轉錄組測序比較分析，旨在探究與茶葉品質相關的差異基因和代謝通路，為后續分子生物學研究提供數據參考。

1 材料與方法

1.1 試驗材料

采摘貴州省都勻市黔南民族師范學院試驗基地長勢一致的福鼎大白茶扦插苗9株，分為3組，每組3盆，采摘嫩葉、嫩根和嫩莖分別作為葉、根和莖處理，每個組織3個生物學重復，共9個樣本，包括葉片L1、L2、L3，根R1、R2、R3，莖S1、S2、S3，放入液氮中進行固樣。將采集的樣本放入干冰中送往北京百邁客生物科技有限公司進行全長轉錄組測序分析。

1.2 RNA提取、文庫構建

利用北京天根生物技術有限公司的RNAprep Pure多糖多酚植物總RNA提取試劑盒（DP441）對都勻福鼎大白茶葉、根和莖的RNA進行提取。采用Nanodrop、Agilent2100及Agilent RNA 6000 Nano Kit對RNA的濃度和完整性進行檢測，檢測合格的RNA樣品，使用oligod（T）磁珠從TotalRNA中純化出poly（A）+RNA。其次利用Superscript Ⅳ reverse transcriptase反轉錄合成cDNA第一鏈，再使用帶barcode的引物及LongAmpTaq2*Master Mix進行PCR擴增。然后使用NEBNext End repair/dA-tailing Module進行末端修復及加A。再次使用ONT SQKLSK109試劑盒及NEBNext Quick Ligationg Module進行測序接頭的連接。最后使用PromethION測序儀及PromethION Flow Cells 9.4進行測序，并將三代轉錄組測序數據提交至NCBI-SRA數據庫，BioProject的編號為PRJNA562747。

1.3 全長轉錄組測序

全長轉錄組測序平臺使用Oxford Nanopore Technologies的PromethION，將原始下機序列中長度小于500 bp、Qscore小于7的低質量序列和核糖體RNA序列過濾掉，根據兩端是否存在引物得到全長序列，對全長序列進行polish獲得一致性序列，然后與安徽農業大學第一版的中國種茶樹[27]基因組或構建的contig進行比對，將identity和coverage的值分別設置為0.9和0.85，去除冗余，得到轉錄本序列，再利用gffcompare v0.9.8將全長轉錄本與基因組已知的轉錄本進行比較，獲得新基因和新轉錄本。隨后進行SSR分析、ORF預測、轉錄因子分析、lncRNA分析，并且將所有的轉錄本、新轉錄本、新基因、開放閱讀框、轉錄因子和LncRNA都上傳到FigShare數據庫，DOI號分別為：10.6084/m9.figshare.13671901；10.6084/m9.figshare.14370011。將Fold Change≥2且FDR＜0.01作為篩選差異表達基因的標準，預測出葉與根、葉與莖、莖與根的差異轉錄本。最后，應用軟件Blast2GO v2.5對差異轉錄本進行GO注釋，并通過將轉錄本的蛋白序列和KEGG數據庫中收錄的蛋白序列進行BLAST比對得到KEGG功能注釋信息，繼而使用軟件blast v2.2.31將轉錄本與kog202101數據庫進行比對，得到轉錄本的KOG注釋信息。葉與根、葉與莖、莖與根的差異轉錄本的GO、KEGG和KOG注釋信息上傳至figshare數據庫，DOI號為：10.6084/m9.figshare.13671901。

1.4 轉錄本表達定量

應用CPM（Counts per million）[28]計算轉錄本表達量。

式中：R：比對到某一轉錄本上的reads數；T：比對到參考轉錄組的片段總數。

1.5 實時熒光定量PCR驗證

為確保轉錄組數據質量準確性高，利用實時熒光定量PCR對隨機選擇的4個基因進行驗證。試驗材料和取樣方法同1.1。利用Aidlab公司反轉錄試劑盒（TUREscript 1st Stand cDNA SYNTHESIS Kit）進行cDNA的合成。使用Primer5.0軟件設計實時定量PCR引物（表1）。以茶樹的GAPDH（GE651107.1 EST1434）基因為內參基因，利用2-△△Ct計算基因相對表達量，試驗設置3個生物學重復。實時熒光定量PCR反應體系及程序按照Fermentas公司SYBR GREEN I說明書進行。

表1 qRT-PCR引物Table 1 qRT-PCR primer

2 結果與分析

通過Nanopore三代全長轉錄組測序對貴州省都勻市種植的福鼎大白茶葉、根和莖進行測序分析。每個樣品測序產出clean data均達到7.93 GB，并且所有樣本的平均質量值都達到Q9，9個樣品得到的全長序列個數介于3768495～5078770（表2），利用經過polish處理的全長序列與第一個版本的中國種茶樹進行minimap2.1.1[29]比對（表2），提取出69379個轉錄本。然后進行融合轉錄本的預測，9個樣品的融合轉錄本個數為188～248。最后，預測出93102個SSR，獲得7556個新基因位點，65795個新轉錄本，45852個ORF，6335個轉錄因子和2229個lncRNA，并完成了58398個新轉錄本的功能注釋。

表2 Clean data數據Table 2 Clean data

2.1 測序數據分析

2.1.1 SSR分析 利用MISA1.0軟件[30]對福鼎大白茶葉、根和莖的全長轉錄組做SSR檢測，將轉錄本的序列長度≥500 bp作為篩選標準，共檢測到50013條序列具有SSR位點，包含7種類型的SSR，共78931個。其中完美單堿基重復SSR個數為53433，數量最多，其次是完美雙堿基重復，為27961個，然后依次是完美三、四、五、六堿基重復，數量分別為10316、586、254、552個。以搜索標準為1～6個堿基基序重復次數分別≥10、6、5、5、5、5在SSR位點中檢測，單堿基重復出現頻率最高是T/A（38139）。雙堿基重復出現最多的是TC/GA（7356），其次是CT/AG（7226）。三堿基重復最多的是GAA/TTC（598）和CCA/TGG（581）。四堿基重復以TTTA/TAAA（70）和TTAT/ATAA（31）占優勢。五堿基和六堿基重復頻率最高的分別是TGTTA/TAACA（13）和GGTGCT/AGCACC（15）。以上研究結果與紫娟茶樹[31]和藤茶[32]的結果基本一致，為未來開展茶樹遺傳圖譜構建、SSR分子標記開發及培育良種提供理論基礎。

2.1.2 新基因編碼區序列預測 利用TransDecoder5.0.0（https://TransDecoder.sourceforge.net）預測出開放閱讀框58355個，其中完整開放閱讀框45852條。預測的完整開放閱讀框編碼蛋白序列長度范圍主要在0～800氨基酸，0～100氨基酸有30455個，占52.19%，100～200個氨基酸的有22949個，占39.33%，200～300氨基酸有3936個，占6.74%（圖1A），與云南金花茶[33]通過ESTScan預測的CDS長度大體趨勢相似，以上結果表明轉錄組的序列質量較高。

2.1.3 轉錄因子預測 使用iTAK1.6軟件[34]預測都勻福鼎大白茶轉錄因子，共預測到轉錄因子6335個，主要分為20類，其中GRAS家族轉錄因子數量最多，其次是MYB-related、RLK-Pelle_DLSV、WRKY和C3H（圖1B），這些轉錄因子家族成員的獲得為后續分子生物學的研究提理論參考。

2.1.4 LncRNA預 測 分 別 應 用Cpc[35]、Cnci[36]、Cpat1.2.2[37]、Pfam1.6[38]預測lncRNA，4種方法取交集共2229個（圖1C）。lncRNA主要分為基因間lncRNA（LincRNA）（1878/84.3%）、反 義lncRNA（Antisense-lncRNA）（110/4.9%）、內 含 子lncRNA（Intronic lncRNA）（34/1.5%）、正義lncRNA（SenselncRNA）（207/9.3%），其中LincRNA最多，以上數據為將來研究lncRNA在茶中的調控機制提供重要參考。

2.2 功能注釋

為了獲得轉錄本的注釋信息（表3），將得到的新轉錄本進行eggNOG、COG、NR、Pfam、Swissprot、KEGG、GO和KOG注釋。

表3 新轉錄本注釋Table 3 New isoform annotation

2.2.1 轉錄本表達量分析 根據CPM計算出轉錄本的表達量。采用CPM箱線圖從整體上對9個樣品的表達量進行比較，發現他們的表達水平基本上一致（圖1D）。

圖1 開放閱讀框、轉錄因子、長鏈非編碼RNA和轉錄本表達量分析Fig. 1 Analyses on ORF, transcription factor, lncRNA, and transcript expression

2.2.2 差異表達轉錄本 利用DESeq1.18.0[39]進行葉、根和莖轉錄本的差異表達分析，篩選條件為Fold Chang≥2且FDR＜0.01，葉和根的差異表達轉錄本最多，莖和根的差異表達轉錄本最少（表4）。對差異表達轉錄本進行數據庫的功能注釋，共注釋了14306個差異表達轉錄本，分別有9649、2825和7349個轉錄本注釋到GO、KEGG和KOG數據庫中。

表4 差異表達轉錄本注釋Table 4 Annotation of DETs

2.2.3 差異表達轉錄本GO注釋 對差異轉錄本進行GO注釋，一級分類主要包含3個類型，分別為生物過程、細胞組分和分子功能；二級分類主要分為52個類別。在葉與根的差異轉錄本中生物學過程涉及21個類別，以代謝過程（總轉錄本/差異轉錄本，29572/5019）、細胞過程（26938/4021）、單生物過程（19040/3581）最多。細胞組分包括16個功能組，其中細胞（27641/4454）、細胞部分（27478/4411）、膜（20620/3584）最多。分子功能分為15個功能類別，催化活性（30151/5408）、結合（24530/3744）、轉運活性（3340/651）最多（圖2）。葉與莖、莖與根的差異轉錄本的GO、KEGG和KOG注釋信息上傳至figshare數據庫，DOI號為：10.6084/m9.figshare.13671901。研究結果豐富了都勻福鼎大白茶的分子生物學信息，可為進一步研究茶葉的品質提供理論參考。

2.2.4 差異表達轉錄本KEGG注釋 KEGG數據庫共注釋到2825個葉與根的差異轉錄本，其中參與碳代謝（270、9.56%）、氨基酸生物合成（192、6.8%）、苯丙素生物合成（182、6.44%）、淀粉和蔗糖代謝（167、5.91%）、植物激素信號轉導（137、4.85%）的差異轉錄本最多（圖3）。參與茶葉品質形成中涉及茶葉滋味相關代謝途徑的有氨基酸生物合成、各種氨基酸代謝、類黃酮生物合成和苯丙素生物合成。與香氣相關的代謝途徑有泛醌和其他萜烯醌生物合成和萜類骨架生物合成。研究結果為后期深入開展提高茶葉品質的研究提供理論基礎。

圖3 葉與根差異表達轉錄本KEGG功能分類Fig. 3 KEGG function annotation of DETs between leaf and root

2.2.5 差異轉錄本的KOG注釋 葉與根的差異轉錄本與KOG數據庫進行比對并根據其功能進行分類，有7349個轉錄本得到注釋。按照功能一共分為25類，注釋到一般功能預測（1531）中的差異轉錄本最多，其次是翻譯后修飾、蛋白質轉換、伴侶蛋白（924），次級代謝生物合成運輸和分解代謝（812）。注釋到KOG數據庫最少的葉與根差異轉錄本數是染色質結構和動力學（30）、真核細胞的細胞外結構（20）、核結構（6）、細胞運動（0）（表5）。

表5 部分葉與根差異轉錄本的KOG分類Table 5 DETs KOG classification between leaf and root

2.3 差異表達轉錄本qRT-PCR分析

以茶樹的GAPDH（GE651107.1 EST1434）基因為內參基因，隨機選取4個差異表達基因（ONT.4041、ONT.8670、TEA019914、TEA007016）進行qRTPCR驗證（圖4），結果表明其表達量與轉錄組測序結果相一致，進一步說明轉錄組測序質量較高。

圖4 轉錄組測序結果與qRT-PCR表達量比對Fig. 4 Comparison between RNA-seq and qRT-PCR

3 討論

ONT第三代測序技術是近年來新興的單分子測序技術[40]，DNA/RNA鏈以一定的速率通過納米孔通道蛋白時，單個堿基會引起不同電學信號的變化，根據電流信號對序列進行實時測定[41-42]。本文利用ONT全長轉錄組測序技術及生物信息學分析獲得了福鼎大白茶葉、根和莖9個樣品的全長轉錄組。通過CPM分析葉、根和莖9個樣品的轉錄本表達量，發現它們的整體表達水平基本上一致。另外，我們共得到93102個SSR，其中完美單堿基重復和完美雙堿基重復最多。單堿基、雙堿基、三堿基和四堿基重復出現頻率最高，分別是T/A、CT/AG、GAA/TTC、TTTA/TAAA，與朱興正等[43]和鞠燁等[44]的研究結果基本一致，為未來開展茶樹遺傳圖譜構建、SSR分子標記開發及培育良種提供理論基礎。預測的45852個完整ORF序列，蛋白序列長度范圍主要在0～800氨基酸，0～300個氨基酸最多，與夏麗飛等[1]和潘敏等[45]的研究結果大體相似，表明轉錄組的序列質量較高。我們還對轉錄因子、LncRNA進行了分析，研究結果為將來研究LncRNA在茶中的調控機制提供重要參考，為后續分子生物學的研究提供理論參考。

根據GO功能顯著性富集分析，上調基因和下調基因在生物過程中，代謝過程（2162/2857），即葉與根中有2162差異轉錄本在生物過程中上調，2857個轉錄本在此過程中下調。單生物過程（1688/1893），細胞過程（1673/2348）、生物調節（1673/633）等顯著富集。在細胞組分中，細胞（1623/2831）、細胞部分（1608/2803）、膜（1597/1987）、膜部分（978/1201）等顯著富集。在分子功能中，催化活性（2496/2912）、結合性（1771/1973）、轉運活性（331/320）等顯著富集。以上結果是葉與根差異轉錄本數最多的生物過程，朱興正[43]等采用PacBio平臺分析了保護品種云茶1號的全長轉錄組，Unigene的GO功能注釋與本研究結果相一致。同時，以上結果顯示在這些生物過程中大部分下調轉錄本數都比上調轉錄本多，即在都勻福鼎大白茶樹生長發育過程中，參與葉生物過程的轉錄本比根多。

在KEGG中，葉與根的差異轉錄本中有182個差異基因參與苯丙素的生物合成，KEGG直系同源基因ID號為ko00940。對該通路的苯丙氨酸、絡氨酸、色氨酸的生物合成到木質素的通路進行深入分析，涉及6個節點，其中4個節點中的基因全部上調，2個節點既有上調又有下調基因。將涉及的所有基因根據基因表達量進行聚類分析，TEA019411、TEA017242、TEA028682、 TEA029025、TEA006839、TEA017067、ONT.2421、TEA008747、TEA031671、TEA005749等基因在葉中高表達，都屬于K00430，并且親緣關系較近，與過氧化物酶有關，參與氧化應激反應。過氧化物酶與光合作用、呼吸作用有關，并能使組織中所含的某些碳水化合物轉化成木質素，增加木質化的程度，是組織老化的一種標志。在氨基酸生物合成中，有192個差異基因參與，其中涉及谷氨酸合成酶的基因有ONT.24127.2和TEA003892.1等，并且這兩個差異基因在根中表達上調。茶氨酸是在茶樹根部由乙胺和谷氨酸在酶的催化下合成的，然后運輸到葉部，參與代謝過程。而茶氨酸的降解與陽光有關，影響茶氨酸向兒茶素轉化，因此，茶氨酸的合成與降解與茶葉品質密切相關。以上結果表明，茶樹葉片通過光合作用等一系列反應形成與茶葉品質相關的化合物。

葉與根上調和下調差異轉錄本KOG分類中，一般功能預測（611/920），上調轉錄本611，下調轉錄本920。翻譯后修飾、蛋白質轉換、伴侶蛋白（349/575），碳水化合物的運輸和代謝（339/405）的分類中，大部分下調轉錄本比上調轉錄本多，再一次說明茶樹生長過程中葉的重要性。本研究的全長轉錄組數據和各項結果豐富了茶樹的分子生物學信息，可為進一步研究茶葉的品質和獨特的芳香及與之相關的基因提供理論基礎。為保障茶葉產量和品質打下良好基礎，為后續良種選育及分子生物學研究提供重要的數據支持。