基于比較轉錄組學的多花黃精黃酮類化合物合成基因表達分析

葉碧歡， 楊 陽， 朱杰麗， 石從廣，陳友吾， 胡傳久， 宋其巖， 李海波*

（1． 浙江省林業(yè)科學研究院，浙江 杭州 310023；2． 浙江農(nóng)林大學 林業(yè)與生物技術學院，浙江 杭州 311300）

黃精屬 Polygonatum Mill. 隸屬百合科（Liliaceae），為多年生草本植物，在世界范圍內(nèi)廣泛分布。 該屬植物在我國有31 種，除南方熱帶以外，南北方各有分布[1]。 在黃精屬植物的化學成分和藥理活性研究方面，多集中于甾體皂苷和多糖二種主要活性成分[2-3]，對黃酮類化合物的報道較少且不深入，大多為總黃酮提取工藝、含量測定以及生物活性評價[4-6]。 通過比較《中國藥典》記載的3 種黃精新鮮塊莖的總黃酮質量分數(shù)顯示， 黃精P. sibiricum Red.、滇黃精P. kingianum Coll. et Hemsl.和多花黃精P. cyrtonema Hua. 的總黃酮質量分數(shù)均比較低，分別為0.018～0.035、0.015～0.030、0.004～0.034 mg/g，且在不同黃精屬種間以及不同產(chǎn)地間存在很大差異[7]。

黃酮類化合物是一類廣泛分布于植物體內(nèi)的次生代謝物，目前已從黃精屬8 種植物中共分離出6 種（高異黃酮類、異黃酮類、黃酮類、查耳酮類、二氫黃酮類、紫檀烷類）不同結構亞型的化合物54 個[8]，其中高異黃酮類（Homoisoflavones）最多，是該屬植物的特征性成分，在自然界頗為少見，有降血糖、抗炎、抗氧化、抗腫瘤、抗衰老等藥理活性，具有重要的研究價值[9-10]。因此，研究黃精屬植物中黃酮類，特別是高異黃酮類化合物的生物合成途徑與調控機制， 可為黃酮類化合物的合成生物學研究奠定基礎，也為通過外源誘導手段促進黃酮類物質的大量積累，以及今后黃酮類藥物和保健品的產(chǎn)業(yè)化研發(fā)開辟新的來源途徑。

近年來，轉錄組測序技術已廣泛用于藥用植物新基因的發(fā)現(xiàn)、分子標記挖掘、代謝途徑的確定等研究，有助于次生代謝相關基因的挖掘與功能研究[11-16]。迄今，黃精屬植物的轉錄組學研究多在黃精和滇黃精2 個物種上開展，且以多糖和甾體皂苷為主[17-21]，在多花黃精物種上僅有多糖合成相關基因的報道[20]，有關黃精屬植物黃酮類化合物生物合成途徑的分析與相關基因鑒定尚屬空白。

作者將對多花黃精組培苗不同發(fā)育期的根莖進行比較轉錄組學研究，擬基于轉錄組數(shù)據(jù)的組裝與注釋，鑒定與黃酮類化合物生物合成相關的關鍵基因，并對在多花黃精幼苗期根莖中表達的轉錄因子（Transcription factor，TF）進行鑒定與分類?；谵D錄組FPKM 數(shù)據(jù)篩選塊莖不同發(fā)育期的差異表達基因， 分析差異表達基因在次生代謝通路的富集，比較黃酮類化合物生物合成關鍵基因在不同發(fā)育期塊莖中的差異表達， 并進行實時熒光定量PCR（RT-qPCR）驗證。 本研究可彌補黃精屬植物黃酮類化合物生物合成途徑基因資源的空白，為后續(xù)研究黃酮類化合物合成與積累的分子調控機制奠定基礎。

1 材料與方法

1.1 材料來源

