基于RNA-Seq分析茶樹新品系‘春綠’生物信息

林鄭和，孔祥瑞，陳常頌,2*，李鑫磊，張亞真，游小妹，陳志輝，單睿陽，鐘秋生

（1. 福建省農(nóng)業(yè)科學(xué)院茶葉研究所/國家茶樹改良中心福建分中心，福建 福州 350013；2. 國家土壤質(zhì)量福安觀測實驗站，福建 福安 355015）

0 引言

【研究意義】伴隨著高通量測序新技術(shù)（RNA-Seq）的廣泛應(yīng)用，越來越多的研究者利用轉(zhuǎn)錄組學(xué)研究茶樹生長發(fā)育、抗性、物質(zhì)代謝的分子調(diào)控機制，探索茶樹基因表達的轉(zhuǎn)錄調(diào)控機制，并試圖通過基因工程技術(shù)培育出優(yōu)良品種?！厩叭搜芯窟M展】昆明植物研究所的研究團隊發(fā)布了茶樹大葉種‘云抗10號’基因組[1]，開啟了茶樹基因組學(xué)的研究序幕，之后安徽農(nóng)業(yè)大學(xué)的茶學(xué)團隊基于三代PacBio測序方法，獲得了茶樹小葉種‘舒茶早’的基因組序列[2]；華南農(nóng)業(yè)大學(xué)也公布了茶樹小葉種‘碧云’的染色體級別基因組序列[3]；同時中國農(nóng)業(yè)科學(xué)院茶葉研究所也公布了‘龍井43’的染色體級別基因組序列[4]；之后‘黃棪’與‘鐵觀音’染色體級別基因組序列的公布[5,6]，探索茶樹基因表達的轉(zhuǎn)錄調(diào)控機制，并試圖通過基因工程技術(shù)培育出優(yōu)良品種，現(xiàn)已成為茶樹研究新興熱點之一。李娜娜等[7]以‘福鼎大白茶’和‘小雪芽’葉片為研究材料，發(fā)現(xiàn)小雪芽新梢白化的可能遺傳機理存在于蛋白質(zhì)翻譯后的修飾過程，而且這些變異主要發(fā)生在細胞內(nèi)。朱興正等[8]采用PacBio平臺進行‘云茶1號’全長轉(zhuǎn)錄組測，發(fā)現(xiàn)云茶1號中216個代謝途徑分支，包括茶葉品質(zhì)、活性物質(zhì)代謝以及抗逆等相關(guān)基因。然而，茶樹龐大的3 Gb大小基因組、超過70%的重復(fù)序列含量以及自交不親和導(dǎo)致的高度雜合特性[7]，遺傳結(jié)構(gòu)十分復(fù)雜，遺傳分析很困難。【本研究切入點】‘春綠’茶樹新品系由福建省農(nóng)業(yè)科學(xué)院茶葉研究所于2012—2022年，以‘福云6號’為母本、‘早春毫’為父本雜交選育而成，小喬木型，樹姿直立，生長勢中等。一芽一葉期3月7日；一芽二葉期3月10日；芽葉綠色，芽葉茸毛中等或較多，一芽三葉長8.2 cm，發(fā)芽較密，一芽三葉百芽重88.6 g。葉片為長橢圓形，葉面微隆，葉片質(zhì)地軟，葉尖漸尖，葉基楔形，葉緣微波，葉齒銳，葉齒密度中。春茶一芽二葉蒸青樣水浸出物含量49.9%、茶多酚含量22.4%、游離氨基酸總量3.4%、咖啡堿含量4.0%。適制綠茶和白茶，在福建、湖北等多地有種植示范，是一個適制性廣、高香且特早生優(yōu)良茶樹品系。【擬解決的關(guān)鍵問題】本研究利用轉(zhuǎn)錄組測序技術(shù)，分析春綠與福云6號的轉(zhuǎn)錄組測序，通過序列分析、基因功能注釋、基因功能分類、差異表達基因篩選等，旨在為進一步挖掘春綠茶樹新品系次生代謝相關(guān)功能基因和抗逆基因奠定基礎(chǔ)，為新品種的選育提供參考依據(jù)。

1 材料與方法

1.1 試驗材料

