多油辣木轉錄組高通量測序及分析

歐奇+李鑫+田洋+曾千春

摘要：為了研究多油辣木（Moringa oleifera Lam.）中黃酮類化合物生物合成的分子基礎，通過Illumina HiSeq 2000高通量測序技術對辣木莖、葉進行轉錄組測序，利用Trinity軟件將數據組裝成Unigene，基于BLAST對所有Unigene進行功能注釋，共獲得49 365個Unigene，平均長度為903 bp。通過GO分類，15 959個Unigene被分成生物學過程、細胞組分和分子功能3個主要類別。通過KOG分類，8 713個Unigene被分為26個種類。通過KEGG分類，8 397個Unigene分屬于130個代謝途徑，其中代謝所含Unigene最多，共3 620個。在代謝途徑中，找到參與黃酮類化合物合成相關的Unigene 45個，其中包括查爾酮合酶、查爾酮異構酶、黃烷酮-3-羥化酶、黃酮醇合成酶和二氫黃酮醇還原酶等。研究結果為挖掘多油辣木黃酮類化合物生物合成關鍵基因提供了基礎數據，并為下一步的資源開發(fā)和利用奠定了基礎。

關鍵詞：多油辣木；轉錄組；黃酮類化合物；高通量測序

中圖分類號： Q522+.6 文獻標志碼： A 文章編號：1002-1302（2017）20-0071-05

多油辣木（Moringa oleifera Lam.）屬辣木科（Moringaceae）辣木屬（M. Adans），是一種具有獨特經濟價值的多年生熱帶植物，現廣泛種植于亞洲、非洲和中美洲[1]。多油辣木中含有豐富的脂肪酸、蛋白質、多糖、維生素和黃酮類化合物[2]。相關研究表明，多油辣木提取液富含黃酮類化合物，具有顯著的抗炎[3]、降血脂[4]、抗腫瘤[5]和抗氧化[6-7]等藥用效果。黃酮類化合物包括黃酮、異黃酮、黃酮醇、查爾酮等，是植物在長期進化過程中產生的一類次生代謝產物[8]，黃酮類化合物廣泛存在于高等植物中，具有抗氧化、抗炎、降血脂、抗腫瘤和擴張血管等藥理活性[9]。黃酮類化合物合成的一般途徑已經明確，首先，丙二酰輔酶A（malonyl-CoA）和4-香豆酰輔酶A（4-coumaroyl-CoA）在查爾酮合酶（chalcone synthase，簡稱CHS）[10]催化下合成查爾酮；接著，查爾酮在查爾酮異構酶（chalcone isomerase，簡稱CHI）[11]的催化下合成黃烷酮；然后，黃烷酮能被黃酮合成酶（flavone synthase，簡稱FNS）[12]催化合成黃酮，可被異黃酮合成酶（isoflavone synthase，簡稱IFS）[13]催化合成異黃酮，也能被黃烷酮-3-羥化酶（flavanone-3-hydroxyalse，簡稱F3H）[14]催化合成二氫黃酮醇；最后，二氫黃酮醇在黃酮醇合成酶（flavonol synthase，簡稱FLS）[15]的催化下合成黃酮醇。黃酮類化合物合成關鍵酶還包括二氫黃酮醇還原酶（dihydroflavonol 4-reductase，簡稱DFR）[16]、類黃酮3-O-葡萄糖基轉移酶（flavonol 3-O-glucosyltransferase，簡稱3GT）[17]等。

轉錄組（transcriptome）是指細胞或組織內全部RNA轉錄本的集合，代表基因在不同生命階段、不同生理狀態(tài)、不同組織類型以及不同環(huán)境條件下的表達水平，包括基因表達的所有轉錄產物（mRNA、Non-coding RNA）的總和。轉錄組測序（RNA sequencing）是利用Illumina HiSeq 2000高通量測序技術進行cDNA測序，能更全面快速地揭示生物個體特定組織和特定時期的基因表達情況[18]。與基因組學研究相比，轉錄組學更偏重基因編碼區(qū)域，它能進一步呈現細胞內生命活動，能為詮釋細胞的功能提供具有價值的參考信息[19]。由于轉錄組學在基因功能研究中的重要作用，近年來，轉錄組的研究對象逐漸從模式生物和重要農作物向其他植物延伸。本研究首次對多油辣木進行轉錄組測序并加以分析，以期從中挖掘具有重要功能的基因。