以多花黃精種子作為外植體，經(jīng)消毒處理后無菌誘導其植株再生，建立組培快繁體系。 組培培養(yǎng)基為1/2 MS（Murashige & Skoog）。 培養(yǎng)條件為溫度25 ℃，濕度50%，光照2500200 lx。 以培養(yǎng)于浙江省林業(yè)科學研究院林木育種組培實驗室的多花黃精Polygonatum cyrtonema 組培苗為試材。 從組培瓶中取出分別培養(yǎng)了3 個月和9 個月的多花黃精幼苗，切去莖葉，取帶芽根莖立即用液氮速凍，存儲于-80 ℃作為轉錄組學研究用材料。 兩種材料均取6 瓶， 其中3 瓶以去除須根后的塊莖作為總黃酮質量分數(shù)分析的3 個生物學重復，另3 瓶作為轉錄組測序的3 個生物學重復。

1.2 總黃酮質量分數(shù)測定

多花黃精塊莖總黃酮的測定參考文獻[22]的方法，以蘆丁作為對照品，利用紫外分光光度計測定各樣品510 nm 處的吸光度（A510nm）。 以蘆丁質量濃度為橫坐標，吸光度為縱坐標，繪制出標準曲線，得到回歸方程。 多花黃精塊莖總黃酮的提取參考文獻[23]的方法并略微修改，即稱取完全干燥粉碎的塊莖10 g，加入甲醇500 mL，于40 ℃下超聲提取1 h，過濾，取濾液60 mL，加入100 mL 酸解液（75 mL 甲醇＋25 mL 鹽酸），水浴回流酸解100 min，加甲醇定容，過濾，取濾液檢測A510nm。 根據(jù)蘆丁標準曲線和回歸方程，計算不同樣品的總黃酮質量分數(shù)。

式中：I 為總黃酮質量分數(shù)，mg/g；C 為樣品提取液的總黃酮質量濃度，mg/mL；N 稀釋倍數(shù)；V 為定容體積，mL；m 為多花黃精塊莖樣品質量，g。

1.3 轉錄組測序與序列組裝和功能注釋

轉錄組測序與序列組裝和功能注釋的具體方法參考作者所在課題組前期研究報告[24]，簡述為：采用Plant Total RNA Miniprep Purification Kit（TR02，GeneMark，Taiwan，China） 提取多花黃精根莖的總RNA。 采用SureSelect Strand-Specific RNA Library Prep for Illumina Multiplexed Sequencing （Agilent Technologies，G9691）構建測序文庫；使用Trinity 軟件（v2012-10-15）對測序得到的高質量clean reads從頭組裝， 再利用CAP3 （r12/21/07） 軟件拼接為unigene； 將unigene 序 列 比 對 到Nr、SwissProt、KEGG、COG/KOG 和Pfam 蛋白數(shù)據(jù)庫以及擬南芥基因組數(shù)據(jù)庫TAIR10（E-value≤1×10-5），得到與給定unigene 具有最高序列相似的蛋白質， 獲得蛋白質功能注釋信息。

1.4 轉錄因子鑒定

利 用Blast 程 序（https：//www.ncbi.nlm.nih.gov/）將unigene 序列與Nr 蛋白數(shù)據(jù)庫中的轉錄因子信息進行相似性比對（閾值：E-value≤1×10-5），得到與給定unigene 具有最高序列相似性的TF， 獲得TF注釋信息。

1.5 差異表達基因的篩選與功能富集分析

采 用 Cuffdiff （v 2.1.1） 軟 件 計 算 FPKM（Fragments per kilobase of transcript per million fragments mapped）值[25]作為衡量轉錄本表達水平的指標，使用錯誤發(fā)現(xiàn)率FDR（False discovery rate）進行可信度檢驗。 差異表達基因的篩選標準為：倍數(shù)差異（Fold change，F(xiàn)C）≥2，F(xiàn)DR＜0.01。 根據(jù)基因在不同樣品中的表達量進行差異表達分析與功能富集分析， 差異表達基因的KEGG 富集分析采用KOBAS（v 2. 0）軟件。 基于KEGG 注釋的代謝途徑統(tǒng)計與次生代謝相關的unigene 數(shù)量， 基于Nr 和TAIR10 注釋信息統(tǒng)計與黃酮生物合成相關的unigene 數(shù)量。根據(jù)注釋結果鑒定參與黃酮類化合物生物合成途徑的unigene 及其編碼的代謝酶。

