基于高通量測序的玄參根部轉錄組學研究及萜類化合物合成相關基因的挖掘

潘媛+陳大霞+宋旭紅+張雪+李隆云

[摘要] 該研究應用新一代測序技術Illumina HiSeq^TM4000在轉錄水平上對藥用植物玄參根部進行測序，結合生物信息學方法開展基因功能注釋和SSR位點搜索。通過測序，共獲得65 602 036條原始序列。利用生物信息學軟件拼接和組裝序列，獲得73 983條unigene，平均長度823 bp。序列同源性比較表明，56 389條unigene與其他物種具有不同程度的同源性。通過Swiss-Prot，GO，KEGG，COG比對注釋，發現520條編碼玄參次生代謝途徑關鍵酶基因和191個相關轉錄因子。利用MISA軟件在所有unigenes中共搜索到11 659個SSR位點，重復類型以二核苷酸為主。該研究所獲得的參與次生代謝的關鍵基因可為研究玄參藥用成分的生物合成和調控機制奠定基礎，獲得的大量SSR位點為后續研究玄參種質鑒定及遺傳多樣性研究提供參考。

[關鍵詞] 玄參；轉錄組；高通量測序；萜類物質

[Abstract] To investigate the profile of gene function and search for SSR， a new technology of high-throughput Solexa / Illumina sequencing was used to generate the root transcriptome of Scrophularia ningpoensis， and 65 602 036 raw reads were obtained. Based on the bioinformatics analysis and Trinity， 73 983 unigenes were obtained with an average length of 823 bp. The comparison of sequence homology in database showed that 56 389 unigenes had different degrees of homology. A total of 520 metabolic pathways related genes and 191 related transcription factors were identified by the Swiss-Prot， GO， KEGG and COG.The 11 659 SSRs were found by MISA and the highest frequency was AG/CT. In this study， we obtained numerous SSRs to provide references for the study of functional gene cloning and genetic diversity of S. ningpoensis. The key genes involved in the secondary metabolism are the basis for the study of biosynthesis and regulatory mechanism of the secondary metabolites.

[Key words] Scrophularia ningpoensis；transcriptome；high throughput sequencing；terpenoids

玄參為玄參科植物玄參Scrophularia ningpoensis Hemsl.的干燥根。玄參為我國常用中藥材，始載于《神農本草經》，列為中品，歷代藥典都有收載。味甘、苦、咸、微寒，具有清熱涼血，滋陰降火，解毒散結等功效[1]。研究發現玄參含有環烯醚萜、苯丙素、多糖等多種化學成分，具有保護心腦血管系統、抗炎、增強免疫等藥理活性[2]。長期以來，由于玄參分子生物學相關研究起步較晚，缺乏玄參生長發育相關的分子標記開發、遺傳圖譜構建以及次生代謝途徑等基礎性研究成果的支撐，玄參分子育種、藥效成分合成研究進展緩慢。高通量測序技術的出現，為研究玄參生長發育及次生代謝的分子機制提供了重要的基因資源，并為開展玄參功能基因組學研究提供了全新的思路和方法[3-4]。

高通量轉錄組測序技術已廣泛應用于生物體轉錄組基因表達分析，采用該技術能全面快速地獲取研究對象在某一狀態下基因轉錄信息，從中挖掘重要功能基因，揭示不同生物學性狀的分子機制[5-7]。開展玄參轉錄組的研究，也可能發現一些與其藥效活性成分生物合成相關的候選基因，為玄參藥效資源的充分利用奠定基礎。本研究擬在轉錄水平上，利用Illumina HiSeq^TM4000測序技術構建玄參根系轉錄組數據庫，獲得玄參轉錄本信息，并進行功能注釋及SSR位點分析，揭示玄參根系轉錄組的整體表達特征，為進一步揭示玄參有效成分的累積、道地性形成等生物學過程的分子生物學研究提供豐富的數據資源。

1 材料與方法

1.1 樣品 藥用植物玄參塊根采自重慶市武隆縣仙女山玄參GAP種植基地，采集時間為2015年8月初（塊根膨大期），經重慶市中藥研究院李隆云研究員鑒定為玄參科玄參屬植物玄參S. ningpoensis。選擇生長健壯無病害的玄參植株，純水洗凈整個塊根，用滅菌后的吸水紙吸干表面水分，迅速將塊根切成約5 mm厚的薄片，立即用液氮速凍，后放入-80 ℃冰箱保存備用。

