利用轉(zhuǎn)錄組數(shù)據(jù)挖掘東紫蘇單萜生物合成相關基因

耿秀文，張愛麗，唐仁華，普春霞

云南中醫(yī)藥大學 中藥材優(yōu)良種苗繁育工程中心，云南 昆明 650500

東紫蘇Elsholtzia bodinieriVaniot 為唇形科荊芥亞科（Nepetoideae）香薷族（Elsholtzieae）香薷屬ElsholtziaL.的植物，又名鳳尾茶、牙刷草、小山茶、野山茶等，主產(chǎn)于云南滇中、滇南地區(qū)，生于海拔1200～3000 m 的山坡草地、稀疏松林中或石山上。東紫蘇全草入藥，氣清香，味辛、微苦，用于感冒、咽喉痛、頭痛、扁桃腺炎、小兒口腔炎、牙痛、眼結膜炎、肝炎、消化不良[1-3]。東紫蘇含有揮發(fā)油、黃酮類、三萜類、酚類等成分[4]?，F(xiàn)代藥理研究表明，東紫蘇具有抑菌[5-6]、抗病毒[6-7]、抗氧化[8-9]、調(diào)血脂[10]的作用。

香薷屬植物的揮發(fā)油被認為具有抗菌、抗氧化等作用[11]。前期課題組對東紫蘇15 個居群葉的揮發(fā)油進行GC-MS 分析[12]，初步確定1,8-桉葉油醇（1, 8-cineole）、α-乙酸松油酯（α-terpinyl acetate）、β-蒎烯（β-pinene）、α-蒎烯（α-pinene）等為東紫蘇葉揮發(fā)油的主成分，其中單萜為 72.469%～91.855%，倍半萜為3.727%～13.121%。付立卓等[13]、程偉賢等[14]的研究結果在主成分上與本課題組的研究結果基本一致[15]。根據(jù)文獻報道，1, 8-桉葉油醇有抗炎和鎮(zhèn)痛作用[16]，對氣道炎癥具有抑制作用[17]，能夠減少鼻竇炎的惡化，預防慢性阻塞性肺病加重和改善哮喘[18-20]，并且還有解痙和抗菌活性[21]，以及一定的止瀉作用[22]。

單萜化合物是由2 個異戊二烯結構單元（C10）組成的鏈狀或環(huán)狀化合物，是一類結構多樣的天然化合物，廣泛存在于植物揮發(fā)油和樹脂中[23]。植物單萜（monoterpenoid）由質(zhì)體內(nèi)的4-磷酸-2-甲基赤蘚糖（2-C-methyl-D-erythritol 4-phosphate，MEP）途徑合成，香葉基二磷酸（geranyl diphosphate，GPP）是所有單萜共同的前體物質(zhì)，GPP 在單萜合酶（monoterpene synthases，mono-TPS）的催化作用下生成單萜，單萜合酶是單萜生物合成的關鍵酶[24]。東紫蘇葉揮發(fā)油的主成分1,8-桉葉油醇應是由GPP在桉葉油醇合酶（cineole synthase）的作用下形成[25]。

轉(zhuǎn)錄組測序數(shù)據(jù)是挖掘基因以及參與各種代謝途徑的酶的有效方法[26]，現(xiàn)已廣泛應用于藥用植物和芳香植物中次生代謝產(chǎn)物生物合成的關鍵基因的研究和鑒定[27]。本研究通過對東紫蘇進行轉(zhuǎn)錄組測序，對測序數(shù)據(jù)進行組裝、拼接及注釋，并初步挖掘了東紫蘇單萜類化合物代謝途徑的相關基因，為東紫蘇揮發(fā)油單萜合酶基因的克隆、次生代謝生物合成途徑解析以及功能驗證的研究提供了基因資源和理論基礎。

1 材料

樣本于2015年11月3日采自云南省昆明市呈貢區(qū)松茂水庫（經(jīng)度102°55′18″，緯度24°53′25″，海拔2114 m），編號SM；于2015年11月7日采自云南省昆明市嵩明縣阿子營村（經(jīng)度102°47′43″，緯度25°21′25″，海拔2086 m），編號AZY；原植物由云南中醫(yī)藥大學普春霞副教授鑒定為東紫蘇E.bodinieriVaniot。樣本材料采集后，移栽至云南中醫(yī)藥大學呈貢校區(qū)，于2017年11月22日采樣后送轉(zhuǎn)錄組測序。

