黃秋葵花和果莢轉錄組測序及類黃酮代謝差異表達分析

姚運法，張少平，練冬梅，賴正鋒，黃慧明，洪建基

(福建省農(nóng)業(yè)科學院 亞熱帶農(nóng)業(yè)研究所，福建漳州 363005)

黃秋葵，學名咖啡黃葵[Abelmoschusesculentus(Linn.)Moench],為錦葵科(Malvaceae)秋葵屬(Abelmoschus)一年或多年生草本植物。原產(chǎn)于非洲，自20世紀90年代初引入中國，現(xiàn)國內各地均有栽培。黃秋葵花和果莢作為重要開發(fā)價值的器官，富含蛋白質、果膠多糖、總黃酮等營養(yǎng)成分[1-4]，其花主要開發(fā)黃秋葵花茶，經(jīng)濟附加值高[5]；果莢具有降血壓、血脂和提高機體抗疲勞等功效，在黃秋葵保健品開發(fā)方面具有巨大潛力[6]。當前，國內對黃秋葵研究主要從栽培技術、營養(yǎng)成分提取、功效分析和產(chǎn)品加工[7-10]等方面，但其基礎理論研究方面還很薄弱， 例如花、果莢中類黃酮物質等主要功能成分及相關代謝途徑等均屬研究空白。

類黃酮(flavonoids)屬于植物重要次生代謝產(chǎn)物之一，是指2個苯環(huán)(A-與B-)通過3個碳原子連接形成具有C6-C3-C6基礎結構的一類化合物[11]，由于其純凈狀態(tài)呈現(xiàn)黃色，故稱黃酮。研究表明，類黃酮物質對人體具有抗癌、抗氧化、抗動脈硬化等功能[12-13]，據(jù)其結構的差異將類黃酮主要分為黃烷酮(flavanols)、異黃酮(isoflavones)、黃酮(flavones)、黃酮醇(flavonols)、二氫黃酮(2H-flavanones)與花色素(anthoyanidins)等6大類[14]。

轉錄組測序技術是連接基因組與代謝組的重要紐帶，本團隊前期對紫色黃秋葵葉片[15]、果莢[16]轉錄組測序和分析，初步探討了黃秋葵次生代謝物質合成的遺傳基礎。本研究利用Illumina Hi-seq 2500高通量測序技術，根據(jù)黃秋葵花、果莢的轉錄組數(shù)據(jù)及其類黃酮代謝途徑關鍵基因分析，探討黃秋葵花、果莢類黃酮合成機理和關鍵差異表達基因，為后續(xù)黃秋葵關鍵基因克隆和功能驗證、遺傳改良和加工利用等提供研究基礎。

1 材料和方法

1.1 試驗材料

試驗材料為黃秋葵品種‘閩秋葵3號’，種植于福建省農(nóng)業(yè)科學院亞熱帶農(nóng)業(yè)研究所試驗農(nóng)場，2017年6月10日種植，8月15日開始采集樣品(花/果莢)。采集同一株黃秋葵植株花朵(上午9:00)和果莢(開花后8 d)，花朵只保留完整花瓣，果莢去除果柄部位?；ê凸v均取3重復，后立即液氮速凍，置于-70 ℃超低溫冰箱備用。

1.2 試驗方法

1.2.1高通量測序及數(shù)據(jù)組裝提取黃秋葵花、果莢總RNA，分別采用Nanodrop、Qubit2.0、Aglient 2100技術檢測RNA樣品的純度、濃度和完整性等，構建cDNA文庫，再分別使用Qubit2.0、Agilent 2100和Q-PCR對文庫的質量進行檢測。合格后，用Illumina HiseqTM2500進行測序。測序讀長為PE125。獲得原生數(shù)據(jù)，進行數(shù)據(jù)過濾，去除Reads中測序引物、接頭等人工序列，用Trinity對有效數(shù)據(jù)(Clean Data)進行組裝。