1.2 RNA的提取與轉錄組測序 采用 Trizol Reagent （Invitrogen）法提取玄參根總RNA，使用Agilen2100生物分析儀和NanoDrop分光光度計對提取的總RNA進行質量檢測?？俁NA質檢合格后，用帶有Oligo （dT）的磁珠富集真核生物mRNA加入fragmentation buffer，將mRNA打斷成短片段，以mRNA為模板，用6堿基隨機引物（random hexamers）反轉錄合成第一條cDNA鏈，然后加入緩沖液、dNTPs、RNase H和DNA polymerase Ⅰ合成雙鏈cDNA鏈，經過QiaQuick PCR試劑盒純化并加EB緩沖液洗脫之后做末端修復、加poly （A）并連接測序接頭，然后用瓊脂糖凝膠電泳進行片段大小篩選，接著進行PCR擴增，構建好的文庫用Illumina HiSeq^TM4000進行測序。

1.3 數據的拼接與組裝 經測序獲得的原始序列（raw reads），去除里面含有帶接頭的、低質量的reads，評估測序數據質量，并對測序數據進行過濾，從而獲得干凈序列（clean reads）。本研究采用Trinity[8]對clean reads進行拼接。該軟件通過序列之間的重疊（overlap）信息組裝得到重疊群（contigs），然后局部組裝得到轉錄本（transcripts），最后用 TGICL和 Phrap 軟件對轉錄本進行同源聚類和拼接得到單基因簇（unigene）。

1.4 功能注釋與分類 通過blastx將拼接所得unigene比對到Nr[10]（Non-redundant protein database，非冗余蛋白數據庫），Nt，Swiss-Prot（SwissProt protein database，蛋白質序列數據庫），GO[11]（Gene Ontology，基因本體論數據庫）、KEGG（Kyoto Encyclopedia of Genes and Genomes，東京基因與基因組百科全書）和COG（Cluster of Orthologous Groups，蛋白質直系同源數據庫） （e-value<10^-5），從而獲得該unigene的功能注釋信息[9]和分類信息，對所有注釋信息進行整理。

1.5 SSR位點篩選 將轉錄組數據用MISA 軟件進行SSR分析。設置參數如下：總重復序列長度不低于20 bp；二核苷酸、三核苷酸、四核苷酸、五核苷酸和六核苷酸至少重復次數分別為10，7，5，4，4 [14]。

2 結果與分析

2.1 轉錄組測序與數據組裝 采用Illumina HiSeq^TM4000高通量測序技術對玄參根系轉錄組進行了測序，共得到6 560萬條raw reads以及6 456萬條clean reads。本研究clean reads Q20為97%（一般為>90%），GC量為44.82%，基本呈正態分布，質量合格。采用Trinity軟件組裝共產生109 260個轉錄本，平均長度為493 nt。一般把所有轉錄本中最長的一個視為unigene，共獲得了73 983個unigene，長度201～15 502 nt，見表1。

2.2 序列功能注釋與分類 使用BLAST程序將組裝得到的unigene與NT，NR，KOG，GO，Swissprot，KEGG數據庫進行比對，進行unigene的序列相似性分析，從而得到該unigene的蛋白質功能注釋信息。其中，匹配到 NR數據庫中的有56 389條，占全部unigene的76.21%，其后依次是Swissprot（56.44%），Nt（55.9%），GO（50.47），KO（31.05%），KOG（21.97%）。對這6種數據庫進行拓撲分析，共有9 494條unigene在所有數據庫中同時標注成功，占總unigene數的12.83%，并且在所有數據庫中至少有1種數據庫注釋成功的unigene有58 948條，占總unigene數的79.67%。

以 NR 數據庫為例進行分析，56 389條unigenes在NR數據庫中可找到相似序列。注釋基因同源序列的物種分布情況見圖1，在相似序列匹配度較高的近緣物種中，芝麻Sesamum indicum所占比例最高，為71.5%；其次是合瓣花Erythranthe guttata所占比例為14.0%，這些物種都為本研究中的序列注釋提供了參考序列。