2 方法

2.1 東紫蘇揮發(fā)油的提取及GC-MS 分析

采用水蒸氣蒸餾法提取SM 和AZY 2 個產(chǎn)地東紫蘇葉片中的揮發(fā)油，利用GC-MS 法分析檢測其主要成分及含量。

2.1.1 氣相條件 色譜柱：Agilent 19091J-115，HP-5（5% Phenyl Methyl Siloxan，50 m×0.32 mm×0.52 μm）；進樣口溫度250 ℃，檢測器（氫火焰離子檢測器，F(xiàn)ID）溫度250 ℃；柱箱升溫程序：50 ℃升到250 ℃，每分鐘升10 ℃；進樣量1 μL；分流比20∶1；進樣口壓力102.28 kPa；體積流量2.1 mL/min；載氣為高純氦氣。

2.1.2 質(zhì)譜條件 電離方式：EI，電子轟擊能量：70 eV，掃描范圍m/z50～550；離子源溫度230 ℃；四極桿溫度150 ℃；定量方法：相對峰面積歸一化法；譜庫：NIST14。

2.2 轉(zhuǎn)錄組測序

采集東紫蘇2 個產(chǎn)地葉的樣品材料，委托北京諾禾致源科技股份有限公司利用Illumina Hiseq 2000 進行高通量轉(zhuǎn)錄組測序。

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

得到原始測序數(shù)據(jù)后，為保證信息分析質(zhì)量，去除帶接頭的、N（N 表示無法確定堿基信息）的比例大于10%的以及低質(zhì)量reads，對raw reads 過濾，得到clean reads，然后采用Trinity 對clean reads進行拼接組裝，最后得到unigenes。

2.4 轉(zhuǎn)錄組功能注釋

將組裝得到的unigenes 序列比對到NCBI 官方的蛋白序列數(shù)據(jù)庫（RefSeq non-redundant proteins，NR）、NCBI 官方的核酸序列數(shù)據(jù)庫（nucleotide sequence database，NT）、京都基因與基因組百科全書（Kyoto encyclopedia of genes and genomes，KEGG）、Swiss-Prot 蛋白質(zhì)數(shù)據(jù)庫（A manually annotated and reviewed protein sequence database，Swiss-prot）、蛋白質(zhì)家族數(shù)據(jù)庫（Protein family，PFAM）、基因本體論數(shù)據(jù)庫（gene ontology，GO）、真核生物蛋白質(zhì)同源簇數(shù)據(jù)庫（clusters of orthologous groups for eukaryotic complete genomes，KOG）這7 個數(shù)據(jù)庫進行基因功能注釋。

2.5 東紫蘇揮發(fā)油單萜合成相關候選基因篩選及克隆

從轉(zhuǎn)錄組數(shù)據(jù)中篩選出東紫蘇揮發(fā)油單萜合成相關候選基因，獲得基因序列和開放閱讀框（ORF），篩選ORF 長度大于700 bp 的基因，分別設計ORF 擴增引物。全部引物由生工生物工程（上海）股份有限公司合成（表1）。

表1 基因與引物序列Table 1 Gene and primer sequences

2.6 東紫蘇單萜合成相關候選基因RT-PCR 驗證

取東紫蘇葉片總RNA 4 μL，按照TaKaRa 公司的PrimeScriptTMII 1st Strand cDNA Synthesis Kit 試劑盒說明書進行反轉(zhuǎn)錄反應，合成cDNA。以cDNA為模板，反應體系（50 μL）：模板2 μL，正反向引物各1 μL，TaKaRa Taq DNA 聚合酶25 μL，滅菌雙蒸水21 μL。PCR 反應程序為98 ℃、30 s；98 ℃、 10 s，58 ℃、20 s，72 ℃、2 min，30 個循環(huán)；72 ℃、8 min，12 ℃保溫。PCR 產(chǎn)物經(jīng)電泳分析后，送生工生物工程（上海）股份有限公司測序。

3 結果與分析