將測序Reads構建K-mer庫，去除錯誤的K-mer；將高頻率的K-mer作為種子向兩端進行擴展，不斷循環(huán)直到耗光K-mer庫；再對Contig進行聚簇，得到片段集合(Component)；對每個Component中的Contig構建De Bruijn圖；De Bruijn圖進行簡化；解開De Bruijn圖，獲得轉錄本序列。

1.2.2基因FPKM值估計與差異表達分析采用Bowtie[17]將測序得到的Reads與單基因數(shù)據(jù)(Unigene)庫進行比對，根據(jù)比對結果，結合RSEM[18]進行表達量水平估計。利用FPKM值(fragments Per Kilobase of transcript per Million mapped reads)[19]表示對應Unigene的表達豐度。使用EBSeq進行差異表達分析[20]。采用Benjamini-Hochberg方法對原假設試驗得到的顯著性P值進行校正，校正后P值，即偽發(fā)現(xiàn)率(false discovery rate)小于0.01且差異倍數(shù)(fold change, FC)≥2作為篩選標準[21]。其中，F(xiàn)C表示兩樣品(組)間FPKM的比值。

1.2.3差異基因的功能注釋與代謝途徑富集分析使用Blast軟件將Unigene序列與NR、Swiss-Pro、GO、COG、KOG、KEGG數(shù)據(jù)庫比對，再用KOBAS2.0得到Unigene在KEGG中的Orthology，預測氨基酸序列，使用HMMER軟件與Pfam數(shù)據(jù)庫分析，獲得Unigene的注釋信息。系統(tǒng)分析基因產(chǎn)物代謝途徑及功能，并將差異表達基因(DEGs)對比到KEGG數(shù)據(jù)庫，得到DEGs的代謝途徑。

1.2.4類黃酮代謝關鍵差異基因分析利用注釋信息檢索法，對黃秋葵花、果莢類黃酮(黃烷酮、黃酮、異黃酮、黃酮醇、黃烷酮和花青素)關鍵詞進行數(shù)據(jù)庫檢索，分析類黃酮物質關鍵基因在KEGG數(shù)據(jù)庫功能注釋及代謝路徑。

1.2.5類黃酮代謝關鍵差異基因驗證取 1 μg黃秋葵花或果莢的總 RNA，利用反轉錄試劑盒反轉錄成 cDNA，設計合成引物(表1)，采用qRT-PCR檢測黃秋葵花、果莢組織部位與類黃酮代謝相關差異表達基因，設置3個重復；統(tǒng)計 8個基因 (含 1 個內參) 在待測樣品中的 Ct 值，計算相對表達量。

2 結果與分析

2.1 數(shù)據(jù)組裝及分析

經(jīng)高通量測序和質量控制，共獲得15.26 Gb有效數(shù)據(jù)，其中黃秋葵花為7.12 Gb，果莢為8.14 Gb，堿基百分比均達到91.0%以上(表2)。數(shù)據(jù)質量良好，適宜后續(xù)分析。對組裝結果進行統(tǒng)計(表3)，轉錄本序列組裝出154 721條轉錄本序列，平均長度為937.5 bp，N50長1 308 bp；由單基因序列(Unigene)組裝出65 436條Unigene，平均長度為782.32 bp，N50長1 215 bp。

2.2 Unigene功能注釋

經(jīng)與NR、Swiss-Prot、KEGG、COG、KOG、GO和Pfam數(shù)據(jù)庫的比對，對Unigene功能注釋進行統(tǒng)計(表4)，在65 436條單基因序列中，共獲得39 245條Unigene的注釋結果，占單基因序列總數(shù)59.97%。其中與COG數(shù)據(jù)庫比對，獲得9 366條同源序列，占單基因序列總數(shù)14.31%；與KOG數(shù)據(jù)庫比對，獲得20 535條同源序列，占單基因序列總數(shù)31.38%；與GO數(shù)據(jù)庫比對，獲得15 931條同源序列，占單基因序列總數(shù)24.35%；與KEGG數(shù)據(jù)庫比對，獲得14 110條同源序列，占單基因序列總數(shù)21.56%；與Pfam數(shù)據(jù)庫比對，獲得22 580條同源序列，占單基因序列總數(shù)34.51%；與Swissprot數(shù)據(jù)庫比對，獲得24 494條同源序列，占單基因序列總數(shù)37.43%；與NR數(shù)據(jù)庫比對，獲得38 905條同源序列，占單基因序列總數(shù)59.46%。