供試茶樹品種福云6號與新品系春綠，種植于福建省農(nóng)業(yè)科學(xué)院茶葉研究所二號山同一地塊，樹齡均為5年，施肥、噴藥等茶園管理均一致。春季統(tǒng)一采摘新梢第一輪一芽二葉，迅速用液氮固樣，送上海歐易生物醫(yī)學(xué)科技有限公司測序。

1.2 表觀性狀觀測及綠茶制作與審評

表觀性狀按陳亮等[9]的標準觀測。采摘2022年春茶第一批新梢（達到一芽一葉標準），綠茶制作參照傳統(tǒng)烘青綠茶加工工藝：鮮葉萎凋→殺青→揉捻→烘干。詳細工藝參數(shù)為：采用自然萎凋方式（水篩攤涼，1.5斤/個水篩），萎凋時間約6 h；采用滾筒殺青，溫度為280℃，投葉量為3斤/鍋（2個水篩的鮮葉），時間約90 s。出鍋后攤涼10～15 min（待熱氣散發(fā)后），放入6CR-15的揉捻機揉捻10～12 min，迅速解塊，毛火120℃約10 min，足火80℃烘箱中烘干至含水率低于5%。

樣品感官審評參照GB/T 23776—2018《茶葉感官審評方法》，審評項目有外形、香氣、湯色、滋味和葉底五項。審評結(jié)果采用加權(quán)評分法，湯色10%、外形占20%、滋味30%、香氣30%、葉底10%。再結(jié)合記錄的評審術(shù)語，綜合評定品質(zhì)。由4名專業(yè)審評人員進行密碼審評。

1.3 RNA的提取與文庫構(gòu)建

茶樹葉片總RNA的提取參考本實驗室改良的Tri-2Reagent法，凝膠電泳檢測RNA完整性；測定260 nm和280 nm波長下的OD值，合格后進行混勻。RNA混勻后使用DNase消化DNA，用帶有Oligo（dT）的磁珠富集mRNA（去除rRNA）；加入打斷試劑將mRNA打斷成短片段，用六堿基隨機引物合成一鏈cDNA，然后配制二鏈合成反應(yīng)體系合成二鏈cDNA，并使用試劑盒純化雙鏈cDNA；純化的雙鏈cDNA再進行末端修復(fù)、加A尾并連接測序接頭，然后進行片段大小選擇，最后進行PCR擴增；構(gòu)建好的文庫用Agilent 2100 Bioanalyzer質(zhì)檢合格后，使用Illumina HiSeqTM 2500測序儀進行測序，產(chǎn)生125 bp或150 bp的雙端數(shù)據(jù)。具體參照林鄭和等[10]的方法。

1.4 Illumina測序、數(shù)據(jù)質(zhì)檢和參考基因組比對

通過Illumina HiSeq 測序平臺，對文庫片段進行雙末端（Paired-end，PE）測序，得到原始下機數(shù)據(jù)fastq 文件。使用FastQC 軟件的默認參數(shù)對fastq 文件進行測序質(zhì)量統(tǒng)計，使用Cutadapt 軟件進行reads的過濾和去除低質(zhì)量序列（overlap≤10 bp，20%的堿基錯誤率），最終得到高質(zhì)量的clean reads。再通過 Bowtie2 和 Tophat2 將 clean reads 比對至茶樹參考基因組[2]，得到SAM/BAM 比對結(jié)果文件，具體方法參照文獻[11]。

1.5 差異表達基因篩選及富集分析

分別以各測序樣本的總表達量為內(nèi)標，對2個測序樣本的RPKM進行Fisher-test差異檢驗。以FDR（錯誤發(fā)現(xiàn)率）≤0.05及Fold-Change（表達差異倍數(shù)）≥2為條件進行差異基因篩選，并進行差異轉(zhuǎn)錄組的GO和KEGG富集分析，以判定差異轉(zhuǎn)錄組主要影響的生物學(xué)功能或者通路。

1.6 差異表達基因的實時熒光定量PCR（qPCR）檢測

