顏仁梁，梁永樞，周國洪，夏 黎，林 勵，周代營*

(1. 廣東食品藥品職業(yè)學(xué)院，廣州 510520； 2. 廣州中醫(yī)藥大學(xué)，廣州 510006)

化橘紅為蕓香科植物化州柚Citrusgrandis‘Tomentosa’或柚Citrusgrandis(L.) Osbeck的未成熟或近成熟的干燥外層果皮，前者習(xí)稱“毛橘紅”，后者習(xí)稱“光七爪”“光五爪”或“光橘紅”，具有燥濕化痰、理氣、消食的功效[1]. 有文獻(xiàn)報道將目標(biāo)起始密碼子多態(tài)性(SCoT)分子標(biāo)記[2]、簡單重復(fù)序列(SSR)分子標(biāo)記技術(shù)[3]、DNA條形碼[4]等基因分析技術(shù)用于不同品種化橘紅的親緣關(guān)系分析. 目前，已知柚皮苷、野漆樹苷、柚皮素、新橙皮苷、芹菜素等黃酮類成分是化橘紅的主要活性成分[5-9]，但毛橘紅與光橘紅的黃酮類成分存在很大差異[10-11]，由于其黃酮類成分的生物合成相關(guān)的基因序列還不清楚[12]，所以不同品種化橘紅的黃酮類成分差異形成的分子基礎(chǔ)有待深入研究.

隨著高通量測序技術(shù)的不斷發(fā)展，測序成本不斷降低，轉(zhuǎn)錄組研究已經(jīng)成為挖掘藥用植物功能基因的重要手段之一. 近年來，國內(nèi)外已經(jīng)開展了三七[13]、茴香[14]、決明子[15]、鐵皮石斛[16]、薯蕷[17]、馬藍(lán)[18]、當(dāng)歸[19]、澤瀉[20]多種藥用植物轉(zhuǎn)錄組的研究；在黃酮類成分的生物合成方面，有文獻(xiàn)報道采用二代測序技術(shù)研究了紫丁香[21]、三葉青[22]、甘草[23]、杭菊[24]、金銀花[25]、紅花[26]、雪蓮[27]等藥用植物體內(nèi)的黃酮類生物合成相關(guān)基因及其表達(dá)差異. 目前，柚及其同屬植物的轉(zhuǎn)錄組也有少量報道，例如：利用轉(zhuǎn)錄組研究低溫對柑橘長時間貯藏后果實品質(zhì)的影響[28]、柚子分子標(biāo)記的開發(fā)[29]、柚子果皮顏色突變相關(guān)差異表達(dá)基因[30]等. 本文采用二代高通量測序技術(shù)，對化州柚和柚進(jìn)行轉(zhuǎn)錄組測序、拼接組裝，建立化州柚和柚葉的轉(zhuǎn)錄組，對黃酮類成分生物合成相關(guān)的差異表達(dá)基因進(jìn)行了功能注釋與比較分析，為闡明化州柚和柚的黃酮類有效成分含量差異的分子基礎(chǔ)研究奠定基礎(chǔ).

1 研究方法

1.1 材料與試劑

化州柚和柚的嫩葉樣本采自廣東食品藥品職業(yè)學(xué)院藥圃. Omega Plant RNA kit 試劑盒(美國Omega)，NEB#7530 試劑盒(新英格蘭Biolabs)，DNA 1000 assay Kit (美國安捷倫科技有限公司).

1.2 儀器

PCR儀(東勝興業(yè)科學(xué)儀器有限公司)，高速離心機(jī)(德國艾本德股份公司)，Nanodrop 2000 (美國賽默飛世爾科技公司)，Illumina HiSeq 2500 (美國Illumina公司).

1.3 樣本采集、總RNA提取及轉(zhuǎn)錄組測定

