基于RNA-Seq技術(shù)的濕加松特定時期產(chǎn)脂差異基因篩選及表達(dá)分析

周晨晨 李義良* 王 哲 趙奮成 鐘歲英 林昌明譚志強(qiáng) 李福明 吳惠姍 郭文冰 廖仿炎 王建革

（1. 廣東省森林培育與保護(hù)利用重點實驗室，廣州 510520；2. 廣東省林業(yè)科學(xué)研究院，廣州 510520；3. 廣東省臺山市紅嶺種子園，臺山 529200；4. 江西農(nóng)業(yè)大學(xué)園林與藝術(shù)學(xué)院，南昌 330045）

濕加松（Pinus elliottii × P. caribaea）是濕地松（P.elliottii）與加勒比松（P.caribaea）的雜交后代，源產(chǎn)于澳大利亞昆士蘭，其在生長量、干形材質(zhì)、抗病蟲害等性狀上均優(yōu)于親本，具有突出的雜種優(yōu)勢［1］。我國于上世紀(jì)90 年代開始濕加松的良種培育工作，培育的濕加松良種干型通直，15年時材積蓄積量可達(dá)225 m3·hm-2［2］。目前，濕加松已在我國廣西、廣東、福建、海南等地廣泛種植，并逐漸北擴(kuò)，成為我國南方培育材脂兼用林和大徑材林的有潛力的重要樹種。濕加松同其他松樹一樣可產(chǎn)生具有防御作用的萜類次生代謝物，人們俗稱為松脂。松脂不僅是松科（Pinaceae）植物抵御昆蟲、病菌、機(jī)械損傷的重要物理及化學(xué)防御手段［3～4］，又是我國林業(yè)經(jīng)濟(jì)的主要商品。我國松脂產(chǎn)量占世界松脂總產(chǎn)量70%左右，年出口創(chuàng)匯達(dá)1 億美元。松脂加工后得到的松香、松節(jié)油可廣泛用于香精香料、合成橡膠、液體生物燃料等方面［5～6］，其關(guān)聯(lián)產(chǎn)業(yè)占國民經(jīng)濟(jì)的1/10。是我國重要的可再生自然資源。然而，近年來由于我國松林面積大幅減少、采脂技術(shù)不規(guī)范、亂采濫采現(xiàn)象嚴(yán)重等原因，導(dǎo)致我國松脂資源極度緊張。為改善我國松脂資源短缺現(xiàn)狀，濕加松產(chǎn)脂優(yōu)良品種的培育工作也逐漸開展，雖然已培育的產(chǎn)脂良種5年連續(xù)割脂，產(chǎn)脂量可達(dá)4.5 萬kg·hm-2［2］，但目前對高產(chǎn)脂良種的研究工作仍處于起步階段，對濕加松個體間出現(xiàn)產(chǎn)脂差異的分子機(jī)制不夠清晰，導(dǎo)致濕加松產(chǎn)脂良種早期選擇困難，短時期內(nèi)很難補(bǔ)缺我國松脂資源短板。

RNA-Seq 技術(shù)已被廣泛使用于林木研究中，可以快速全面反映植物組織在不同方式處理下或不同生長發(fā)育期內(nèi)基因的表達(dá)情況，通過對獲得的轉(zhuǎn)錄本信息進(jìn)行處理可獲得一些重要的功能基因，目前已廣泛使用在植物次生代謝產(chǎn)物、抗病蟲、非生物抗逆等研究中。劉彬等人［7］用接種松材線蟲后表現(xiàn)出易感病及抗病的馬尾松轉(zhuǎn)錄測序，通過數(shù)據(jù)分析得到一系列差異基因并篩選出參與抗蟲的相關(guān)基因。Fox H 等人［8］基于RNA-Seq 技術(shù)，揭示了地中海松（P. halepensis）應(yīng)對干旱的響應(yīng)機(jī)制以及從干旱至恢復(fù)過程中的動態(tài)變化。本研究利用松脂產(chǎn)量受中度至高度的遺傳控制這一特征［9～11］，選取生長量一致但產(chǎn)脂量顯著差異的濕加松個體為材料，通過對特定時間下的樣品進(jìn)行轉(zhuǎn)錄組分析，篩選調(diào)控松脂產(chǎn)量差異的相關(guān)候選基因，利用qRT-PCR 技術(shù)對篩選的基因進(jìn)行表達(dá)量驗證分析，為后續(xù)分子標(biāo)記、產(chǎn)脂功能基因克隆及驗證奠定基礎(chǔ)。