為了驗證RNA-Seq數(shù)據(jù)的結(jié)果，隨機選取14個ungene進行RT-qPCR分析。RT-qPCR引物對采用 primer Premier 5.0（Premier Biosoft，Palo Alto，CA，USA）設(shè)計，內(nèi)參基因為GAPDH（CSA024857.1），見表1。用RNA純化試劑盒（中國天根）純化總RNA，按照說明書使用TaKaRa（中國大連）的 PrimeScript II 1st Strand cDNA synthesis kit進行第一鏈cDNA的合成。使用 SYBR Premix Ex Taq II（Tli RnaseH Plus）試劑盒（TaKaRa，日本）在Lightcycler 480 Real-Time PCR檢測系統(tǒng)（Roche，德國）中進行反應(yīng)。采用3 μL 模板 cDNA、1 μL 正向引物（5 pmol）、1 μL反向引物（5 pmol）和 5 μL SYBR Green混合物（Qiagen，Hilden，Germany） 進 行Real-time qPCR。PCR過程包括95℃預(yù)熱5 min，95℃ 10 s，60℃ 20 s，45個循環(huán)，然后生成解離曲線（95℃15 s，60℃ 60 s，95℃ 15 s）。每個實驗設(shè) 3個重復(fù)，用2-ΔΔCt法計算各基因的相對表達量。對不同樣品進行歸一化處理，以GAPDH基因（accession number：CSA024857.1）為內(nèi)標，以福云6號（CK）的葉片為參比樣本，設(shè)置其基因表達量為1。具體的方法參照文獻[10]。

表1 引物序列及特征Table 1 Primer pair used for RT-qPCR analysis

1.7 數(shù)據(jù)處理

利用軟件SIMCA-P 11.5對所測定的數(shù)據(jù)進行分析，篩選出差異基因，并對兩組樣本進行主成分分析（PCA）。

2 結(jié)果與分析

2.1 表觀性狀與品質(zhì)差異分析

從表2可以發(fā)現(xiàn)，春綠葉片呈稍上斜，福云6號的呈水平或者稍下垂，二者有明顯的區(qū)別；葉身內(nèi)折也較福云6號大，葉齒銳深與福云6號的鈍淺有明顯區(qū)別，嫩梢茸毛較福云6號少，萌發(fā)期大約晚3 d。從表3可看出，新品系春綠香氣與滋味得分明顯高于福云6號。春綠品質(zhì)明顯優(yōu)于福云6號。

表2 表觀性狀分析Table 2 Apparent trait analysis

表3 感官品質(zhì)審評表Table 3 Sensory evaluation form

2.2 測序原始序列數(shù)量及質(zhì)量

測序共得到春綠與福云6號品種（系）葉片轉(zhuǎn)錄組原始序列，通過計算每千個堿基的轉(zhuǎn)錄每百萬映射RPKM（Reads Per Kilobase of exon model per Million mapped reads）的讀取量來分析RNASeq數(shù)據(jù)的轉(zhuǎn)錄譜。構(gòu)建了12個cDNA文庫，包括福云6號葉片（FY6H1-6）與春綠（CL1-6）各6個生物重復(fù)，共獲得81.51 G的Clean data，各樣本的有效數(shù)據(jù)量分布在6.59～7.2 G，Q30 堿基分布在89.6%～90.86%，從每個庫生成的原始讀取數(shù)從47 385 674到51 271 456。Clean reads和Q20（測序錯誤率低于1%）的百分比分別超過97%和95%（表4），堿基G和C數(shù)占總堿基數(shù)的44%以上，平均GC含量為45.00%。通過將reads比對到參考基因組上，得到各個樣本的基因組比對情況，比對率為90.91%～91.52%，不能測序的核苷酸N為0，說明轉(zhuǎn)錄組測序質(zhì)量整體較高。

表4 測序數(shù)據(jù)統(tǒng)計結(jié)果Table 4 Summary of RNA-Seq data on Fuyun No. 6 (FY6H) and Chunlv (CL)

2.3 樣本的PCA分析與差異基因分析

從圖1-A可以看出，該組原始試驗數(shù)據(jù)所得在兩個主成分PC1、PC2中良好地呈現(xiàn)，第一主成分的貢獻率為99.26%，第二主成分的貢獻率為0.22%，兩種貢獻率和為99.46%。通過樣品分布散點圖發(fā)現(xiàn)，2個品種的信息可以得到很好的區(qū)分，說明擬合性良好，能夠直觀地反映樣本間的差異性。通過對2個品種的基因測序，以福云6號為對照，在錯誤發(fā)現(xiàn)率（false discovery rate，F(xiàn)DR）＜0.05且|log2 fold change（FC）|＞1的篩選條件下，發(fā)現(xiàn)春綠中有7 042個差異基因（DEGs），其中3 273個表達上調(diào)，3 769個表達下調(diào)，見圖1-B。

