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

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

0 引言

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

1 材料與方法

1.1 試驗(yàn)材料

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

1.2 表觀性狀觀測(cè)及綠茶制作與審評(píng)

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

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

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

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

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

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

1.5 差異表達(dá)基因篩選及富集分析

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

1.6 差異表達(dá)基因的實(shí)時(shí)熒光定量PCR（qPCR）檢測(cè)

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

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

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

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

2 結(jié)果與分析

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

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

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

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

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

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

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

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

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

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

2.4 GO功能富集分析

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

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

2.5 KEGG代謝通路分析

對(duì)春綠與福云6號(hào)樣本進(jìn)行KEGG代謝通路分析（圖3），有8 657個(gè)基因注釋到191個(gè)代謝通路中，這些基因中有1 277個(gè)顯著差異基因。其中Unigenes排名前20代謝通路有：40條DEGs參與卵母細(xì)胞減數(shù)分裂（ko04114），45條DEGs參與RNA降解（ko03018），29條DEGs參與脂肪酸代謝（ko01212），17條DEGs參與不飽和脂肪酸的生物合成（ko01040），40條DEGs參與藥物代謝-細(xì)胞色素P450（ko00982），40條DEGs參與細(xì)胞色素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號(hào)KEGG代謝通路富集顯著差異基因表達(dá)模式Fig. 3 Top enriched KEGG pathways of DEGs in CL and FY6H

2.6 差異表達(dá)基因的qPCR檢測(cè)

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

圖4 差異基因的表達(dá)驗(yàn)證Fig. 4 Expression verification of differential genes

3 討論與結(jié)論

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

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

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

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