表1 實時熒光定量PCR引物

表2 有效數(shù)據(jù)評估統(tǒng)計

2.3 DEGs的篩選與功能注釋

2.3.1DEGs篩選通過DEGs的篩選，獲得黃秋葵花與果莢DEGs 1 336個，其中表達量上調的基因319個，下調基因1 017個；將DEGs單基因序列分別注釋到COG、GO、KEGG、KOG、Pfam、Swiss-Prot、NR7大數(shù)據(jù)庫，共有1 131個基因獲得功能注釋，其中GO數(shù)據(jù)庫455個，COG數(shù)據(jù)庫281個，KOG數(shù)據(jù)庫472個，KEGG數(shù)據(jù)庫372個，Pfam數(shù)據(jù)庫844個，Swiss-Prot數(shù)據(jù)庫807個，NR數(shù)據(jù)庫1 123個，其中NR數(shù)據(jù)庫注釋比高最高，達99.29%。

2.3.2GO功能注釋GO數(shù)據(jù)庫分為3大類，分別為B生物學過程(biological process)，C細胞組成(cellular component)和M分子功能(molecular function)，分別用來描述基因編碼產(chǎn)物所參加的生物過程、所具有的分子功能和所處的細胞環(huán)境等[22]。對黃秋葵花、果莢進行GO分類統(tǒng)計顯示，455個DEGs被歸到41個功能小類(表5)。生物學過程中DEGs“代謝過程”、“單一生物過程”和“細胞過程”3個功能小類占比最高；細胞組分過程中DEGs在“細胞組分”、“細胞”和“細胞器”3個功能小類占比最高；分子功能過程中DEGs在“催化活性”和“結合活性”2個功能小類占比最高。

2.3.3KOG功能注釋將注釋到KOG數(shù)據(jù)庫的472個DEGs進行直系同源分類，獲得23個功能分類，其中功能類別為R(一般功能預測)，獲得96個注釋結果，占比20.33%；O(次生代謝產(chǎn)物生物合成、轉運和代謝)，獲得67個注釋結果，占比14.20%；T(信號轉導機制)獲得52個注釋結果，占比11.01%；Q(次生代謝產(chǎn)物的合成、轉運和代謝)獲得44個注釋結果，占比9.32%。另外在G(碳水化合物轉運與代謝)、K(轉錄)也分別獲得44和39個注釋結果，分別占比9.32%和8.26%(表6)。

表3 組裝結果統(tǒng)計分析

表5 差異表達基因GO功能注釋

2.3.4KEGG功能注釋將DEGs通過KEGG數(shù)據(jù)庫比對(表7)，有372條基因得到注釋，分別富集在80條代謝通路，包括植物激素信號轉導(39個)、淀粉和蔗糖代謝(33個)、戊糖、葡萄糖醛酸轉換(25個)、氨基酸的生物合成(16個)、類苯基丙烷生物合成(14個)和碳代謝(14個)等22條代謝通路。本研究著重選擇類苯基丙烷生物合成代謝通路，尋找類黃酮代謝差異發(fā)生的關鍵基因。

2.4 類黃酮代謝關鍵差異基因分析