1 材料和方法

1.1 植物材料

試驗材料來源于廣東省臺山市紅嶺種子園1996 年培育的濕加松測定林。參照《松脂采集技術(shù)規(guī)程（LY/T 1694-2007）》，對測定林內(nèi)的155 株濕加松進(jìn)行連續(xù)多次產(chǎn)脂量測定。利用SPSS 軟件分析樹高、胸徑、材積與產(chǎn)脂量的相關(guān)性發(fā)現(xiàn)，四者之間存在較強(qiáng)的正相關(guān)。15 年時樹高、胸徑與產(chǎn)脂量的相關(guān)性分別為0.204、0.389、0.425（見表1）。因此為降低樹體生長對松脂產(chǎn)量產(chǎn)生的影響，最終以材積為控制指標(biāo)，選取材積無顯著差異但產(chǎn)脂量連續(xù)多年穩(wěn)定且有顯著差異的高產(chǎn)脂（月產(chǎn)脂量＞1 350 g）與低產(chǎn)脂（月產(chǎn)脂量＜615 g）濕加松作為試驗材料。

表1 樹高、胸徑、材積與產(chǎn)脂量的相關(guān)性Table 1 Correlation between tree height，DBH and yield of oleoresin

表2 材積、產(chǎn)脂量方差分析Table 2 Analysis of Variance of Volume and yield of oleoresin

為了減少不同基因型導(dǎo)致與產(chǎn)脂性狀無關(guān)的差異，本試驗共設(shè)置高產(chǎn)松脂組、低產(chǎn)松脂組共兩組濕加松為材料，每組設(shè)置6個重復(fù)。高產(chǎn)松脂組各株編號分別為H1、H2、H3、H4、H5、H6；低產(chǎn)松脂組中各株編號為：L1、L2、L3、L4、L5、L6。兩組濕加松材積無差異但產(chǎn)脂量差異顯著（見表2）。選擇樹體產(chǎn)脂高峰的9 月中旬進(jìn)行采樣。由于樹體在受到傷害后需經(jīng)過一段時間相關(guān)基因才會啟動表達(dá)，所以在割脂1 h 后再進(jìn)行取樣。在距離地面1 m 位置取剖面為10 cm×10 cm 的正方形，采用與割脂相同的下降式采脂法，采集樹干次生木質(zhì)部組織并立即放于液氮中保存，凍存后置于-80℃冰箱備用。

1.2 試驗方法

1.2.1 RNA提取

參照TIANGEN 公司提供的多糖多酚植物總RNA（RNAprep Pure Plant Plus Kit）提取試劑盒說明書提取試驗材料的總RNA，用0.5%瓊脂糖凝膠電泳檢測RNA 是否降解及污染，用超微量蛋白核酸分析儀—柏精（柏精-BioDrop μlite）檢測RNA 的純度（OD260/280比值在1.8～2.0）。

1.2.2 文庫構(gòu)建及高通量測序

將提取的RNA 純化并檢測質(zhì)量，對符合測序要求的RNA 構(gòu)建cDNA 文庫。為獲取更加完整的濕加松次生木質(zhì)部信息，利用Illumina HiSeq 2500測序平臺對12 個樣品進(jìn)行RNA 測序并按照高產(chǎn)脂、低產(chǎn)脂樣品劃分混合測序結(jié)果。

1.2.3 基因功能注釋及表達(dá)分析