3.1 東紫蘇揮發(fā)油GC-MS 成分測定

按設定的GC-MS 條件分析東紫蘇揮發(fā)油的化學成分，經(jīng)NIST14 譜庫檢索、質(zhì)譜分析等確定揮發(fā)油的化學成分，并用面積歸一化法計算各組分的相對含量，其中1,8-桉葉油醇和α-乙酸松油酯均為含量較高的成分。從圖1 可以看出，產(chǎn)自松茂的東紫蘇葉揮發(fā)油含量最高的成分為α-乙酸松油酯，其次是β-蒎烯；而產(chǎn)自阿子營的東紫蘇葉揮發(fā)油含量最高的成分為1,8-桉葉油醇，其次是β-蒎烯。

圖1 SM (A) 和AZY (B) 產(chǎn)東紫蘇揮發(fā)油GC-MSFig.1 GC-MS charts of essential oil from E.bodinieri from SM (A) and AZY (B)

3.2 東紫蘇轉(zhuǎn)錄組數(shù)據(jù)組裝與質(zhì)量分析

東紫蘇轉(zhuǎn)錄組測序總計得到93 531 502 條raw reads，過濾產(chǎn)生了91 169 356 條高質(zhì)量clean reads，共計13.67 Gb 的有效數(shù)據(jù)，Q20（Phred 數(shù)值大于20的堿基占總體堿基的百分比）SM 為94.67%，AZY為94.62%；Q30（Phred 數(shù)值大于30 的堿基占總體堿基的百分比）SM 為86.81%，AZY 為86.76%；GC 含量（clean reads 中G 與C 占4 種堿基的百分比）SM 為47.99%，AZY 為47.17%。測序后得到的東紫蘇轉(zhuǎn)錄組數(shù)據(jù)產(chǎn)出質(zhì)量見表2。運用Trinity 軟件將質(zhì)控后得到的高質(zhì)量序列進行組裝，得到153 657 個transcripts，獲得93 327 個unigenes，其中unigenes 的總長度為163 878 269 bp，平均長度為1756 bp，N50為2333 bp。unigenes 長度分布圖（圖2）顯示，63 142條unigenes 長度超過1000 bp，31 407 條unigenes 長度超過2000 bp。以上結果表明本研究轉(zhuǎn)錄組測序及組裝完整性較高，能夠進行后續(xù)的注釋分析。轉(zhuǎn)錄組數(shù)據(jù)已上傳至 NCBI 公共數(shù)據(jù)庫，登錄號為SRR14567298、SRR14567299。

3.3 東紫蘇轉(zhuǎn)錄組unigenes 的功能注釋

使用BLAST 將拼裝所得的unigenes 分別與NR、NT、KEGG、Swiss-prot、PFAM、GO、KOG這7 個數(shù)據(jù)庫進行比對分析，對各數(shù)據(jù)庫注釋的unigenes 數(shù)目進行統(tǒng)計，進而獲得東紫蘇轉(zhuǎn)錄組unigenes 的功能注釋信息，比對結果見表3。結果表明，80 420 條unigenes（86.17%）在NR 數(shù)據(jù)庫比對成功得到注釋，在KEGG、Swiss-port、GO、KOG 等數(shù)據(jù)庫獲得注釋的unigenes 數(shù)目依次為34 354 條（36.81%）、64 873 條（69.51%）、61 194條（65.56%）、25 391 條（27.2%）。15 115 條unigenes同時在所有數(shù)據(jù)庫中注釋，至少在1 個數(shù)據(jù)庫注釋成功的unigenes 有83 704 條（89.68%）。

表2 東紫蘇轉(zhuǎn)錄組數(shù)據(jù)產(chǎn)出質(zhì)量Table 2 Quality of output of E.bodinieri transcriptome data

圖2 東紫蘇轉(zhuǎn)錄組unigenes 長度分布圖Fig.2 Length distribution of unigenes in transcriptome of E.bodinieri

表3 東紫蘇轉(zhuǎn)錄組unigenes 在各數(shù)據(jù)庫注釋情況Table 3 E.bodinieri transcriptome unigenes annotated in each database