1 材料與方法

1.1 轉錄組測序

多油辣木樣品采自云南昆明，選取同一生活環(huán)境下的一年生辣木莖、葉作為混合研究材料，于-80 ℃超低溫冰凍保存?zhèn)溆?。用TRIzol法[20]提取材料總RNA，構建測序文庫，使用Illumina HiSeq 2000進行測序。

1.2 測序序列拼裝

通過Base Calling將測序得到的原始圖像轉化為序列數據，這些數據稱為原始測序序列（Raw Reads）。將獲得的原始測序序列進行過濾得到處理序列（Clean Reads）。然后利用Trinity軟件[21]對處理序列進行拼接。取每條基因中最長的轉錄本（transcript）作為Unigene，以此進行后續(xù)分析。

1.3 Unigene功能注釋、GO分類和代謝路徑分析

通過BLAST將Unigene與Nr、Nt、Pfam、KOG、Swiss-Prot、KEGG、GO七大數據庫進行比對（E值<0.000 01），獲得序列相似性最高的蛋白，從而將此蛋白信息用作該Unigene的蛋白功能注釋信息。GO是一套國際標準化的基因功能描述的分類系統(tǒng)，分為生物過程（biological process）、分子功能（molecular function）和細胞組分（cellular component）三大類。它們分別用來描述基因編碼的產物所參與的生物過程、所具有的分子功能及所處的細胞環(huán)境。通過Blast2GO軟件[22]對Unigene進行GO注釋，然后用WEGO軟件[23]對其進行GO功能分類統(tǒng)計，從宏觀上了解辣木基因功能分布特征。KEGG是系統(tǒng)分析基因產物和化合物在細胞中的代謝途徑以及這些基因產物功能的數據庫[24]。根據KEGG注釋進一步得到Unigene的代謝路徑。

2 結果與分析endprint

2.1 轉錄組序列的產出、組裝和基因表達

通過Illumina Hiseq 2000高通量測序共獲得了 57 818 454 條高質量短讀序片段，總長度為9 309 341 614 bp，Q20值為95.89%，GC含量占比為46.92%。由Trinity軟件組裝共獲得72 527個Contig，總長度為88 344 692 bp，平均長度為1 218 bp，N50長度為2 144 bp。共組裝49 365個Unigene，總長度為44 565 055 bp，平均長度為903 bp，N50長度為 1 771 bp（表1）。對Contig、Unigene的長度分布特征進行分析可知，在總Contig中，0～300 bp區(qū)段所含的Contig比例最高，為23.34%，共16 931個；1 001～2 000 bp區(qū)段的Contig占21.47%；>2 000 bp區(qū)段的Contig占20.73%。在Unigene中，同樣以0～300 bp區(qū)段所含Unigene比例最高，為 32.38%，共15 983個；1 001～2 000 bp區(qū)段的Unigene占 15.14%；大于2 000 bp區(qū)段的Unigene占12.75%（表2）。

2.2 功能注釋

使用BLAST將所有Unigene與7大數據庫進行比對，由表3可以看出，在49 365個Unigene中，共有23 810個Unigene獲得注釋信息，占總Unigene數的48.23%（表3）。在Nr數據庫中注釋Unigene數量最多，為21 442個，占總數的43.43%。在KEGG、KOG數據庫中注釋的Unigene數量都在10 000個以下，分別為8 397個（17.01%）、8 713個（17.65%）；其他數據庫中注釋的Unigene數量都在10 000～20 000個之間，但仍有25 555個（51.77%）功能未知的Unigene。