將玄參Unigene與KOG數據庫進行對比，可預測 unigene功能并進行分類統計。結果表明，共有16 126條 unigene（14.66%）被注釋到26種KOG分類中，見圖2。從圖中可以看出unigene涉及的KOG功能類別比較全面，涉及了大多數的生命活動。如RNA加工與修飾、能量的合成與運輸、氨基酸轉運與代謝、染色體結構和動力學等。其中，“翻譯后修飾，蛋白質轉運”是最大類別，包含2 233條unigene，結合GO數據庫對玄參根系的unigene進行功能分類，可從宏觀上認識玄參根系表達基因的功能分布特征。試驗結果表明，有 37 346條unigene被注釋到GO分類，其中參與生物學過程（biological process）分類中主要聚集于細胞過程（cellular process，21 126個），代謝過程（metabolic process，19 743個）和生物調節（biological regulation，7 192個）；在細胞組分（cellular component） 主要聚集于細胞 （cell，9 369個）、細胞成分（cell part，9 364個）和細胞器（organelle，7 780個）；在分子功能（molecular function）分類中主要聚集于結合蛋白（binding，21 810個）和催化活性（catalytic activity，16 641個），見圖3。

2.3 序列代謝通路分析 根據KEGG數據庫的注釋信息能進一步得到unigene的代謝通路注釋。本研究將unigene根據參與的KEGG代謝通路分為5個分支：細胞過程（A），環境信息處理（B），遺傳信息處理（C），代謝（D）和有機系統（E），其中涉及較多的有遺傳信息處理中的翻譯（2 096個）、折疊、分類和降解（1 897條），涉及最少的是環境信息處理中的膜轉運（91條），見圖4。

結合KEGG數據庫，對玄參根系的 unigene 可能參與或涉及的代謝途徑進行了統計分析。結果表明，22 972條unigene參與到129個代謝通路中，與玄參次生代謝相關的unigene有782條。主要代謝產物有16種，這些代謝產物分別為花青素（anthocyanin）、咖啡因（caffeine）、黃酮和黃酮醇（flavone and flavonol）、類黃酮（flavonoid）、芥子油苷（glucosinolate）、異黃酮（isoflavonoid）、異喹啉類生物堿（isoquinoline alkaloid）、苯丙素（phenylpropanoid）、類固醇（steroid）、生物素（biotin）、油菜素內酯（brassinosteroid）、類胡蘿卜素（carotenoid）、萜類化合物（terpenoid）、檸檬烯和蒎烯（limonene and pinene）和玉米素（zeatin）。玄參藥用成分主要有環烯醚萜類、苯丙素、多糖、部分黃酮類等，其中注釋到萜類、苯丙素類、黃酮類物質生物合成與代謝途徑的unigene分別有56，249，52條，見圖5。

2.4 玄參次生代謝途徑相關基因的挖掘 環烯醚萜類、苯丙素類是玄參的主要藥用成分，它們的生物合成和代謝涉及到細胞色素P450、DXR-1-脫氧-D-木酮糖-5-磷酸還原異構酶（1-deoxy-D-xylulose 5-phosphate reductoisomerase）、FPPS-法呢基焦磷酸合成酶（farnesyl pyrophosphate synthase ）、HMGS-3-羥基-3-甲基戊二酰輔酶A（3-hydroxy-3-methyglutaryl-CoA）及HMGR-HMG-CoA還原酶等酶的作用[15]，以上提到的酶都存在于玄參根中，其中編碼細胞色素P450家族相關酶的unigene共搜索到504條，編碼1-脫氧-D-木酮糖-5-磷酸還原異構酶的unigene共搜索到9條，法呢基焦磷酸合成酶共搜索到1條，3-羥基-3-甲基戊二酰輔酶A 和HMG-CoA還原酶各搜索到3條，見表2。

轉錄因子也稱反式作用因子，是能夠與真核基因啟動子區域中順式作用原件發生特異性相互作用的DNA結合蛋白，通過他們之間以及與其他相關蛋白之間的相互作用，激活或抑制轉錄。本研究使用iTAK軟件對玄參轉錄組序列信息進行轉錄因子預測，發現有3 919條unigene分屬于72個轉錄因子家族。目前發現的植物萜類轉錄因子主要包括AP2/ERF類、WRKY類、鋅指類、bZIP類、bHLH類等[15]。在玄參轉錄組信息中與萜類合成相關的AP2/ERF類轉錄因子的表達豐度最高，涉及到的unigene有191條，見圖6。AP2/ERF類轉錄因子是植物特有的一類轉錄因子，AP2/ERF家族成員在結構上含有一個或多個AP2/ERF結構域。每個AP2/ERF結合域有2個保守序列塊—YRG原件和RAYD原件[16-18]。該轉錄因子已從擬南芥、煙草、水稻、玉米等多種植物中分離獲得，他們在植物的生長、發育、各種生物和非生物脅迫以及多種生理生化反應中發揮重要作用。此外，WRKY類轉錄因子的表達豐度也較高，它是近年來新發現的植物特有的鋅指型轉錄調控因子，能夠調控植物信號轉導和生理生化過程，調控植物次生代謝途徑中編碼關鍵酶基因的活性，并在植物抗病及免疫方面具有重要作用[19-20]。這些轉錄因子的發現將有助于玄參次生代謝成分生物合成途徑的進一步研究。