據(jù)KEGG類黃酮代謝途徑和表8數(shù)據(jù)分析，黃秋葵類黃酮合成代謝途徑中，查爾酮合酶(chalcone synthase, CHS)和查爾酮-黃烷酮異構酶(chalcone isomaerase, CHI)基因在花、果莢中均無差異表達；3-黃烷酮羥化酶(flavanone 3-hydroxylase, F3H)、二氫黃酮醇4-還原酶(Dihydroflavonol-4-reductase, DFR)在黃秋葵花中具有顯著表達優(yōu)勢；類黃酮3′,5′羧化酶(flavanone 3′,5′-hydroxylase, F3′5′H)、花青素還原酶(anthocyanidin reductase, ANR)、無色花色素還原酶基因(leucoanthocyanidin reductase, LAR)則在黃秋葵果莢中具有顯著表達優(yōu)勢；而花青素苷合成酶(anthocyanidin, ANS)、葡萄糖基轉移酶(glycosyl transferases, GT)基因則分別在花、果莢中具有顯著表達。黃秋葵花、果莢中類黃酮代謝表現(xiàn)為：P-香豆素-CoA(P-coumaroyl-CoA)和3分子丙二醛-CoA(malnoyl-CoA)，在CHS催化下，生成查爾酮(chalcone)，查爾酮又在CHI作用下形成柚皮素(naringenin, NAR)，此過程花與果莢均無顯著差異。黃秋葵花中NAR在F3H催化下生產(chǎn)二氫山奈酚(dihydokaempferol, DHK)，DHK在DFR作用下，生成無色天竺葵苷元，后在ANS(c83240)作用下生成花青素苷元，花青素苷元分別在GT(c44799/c82004)等轉移酶的作用下，生成穩(wěn)定的花青素苷；黃秋葵果莢則在F3′5′H作用下將NAR生成二氫楊梅素(dihydromyricetin, DHM)，后在FLS催化下，進入黃酮醇代謝途徑，部分DHM在DFR、ANS (c94125)、GT(c44799)作用下生成飛燕草素苷元(delphinidin)，飛燕草素苷元在ANR作用下，進入原花青素代謝途徑，無色飛燕草素苷元在LAR催化下也進入原花青素代謝途徑。另外，GT、ANS在黃秋葵花和果莢中分別有3個和2個差異表達拷貝，且表達量在不同組織中具有顯著互補性。

表6 差異表達基因KOG功能注釋

表7 差異表達基因KEGG功能注釋

表8 黃秋葵花、果莢類黃酮代謝關鍵基因差異表達情況

2.5 類黃酮代謝關鍵差異基因驗證分析

將表10中部分DEGs進行熒光實時定量(RT-PCR)，其中GT和ANS基因隨機選取多拷貝中1個，共7個。差異表達基因中c82420、c78393和c44799 3個基因表達量上調；c55984、c76561、c89149和c83240 4個基因表達量下調，以黃秋葵c84471基因為內參基因, 進行qRT-PCR驗證。由圖1可知，7個基因qRT-PCR分析得到的相對表達量與轉錄組表達譜分析趨勢完全一致, 但表達的變化大小存在一定差異，說明基因表達譜的分析結果基本可靠，其中DFR和ANR基因分別特異性在黃秋葵花、果莢中表達，且表達量差異極大。

1.c55984；2.c82420；3.c76561；4.c78393；5.c44799；6.c89149；7.c94125圖1 差異表達基因qRT-PCR相對表達Fig.1 The relative expression levels of DEGs by qRT-PCR

3 討 論

通過Illumina HiSeqTM2500測序技術構建黃秋葵花、果莢轉錄組數(shù)據(jù)庫，花和果莢分別獲得7.12 Gb和8.14 Gb 有效數(shù)據(jù)，堿基百分比(Q30)均達到91.0%以上。對Unigene進行功能注釋，在65 436條單基因序列中，共獲得39 245條Unigene的注釋結果，占單基因序列總數(shù)59.97%。通過DEGs分析，獲得差異基因1 336個，其中上調基因319個，下調基因1 017個。獲得功能注釋基因有1 131個，GO數(shù)據(jù)庫將455個DEGs歸到41個功能小類，代謝過程、催化活性、單一生物過程和細胞過程等4個功能小類占比最高；KOG數(shù)據(jù)庫將472個DEGs進行直系同源分類，獲得23個功能分類，其中與次生代謝直接相關過程O和Q類別獲得111個注釋結果，合計占比23.52%；通過KEGG數(shù)據(jù)庫比對，將372個DEGs注釋到80條代謝通路上，其中富集在類苯基丙烷生物合成途徑有14個DEGs，再通過類黃酮關鍵基因分析，共獲得10個關鍵差異基因。