2.3 GO注釋

通過GO分析，15 959個Unigene被分成了生物學過程、細胞組分和分子功能3個主要類別（圖1）。在“生物學過程”類別中，細胞過程（9 229個Unigene）功能組所含Unigene數量最多，細胞聚集（1個Unigene）所含Unigene數量最少。在“細胞組分”類別中，細胞（5 075個Unigene）和細胞部分（5 075 個Unigene）功能組所含Unigene數量最多，共質體（1個Unigene）所含Unigene數量最少。在“分子功能”類別中，結合功能（9 009個Unigene）所含Unigene數量最多，金屬伴侶（3個Unigene）所含Unigene數量最少。

2.4 KOG注釋

為了進一步評估Unigene的完整性和注釋的有效性，對 8 713 個Unigene進行KOG注釋和分類，共獲得26個種類（圖2）。在26個KOG種類中，一般功能預測（R）為最大類，共有1 492個Unigene；然后是翻譯后修飾、蛋白轉換和分子伴侶（O），有1 175個Unigene；未知功能蛋白（X）、細胞運動（N）為最小類，分別只有1、3個Unigene。

2.5 KEGG注釋

為了識別辣木中活性高的代謝通路，對8 397個Unigene進行KEGG代謝通路分析，將其分為代謝（3 620個Unigene）、遺傳信息調控（1 796個Unigene）、環(huán)境信息調控（337個Unigene）、有機系統(tǒng)（346個Unigene）和細胞過程（382個Unigene）5大類130個代謝途徑（圖3）。其中包含Unigene數量最多的是遺傳信息調控類中的翻譯，共有749個；代謝類中的糖代謝有732個Unigene；細胞過程類中的運輸和代謝有382個Unigene；有機系統(tǒng)類中的環(huán)境適應性有346個Unigene；環(huán)境信息調控類中的膜運輸包含Unigene數量最少，僅有55個。

2.6 黃酮類化合物合成相關Unigene

基于數據庫中已知的基因信息，找到與黃酮類化合物合成相關的Unigene 45個，其中與查爾酮合酶有關的Unigene有3個，每100萬片段中來自某一基因每1 000堿基長度的片段數（expected number of fragments per kilobase of transcript sequence per millions base pairs sequenced，簡稱FPKM，FPKM同時考慮了測序深度、基因長度對片段計數的影響，是目前最為常用的基因表達水平估算方法）[25]最高的Unigene基因號為c32808_g1；與查爾酮異構酶有關的Unigene有5個，FPKM值最高的Unigene基因號為c20602_g1；與黃烷酮-3-羥化酶有關的Unigene有4個， FPKM值最高的Unigene基因號為c32320_g1；與異黃酮合成酶有關的Unigene有5個，FPKM值最高的Unigene基因號為c11721_g1；與黃酮醇合成酶有關的Unigene有4個，FPKM值最高的Unigene基因號為c16354_g1；與二氫黃酮醇還原酶有關的Unigene有2個，FPKM值最高的Unigene基因號為c18805_g2；與黃酮3-O-葡萄糖基轉移酶有關的Unigene有4個，FPKM值最高的Unigene基因號為c13219_g1（表4）。

3 討論與結論

本研究首次采用Illumina HiSeq 2000高通量測序技術對多油辣木進行轉錄組測序，共獲得57 818 454條高質量短讀序，經拼裝后得到Unigene 49 365個，平均長度為903 bp，N50長度為1 771 bp。在Unigene中，0～300 bp區(qū)段所含比例最高，為32.38%；301～1 000 bp區(qū)段占39.74%；1 001～2 000 bp 區(qū)段占15.14%，>2 000 bp區(qū)段占12.75%。由此可見Unigene的整體長度分布均勻，組裝的質量和長度都能滿足轉錄組分析的基本要求。測序結果的好壞除了與組裝質量和長度有關外，還與測序數據量密切相關，數據量越大得到的基因序列，越長越準確。相關研究表明，辣木基因組總大小為289 Mb[26]，因此本研究8Gb的測序量足以保證獲得足夠長的組裝序列和足夠多的Unigene。在與7大數據庫比對時，發(fā)現有25 555（51.77%）個功能未知的Unigene。對這些Unigene進行深入研究，極有可能獲得與多油辣木特性有關的新基因。endprint