1.6 RT-qPCR 分析基因表達

取上述1.3 中用于轉錄組測序的RNA 樣品進行RT-qPCR 分析。 以RNA 為模板使用Prime ScriptTMRT reagent Kit with gDNA Eraser Kit（Perfect Real Time）逆轉錄合成cDNA 第一鏈，逆轉錄方法參照試劑盒說明。 用于RT-qPCR 分析的內(nèi)參基因為UBQ-E2-10 和EF-1α2。 擴增目的基因和內(nèi)參基因的引物序列見2.1 部分。 使用Oligo 7.57 軟件設計特異性引物，引物合成由杭州有康生物科技有限公司完成。

RT-qPCR 分析參考作者所在課題組前期研究報告[26]，具體為：反應在LineGene 9600 Plus 上進行。 反應體系為（10 μL）：cDNA 模板1 μL，上、下游引 物 均 為0.2 μL，2×SYBR Green Mix （BioEasy Master Mix）5 μL，ddH2O 3.6 μL。 擴增程序為：恒溫段95 ℃2 min；循環(huán)段95 ℃15 s，60 ℃20 s，40 個循環(huán)；溶解段95 ℃15 s，60 ℃1 min，95 ℃15 s。 每個實驗設3 個技術重復。

采用2-△△Ct法[27]計算RT-qPCR 實驗中的基因相對表達量。

2 結果與分析

2.1 多花黃精轉錄組測序結果

基于Illumina HiSeq 2500 測序平臺獲得多花黃精3 個月和9 個月幼苗期根莖的轉錄組數(shù)據(jù)。 經(jīng)對原始數(shù)據(jù)進行統(tǒng)計， 去除raw data 中的接頭序列及低質量reads，共獲得約12 GB 的高質量clean data。通過對2 個時期轉錄本數(shù)據(jù)的合并組裝，經(jīng)拼接后共得到73218 條unigenes 序列（≥300 bp），其中長度在1 kb 以上的有17249 條。 所有unigenes 序列與各大數(shù)據(jù)庫進行注釋比對， 結果顯示， 共有35 511 條unigenes 獲得基因注釋信息， 占全部unigenes 的48.5%， 其中長度大于1 kb 以上的有16640 條（見表1）。

表1 數(shù)據(jù)庫注釋的unigene 數(shù)量Table 1 Annotated unigene numbers in databases

2.2 黃酮類化合物生物合成相關酶基因鑒定

根據(jù)對2 個不同發(fā)育期根莖樣品在各大數(shù)據(jù)庫的整體注釋，共鑒定出與黃酮類化合物生物合成相關的unigenes 83 條，編碼18 個代謝酶（見表2）。目前，異黃酮類化合物（Isoflavones）的生物合成途徑比較明確， 涉及的關鍵酶包括PAL、CHS、CHI 和異黃酮合成酶（IFS）[28-31]。 本研究在鑒定出的83 個酶基因中包括了PAL、CHS、CHI、F3H、FLS、DFR、F3′，5′H 等一批黃酮類化合物生物合成的關鍵性酶基因，但IFS 酶基因沒有得到注釋，這可能是由于在培養(yǎng)了3～9 個月的多花黃精幼苗期，其根莖中IFS 基因尚未充分表達，也可能與unigenes 的組裝拼接質量有關。

表2 多花黃精轉錄組黃酮類化合物合成途徑中的酶基因Table 2 Genes involved in the synthetic pathway of flavonoids in P. cyrtonema transcriptome

2.3 轉錄因子分析