濕加松為無參基因組，為獲得更全面的數(shù)據(jù)信息，將測序所得raw reads 數(shù)據(jù)質(zhì)控后得到的clean reads 利用Trinity 軟件進(jìn)行拼接組裝。組裝后得到的unigenes 經(jīng)NCBI Blast+與Nr，NT，Pfam，KOG/COG，Swiss-prot，KEGG，GO 共7 個數(shù)據(jù)庫比對，得到相應(yīng)的功能注釋，同時將Trinity 拼接得到的轉(zhuǎn)錄本作為參考序列，采用RSEM 軟件計算各個樣品的基因表達(dá)水平。采用DESeq2 進(jìn)行分析，篩選差異表達(dá)基因條件為閾值padj＜0.05 且|log2FoldChange|＞1。若log2FoldChange＞0 則該基因上調(diào)，反之該差異基因下調(diào)。將得到的差異基因進(jìn)行GO功能分析及KEGG pathway 通路分析。

1.3 實時熒光定量PCR

用天根FastKing gDNA Dispelling RT Super-Mix FastKing 一步法除基因組cDNA 試劑盒對1 μg的總RNA 進(jìn)行反轉(zhuǎn)錄，轉(zhuǎn)錄體系為20 μL。將反轉(zhuǎn)錄得到的產(chǎn)物進(jìn)行10 倍稀釋，使用Applied Biosystems 7500 儀器，熒光染料為ChamQTMUniversal SYBR qPCR Master Mix（Vazyme，China），進(jìn)行實時熒光定量PCR 試驗，各試劑加量及擴(kuò)增程序參照熒光染料說明書，內(nèi)參基因為TUB 基因，其余基因采用2-ΔΔCt算法計算在各個樣品中的相對表達(dá)量。

2 結(jié)果與分析

2.1 RNA-Seq測序、拼接及功能注釋

對高產(chǎn)脂與低產(chǎn)脂植株的總RNA 進(jìn)行測序，測序所得序列經(jīng)過過濾后得到594 108 382 條clean reads，總堿基為89.12 Gb。其中Q30 值達(dá)到了92%以上（見表3），clean reads 進(jìn)行數(shù)據(jù)整合和拼接后，得到270 815 個transcripts 平均長度為934 bp，以及190 514個unigenes，平均長度為1 208 bp，有78 869個unigenes長度大于1 kB，200～500 bp的Genes 數(shù)量最多，N50 長度為1 886 bp，為后續(xù)分析提供了高質(zhì)量的數(shù)據(jù)（見表4～5）。

2.2 基因功能注釋

為更加直觀地研究濕加松轉(zhuǎn)錄本信息，本研究將得到的190 514 條unigenes 在7 大數(shù)據(jù)庫中進(jìn)行注釋，有144 023 條unigenes 至少被1 個數(shù)據(jù)庫注釋，占總量的75.59%，有10 265 條unigenes 在所有數(shù)據(jù)庫中均有注釋，占總量的5.38%（見表6）。其中有87 608 條unigenes 注釋到NR 數(shù)據(jù)庫，43.78%的unigenes 被注釋到北美云杉（Picea sitchensis）中，其次是無油樟（Amborella trichopoda）、水芙蓉（Nelumbo nucifera）等，有36.62%的基因未得到明確的功能注釋（見圖1）。有30 323條unigenes注釋到19 條KEGG 代謝通路中，主要富集的通路有“折疊、分類和降解（Folding，sorting and degradation）”，“碳水化合物代謝（Carbohydrate metabolism）”等。其中有770條unigenes富集到萜類物質(zhì)合成“（Metabolism of terpenoids and poly ketides）”通路上（見圖2）其中有304條unigenes直接參與到萜類物質(zhì)合成過程中，涵蓋了萜類物質(zhì)合成的全過程。

表3 高產(chǎn)脂、低產(chǎn)脂轉(zhuǎn)錄組數(shù)據(jù)信息Table 3 High-yield and Low-yield oleoresin transcriptome data information

表4 拼接長度分布頻數(shù)Table 4 Distribution frequency of splicing length

表5 濕加松All-unigenes長度分布Table 5 Distribution of All-unigenes length of P.elliottii×P.caribaea

表6 濕加松All-unigenes功能注釋Table 6 All-unigenes function annotation of P. elliottii×P.caribaea

2.3 差異基因的篩選及分析