通過與NR 數(shù)據(jù)庫進行比對注釋，可以獲取本物種基因序列與近緣物種基因序列的相似性，以及 本物種基因的功能信息。從圖3 物種分布圖可以看出與其他物種序列的相似性，在比對到NR 數(shù)據(jù)庫中的unigenes 中，有65.2%和芝麻Sesamum indicumL.相匹配，其次是黃色猴面花Erythranthe guttata(DC.) G.L.Nesom.，匹配度為13.5%，丹參Salvia miltiorrhizaBunge、中?？Х菴offea canephoraPierre ex Froehn.、葡萄Vitis viniferaL.的匹配度分別為1.6%、1.3%、1.2%，其他物種為17.3%。

圖3 東紫蘇轉(zhuǎn)錄組unigenes 的物種分布圖Fig.3 Species classification of unigenes in transcriptome of E.bodinieri

對基因進行GO 注釋之后，共有61 194 條（65.56%）unigenes 注釋到GO 數(shù)據(jù)庫，并根據(jù)GO 功能分別將它們注釋到生物過程（biological process）、細胞組分（cellular component）和分子功能（molecular function）3 大類中，共分為56個小組，注釋結果見圖4。其中生物過程主要聚集在細胞過程（cellular process）、代謝過程（metabolic process），涉及的基因分別有36 783、34 872 條；其次是單組織過程（single-organism process）和生物調(diào)節(jié)（biological regulation），基因數(shù)量分別為27 213、12 627 條。細胞組分中細胞（cell）和細胞部分（cell part）相關基因數(shù)量較多，有19 561 和19 559 條，其次是細胞器（organelle），有13 142 條。分子功能中具有結合功能（binding）和催化活性（catalytic activity）的基因數(shù)量較多，分別為37 797 和30 440 條，其他類別的基因數(shù)目普遍較少。

圖4 東紫蘇轉(zhuǎn)錄組unigenes 的GO 功能分類注釋Fig.4 GO functional classification annotation of unigenes of E.bodinieri transcriptome

為了進一步分析東紫蘇unigenes 的功能，進行了KOG 功能分類，共有25 391 條（27.2%）unigenes注釋到KOG 數(shù)據(jù)庫，得到25 個不同的KOG 功能類群（圖5）。其中翻譯后修飾、蛋白反轉(zhuǎn)和伴侶（ posttranslational modification, protein turnover, chaperones）的數(shù)量最多，有3509 條；其次是一般功能預測（general function prediction only）的基因，數(shù)量為3458 條；而數(shù)量最少的為細胞運動性（cell motility），僅有23 條。 將東紫蘇unigenes 與KEGG 數(shù)據(jù)庫比對，有34 354條（36.81%）unigenes 注釋到該數(shù)據(jù)庫中，分布于129 條代謝通路中。KEGG 代謝通路（圖6）分為5大分支：細胞過程（cellular processes，A）1545 條、環(huán)境信息處理（ environmental information processing，B）1299 條、遺傳信息處理（genetic information processing ， C ） 7448 條、 代 謝（metabolism，D）13 643 條和有機系統(tǒng)E（organismal systems）1322 條。其中注釋數(shù)量最多的是碳水化合物代謝（carbohydrate metabolism），為2983 條，其次是折疊、分類和降解（folding，sorting and degardation）為2687 條，數(shù)量最少的是膜運輸（membrane transport），為135 條。

圖5 東紫蘇轉(zhuǎn)錄組unigenes 的KOG 功能分類注釋Fig.5 KOG functional classification annotation of unigenes of E.bodinieri transcriptome

圖6 東紫蘇轉(zhuǎn)錄組unigenes 的KEGG 分類注釋Fig.6 KEGG classification annotation of unigenes of E.bodinieri transcriptome

3.4 東紫蘇單萜代謝途徑與相關基因

通過KEGG 代謝通路分析，發(fā)現(xiàn)與東紫蘇揮發(fā)油萜類生物合成相關的代謝通路有4 條，分別是：萜類骨架生物合成，編號為ko00900，相關基因有218 條；單萜類生物合成，編號為ko00902，相關基因有30 條；檸檬烯和蒎烯降解，編號為ko00903，相關基因有38 條；倍半萜和三萜生物合成，編號為ko00909，相關基因有95 條。其中與東紫蘇單萜合成相關的代謝通路是萜類骨架生物合成和單萜類生物合成。