采集化州柚和柚的嫩葉適量，在液氮中速凍15 min后，保存于-85 ℃超低溫冰箱備用. 按照Qiagen RNeasy Mini Kit提取試劑盒說明書的方法提取各樣本總RNA，用Nanodrop 2000 與Agilent 2100 Bioanalyzer檢測所提取總RNA的濃度與純度，A260/A280比值在1.8～2.0，RNA 分子完整數(shù)(RIN，RNA Integrity Number)大于8，則符合轉(zhuǎn)錄組測序樣本要求，轉(zhuǎn)錄組測定由廣州基迪奧生物科技有限公司完成.

1.4 轉(zhuǎn)錄組數(shù)據(jù)組裝

測序產(chǎn)生的原始數(shù)據(jù)(Raw Reads)，經(jīng)過初步過濾后得到待分析數(shù)據(jù)(Clean Reads)，去除有帶接頭的Reads、N的比例大于10%的Reads以及質(zhì)量很低的 Reads 后得到高質(zhì)量的Clean Reads. 使用短序列組裝軟件 Trinity將高質(zhì)量待分析數(shù)據(jù)從頭組裝得到單基因(Unigene)，用N50數(shù)值來評估組裝結(jié)果質(zhì)量.N50越長，數(shù)量越少，說明組裝質(zhì)量越好. Unigene 表達(dá)量的計算使用 RPKM 法(Reads Per kb Per Million Reads)，其計算公式為：RPKM=(1 000 000×C)/(N×L/1 000)，假設(shè)RPKM為某Unigene 的表達(dá)量，則C為比對到某Unigene的Reads數(shù)，N為比對到所有Unigene的總Reads數(shù)，L為某Unigene的堿基數(shù).

1.5 Unigene基本功能注釋

首先，通過比對軟件Blastx將Unigene序列比對到美國國家生物信息中心(NCBI)非冗余蛋白序列數(shù)據(jù)庫(Nr)、由歐洲生物信息學(xué)研究所(EBI)維護(hù)的經(jīng)過注釋的蛋白質(zhì)序列數(shù)據(jù)庫(Swiss-Prot)、京都基因與基因組百科全書數(shù)據(jù)庫(KEGG)和蛋白相鄰類的聚簇數(shù)據(jù)庫(COG/KOG)，得到與給定 Unigene 具有最高序列相似性的蛋白，從而得到該 Unigene 的蛋白功能注釋信息.

1.6 差異表達(dá)基因分析

先用edgeR軟件包進(jìn)行化州柚與柚的轉(zhuǎn)錄組差異統(tǒng)計分析，利用錯誤發(fā)現(xiàn)率(FDR)與log 2FC(FC為兩樣品間表達(dá)量的比值)來篩選差異基因，篩選條件為 FDR<0.05 且 |log 2FC|>1，得到差異表達(dá)基因之后，對差異表達(dá)基因再進(jìn)一步做基因本體(GO)功能分析和KEGG Pathway分析. GO功能分析給出差異表達(dá)基因的GO功能分類注釋和GO功能顯著性富集分析. 通過Pathway顯著性富集可以確定化州柚和柚轉(zhuǎn)錄組的差異表達(dá)基因參與的最主要生化代謝途徑和信號轉(zhuǎn)導(dǎo)途徑.

2 結(jié)果

2.1 測序與組裝結(jié)果

本次測定共獲得86 298 757 bp 堿基，組裝得到116 202個單基因(Unigene)，N50為1 543，GC%為42.077 5，最長基因長度20 219 nt，最短基因長度201 nt，基因平均長度742 nt. Unigene 長度分布見圖1.

圖1 化州柚和柚轉(zhuǎn)錄組Unigene基因的長度分布圖

2.2 Unigene基本注釋

運(yùn)用BLAST比對軟件將組裝得到的116 202個Unigene與已知的4個公共數(shù)據(jù)庫Nr、Swiss-Prot、KOG、KEGG 進(jìn)行比對，其中Nr注釋了53 738個基因，Swiss-Prot注釋了52 965個基因，KOG注釋了40 293個基因，KEGG注釋了24 891個基因，總共獲得同源比對信息的Unigene有68 923條，占 59.31%，沒有獲得同源比對信息的基因有47 279個，占40.69%(圖2).

圖2 化州柚與柚轉(zhuǎn)錄組在4個數(shù)據(jù)庫基因注釋的維恩圖