為探究濕加松高、低產(chǎn)脂基因差異表達(dá)的分子機(jī)理，采用DESeq2進(jìn)行分析，共篩選出414個差異基因，其中183 個基因上調(diào)，231 個基因下調(diào)。將414 個差異基因進(jìn)行GO 富集分析。261 條差異表達(dá)unigenes 在GO 數(shù)據(jù)庫中具有功能注釋，顯著性富集的條目群主要是“氧化磷酸化過程”（oxidoreductase activity，acting on）及“雙加氧酶活性”（dioxygenase activity）等。

在生物體內(nèi)，不同基因相互協(xié)調(diào)行使其生物學(xué)功能。為了探究濕加松差異表達(dá)基因主要富集的代謝通道。使用KEGG 數(shù)據(jù)庫，對差異基因進(jìn)行顯著性富集分析，得到97 個差異表達(dá)unigenes的功能定義，被映射到54 條KEGG Pathways 中。差異基因顯著富集的20條代謝通路有內(nèi)質(zhì)網(wǎng)蛋白質(zhì)加工（Protein processing in endoplasmic reticulum）、亞油酸（Linoleic acid metabolism）、苯丙氨酸、酪氨酸和色氨酸的生物合成（Phenylalanine，tyrosine and tryptophan biosynthesis）、剪接體（Spliceosome）、植物病原相互作用（Plant-pathogen interaction）等（見表7）。

2.4 轉(zhuǎn)運與修飾過程對萜類物質(zhì)的調(diào)節(jié)作用

根 據(jù) 差 異 基 因 在Nr，Nt，Pfam，KOG/COG，Swiss-prot，KEGG，GO 七大數(shù)據(jù)庫中的功能釋義發(fā)現(xiàn)，與萜類物質(zhì)合成有關(guān)的基因多表現(xiàn)為轉(zhuǎn)運與修飾過程。

表7 差異表達(dá)基因KEGG功能注釋Table 7 Differential expression genes KEGG function annotation

差異基因中編碼ABC 轉(zhuǎn)運蛋白的基因共6 個（Cluster-37595.96613、Cluster-37595.128007、Cluster-37595.77181、 Cluster-37595.89400、 Cluster-37595.20113、Cluster-37595.127923），其中Cluster-37595.96613、 Cluster-37595.77181、 Cluster-37595.89400、Cluster-37595.20113 基 因 在 高 脂 樹中的表達(dá)量顯著高于低脂樹中的表達(dá)量，Cluster-37595.128007、Cluster-37595.127923 在 低 產(chǎn) 脂 樹中的表達(dá)量更高。編碼細(xì)胞色素P450（Cytochrome P450）蛋白共有3 個基因（Cluster-36047.0、Cluster-31181.0、Cluster-14555.0），3 個基因在高產(chǎn)脂樹中的表達(dá)量均高于低產(chǎn)脂樹（見表8）。

2.5 抗逆蛋白在松脂產(chǎn)量中的調(diào)節(jié)作用

為更加全面探討濕加松產(chǎn)脂差異的形成原因，進(jìn)一步從植物防御系統(tǒng)方面對差異基因進(jìn)行分析。差異基因注釋表明，基因中編碼植物抗病的TIR-NBS-LRR 蛋白有2 個（Cluster-2245.0、Cluster-37595.30441），均在高產(chǎn)脂植株中呈現(xiàn)高表達(dá)。編碼雙功能醇落葉松醇還原酶（PLR）蛋白1 個（Cluster-37595.73116），其在高產(chǎn)脂單株中表達(dá)量更高（見表9）。

表9 TIR-NBS-LRR、PLR的表達(dá)Table 9 The expression of TIR-NBS-LRR and PLR

2.6 qRT-PCR驗證與產(chǎn)脂相關(guān)的候選基因