2.5 SSR位點分析 SSR，簡單重復序列標記（simple sequence repeats），又稱為短串聯重復序列或微衛星標記，是一類由幾個核苷酸（1～6個）為重復單位組成的長達幾十個核苷酸的重復序列，長度較短，且廣泛均勻分布于真核生物基因組中。由于重復單位的核苷酸不同以及重復次數不完全相同，造成了SSR長度的高度變異性，其中最常見的雙核苷酸重復類型，如（CA）n。一般采用SSR分子標記法對物種種質資源進行遺傳多樣性分析。本實驗利用MISA軟件在玄參根系的73 983條unigenes中共搜索到11 659個SSR位點，其中10 022條序列都存在SSR位點。SSR 的類型豐富，單核苷酸至六核苷酸重復類型均存在，所占比例變化較大，見表3。其中，二核苷酸重復所占比例最高，達到了40.13%；比例最低的是五核酸重復，僅為 0.20%；單核苷酸重復和三核苷酸重復所占比例大致相當，分別為30.46%，27.87%。在檢測結果中，共出現61種基序類型，出現頻率最高的6類基序為：AG/CT（2475），AT/AT（1316），AC/GT（885），ATC/ATG（692），AAG/CTT（590）和ACC/GGT（572）。上述 SSR 特征分析，有助于開展玄參及其同屬物種的基因組差異分析、分子標記開發和遺傳連鎖圖譜構建的研究。

3 討論

目前，高通量轉錄組測序技術已經廣泛應用于藥用植物轉錄組分析中。本研究首次采用高通量測序技術對玄參根進行轉錄組測序和功能分析，深一步挖掘其次生代謝相關基因，填補了玄參轉錄組信息的空白。測序數據采用Trinity軟件共拼接得到73 983條unigene，平均序列長度823 bp，約73%的reads參與了拼接，拼接的N50長度為1 546 bp，所測得的unigene數量基本涵蓋了全部轉錄組信息。測序數據質控合格，測序質量良好。獲得如此大的序列信息量，表明高通量測序技術是批量發現玄參功能基因的有效手段。本研究利用生物信息學方法對拼接序列進行注釋和功能分類，其中56 389條unigene在Blast、同源性搜索中得到注釋，注釋率達76.2%，剩下的17 594條unigene可能是由于長度較短而未與公共數據庫中的序列比對上，也可能是非編碼序列或者是新的基因[21]。

本研究通過同源搜索，共發現520條編碼玄參次生代謝途徑關鍵酶的相關基因和191個相關轉錄因子。這些基因的發現，為后續開展的玄參次生代謝物合成關鍵基因的鑒定和克隆提供了基礎數據。眾所周知，萜類物質結構復雜，化學合成較困難，目前主要以原植物提取獲得。因此，開展玄參次生代謝物合成關鍵酶基因及轉錄因子的表達調控分子機制尤為重要，隨著后基因組工作的深入，這些關鍵基因將作為改造植物代謝途徑的有力工具，人為控制次生代謝物的合成量。本研究所獲得的轉錄組信息不光為玄參次生代謝物生物合成研究提供基礎數據，同時也為進一步開展玄參生長發育、抗病抗逆等相關分子機制研究提供可靠信息。

此外，與傳統測序方法相比，高通量測序技術操作簡單，能夠挖掘出大量的SSR位點信息。本研究發現玄參根SSR位點11 659個，重復類型以二核苷酸為主，占全部SSR的40.13%。這些SSR位點的發現可為玄參分子標記的開發、群體遺傳多樣性分析、種質鑒定、標記輔助選擇、基因定位、親緣鑒定等方面的研究提供依據。