圖1 春綠與福云6號PCA分析（A）與差異基因分析（B）Fig. 1 Principal component and differential gene analyses on CL and FY6H

2.4 GO功能富集分析

春綠與福云6號樣本的DEGs主要分為三類，包括23個基于生物過程的GO組、20個基于細胞成分的GO組和21個基于分子功能的GO組（圖2）。通過GO（Gene Ontology）富集分析表明，差異顯著的生物過程功能基因中（前10）有68條基因顯著上調(diào)、268條顯著下調(diào)；參與細胞功能組成基因中有67條顯著差異基因上調(diào)、1 444條顯著差異下調(diào)；參與分子功能基因中發(fā)現(xiàn)220條顯著差異基因上調(diào)、623顯著差異基因下調(diào)。以上發(fā)現(xiàn)，參與細胞功能組成差異基因最多。GO項富集分析顯示，20個DEGs在生物堿代謝過程（GO：0009820）富集，26個DEGs在營養(yǎng)活動（GO：0045735）中富集，393個DEGs在ATP結(jié)合中富集（GO：0005524），13個DEGs富集在次生代謝過程（GO：0019748），3個DEGs富集在精胺酸分解過程(GO：0046208)，28個DEGs富集在葡糖基轉(zhuǎn)移酶活性（GO：0035251），3個DEGs富集在苯甲酸葡萄糖基轉(zhuǎn)移酶活性(GO：0052641)，8個DEGs富集在水楊酸糖苷轉(zhuǎn)移酶(GO：0052640)。28個DEGs富集在葡糖基轉(zhuǎn)移酶活性。

圖2 春綠與福云6號差異顯著富集GO分類Fig. 2 Significantly enriched GO groups in CL and FY6H

2.5 KEGG代謝通路分析

對春綠與福云6號樣本進行KEGG代謝通路分析（圖3），有8 657個基因注釋到191個代謝通路中，這些基因中有1 277個顯著差異基因。其中Unigenes排名前20代謝通路有：40條DEGs參與卵母細胞減數(shù)分裂（ko04114），45條DEGs參與RNA降解（ko03018），29條DEGs參與脂肪酸代謝（ko01212），17條DEGs參與不飽和脂肪酸的生物合成（ko01040），40條DEGs參與藥物代謝-細胞色素P450（ko00982），40條DEGs參與細胞色素P450外源生物代謝（ko00980），27條DEGs參與氨?；?tRNA生物合成（ko00970），66條DEGs參與苯丙素的生物合成（ko00940），13條DEGs參與氮代謝（ko00910），13條DEGs參與類胡蘿卜素生物合成（ko00906），11條DEGs參與油菜素甾醇生物合成（ko00905），10條DEGs參與生物素代謝（ko00780），37條DEGs參與光合作用生物的碳固定（ko00710），27條DEGs參與甲烷代謝（ko00680），8條DEGs參與苯乙烯退化（ko00643），5條DEGs參與氨基苯甲酸酯退化代謝（ko00627），53條DEGs參與谷胱甘肽代謝（ko00480），14條DEGs參與氨酸代謝（ko00380），23條DEGs參與精氨酸和脯氨酸代謝（ko00330），14條DEGs參與角質(zhì)、亞伯堿和蠟的生物合成（ko00073）。

圖3 春綠與福云6號KEGG代謝通路富集顯著差異基因表達模式Fig. 3 Top enriched KEGG pathways of DEGs in CL and FY6H

2.6 差異表達基因的qPCR檢測

為驗證轉(zhuǎn)錄組測序結(jié)果的可靠性，從轉(zhuǎn)錄組數(shù)據(jù)鑒定到的差異基因中選取了14個基因進行實時熒光定量PCR分析。其將RT-qPCR得到的基因表達趨勢與轉(zhuǎn)錄組測序結(jié)果進行比較，發(fā)現(xiàn)基因的表達趨勢基本是一致的（圖4），說明轉(zhuǎn)錄組數(shù)據(jù)是可靠的。

圖4 差異基因的表達驗證Fig. 4 Expression verification of differential genes

3 討論與結(jié)論