為了驗證轉(zhuǎn)錄組數(shù)據(jù)結(jié)果可靠性以及與產(chǎn)脂有關(guān)的候選基因，我們選擇與產(chǎn)脂有關(guān)的4個差異表達(dá)基因進(jìn)行qRT-PCR 分析（見圖3）。發(fā)現(xiàn)同一基因在高、低產(chǎn)脂不同單株中表現(xiàn)出差異性。CYP450 基因在H1、H4、H5 的表達(dá)量高于H6、H2、H3，在低產(chǎn)脂植株中趨于穩(wěn)定；PLR 基因在低脂樹中表達(dá)量較為均一，但高脂樹中H1、H2、H6 明顯高于H3、H4、H5。ABC transporter 在高脂樹中H1、H2、H5 的表達(dá)量高于H3、H4、H6，其中H5 與H6的表達(dá)量差異顯著；異戊烯基轉(zhuǎn)移酶在L5 中的表達(dá)量最高，其次是L1、L3 的表達(dá)量最低，但也顯著高于高脂樹中的表達(dá)。GGPS 基因在高產(chǎn)脂、低產(chǎn)脂樹中的表達(dá)趨于一致性，高產(chǎn)脂植株中H1、H2、H3 的表達(dá)顯著高于H4、H5、H6，低脂樹中L5、L6的表達(dá)量最高，L2 中的表達(dá)量最低。總體來看CYP450、ABC transporter、PLR 基因在高脂樹中的表達(dá)量顯著高于低脂樹中的表達(dá)量，異戊烯基轉(zhuǎn)移酶在低脂樹中的表達(dá)量更高，GGPS 基因在高脂樹與低脂樹的表達(dá)量無差異，此結(jié)果與轉(zhuǎn)錄組測序結(jié)果一致。

3 討論

松脂中95%以上是由萜類物質(zhì)組成，因此人們認(rèn)為萜類物質(zhì)合成途徑是松脂合成的主要通路。異戊烯基轉(zhuǎn)移酶是萜類物質(zhì)合成中的一類關(guān)鍵酶。根據(jù)產(chǎn)物不同，異戊烯基轉(zhuǎn)移酶可分為香葉基二磷酸合酶（GGPS）、法尼基二磷酸合酶（FPS）、香葉基香葉基二磷酸合酶（GGPPS）［12］。思茅松（Pinus kesiya var.langbia nensis）中GGPS 基因表達(dá)量與產(chǎn)脂量呈顯著正相關(guān)［13］。劉美佳等基于珠子參（Panax japonicus）轉(zhuǎn)錄組結(jié)果，克隆出一條FPS 的cDNA 序列，并構(gòu)建珠子參過表達(dá)載體導(dǎo)入珠子參細(xì)胞中，發(fā)現(xiàn)轉(zhuǎn)基因植物的甾醇及三萜含量都有所升高［14］。楊銳等在丹參（Salvia miltiorrhiza Bge.）中過表達(dá)GGPPS 基因，發(fā)現(xiàn)丹參中丹參酮的含量顯著高于對照組［15］。本實驗中異戊烯基轉(zhuǎn)移酶在低脂樹中的表達(dá)量更高，此研究結(jié)果與Bai Q S 等人［16］的研究結(jié)果相似。當(dāng)樹體遭遇傷害后，高脂樹本身儲存的松脂大量流出可隔絕傷口與外界的接觸，而低脂樹樹體內(nèi)的松脂流動性較弱，樹體儲存的松脂較少，為保護(hù)自身免受外界的進(jìn)一步傷害而增加萜類物質(zhì)含量以封閉傷口形成物理屏障，因此異戊烯基轉(zhuǎn)移酶在低脂樹中的表達(dá)量更高。

萜類物質(zhì)是目前植物體內(nèi)種類最多的一類次級代謝產(chǎn)物，其種類多樣離不開各種修飾酶的作用，其中細(xì)胞色素P450 是修飾萜類物質(zhì)種類多樣的一類關(guān)鍵酶，97%的萜類物質(zhì)的多樣性是通過細(xì)胞色素P450（CYPs）氧化催化完成［17］。CYP450是一類家族基因，針葉樹種中，CYP450 多參與樹脂酸的修飾與合成［18～20］，本次測序差異基因中細(xì)胞色素P450 上調(diào)表達(dá)，表明高脂樹的CYP450 表達(dá)量更高。松脂由樹脂道上皮細(xì)胞合成并在樹脂道中儲存。松脂從上皮細(xì)胞轉(zhuǎn)移到樹脂道的過程中需ABC 轉(zhuǎn)運蛋白的參與。Westbrook 等人［21］發(fā)現(xiàn)ABC 轉(zhuǎn)運蛋白與松脂的跨位點流動有關(guān)，本研究中轉(zhuǎn)錄組與qRT-PCR 數(shù)據(jù)表明，參與轉(zhuǎn)運的ABC 轉(zhuǎn)運蛋白上調(diào)表達(dá)，與Liu［22］、陳曉明［23］等人的研究結(jié)果一致。這些研究結(jié)果表明濕加松松脂產(chǎn)量差異主要表現(xiàn)在萜類物質(zhì)的修飾與轉(zhuǎn)運過程。