KOG注釋結(jié)果顯示：化州柚和柚的40 293個Unigene分布于25個KOG分類中，其中聚類到一般性功能預(yù)測、信號轉(zhuǎn)導(dǎo)機(jī)制、翻譯后的修飾、蛋白轉(zhuǎn)換、分子伴侶等功能基因分布較多(圖3).

化州柚與柚轉(zhuǎn)錄組GO功能注釋的具體亞類與基因數(shù)量見圖4. 圖4顯示：分子功能(Molecular Function)、細(xì)胞組分(Cellular Component)、生物過程(Biological Process)又可細(xì)分為 50個功能亞類，其中生物學(xué)過程分為23個亞類，細(xì)胞組分分為16個亞類，分子功能分為11個亞類. 在生物學(xué)過程中，代謝過程(Metabolic Process)、細(xì)胞過程(Cellular Process)和單有機(jī)體過程(Single-Organism Process)所占基因數(shù)較多；在細(xì)胞組分中，涉及細(xì)胞(Cell)、細(xì)胞組分(Cell Part)和細(xì)胞器(Organelle)的基因較多；在分子功能中，涉及催化活性(Catalytic Activity)和結(jié)合(Binding)基因較多.

圖3 化州柚與柚轉(zhuǎn)錄組Unigene基因的KOG功能分類圖

圖4 化州柚與柚轉(zhuǎn)錄組Unigene基因的GO功能分類圖

KEGG 數(shù)據(jù)庫注釋結(jié)果(圖5)顯示：具有通路注釋的基因共14 833個，涉及 136個代謝通路. 其中代謝途徑相關(guān)基因6 430個(43.35%);次級代謝產(chǎn)物的生物合成相關(guān)基因3 529個(23.79%);核糖體相關(guān)基因1 665個(11.22%)，碳代謝相關(guān)基因1 423個(9.59%)，氨基酸的生物合成相關(guān)基因1 032個(6.96%)，植物-病原體相互作用相關(guān)基因743個(5.01%)，內(nèi)質(zhì)網(wǎng)中的蛋白質(zhì)加工相關(guān)基因684個(4.61%)，氧化磷酸化相關(guān)基因671個(4.52%)，糖酵解/糖異生相關(guān)基因661個(4.46%)，嘌呤代謝相關(guān)基因632個(4.26%).

圖5 化州柚與柚轉(zhuǎn)錄組Unigene基因的KEGG通路注釋

2.3 基因差異性表達(dá)分析

本研究共篩選得到化州柚與柚葉轉(zhuǎn)錄組的差異表達(dá)基因共 6 419個，其中化州柚相對于柚上調(diào)基因有3 799個，占差異表達(dá)基因的59.18%，下調(diào)基因2 620個，占40.82%. 獲得 GO 注釋的差異表達(dá)基因有4 826個，1 346 條GO注釋中，937個注釋為生物學(xué)過程相關(guān)，111個注釋為細(xì)胞組分相關(guān)，298個注釋為分子功能相關(guān).

與KEGG數(shù)據(jù)庫進(jìn)行比對，共有771個基因獲得注釋，涉及125條代謝通路，其中顯著差異表達(dá)富集通路的基因有：植物-病原體相互作用相關(guān)基因134個(17.38%);MAPK信號傳導(dǎo)途徑-植物相關(guān)基因59個(7.65%)；植物激素信號轉(zhuǎn)導(dǎo)40個(5.19%)；苯丙烷生物合成相關(guān)基因37個(4.80%)、倍半萜類化合物和三萜類化合物生物合成相關(guān)基因10個(1.30%)；胡蘿卜素生物合成相關(guān)基因17個(2.20%)；油菜素類固醇生物合成相關(guān)基因7個(0.91%)；泛醌和其他萜類化合物-醌生物合成相關(guān)基因相關(guān)基因14個(1.82%)；類黃酮生物合成相關(guān)基因11個(1.43%)；氨基糖和核苷酸糖代謝相關(guān)基因36(4.67%)；次級代謝產(chǎn)物的生物合成相關(guān)基因223個(28.92%)；芪類化合物、二苯基庚烷類和姜辣素生物合成相關(guān)基因8個(1.04%)；黃酮和黃酮醇生物合成相關(guān)基因3個(0.39%);谷胱甘肽代謝相關(guān)基因29個(3.76%)等.