轉錄因子是在轉錄水平上調控結構基因表達的關鍵因素之一，在調節(jié)植物生長發(fā)育、活性成分合成和環(huán)境脅迫中的基因表達發(fā)揮了主要作用。 許多轉錄因子家族由不同的基因家族組成， 包括bHLH、ERF、MYB-relaed、C2H2、NAC（NAM，ATAF1/2，CUC1/2）等。轉錄組數(shù)據(jù)顯示，在多花黃精幼苗期的根莖中，共有1359 個編碼轉錄因子的基因表達，可以分類為56 個TF 家族（見圖1），以bHLH 類型最豐富（94 個，占比6.92%），其次是ERF（82 個，占比6.03%），MYB-relaed（79 個，占比5.81%）、C2H2（67個，占比4.93%）和NAC（58 個，占比4.27%）。MYBrelaed 和bHLH 轉錄因子單獨或協(xié)作調節(jié)類黃酮次生代謝生物合成途徑中結構基因的表達，從而控制植物體內(nèi)類黃酮的生物合成[32]。 利用bHLH 轉錄因子調節(jié)藥用植物活性物質生物合成的研究備受關注，其對黃酮類物質（花青素等）、生物堿（尼古丁等）、 萜類等活性成分生物合成的調控作用已被證實[33]。 本研究中鑒定的轉錄因子為今后研究多花黃精活性成分（多糖、皂苷和黃酮類化合物等）的積累與外源調控機制提供了豐富的基礎數(shù)據(jù)。

圖1 多花黃精轉錄組中轉錄因子的分布Fig. 1 Distribution of transcription factors in P.cyrtonema transcriptome

2.4 不同發(fā)育期差異表達基因的篩選與分析

篩選多花黃精2 個不同發(fā)育期根莖樣品間表達水平顯著差異的基因，形成差異基因數(shù)據(jù)庫。 通過分析差異基因數(shù)據(jù)庫共獲得2602 個差異基因，占總被注釋unigene 的7.33%，其中上調基因1866個，下調基因736 個（見圖2）。這些差異表達基因為后續(xù)開展次生代謝通路分析提供了基礎數(shù)據(jù)。

圖2 多花黃精不同發(fā)育期塊莖的差異表達基因Fig. 2 Differentially expressed genes in P. cyrtonema rhizome at different development stages

2.5 差異表達基因的KEGG 通路富集分析

為了進一步了解2 個不同發(fā)育期根莖樣品間差異表達基因參與的代謝通路，對這些差異表達基因進一步利用KOBAS 進行KEGG 富集分析， 發(fā)現(xiàn)有497 條unigenes 注釋到了109 個代謝通路中，且主要富集在代謝組（Metabolism）和遺傳信息處理組（Genetic information processing）。 其中富集的前20條通路統(tǒng)計如圖3 所示。

圖3 差異表達基因的KEGG 注釋分類Fig. 3 KEGG classification for differentially expressed genes in P. cyrtonema rhizome

在497 條unigenes 富 集 的109 個KEGG 代 謝通路中， 其中186 條定位到了19 個次生代謝物生物合成通路中，數(shù)量最多的是參與次生代謝物合成（Biosynthesis of secondary metabolites）通路，共118條，最少的是咖啡堿代謝（Caffeine metabolism）、吲哚類生物堿合成（Indole alkaloid biosynthesis）、黃酮和黃酮醇類合成（Flavone and flavonol biosynthesis）和硫代葡萄糖苷合成（Glucosinolate biosynthesis）通路，都僅只有1 條（見表3）。 此外，這些差異表達的酶基因也參與了像二苯乙烯類（Stilbenoids）、苯丙素（Phenylpropanoids）、萜 類（Terpenoids）、托 品 烷 類（Tropanes）、類胡蘿卜素（Carotenoids）、油菜素甾醇（Brassinosteroids）、 黃 酮 類 （Flavonoids） 和 甾 體（Steroids）等次生代謝物的合成。

表3 與多花黃精次生代謝相關的unigenesTable 3 Unigenes related to secondary metabolism of P. cyrtonema

2.6 黃酮類化合物合成途徑中基因表達差異的初步分析