對植物體而言，松脂是其重要的物理及化學(xué)防御手段。TIR-NBS-LRR 蛋白屬于R 基因的一種，主要參與植物體對病蟲的防御過程。李海霞等人［24］從毛白楊（Populus tomentosa）中分離出的一條TIR-NBS-LRR 抗病基因（PtDRG01）并將其導(dǎo)入到擬南芥（Arabidopsis thaliana）中，發(fā)現(xiàn)轉(zhuǎn)基因煙草株系的發(fā)病程度明顯低于非轉(zhuǎn)基因煙草，同時轉(zhuǎn)基因株系中PtDRG01的基因表達(dá)水平顯著高于非轉(zhuǎn)基因株系。張穎等［25］從野生刺葡萄中克隆出10 條NBS-LRR 類抗病基因片段，同時證實了NB6、NB7的表達(dá)受白腐病的調(diào)控。王猛［26］也提出TIR-NBS-LRR 的表達(dá)在馬尾松抗松材線蟲中發(fā)揮作用。這些結(jié)果說明TIR-NBS-LRR 蛋白參與了植物體抗病蟲的防御過程。木脂素是植物體重要的防御物質(zhì)，其合成受雙功能松脂醇—落葉松醇還原酶（PLR）和PR基因的調(diào)控。本試驗中雙功能松脂醇—落葉松醇還原酶在高脂樹中表達(dá)量顯著高于低脂樹中的含量。野生金蕎麥（Fagopyrum cymosum.）的研究中發(fā)現(xiàn)，PLR 可催化松脂醇—落葉松脂醇還原開環(huán)異構(gòu)落葉松脂并參與植物體的防御［27］。本研究中編碼TIR-NBS-LRR 蛋白的基因及PLR 基因在高產(chǎn)脂樹中的表達(dá)量顯著高于低產(chǎn)脂植株，說明高產(chǎn)脂植株對抵御病蟲傷害時表現(xiàn)出更強(qiáng)的抗性。研究表明輻照西洋參中隨著ATP 濃度會影響其次生產(chǎn)物皂苷的合成［28］。馬尾松中高產(chǎn)樹脂的ATP 含量要比低產(chǎn)樹脂高［23］。本研究發(fā)現(xiàn)GO 項主要富集在氧化磷酸過程及雙加氧酶過程中，氧化磷酸化可為植物體生命活動提供必要的ATP，ATP 是細(xì)胞進(jìn)行代謝活動最直接的能量來源，對代謝產(chǎn)物的生成有調(diào)控作用，高脂樹的代謝活動能力更強(qiáng)。

濕加松松脂合成途徑中，部分相關(guān)調(diào)控基因在高脂樹與低脂樹中的表達(dá)量不同。調(diào)控萜類物質(zhì)轉(zhuǎn)移的細(xì)胞色素P450 等基因、ABC 轉(zhuǎn)運體家族基因、參與防御系統(tǒng)的PLR 基因、TIR-NBS-LRR 基因上調(diào)表達(dá)直接或間接促進(jìn)了松脂總量的增加，異戊烯基轉(zhuǎn)移酶在低產(chǎn)脂植株中表達(dá)量更高。本試驗可以作為從轉(zhuǎn)錄層面探討濕加松產(chǎn)脂差異的形成機(jī)制，為后續(xù)相關(guān)分子標(biāo)記和功能基因克隆奠定重要基礎(chǔ)。