3個與黃酮和黃酮醇生物合成相關(guān)基因注釋為黃烷酮7-O-葡萄糖苷2″-O-β-L-鼠李糖基轉(zhuǎn)移酶(C12RT1)、黃酮醇-3-O-葡萄糖苷L-鼠李糖基轉(zhuǎn)移酶(RT)、黃酮醇-3-O-葡糖苷/半乳糖苷葡糖基轉(zhuǎn)移酶(3GGT)；11個與類黃酮生物合成相關(guān)的差異表達(dá)基因注釋為反式肉桂酸4-單加氧酶(CYP73A12)、咖啡酰-CoAO-甲基轉(zhuǎn)移酶(TSM1)、查爾酮合成酶(CHS1)、黃酮醇合成酶(FLS)、莽草酸O-羥基肉桂酰轉(zhuǎn)移酶(HCT)(5個)、黃烷酮7-O-葡萄糖苷2″-O-β-L-鼠李糖基轉(zhuǎn)移酶(C12RT1)、白細(xì)胞花青素還原酶(LAR)，這14個黃酮類生物合成相關(guān)差異表達(dá)基因表達(dá)量rpkm值見表1.

表1 化州柚與柚轉(zhuǎn)錄組中與黃酮類成分生物合成相關(guān)的差異表達(dá)基因Table 1 The differentially expressed genes related to flavonoids biosynthesis in the transcriptome of CGT and CGO

3 討論

本次研究發(fā)現(xiàn)化州柚和柚的轉(zhuǎn)錄組有6 419個差異性表達(dá)基因，Pathway富集分析結(jié)果顯示：與次級代謝產(chǎn)物的生物合成相關(guān)基因最多(占28.92%)，這為化州柚和柚的化學(xué)成分及含量差異分子基礎(chǔ)研究奠定了基礎(chǔ)；與植物-病原體相互作用相關(guān)基因也較多(占17.38%)，為闡明二者植物生態(tài)適應(yīng)性差異奠定了基礎(chǔ). 另外，本研究結(jié)果顯示：與黃酮和黃酮醇生物合成相關(guān)差異表達(dá)基因僅3個，與類黃酮生物合成相關(guān)的差異表達(dá)基因有11個，說明化州柚和柚在黃酮和黃酮醇生物合成方面很類似，但在類黃酮生物合成方面差異較大. 由于本次研究采集的是化州柚和柚的嫩葉樣本，而化橘紅藥用部位是未成熟或近成熟的干燥外層果皮，所以本次研究為不同品種來源化橘紅黃酮類成分差異奠定了基因數(shù)據(jù)基礎(chǔ)，其生物合成的相關(guān)合成酶差異表達(dá)還有待于下一步深入研究. 尚有47 279個(約40%)Unigene未獲得功能注釋，與二代測序技術(shù)所獲得的基因序列較短，不易獲得同源基因信息有關(guān)；本次研究是無參考基因組的轉(zhuǎn)錄組測序，由于目前還未見化州柚和柚的全基因組信息，只能參考相關(guān)物種的基因組數(shù)據(jù)進(jìn)行分析，故所得注釋信息有所欠缺，這也是目前藥用植物轉(zhuǎn)錄組數(shù)據(jù)分析存在的共性問題[31-35].

4 結(jié)論

本研究首次獲得了化州柚和柚的轉(zhuǎn)錄組數(shù)據(jù)，共得到116 202個Unigene，平均長度為742 nt，組裝效果較好. 通過與Nr、Swiss-Port、KOG、KEGG等數(shù)據(jù)庫比對獲得的生物信息學(xué)注釋結(jié)果，揭示了化州柚和柚轉(zhuǎn)錄組主要的整體表達(dá)特征(約60%).