滇黃精轉錄組測序及類黃酮合成相關基因的分析

肖韻錚， 韓世明， 秦昭， 李春奇

(1.河南農業大學生命科學學院，河南 鄭州 450002；2.六盤水師范學院生物科學與技術學院，貴州 六盤水 553000)

黃精(Polygonati rhizoma)是著名的藥食兩用植物,包括黃精(PolygonatumsiricumRed.)、滇黃精(PolygonatumkingianumColl.et Hemsl.)等品種。滇黃精(PolygonatumkingianumColl.et Hemsl)為百合科黃精屬多年生草本植物，主要分布于中國西南地區,主產于云南省,以根莖入藥，是中國常用的中草藥[1]。滇黃精富含黃酮類化合物(flavonoids)。黃酮類化合物又稱類黃酮，是植物中重要的次生代謝產物，具有抗衰老、降血糖、改善血液循環等保健功能，具有很高的經濟價值和藥用價值[2-3]。目前，滇黃精和其次生代謝產物合成機理鮮有報道，尤其對滇黃精類黃酮合成的分子機理研究較少。因此，本試驗以滇黃精為研究對象進行轉錄水平分析，并對影響滇黃精類黃酮含量的關鍵基因進行了系統研究，為解析類黃酮生物合成分子機制奠定基礎。

轉錄組測序(RNA-Seq)技術是利用高通量測序技術對生物進行轉錄組分析，目前已廣泛應用于生物基因組基因表達的研究[4-7]，包括人參[8]、柴胡[9]以及甘草[10]等藥用植物已通過該技術進行了轉錄組測序。RNA-Seq技術不僅可以用于檢測和分析有基因組序列的生物體中的轉錄本，還可以用于組裝短片段，在無參考基因組的生物體中進行基因發現和基因表達分析[11]。鑒于此，本研究以滇黃精根、莖和葉為對象進行轉錄組測序，通過數據分析挖掘與類黃酮生物合成途徑相關的功能基因，為滇黃精功能基因鑒定和遺傳改良奠定研究基礎。

1 材料和方法

1.1 試驗材料

滇黃精植株由六盤水師范學院提供。選擇9株生長狀態良好且長勢相同的6個月盆栽苗。每3株為1組，共3組。第1組取根、莖和葉分別放入 3 個離心管中，其余組取材時選取發育階段相似、部位相同的根、莖和葉，共計9個樣本，3次生物學重復。取樣后立即置于液氮速凍保存于-80 ℃備用。

1.2 RNA提取、文庫的制備和測序

從滇黃精器官樣品(根、莖和葉)中提取總RNA，RNA提取使用植物總RNA提取試劑盒(Quick RNA isolation Kit，北京華越洋生物科技有限公司)。使用微量紫外分光光度計(Thermo Scientific, NANO, DROP2000)對總RNA進行定量，并通過測量A260/A280和A260/A230檢測總RNA的純度，將純化的RNA溶解在RNasefree水中并儲存在-80 ℃。測得樣品質量濃度不小于400 mg·L-1，OD260/OD280=1.8～2.2，28S∶18S≥1.5∶1.0。將總RNA反轉錄成cDNA，并通過PCR擴增富集得到最終的cDNA文庫，完成文庫制備。庫檢合格后，使用Illumina HiSeq 2500測序平臺(北京百邁客生物科技有限公司)進行測序。

1.3 質量控制、序列組裝和功能注釋