由于玄參未開展全基因組測序，可供參考的遺傳信息非常少，因此對玄參根轉錄組的特性分析還有待于進一步的深入研究。本研究所獲得的玄參根轉錄組信息，一方面獲得大量SSR位點，為后續研究玄參的功能基因克隆及遺傳多樣性研究提供參考；另一方面獲得了豐富的參與次生代謝的關鍵基因，也為玄參藥用成分的生物合成和調控機制奠定基礎。

[參考文獻]

[1] 中國藥典.一部[S].2015：108.

[2] Qian J， Hunkler D， Safayhi H， et al. New iridoid-related constituents and the anti-inflammatory activity of Scrophularia ningpoensis[J]. Planta Med， 1991， 57： 56.

[3] Grabherr M G， Haas B J， Yassour M， et al..Full-length transcriptome assembly from RNA-Seq data without a reference genome[J]. Nat Biotechnol， 2011， 29， 644.

[4] Fullwood M J， Wei C L， Liu E T，et al. Next-generation DNAsequencing of paired-end tags for transcriptom and genome analysis[J]. Genome Res， 2009， 19（4）： 521.

[5] Dassanayake M，Haas J S，Bohnert H J，et al. Shedding light on an extremophile life style through transcriptomics[J]. New Phytologist，2010， 183 （3）： 764.

[6] Lu T T，Lu G J，Fan D L，et al. Function annotation of the rice transcriptome at single-nucleotide resolution by RNA-Seq[J]. Genome Res，2010， 20 （1）：1238.

[7] Sangwan R S，Tripathi S，Singh J，et al. De novo sequencing and assembly of Centella asiatica leaf transcriptome for mapping of structural，functional and regulatory genes with special reference to secondary metabolism[J]. Gene，2013， 525 （2）： 58.

[8] Grabherr M G， Haas B J， Yassour M， et al.Full-length transcriptome assembly from RNA-Seq data without a reference genome[J]. Nat Biotechnol， 2011， 29： 644.

[9] Conesa A， Gtz S， García-Gómez J M， et al. Blast GO： a universal tool for annotation， visualization and analysis in functional genomics research[J]. Bioinformatics， 2005， 21（18）：3674.

[10] Altschul S F， Madden T L， Schffer A A， et al. Gapped BLAST and PSI-BLAST： a new generation of protein database search programs[J]. Nucleic Acids Res， 1997， 25：3389.

[11] Young M D， Wakefield M J， Smyth G K， et al. Gene ontology analysis for RNA-seq： accounting for selection bias[J]. Genome Biol， 2010， 11 （2） ：R14.

[12] Gtz S， García-Gómez J M， Terol J， et al. High-throughput functional annotation and data mining with the Blast2 GO suite[J]. Nucleic Acids Res， 2008， 36： 3420.

[13] Kanehisa M， Araki M， Goto S， et al. KEGG for linking genomes to life and the environment[J]. Nucleic Acids Res， 2008， 36：480.

[14] Simbaqueba J， Sanchez P， Sanchez E， et al. Development and characterization of microsatellite markers for the cape gooseberry physalisperuviana[J]. PLoS ONE， 2011， 6（10）： e26719.

[15] 趙恒偉，葛峰，孫穎，等. 植物萜類物質生物合成的相關轉錄因子及其應用前景[J].中草藥，2012，10（43）：2512.

[16] Iwase A， Mitsuda N， Koyama T， et al. The AP2/ERF transcription factor WIND1 controls cell dedifferentiation in Arabidopsis [J]. Curr Biol， 2011， 21（6）： 508.

[17] Qi W W， Sun F， Wang Q J， et al. Rice ethylene-response AP2/ERF factor OsEATB restricts internode elongation by down-regulating a gibberellin biosynthetic gene [J]. Plant Physiol， 2011， 157（1）： 216.

[18] Yang C Y， Hu F C， Li J P， et al. The AP2/ERF transcription factor AtERF73/HRE1 modulates ethylene responses during hypoxia in Arabidopsis [J]. Plant Physiol， 2011， 156（1）： 202.

[19] Rushton P J， Somssich I E， Ringler P， et al. WRKY transcription factors [J]. Trends Plant Sci， 2010， 15（5）： 1360.

[20] Tripathi P， Rabara R C， Langum T J， et al. The WRKY transcription factor family in Brachypodium distachyon[J]. BMC Genomics， 2012， 13（270）： 1.

[21] 李瀅，孫超，羅紅梅，等. 基于高通量測序 454 GS FLX的丹參轉錄組學研究[J]. 藥學學報，2010，45（4）：524.

[責任編輯 呂冬梅]