與東紫蘇單萜合成相關的酶分別有香葉基二磷酸合酶（geranyl diphosphate synthase，GPS），相關 unigenes 有11條；新薄荷醇脫氫酶［(+)-neomenthol dehydrogenase，NDS］，相關unigenes 數(shù)量最多，有 22 條；α-松油醇合酶［(-)-alpha-terpineol synthase，ATS］，芳樟醇合酶［(3S)-linalool synthase，LIS］，月桂烯/羅勒烯合酶（myrcene/ocimene synthase，MS/OS），這3 個酶相關unigenes 數(shù)量較少，共有8 條，具體見表4。另外，通過比較這些unigenes 的FPKM 值可以看出，NDS 基因平均表達量是最高的，其中SM 為8.97，AZY 為17.06；其次是ATS 基因，平均表達量SM為7.16，AZY 為3.92；GPS 基因平均表達量SM 為2.29，AZY 為1.64；LIS 基因和MS/OS 基因平均表達量均偏低。

篩選得到的東紫蘇11 個單萜合成相關候選基因在SM 和AZY 2 個產(chǎn)地的FPKM 值（圖7）。從圖7中可以看到，基因2（Cluster-10292.10951）注釋為NDS，在2 個產(chǎn)地中的表達量均為最高，其中SM 為25.66，AZY 為29.56；而基因1（Cluster-11595.0）注釋為ATS，在2 個產(chǎn)地中的表達量均較低，其中SM為0.18，AZY 則為0；基因9（Cluster-10292.51149）注釋為NDS，基因10（Cluster-10292.60790）和基因11（Cluster-10292.60789）均注釋為GPS，這3個基因在2 個產(chǎn)地中的表達量基本持平。

表4 東紫蘇揮發(fā)油單萜代謝途徑的酶與相關基因Table 4 Enzymes and related genes of monoterpene metabolic pathway of essential oil from E.bodinieri

圖7 東紫蘇單萜合成相關候選基因的表達量Fig.7 Expression of candidate genes related to monoterpene synthesis of E.bodinieri

3.5 東紫蘇揮發(fā)油單萜合成相關候選基因全長unigenes 的RT-PCR 驗證

利用表1 中的ORF 擴增引物進行PCR 擴增，并進一步對PCR 產(chǎn)物進行測序分析，得到以下結果：基因 2（Cluster-10292.10951）、基因 4（ Cluster-10292.60008 ）、 基因 6 （ Cluster- 10292.56890）、基因7（Cluster-5267.0）、基因10（ Cluster-10292.60790 ）、 基因11 （ Cluster- 10292.60789）這6 個基因擴增后，經(jīng)電泳分析，均產(chǎn)生明亮的目標條帶（圖8），經(jīng)測序分析后表明這6 個基因與原unigenes ORF 一致，這些候選unigenes為全長基因；另外基因1（Cluster-11595.0）目前擴增出后半段，基因 5（Cluster-10292.9722）、8（Cluster-1426.0）這2 個基因擴增出中間部分，下 一步將繼續(xù)嘗試全長擴增；基因 3（Cluster- 10292.50707）、基因9（Cluster-10292.51149）這2個基因擴增失敗。

圖8 東紫蘇揮發(fā)油6 個單萜合成相關候選基因RT-PCR 凝膠電泳圖Fig.8 RT-PCR gel electrophoresis of six candidate genes related to monoterpene synthesis in essential oil from E.bodinieri

4 討論

東紫蘇作為云南特產(chǎn)物種，在云南主要用作民間藥物及保健飲品[28]，其作為藥用植物的基因資源部分一直是空白。本研究采用Illumina Hiseq 2000高通量測序平臺，首次對東紫蘇進行轉(zhuǎn)錄組測序分 析，填補了東紫蘇轉(zhuǎn)錄組數(shù)據(jù)信息的空白。通過測序共獲得13.67 Gb 數(shù)據(jù)，91 169 356 條高質(zhì)量序列（clean reads），Trinity 組裝獲得93 327 條unigenes，平均長度為1756 bp。