采用Illumina HiSeq 2500高通量測序平臺對cDNA文庫進行測序，去除reads中的測序接頭、引物序列以及低質量reads(質量值Q<30)，對于得到的高質量reads(clean reads)用Trinity軟件進行組裝。從獲得的基因功能信息中，將組裝得到的滇黃精的轉錄本數據與公共數據庫進行BLAST比對，并使用HMMER軟件與Pfam數據庫比對獲得注釋信息。公共數據庫為非冗余蛋白質數據庫(Nr)(ftp://ftp.ncbi.nih.gov/blast/db/)、蛋白質家族(Pfam)(http://pfam.xfam.org/)、直系同源蛋白質簇(KOG/COG/eggNOG)(http://www.ncbi.nlm.nih.gov/KOG/; http://www.ncbi.nlm.nih.gov/COG/; http://eggnogdb.embl.de/)、Swiss-Prot 蛋白序列數據庫(Swiss-Prot)(http://www.uniprot.org/)、基因本體數據庫(GO)(http://www.geneontology.org/)和京都基因與基因組百科全書(KEGG)(http://www.genome.jp/kegg/)。BLAST參數E-value≤10-5,HMMER軟件參數E-value≤10-10，通過序列相似性進行功能注釋。

1.4 差異表達基因的篩選

采用Bowtie軟件[12]將測序獲得的 reads與unigenes庫比對，并使用 RSEM軟件[13]根據比對結果對表達量水平進行估計。使用每百萬reads中比對到某一基因每千堿基長度的reads數目(fragments per kilobase of transcript per million mapped reads, FPKM)表示對應unigenes的表達豐度。然后用DESeq方法[14]對樣品進行差異表達分析，錯誤發現率(false discovery rate，FDR)<0.05且倍數差異(fold change，FC)>2作為篩選標準。GO功能富集采用topGO軟件(http://www.bioconductor.org/packages/release/bioc/html/topGO.html)進行分析，使用KOBAS2.0(http://kobas.cbi.pku.edu.cn/help.do)對差異表達基因(differentially expressed genes, DEGs)進行KEGG富集分析。

1.5 類黃酮含量的測定

以蘆丁標準品(生工生物工程上海股份有限公司，分析純)來配制蘆丁標準溶液(質量分數80%乙醇為溶劑)，以質量分數80%乙醇作為對照，利用紫外分光光度計測定510 nm 處的吸光值。以蘆丁質量濃度為橫坐標，吸光度值為縱坐標，繪制出標準曲線。

測定滇黃精不同器官中類黃酮的含量，參照劉清華等[15]的方法稍作修改。稱取滇黃精根、莖和葉樣品各0.2 g，按照料液比m(樣品)∶V(質量分數80%乙醇)=1∶15分別置于10 mL試管中，超聲提取50 min，后10 000 r·min-1離心10 min，取1 mL上清液用質量分數80%乙醇溶液稀釋到5 mL，吸取2 mL定容至10 mL，測定其在510 nm處吸光值，以質量分數80%乙醇作為對照，查詢標準曲線計算提取液中類黃酮質量濃度，再按下式計算不同器官中類黃酮含量：

M=(c×v×d)/m

(1)

式中：M表示類黃酮含量；c表示提取液中類黃酮質量濃度；v表示待測提取液總體積；d表示總稀釋倍數；m表示滇黃精樣品質量。

1.6 差異表達基因的qRT-PCR

總RNA提取使用植物總RNA提取試劑盒(Quick RNA isolation Kit，北京華越洋生物科技有限公司)。使用FSQ-101反轉錄試劑盒(東洋紡(上海)生物科技有限公司)將mRNA反轉錄為cDNA。依據轉錄組測序結果，篩選得到10個與類黃酮合成相關的差異基因，通過實時熒光定量PCR儀(LC480，瑞士羅氏集團)以及SYBR qPCR Mix試劑盒(北京艾德萊生物科技有限公司)進行qRT-PCR實驗。qPCR的引物由NCBI/Primer-BLAST界面進行設計，引物由生工生物工程(上海)股份有限公司合成，如表1所示。引物序列長度范圍18～25 bp，基因擴增產物長度限制100～350 bp。以18S rRNA為內參基因，上游引物為CTTGCGACCTGGAGTTATTGAT，下游引物為CTGCAACACTGCTCTTGCTC[16]。PCR反應體系和反應程序參照崔雯等[17]方法。

表1 熒光定量PCR引物

2 結果與分析

2.1 滇黃精轉錄組數據組裝分析結果

采用Illumina HiSeq 2500測序平臺技術得到滇黃精轉錄組數據。對獲得的原始數據進行統計，完成滇黃精樣品的轉錄組測序，去除raw data中的接頭序列及低質量reads獲得高質量的clean data 40.92 Gb，各樣品clean data均達到5.93 Gb，Q30堿基比例在91.94%以上，通過組裝后共產生了390 456條轉錄本，平均長度為770.38 bp，轉錄本N50長度為1 350 bp。對轉錄本序列進一步組裝拼接，共獲得了219 769條unigenes序列，平均長度為502.31 bp，N50長度為627 bp。其中，轉錄本在200～300 bp之間有139 323條，300～500 bp 之間有83 204條，大于2 000 bp的有32 976條，如圖1所示。

2.2 滇黃精unigenes序列功能注釋

將獲得的unigenes進行基因功能注釋以得到全面的基因功能信息，數據庫分別為Nr、Pfam、KOG/COG、Swiss-Prot、eggNOG、KEGG以及GO。經過BLAST程序對比，在219 769條滇黃精unigenes中，共有81 853條unigenes被注釋到，覆蓋率為37.2%。如表2所示，以數據庫Nr注釋成功率最高(34.7%)，其次為eggNOG(30.5%)、Pfam(20.7%)、KOG(19.7%)、GO(18.0%)、Swiss-Prot(16.1%)、KEGG(11.9%)和COG(11.2%)。

圖1 滇黃精轉錄本數據和unigenes的長度分布情況

2.3 滇黃精unigenes的GO和KOG功能分類

2.3.1 滇黃精unigenes的GO分類 根據unigenes的GO功能分析，滇黃精unigenes的GO注釋結果如圖2所示。在219 769條unigenes中，共有39 479條unigenes(18.0%)被注釋到GO條目中，并被分為3個主要類別，即生物過程(99 866條)、分子功能(48 515條)和細胞組分(78 215條)。在20個生物過程中，unigenes最多的為參與代謝過程(28 018條)和細胞過程(23 119條)，其次是單一有機體過程(17 934條)；在18個細胞組分中，unigenes最多的是細胞(18 352條)和細胞組分(18 352條)，其次是細胞器(14 245條)和細胞膜(8 512條)；在16個分子功能中，unigenes最多的是催化活性(21 212條)和綁定功能(19 999條)。

表2 Unigenes注釋的成功率統計

注： 1.胞外區； 2.膠原三聚體；3.細胞； 4.擬核； 5.細胞膜； 6.病毒粒子；7.細胞連接；8.胞外基質； 9.膜內腔；10.大分子復合物； 11.細胞器 12.細胞外基質部分；13.胞外區組件；14.細胞器組件；15.病毒粒子部分；16.細胞膜組件；17.細胞組件；18.共質體；19.蛋白質結合轉錄因子活性；20.核酸結合轉錄因子活性；21.催化活性； 22.受體活性； 23.鳥苷酸交換因子活性；24.結構分子活性；25.轉運器活性；26.綁定； 27.電子載體活性； 28.抗氧化活性； 29.金屬伴侶活性 30.酶調節劑活性； 31.蛋白標簽 32.翻譯調節器活性；33.營養儲層活動； 34.分子傳感器活性；35.生殖；36.細胞殺傷；37.免疫系統進程 ；38.代謝進程；39.細胞進程； 40.生殖進程；41.生物黏附 ； 42.信號；43.多細胞有機體進程；44.發育進程 ；45.生長； 46.運動力； 47.單一有機體進程；48.生物相；49.節律進程；50.應激反應； 51.定位；52.多有機體進程；53.生物調節； 54.細胞成分組織或生物合成。

2.3.2 滇黃精unigenes的KOG分類 將組裝得到的unigenes在KOG數據庫中進行比對，共有43 295條unigenes具有明顯的序列相似性且被注釋到25個分類中，如圖3所示。在滇黃精器官中，一般功能預測類的unigenes最多為10 244條，其次是翻譯后修飾為4 904條；而細胞運動是最少的類別為44條。

注: A: RNA加工和修飾; B: 染色質結構與動力; C: 能量的產生及轉換; D: 細胞周期控制, 細胞分裂, 染色體分區; E: 氨基酸運輸與代謝; F: 核苷酸運輸與代謝; G: 碳水化合物的運輸與代謝; H: 輔酶運輸與代謝; I: 脂質運輸與代謝; J: 翻譯, 核糖體結構和生物轉化; K: 轉錄; L: 復制, 重組和修復; M: 細胞壁, 細胞膜的發生; N: 細胞移動; O:蛋白翻譯后修飾, 蛋白質轉換,分子伴侶; P: 無機離子運輸與代謝; Q: 次生代謝產物生物合成, 運輸和分解代謝; R: 一般功能基因; S:功能未知; T: 信號轉導機制; U: 細胞內運輸, 分泌和膜泡運輸; V: 防御機制; W: 細胞外結構; Y: 細胞核結構; Z: 細胞骨架。

2.4 滇黃精差異基因的表達分析

滇黃精3個器官中共同表達的 unigenes有20 573條，根和莖、根和葉以及莖和葉中分別有930、 256、 1 171條unigenes特異表達，如圖4所示。以莖為對照，根中有2 361條DEGs，包括921條上調和1 440條下調；以葉為對照，根中有2 862條DEGs，包括1 340條上調和1 522條下調；以葉為對照，莖中有100條DEGs，包括27條上調和73條下調，如圖5所示。

圖4 滇黃精Unigenes的維恩圖

圖5 滇黃精差異表達基因數目

2.5 滇黃精差異基因功能注釋和富集分析

在7個數據庫中，本研究對差異表達基因進行功能注釋，并且繪制了不同器官之間(根和莖、根和葉、莖和葉)差異表達基因的GO功能分類圖，如表3和圖6所示，圖中橫坐標表示GO的3個模塊(生物學過程、 細胞組分、分子功能)，縱坐標左側表示差異表達基因數量，右側表示所有基因數量。利用KEGG代謝途徑數據庫對差異表達基因KEGG的注釋結果進行分類，如圖7所示，圖中橫坐標表示的是注釋到該通路下的基因數量及其與被注釋上的基因總數的比例，縱坐標表示的是KEGG代謝通路的名稱。根據注釋結果，共有26 259條(11.9%)unigenes注釋到KEGG數據庫中，和通路相關的為16 574條。在根和莖中，和通路相關的差異表達基因總數為483條，分別注釋到105條代謝通路；在根和葉中，和通路相關的差異表達基因總數為560條，分別注釋到108條代謝通路；在莖和葉中，和通路相關的差異表達基因總數為18條，分別注釋到22條通路中。富集通路主要集中在核糖體(ko03010)，另外是光合作用(ko00195)、碳代謝(ko01200)、苯丙烷類生物合成(ko00940)、淀粉和蔗糖代謝(ko00500)，這表明與類黃酮合成相關通路是差異基因主要富集通路之一。滇黃精unigenes代謝途徑可分成5大類，涉及unigenes數量最多的前3個代謝途徑分別為核糖體、碳代謝和光合作用。不同組間代謝通路比例最高，而有機系統通路比例最低。與核糖體相關的unigenes為2 254條，其次是碳代謝相關的unigenes為1 418條。對KEGG中的pathway通路進行分析，結合unigenes的注釋情況，發現滇黃精根中高表達基因的功能主要與自身能量儲藏代謝相關，莖和葉高表達基因功能主要與葉綠素代謝相關。

表3 滇黃精器官間注釋的差異表達基因數(DEGs)

注：1.代謝過程； 2.細胞過程； 3.單一有機體進程； 4.定位；5.應激反應； 6.生物調節；7.細胞成分組織或生物合成； 8.發育過程； 9.多細胞生物過程； 10.信號；11.生殖過程； 12.多組織過程； 13.生殖； 14.生長； 15.免疫系統過程； 16.生物附著；17.生物相； 18.節律性過程；19.運動；20.細胞殺傷性；21.細胞；22.細胞組件； 23.細胞器；24.細胞膜；25.大分子復合物；26.細胞器組件；27.細胞膜組件；28.膜內腔； 29.胞外區； 30.細胞連接；31.共質體；32.擬核； 33.胞外區組件； 34.細胞外基質；35.細胞外基質部分36.病毒粒子； 37.病毒粒子組件；38.膠原三聚體; 39.催化活性；40.結合；41.轉運活性； 42.結構分子活性； 43.電子載體活性；44.核酸結合轉錄因子活性；45.抗氧化活性；46.酶調節劑活性47.分子轉導活性；48.受體活性； 49.鳥苷酸交換因子活性 50.轉錄因子活性，蛋白結合；51.營養庫活性；52.金屬伴侶活性；53.蛋白標簽； 54.翻譯調控活性。橫坐標表示GO的3個模塊(生物學過程、 細胞組分、分子功能)。

圖7 滇黃精器官間差異表達基因的KEGG代謝途徑分析圖

2.6 滇黃精類黃酮合成基因的篩選

2.6.1 滇黃精類黃酮含量測定結果 滇黃精根、莖和葉中類黃酮含量測定結果如表4所示。滇黃精中不同器官的類黃酮的含量差異較大，從大到小排列的順序依次是：葉>莖>根，表現出明顯的器官特異性。結合根、莖和葉中差異基因的富集情況，發現莖與葉的差異基因數量較少且富集通路分布規律相同，而根與葉和莖的差異基因數量較多，差別較大，推測可能與類黃酮代謝通路中關鍵酶的合成有關。

表4 滇黃精根、莖、葉中類黃酮的含量

2.6.2 滇黃精類黃酮途徑上差異基因的表達 經過數據挖掘，篩選出與類黃酮生物合成相關的10個差異表達基因進行qRT-PCR，發現這些差異基因在器官中均有表達，但規律不盡一致，如圖8所示。CHI和CHS在葉中呈現高水平表達，而在根和莖中的表達量較低。DFR、F3’H、C3’H、LAR、CYP73A、CCOMT在根中呈現高水平表達，莖和葉中的表達量較低。ANR在莖中呈現高水平表達，但根和葉中表達量較低。FLS在根、莖和葉中均有表達，器官之間無顯著性差異。對滇黃精根、莖和葉中類黃酮的含量與類黃酮途徑中10個結構基因進行相關性分析，如表5所示。CHI和CHS與滇黃精器官中類黃酮含量呈顯著正相關，說明CHI和CHS是類黃酮生物合成途徑中的關鍵基因，其基因的表達量決定了類黃酮的合成量。

3 結論與討論

滇黃精是中國著名的藥食兩用植物，具有健脾補氣、降血糖和抗衰老等功效，黃酮類化合物是其主要的活性成分[18]。目前，關于滇黃精的研究主要集中在多糖和甾體皂苷類等途徑[19-22]，對類黃酮代謝途徑的研究卻鮮有報道。為了挖掘參與類黃酮生物合成的功能基因，本試驗通過高通量測序技術，對滇黃精根、莖和葉進行 RNA-Seq轉錄組測序并進行組裝，最終獲得219 769 條unigenes，平均長度為502.31 bp，N50為627 bp，Q30值分別為92.70%、 93.27%和92.54%，均大于90%；獲得390 456條轉錄本，其中在200～1 000 bp范圍內的有91.55%，32 976條轉錄本長度超過2.0 kb。這些結果表明滇黃精轉錄組序列組裝質量較高。將組裝后的序列在多個公共數據庫中進行注釋，共有81 853條unigenes至少在1個數據庫中獲得功能注釋，但仍有137 916條unigenes未被注釋，推測與序列片段較短及滇黃精缺少保守的核心序列和基因組信息有關。

本研究通過pathway分析，發現有26 259條(11.9%)unigenes注釋到KEGG數據庫中，和通路相關的所有基因總數為16 574條。在根和莖中，和通路相關的差異表達基因總數為483條，在根和葉中，和通路相關的差異表達基因總數為560條，在莖和葉中，和通路相關的差異表達基因總數為18條。結合篩選到的 unigenes 在相關數據庫中予以的功能注釋，進一步確保了所獲得基因的可靠性。通過不同器官比較，發現根中的差異基因數量較多，而大多數基因在莖和葉中的表達模式比根中的表達模式更穩定。對差異基因的注釋和富集情況進一步分析發現，滇黃精根中高表達基因的功能主要與自身能量儲藏代謝相關，莖和葉高表達基因功能主要與葉綠素代謝相關，這一研究結果與器官功能特性相符，說明根的主要功能與物質的吸收和能量儲藏有關，莖和葉的主要功能與物質運輸、光合作用和呼吸作用等有關。同時，鑒定到與類黃酮生物合成相關的unigenes有46條，以及大量與萜類物質合成相關的unigenes，這也為解析滇黃精黃酮類和萜類物質的合成途徑提供了數據支撐。

注：A～J為CHS、CHI、DFR、F3’H、C3’H、CCOMT、LAR、CYP73A、ANR和FLS的表達分析。不同小寫字母表示差異性顯著，顯著性水平P<0.05。

表5 滇黃精不同器官中類黃酮含量與類黃酮途徑差異表達基因相關性分析

為了進一步解析滇黃精類黃酮合成的分子機制，本研究對滇黃精不同器官中的類黃酮的總含量進行了測定，結果表明不同器官中類黃酮含量存在顯著差異，其中葉中類黃酮含量最高，莖次之，根最少，這與北京檸檬[23]、溪黃草[24]、棉花[25]等物種研究結果相一致，表現出明顯的器官特異性?；谵D錄組測序分析，本試驗共篩選出10個與滇黃精類黃酮生物合成相關的差異基因，通過qRT-PCR 檢測了其在滇黃精不同器官中的表達水平，并與類黃酮含量進行相關性分析，發現CHS和CHI基因表達水平與不同器官中類黃酮含量呈顯著正相關關系。在柑橘[23]、銀杏[26]等物種中，CHS基因的表達與類黃酮含量變化趨勢基本一致。在藥用植物鹽膚木[27]中，CHI基因的表達與不同器官類黃酮含量也呈良好的正相關關系。這些研究與本研究結果相符，表明CHI和CHS是滇黃精類黃酮生物合成的關鍵基因，且基因表達水平是導致類黃酮含量差異的重要因素。本研究為滇黃精提供了根、莖和葉的轉錄組數據，通過對滇黃精類黃酮代謝通路的分析，為挖掘類黃酮生物合成相關基因提供了重要理論依據，也為探究該植物次生代謝途徑和更多分子標記奠定基礎，有利于將滇黃精作為藥用植物資源的開發和利用。