黃酮類化合物廣泛存在于高等植物中[27-28]，本研究找到多油辣木中與黃酮類化合物合成相關的Unigene 45個，包括CHS、CHI、F3H、FLS和DFR等關鍵合成基因，與其他植物中參與生物合成的關鍵基因基本一致[29]。多油辣木中與黃酮、黃酮醇合成的相關Unigene有8個，與異黃酮生物合成相關的Unigene有6個，與花青素生物合成有關的Unigene有5個，與其他黃酮類化合物合成有關的Unigene有26個。這些基因與淫羊藿屬（Epimedium L.）[30]植物存在一定差異，造成差異的原因可能是不同物種所含的黃酮類化合物存在差異、特定時間和特定組織內基因表達存在差異、轉錄組測序過程中提取和拼接存在誤差。

田洋等完成了多油辣木基因組草圖，測序結果表明，辣木基因組289 Mb，共注釋出19 465個基因，每個基因的平均長度為3 354.22 bp。這19 465個基因在GO數據庫中共注釋基因有10 476個，占53.82%；在KEGG中注釋了10 936個基因，占56.18%；未知功能的基因有1 166個，占5.99%[26]。該團隊還挖掘了熱休克蛋白（HSP）基因、BAK1（BRI1 associated receptor kinase 1）基因、γ-氨基丁酸和谷甾醇的合成途徑相關基因。通過比較可知，本研究獲得的Unigene （49 365個）和在GO數據庫中注釋的Unigene（15 959個）均較多，但在KEGG數據庫中注釋的Unigene（8 397個）較少。本研究彌補了多油辣木轉錄組數據的缺陷，獲得了多油辣木黃酮類化合物生物合成相關的候選基因，為多油辣木資源開發(fā)和利用提供了重要的理論基礎。

參考文獻：

[1]劉昌芬，李國華. 辣木的研究現狀及其開發(fā)前景[J]. 云南熱作科技，2002，25（3）：20-24.

[2]劉鳳霞，王苗苗，趙有為，等. 辣木中功能性成分提取及產品開發(fā)的研究進展[J]. 食品科學，2015，36（19）：282-286.

[3]Waterman C，Cheng D M，Rojas-Silva P，et al. Stable，water extractable isothiocyanates from Moringa oleifera leaves attenuate inflammation in vitro[J]. Phytochemistry，2014，103（10）：114-122.

[4]Ghasi S，Nwobodo E，Ofili J O. Hypocholesterolemic effects of crude extract of leaf of Moringa oleifera Lam in high-fat diet fed wistar rats[J]. Journal of Ethnopharmacology，2000，69（1）：21-25.

[5]Jung I L. Soluble extract from Moringa oleifera leaves with a new anticancer activity[J]. PLoS One，2014，9（4）：e95492.

[6]Vongsak B，Sithisarn P，Gritsanapan W. Bioactive contents and free radical scavenging activity of Moringa oleifera leaf extract under different storage conditions[J]. Industrial Crops and Products，2013，49（4）：419-421.

[7]Verma A R，Vijayakumar M，Mathela C S，et al. In vitro and in vivo antioxidant properties of different fractions of Moringa oleifera leaves[J]. Food and Chemical Toxicology，2009，47（9）：2196-2201.

[8]諸 姮，胡宏友，盧昌義，等. 植物體內的黃酮類化合物代謝及其調控研究進展[J]. 廈門大學學報（自然科學版），2007，46（增刊1）：136-143.

[9]Chen Z，Hu Y，Wu H，et al. Synthesis and biological evaluation of flavonoids as vasorelaxant agents[J]. Bioorganic & Medicinal Chemistry Letters，2004，14（15）：3949-3952.