結合生物信息學分析方法對東紫蘇轉(zhuǎn)錄組數(shù)據(jù)進行序列相似性、基因注釋和功能分類，將所有unigenes 分別與NR、NT、KEGG、Swiss-prot、PFAM，GO、KOG 等數(shù)據(jù)庫進行比對，經(jīng)過BLAST 比對分析后發(fā)現(xiàn)，有80 420 條unigenes 與其他近緣物種的已知基因具有相似性，與芝麻相似度最高，為65.2%。在KEGG 數(shù)據(jù)庫比對中，有34 354 條（36.81%）unigenes 注釋到該數(shù)據(jù)庫中，分布于129條代謝通路中，并發(fā)現(xiàn)4 條萜類化合物代謝相關的途徑。

RNA-Seq 技術的發(fā)展為轉(zhuǎn)錄組研究提供了機會，已成為發(fā)現(xiàn)與藥用植物各種次生代謝途徑生物合成相關基因的有力工具，即使是在沒有參考基因組的物種中也是如此[29-30]。本研究是通過已注釋基因的挖掘，在無參轉(zhuǎn)錄組中篩選目的基因，因與芝麻同源性匹配度最高，因此以芝麻作為參照來尋找目的基因。

前期課題組對東紫蘇15 個居群葉的揮發(fā)油進行GC-MS 分析[12]，并根據(jù)主成分的差異，將其劃分為不同的化學型，其中SM 屬于T 型（1,8-桉葉油醇＜5%，5%＜α-乙酸松油酯＜20%），阿子營（AZY）屬于C 型（1,8-桉葉油醇>50%，α-乙酸松油酯＜5%），與之呼應的是2 個個體在單萜合成相關候選基因表達量上的明顯差異。已有研究表明與植物揮發(fā)油中高含量成分對應的單萜合酶基因其克隆成功率更高，主要原因可能在于基因表達量與其產(chǎn)物的積累量是一致的[31-34]，由于2 個產(chǎn)地的東紫蘇葉揮發(fā)油主成分差異明顯，具有代表性，因此分析這2 個產(chǎn)地的東紫蘇葉在轉(zhuǎn)錄組數(shù)據(jù)中相關基因的表達量，有可能提高單萜合酶基因克隆的成功率?；诖?，本研究在明確東紫蘇揮發(fā)油主要成分的基礎上，篩選出與單萜合成相關的代謝途徑2 條，單萜合成相關候選基因41 個，通過RT-PCR 驗證其中的6 個unigenes為全長基因，這為下一步研究東紫蘇單萜生物合成及基因功能驗證提供數(shù)據(jù)支持。

本研究篩選出的單萜合酶基因克隆成功率低，并且轉(zhuǎn)錄組注釋的功能信息與東紫蘇揮發(fā)油GC-MS 數(shù)據(jù)分析結果不吻合，即注釋的功能與揮發(fā)油的成分并不匹配，推測是受測序技術所限，以及數(shù)據(jù)庫可供比對的基因資源有限，增加了功能注釋的難度，致使一些東紫蘇揮發(fā)油單萜合成相關的關鍵酶基因未能得到精準拼接、組裝和注釋。此外，由于相同物種來源的單萜合酶的相似性高于不同物種來源而具有相同功能的單萜合酶，序列相關性并不直接轉(zhuǎn)化為產(chǎn)物相關性[23,33,35-38]，因此基于序列相似性僅可以判斷其是否為單萜合酶，卻不能預測其產(chǎn)物，所以已經(jīng)克隆成功的6 個單萜合成相關候選基因需要進行功能驗證，確認其在東紫蘇揮發(fā)油單萜生物合成過程中的催化功能，才能確定該酶屬于哪一種單萜合酶，以便最終得到東紫蘇1,8-桉葉油醇合酶基因。

本研究首次采用高通量測序技術建立了東紫蘇轉(zhuǎn)錄組數(shù)據(jù)庫，獲得了大量的轉(zhuǎn)錄本信息，并挖掘出單萜生物合成途徑的關鍵酶基因，填補了東紫蘇揮發(fā)油單萜合成途徑分子研究的空白，為東紫蘇揮發(fā)油單萜合酶基因的克隆、功能驗證以及品種選育提供了基因資源及理論指導基礎。

利益沖突所有作者均聲明不存在利益沖突