植物類黃酮合成代謝是目前研究最清楚的此生代謝路徑之一[23]，該通路上存在2個重要基因群，即上游基因群(early biosynthetic genes, EBGs)與下游基因群(late biosynthetic genes, LBGs)。EBGs主要包括 CHS、 CHI、F3H、F3′H、F3′5′H等基因[24-25]；LBGs主要包括 DFR、 FLS、 ANS、ANR、 GT、?；D移酶(Acyl transferase, AT)和甲基轉移酶(Aethyl transferase, MT)等[26]。本研究發(fā)現(xiàn)，前體合成階段關鍵酶(CHS、CHI)均無顯著差異表達，這與CHS 和CHI基因編碼區(qū)和結構都十分保守，在不同科植物間、不同組織部位上均具有較高的保守性相一致[27]。黃秋葵花、果莢類黃酮代謝路徑中發(fā)現(xiàn)F3H、F3′5′H、DFR、ANR、ANS、LAR和GT存在差異表達，其中F3H、DFR在黃秋葵花中表現(xiàn)上調效應，F(xiàn)3′5′H、ANR、LAR在黃秋葵果莢中表現(xiàn)顯著上調效應，ANS、GT則分別在花和果莢中均有上調或下調效應。NAR作為類黃酮合成的關鍵分支點，黃秋葵花通過F3H、DFR、ANS(c83240)、GT(c82265/c82004)途徑，生成紫紅色花葵素-3-葡萄糖苷(pelagonidin-3-glucoside)[28]，這與黃秋葵花喉深紅色表現(xiàn)可能有關。黃秋葵果莢則通過F3′5′H、FLS催化下，將NAR生成楊梅素(Myricetin)[29]等；部分NAR在F3H、DFR、ANS、LAR作用下，進入花青素苷元(表焙兒茶素)和原花青素(沒食子兒茶素)合成途徑。初步推斷，黃秋葵花中富含花青素苷(葵素-3-葡萄糖苷)成分，果莢富含黃酮醇(楊梅素等)、花青素苷元(表焙兒茶素)原花青素(沒食子兒茶素等)等類黃酮物質。另外ANS和GT基因在黃秋葵花、果莢中均存在多拷貝和底物特異性現(xiàn)象[30]，推斷GT(c44799)、GT(c82004)在黃秋葵花中特異性的催化無色天竺葵素苷元，GT(c44799)在黃秋葵果莢中特異性的催化無色飛燕草苷元；ANS分別在花和果莢里高表達，ANS(c83240)可能與黃秋葵花瓣顏色(黃色)有關，ANS (c94125)可能與果莢顏色(濃綠色)有關，具體GT、ANS特異性底物種類、產(chǎn)物及最終呈色關系還需要進一步功能驗證。

利用高通量轉錄組測序技術，深入挖掘黃秋葵不同組織部位DEGs，從宏觀上理清黃秋葵花、果莢差異基因及調控機理等。近年來，黃秋葵營養(yǎng)成分研究在國內越來越受到重視，特別是黃酮提取[31]和功能分析[32]逐漸成為研究熱點，備受食品、醫(yī)藥等行業(yè)青睞。本研究豐富了黃秋葵花、果莢轉錄組水平生物數(shù)據(jù)信息，為類黃酮等關鍵基因克隆和功能驗證等提供遺傳基礎。