[10]Koes R E，Quattrocchio F，Mol J N M. The flavonoid biosynthetic pathway in plants：function and evolution[J]. Bioessays，1994，16（2）：123-132.

[11]Jez J M，Bowman M E，Dixon R A，et al. Structure and mechanism of the evolutionarily unique plant enzyme chalcone isomerase[J]. Nature Structural Biology，2000，7（9）：786-791.

[12]Martens S，Forkmann G，Matern U，et al. Cloning of parsley flavone synthase I[J]. Phytochemistry，2001，58（1）：43-46.endprint

[13]Jung W，Yu O，Lau S M，et al. Identification and expression of isoflavone synthase，the key enzyme for biosynthesis of isoflavones in legumes[J]. Nature Biotechnology，2000，18（2）：208-212.

[14]Jin Z，Grotewold E，Qu W，et al. Cloning and characterization of a flavanone 3-hydroxylase gene from Saussurea medusa[J]. DNA Sequence，2005，16（2）：121-129.

[15]Forkmann G，Martens S. Metabolic engineering and applications of flavonoids[J]. Current Opinion in Biotechnology，2001，12（2）：155-160.

[16]Dick C A，Buenrostro J，Butler T，et al. Arctic mustard flower color polymorphism controlled by petal-specific downregulation at the threshold of the anthocyanin biosynthetic pathway[J]. PLoS One，2011，6（4）：e18230.

[17]Owens D K，Mcintosh C A. Identification，recombinant expression，and biochemical characterization of a flavonol 3-O-glucosyltransferase clone from Citrus paradisi[J]. Phytochemistry，2009，70（11/12）：1382-1391.

[18]Ajay S S，Parker S C，Abaan H O，et al. Accurate and comprehensive sequencing of personal genomes[J]. Genome Research，2011，21（9）：1498-1505.

[19]吳 瓊，孫 超，陳士林，等. 轉錄組學在藥用植物研究中的應用[J]. 世界科學技術：中醫(yī)藥現代化，2010，12（3）：457-462.

[20]姚寧濤，祝建波，鄧福軍. 改良Trizol法快速提取棉葉片總RNA[J]. 生物技術通報，2010（7）：125-127.

[21]Grabherr M G，Haas B J，Yassour M，et al. Full-length transcriptome assembly from RNA-Seq data without a reference genome[J]. Nature Biotechnology，2011，29（7）：644-652.

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

[23]Ye J，Fang L，Zheng H，et al. WEGO：a web tool for plotting GO annotations[J]. Nucleic Acids Research，2006，34（Web Server issue）：W293-W297.

[24]Kanehisa M，Araki M，Goto S，et al. KEGG for linking genomes to life and the environment[J]. Nucleic Acids Research，2008，36（Database issue）：D480-D484.

[25]Trapnell C，Williams B A，Pertea G，et al. Transcript assembly and abundance estimation from RNA-Seq reveals thousands of new transcripts and switching among isoforms[J]. Nature Biotechnology，2010，28（5）：511-515.

[26]田 洋，曾 嚴，張 靜，等. 辣木（Moringa oleifera Lam.）的高質量參考基因組[J]. 中國科學：生命科學，2015（5）：488-497.

[27]王興娜，汪 晶，黃午陽，等. 青梅果實不同極性組分的抗氧化活性[J]. 江蘇農業(yè)學報，2016，32（1）：211-215.

[28]朱海軍，生靜雅，張普娟，等. 貯藏溫度對薄殼山核桃抗氧化功能及品質的影響[J]. 江蘇農業(yè)學報，2015，31（2）：449-453.

[29]康亞蘭，裴 瑾，蔡文龍，等. 藥用植物黃酮類化合物代謝合成途徑及相關功能基因的研究進展[J]. 中草藥，2014，45（9）：1336-1341.

[30]張華峰，王 瑛，黃宏文.endprint