近年來，隨著基因組測序技術(shù)的進步，茶樹研究已經(jīng)步入基因組時代。首先是大葉種茶樹云抗10號的基因組序列，拉開了茶樹基因組學(xué)研究的帷幕[1]，之后基于三代測序補洞的方法測序獲得了茶樹小葉種品種舒茶早的基因組序列[2]，接著烏龍茶品種黃棪[5]和鐵觀音[6]的染色體級別基因組發(fā)布等。傳統(tǒng)的茶樹育種和改良方式（新品種選育需15～20年）已很難滿足目前市場多元化的需求[12]，分子生物學(xué)技術(shù)介入茶樹育種，能加速理想性狀的茶樹種質(zhì)的分子育種。本研究通過對茶樹（春綠與福云6號）進行轉(zhuǎn)錄組測序，獲得了大量轉(zhuǎn)錄因子Unigene，可為后期茶樹新品系（春綠）抗性研究、次生代謝研究、新品種選育等提供理論依據(jù)。

本研究通過測序分別獲得春綠為47 367 063.33條Clean reads，福云6號為48 130 043.67條Clean reads，福云6號比春綠有著更長的序列，通過將reads 比對到參考基因組上，得到各個樣本的基因組比對情況，比對率為90.91%～91.52%，說明轉(zhuǎn)錄組測序質(zhì)量整體較高。通過對差異表達基因的GO富集分析表明，新品系春綠與福云6號在生物過程、細胞組分和分子功能三大類別上均分布有差異基因。其中393個DEGs在ATP結(jié)合中富集（GO：0005524），28個DEGs富集在葡糖基轉(zhuǎn)移酶活性（GO：0035251），26個DEGs在營養(yǎng)活動(GO：0045735)中富集，20個DEGs在生物堿代謝過程（GO：0009820）富集。進一步分析發(fā)現(xiàn)，新品系春綠中66條DEGs參與苯丙素的生物合成（ko00940），其中參與苯丙氨酸與酪氨酸代謝中的大部分基因表達都增加，其中苯基丙酸合成基因（LOC114255969）上調(diào)2.49倍、甘露醇脫氫酶（LOC114256196）上調(diào)3.2倍、β-葡萄糖苷酶（LOC114256477）上調(diào)4.1倍。說明新品系春綠在氨基酸合成、香氣合成等方面明顯強于福云6號。這也進一步說明了新品系春綠制綠茶品質(zhì)明顯優(yōu)于福云6號。

脂肪酸代謝由涉及脂肪酸或與脂肪酸密切相關(guān)的各種代謝過程組成，脂肪酸是屬于脂質(zhì)常量營養(yǎng)素類別的分子家族。這些過程主要可分為：（1）產(chǎn)生能量的分解代謝過程；（2）合成代謝過程，它們作為其他化合物的構(gòu)建塊。本研究發(fā)現(xiàn)，29條DEGs參與脂肪酸代謝（ko01212），大部分都顯著上調(diào)，其中?；o酶A氧化酶（LOC114262883）上調(diào)3.1倍；ACP脫氫酶（LOC114265912）上調(diào)3.6倍。結(jié)合GO富集分析，發(fā)現(xiàn)393個DEGs在ATP結(jié)合中富集（GO：0005524），說明新品系春綠能量代謝明顯比福云6號強。KEGG分析還發(fā)現(xiàn)，14條DEGs參與角質(zhì)、亞伯堿和蠟的生物合成（ko00073），其中蠟粉基因CER1（LOC114261760）上調(diào)3.2倍，有研究發(fā)現(xiàn)CER1基因編碼醛脫羰酶，在烷烴生物合成通路中催化醛脫羰形成烷烴，其突變體中烷烴的含量顯著減少，而過表達CER1基因時，烷烴的含量會增加，器官呈現(xiàn)蠟粉合成減少[13]。還發(fā)現(xiàn)羥基棕櫚酸酯-O-阿魏酰轉(zhuǎn)移酶（LOC114285233、LOC11 4287670）表達下調(diào)，羥基棕櫚酸酯-O-阿魏酰轉(zhuǎn)移酶主要與萌芽相關(guān)，福云6號萌發(fā)期相對較早，這與表觀性狀觀測結(jié)果基本一致。

本研究通過差異轉(zhuǎn)錄因子的聚類分析，有助于從分子水平進一步研究新品系春綠與福云6號的品質(zhì)、性狀差異，建立自己優(yōu)選品種的基因表達信息，為茶樹的本土育種和改良建立分子信息庫，這對提高茶葉品質(zhì)和產(chǎn)量具有重要意義。