對多花黃精組培苗在3 個月和9 個月發(fā)育期塊莖的總黃酮質量分數(shù)分析顯示，3 個月發(fā)育期塊莖中總黃酮的平均質量分數(shù)為0.008 mg/g，9 個月為0.012 mg/g，表明隨著根莖的不斷發(fā)育，其總黃酮質量分數(shù)也在不斷積累。 基于FPKM 值比較在3 個月和9 個月發(fā)育期塊莖的黃酮類化合物合成途徑中關鍵酶基因的表達水平， 發(fā)現(xiàn)9 個基因顯著上調， 分別是PAL、4CL、DFR、ANS、FLS、UGT、F3′H、COMT 和F3′，5′MT 基因， 上調倍數(shù)為1.9～11.2；3個基因顯著下調， 分別是CHS、CHI 和F3′，5′H 基因， 下調倍數(shù)為0.3～0.5；6 個基因的表達無明顯表達差異， 分別是CYP73A、F3H、LAR、CCOMT、HCT和C3′H 基因，這些基因的表達變化倍數(shù)為0.9～1.3（見表4）。 這表明隨著塊莖中總黃酮質量分數(shù)的不斷積累，其黃酮類化合物合成途徑中基因的差異表達并不一致，這也預示了多花黃精復雜的黃酮類化合物合成機制。

表4 黃酮類化合物合成途徑中的基因表達差異Table 4 Different gene expression involved in flavonoids biosynthesis

2.7 RT-qPCR 驗證黃酮類化合物生物合成相關酶基因的表達

為驗證18 個黃酮類化合物合成酶基因在不同發(fā)育期根莖中基于FPKM 值表達變化的可靠性，選取其中的7 個酶基因進一步進行RT-qPCR 分析，并繪制標準曲線。 結果顯示，所有酶基因和內(nèi)參基因（UBQ-E2-10 和EF-1α2）標準曲線的相關系數(shù)（R2）均≥0.97，擴增效率均在90%～110%（見表5），表明所有RT-qPCR 引物的擴增具有高特異性， 且擴增效率一致，滿足利用2-△△Ct法進行相對定量分析。基因表達結果顯示（見圖4），這7 個酶基因在不同發(fā)育期根莖中的差異表達與轉錄組FPKM 結果基本吻合，只是這兩種分析結果的差異表達倍數(shù)略有不同。

表5 用于RT-qPCR 分析的多花黃精黃酮類化合物生物合成相關酶基因Table 5 Genes involved in flavonoid biosynthesis in P. cyrtonema used for RT-qPCR analysis

圖4 多花黃精7 個黃酮類化合物生物合成相關酶基因的表達Fig. 4 Expression of 7 genes involved in flavonoid biosynthesis in P. cyrtonema

3 結語

利用二代高通量測序技術，對黃精屬原生藥材多花黃精幼苗期根莖進行轉錄組測序，經(jīng)組裝拼接共產(chǎn)生了73218 條unigenes 序列， 其中35511 條被注釋；共鑒定出83 條unigenes 序列，分別編碼18個參與黃酮類化合物生物合成的關鍵酶；在幼苗期根莖中表達的1359 個轉錄因子中， 以bHLH 最為豐富，其次是ERF、MYB-relaed、C2H2 和NAC?；贔PKM 表達分析顯示， 隨著多花黃精塊莖中黃酮化合物的不斷積累，其黃酮類化合物合成途徑中基因的表達并不一致，18 個酶基因中的9 個顯著上調，3個顯著下調，6 個無明顯表達差異。與前期報道的藥用植物黃酮類化合物轉錄組學研究[34-38]相比，本研究鑒定出的黃酮類化合物合成酶的基因數(shù)量和種類有些許不同，例如細胞色素CYP450 基因、異黃酮合成酶（IFS）基因等，這可能是由于藥用植物在幼苗期根莖發(fā)育過程中黃酮類化合物生物合成相關酶基因的表達與葉、花、莖等器官中的表達水平不同而致，也可能與多花黃精的物種特異性有關，這也預示了不同物種間黃酮類化合物的生物合成機理存在一定差異，百合科的黃精屬植物可能具有更為復雜的黃酮類化合物生物合成機制。 本研究為后續(xù)解析黃精屬植物黃酮類活性成分的合成、積累與外源調控機制奠定了基礎。