基于轉錄組的瑤藥半楓荷根和葉的差異分析

周 云,韋玉丹,黃日濤

(廣西醫科大學 藥學院,南寧 530021)

半楓荷始載于《嶺南采藥錄》[1],以半楓荷同名同用途的植物有5科7屬14種[2]?！稄V西壯族自治區瑤藥材質量標準》記載的半楓荷(或半荷風、扁荷崩)是半楓荷(SemiliquidambarcathayensisChang)的干燥地上部分,隸屬于蕈樹科(Altingiaceae)半楓荷屬(Semiliquidambar)[3-4]或金縷梅科(Hamamelidaceae)楓香樹亞科(Subfam.Liquidambaroideae)半楓荷屬[5],因富含黃酮類、萜類等活性成分而具有祛風除濕、活血舒筋等療效,且治療效果優于其他來源植物[6-9]。半楓荷的根、莖、樹皮、枝葉甚至花蜜都具有祛風除濕、活血消腫的獨特療效[2,10]?，F有的半楓荷化學成分研究表明其根和葉主要差異物質為黃酮類、萜類和多糖類等活性成分[11-13],而藥理學研究表明半楓荷中的黃酮類、三萜類、環烯醚萜類、單萜類、酚類和氨基酸等化合物是消炎、鎮痛、抗風濕的重要藥效物質[7-8]。由于半楓荷的藥理功能和來源植物多樣性,目前對于半楓荷的研究主要集中在組培[14-16]、分類系統[4,17-25]、化學成分和藥理作用[6,9,11,26-28]等方面,但對于半楓荷根及其葉的有效成分生物合成途徑及其相關酶基因表達調控的分子機制研究相對較少[29-30]。

在缺少植物參考基因組信息時,高通量轉錄組測序技術能夠針對植物的基因序列和轉錄本進行分析,目前廣泛應用于人參[31]、黃芪[32]、當歸[33]、防風[34]等。通過對藥用植物的不同組織或器官的轉錄組測序,分析其差異表達基因(differentially expressed genes,DEGs),從而為其次生代謝產物生物合成調控機制、種質資源保護和開發等研究提供理論基礎[35]。半楓荷葉片的轉錄組數據表明半楓荷活性成分基因調控網絡涵蓋了生物堿、萜類、黃酮類等多個關鍵次生代謝產物,推測半楓荷根相關的次生代謝產物更豐富[30],但仍需進一步分析。該研究選擇半楓荷幼苗的根和葉進行轉錄組分析,以期挖掘半楓荷根和葉的差異基因表達變化,為該植物的黃酮類、萜類關鍵化合物的次生代謝途徑解析和調控提供基礎數據。

1 材料和方法

1.1 RNA提取與測序

該研究的半楓荷植株為蕈樹科(Altingiaceae)半楓荷屬(Semiliquidambar)半楓荷(SemiliquidambarcathayensisChang),于2020年9月采自廣西桂林的廣西植物研究所,標本存于廣西醫科大學藥學院。將樣品的根和葉片放入液氮冷凍后寄往生工生物工程(上海)股份有限公司進行測序。利用Total RNA Extractor(Trizol)提取半楓荷根和葉的組織樣本RNA,并利用Qubit2.0 RNA 檢測試劑盒(Life)檢測RNA 濃度,以期獲得滿足后續文庫構建的總RNA 樣本,最終利用Illumina平臺進行高通量測序分析。

1.2 轉錄組序列的組裝和功能注釋

利用FastQC(http://www.bioinformatics.babraham.ac.uk/projects/fastqc/)對測序的原始數據(raw reads)進行質量評估,并利用Trimmomatic[36]進行數據過濾,從而獲得有效數據或高質量數據(clean reads)。利用Trinity[37]將原始數據進行混合拼接獲得Unigenes,參數min_kmer_cov 2,其余默認。采用Blast將單基因Unigenes序列與CDD[38]、KOG、COG[39]、NR、NT 等數據庫進行比對,得到其功能注釋信息。根據Unigenes序列與Swissprot、TrEMBL 的注釋結果得到GO功能注釋信息。利用KAAS[40]得到Unigenes序列的KEGG[41]注釋信息。利用HMMER 軟件[42]將Unigenes序列與PFAM 數據庫進行比對。

1.3 基因結構分析

根據Mapping結果使用BCFtools找出可能的SNP位點,抽提各樣本中的SNP/InDel,并根據以下條件對其進行過濾:(1)質量值大于20;(2)覆蓋度大于8?；谄唇铀肬nigenes序列信息采用MISA 進行SSR 分析,并統計相關信息。

1.4 表達差異分析及富集分析

使用Salmon計算基因表達量,利用TPM(transcripts per million)估算基因表達水平。利用DESeq2進行基因表達差異分析,差異表達基因(differentially expressed genes,DEGs)篩選條件設為:P＜0.05且差異倍數|αFC|＞2,隨后對差異表達基因做GO 功能富集分析、KOG 注釋和KEGG 注釋與分類。應用超幾何檢驗對KEGG 中每個通路進行富集分析,以q＜0.05為閾值,滿足此條件的通路定義為在差異表達基因中顯著富集。

2 結果與分析

2.1 測序數據與組裝結果

分別獲得半楓荷的葉和根44 482 222,65 361 772條原始數據,質控后獲得41 805 222條和62 109 246條有效數據,且各樣品的Q20均大于97.20%,Q30均大于90.21%,說明測序數據控制良好。組裝后獲得73 957條Unigenes序列,N50為1 335 bp,其中Unigenes序列最長為15 692 bp,最短為201 bp,平均長度688.15 bp,長度在1 000 bp以上的Unigenes有15 087個,占總數的20.40%。

2.2 半楓荷轉錄本的注釋

用BLAST 將所獲得的Unigenes比對到CDD、PFAM 等數據庫,并統計所注釋到的Unigenes數目和功能信息。結果顯示,共有39 822條Unigenes被注釋,占總數的53.84%,Unigenes 在CDD、PFAM、KEGG、eggnog、KOG、GO、NR、NT 數據庫得到注釋分別為12 733(17.22%),16 910(22.86%),9 879(13.36%),37 137(50.21%),15 987(21.62%),19 559(26.45%),35 136(47.51%),20 069(27.14%)條。以eggnog和NR數據庫注釋較多。共有37 137條Unigenes被注釋到24種eggnog分類中,其中未知功能有9 125個,是注釋最多的類群,核結構功能基因最少,有9 個。其他分類基因豐度參差不齊。Unigenes注釋同源基因的物種中,匹配度最高的同源物種是葡萄(Vitisvinifera),為4 875(6.59%),其次為藍果樹(Nyssasinensis),為3 895(5.27%)。

2.3 SNP和SSRs特征分析

檢測到半楓荷葉片的SNP 位點數為98 726,InDel為9 153;根的SNP 位點數為112 986,InDel為10 269。葉片的轉換(A-G,C-T)次數為31 967,顛換(A-C,A-T,C-G,G-T)次數為17 448;根的轉換(A-G,C-T)次數為36 764,顛換(A-C,A-T,C-G,G-T)次數為19 984;均表明轉換的次數多于顛換次數。10 830個SSRs中,二堿基重復SSRs豐度最高(5 922),其次為復雜重復類型(4 293),三、四、五和六堿基分別為2 366,173,44,88個。

2.4 基因差異表達分析

通過TPM 密度分析表明半楓荷根的轉錄組數據中TPM 在0～1 的基因數目多于葉的,而其TPM 在4～6.5 的基因數目少于葉的。從整體分析,半楓荷根和葉中基因表達存在明顯差異。17 515條基因發生顯著差異表達,其中上調表達基因數為9 421,占差異表達基因個數的53.79%;下調表達基因數為8 094,占差異表達基因數的46.21%(圖1)。

圖1 差異表達基因火山圖每個點代表1個轉錄本,越接近于原點的點表達量越低。其中紅色表示上調基因,綠色表示下調基因,黑色表示非差異基因。Fig.1 Volcano map of the differentially expressed genesEach point represents a transcript,and the closer the point is to the origin,the lower the expression is.Red represents up-regulated genes,green represents down-regulated genes,and black represents non-differential genes.

有59 378條差異表達基因成功進行了GO 分類(圖2)。在GO 分類的3個大類中有30 037條差異表達基因與生物學過程有關,占比為50.59%,其中富集最多的是細胞過程亞類,共比對到該功能類別總數的14.78%;與細胞成分有關的差異表達基因有22 751條,占比38.31%,富集最多的是細胞亞類,共比對到該功能類別總數的22.26%;有6 590條差異表達基因與分子功能相關,占比11.10%,其中富集最多的是催化活性亞類,共比對到該功能類別總數38.18%。

圖2 差異表達基因GO 注釋分類橫軸為GO 的二級分類,統計基因在生物過程,細胞組分,分子功能3個類別的各GO 項目?？v軸為該分類內基因個數(右)及其占被注釋上基因總數的百分比(左)?？v軸(右)黑色數字注釋到GO 的基因個數,灰色數字代表差異基因。深藍、深綠和深橘色代表根樣本的基因個數,淺色則是代表葉樣本的基因個數。Fig.2 GO annotationclassification ofthe differentially expressed genesThe horizontal axis represent the secondary classification of GO,including the GO items of the three categories of genes in biological process,cell component and molecular function.The vertical axis shows the number of genes in the classification(right) and their percentage to the total number of annotated genes (left).Black numbers (right) represent the number of annotate genes,and gray represent differential expressed genes.Dark blue,dark green,and dark orange represent the number of genes in the root of S.cathayensis,and light colors represent the number of genes in the leaf of S.cathayensis.

在半楓荷的葉和根中有626個注釋基因注釋到24 個KOG 分類中,其中有185 個DEGs,76 個DEGs表達量上調,109個DEGs表達量下調,其中一般功能預測有310個注釋基因,有20個表達量上調,70個表達量下調(表1)。

表1 差異表達基因的KOG 功能注釋Table 1 KOG annotation of the differentially expressed genes

半楓荷葉片和根轉錄組中有17 636個差異表達基因注釋到KEGG 數據庫中(圖3),其中注釋到細胞過程的1 883條(10.68%)、環境信息過程的2 140條(12.13%),遺傳信息過程的3 027條(17.16%),新陳代謝的7 132 條(40.44%),生物系統的3 454 條(19.58%)。進一步在KEGG 代謝通路富集分析中,有3 032個差異表達基因富集到290個KEGG通路中,其中上調表達的有2 099個,下調表達的有933個。差異表達基因在倍半萜類和三萜類、苯丙烷類、黃酮類、二萜類等通路中顯著富集。

圖3 半楓荷葉和根的轉錄組中差異基因的KEGG 分類Fig.3 KEGG pathways of the differentially expressed genes in leaves and roots of S.cathayensis

2.5 關鍵次生代謝產物合成途徑相關酶的鑒定

2.5.1 萜類生物合成通路差異表達基因分析

半楓荷葉片和根中與萜類合成相關的代謝通路有5條,包括倍半萜類和三萜類、萜類化合物骨架生物合成、二萜類、單萜類、其他萜類化合物生物合成,前兩條通路為顯著富集。5 條通路中共注釋到86個DEGs。

萜類化合物由甲戊酸(MAV)和甲基赤蘚糖醇磷酸(MEP)2個途徑合成(圖4)。其中有10條基因編碼MVA 途徑的3個關鍵酶,包括羥甲基戊二酰輔酶A 合酶(HMGS,)、羥甲基戊二酰輔酶A 還原酶(HMGR)、甲羥戊酸焦磷酸脫羧酶(MVD);有7條基因編碼MEP途徑的5個關鍵酶,包括脫氧木酮糖-5-磷酸合酶(DXS)、4-二磷酸胞苷-2-C-甲基-D-赤蘚醇合酶(CMS)、2-甲基赤蘚糖-2,4-環二磷酸合酶(MCS)、1-羥基-2-甲基-2-(E)-丁烯基4-二磷酸合酶(HDS)、1-羥基-2-甲基-2-(E)-丁烯基-4-二磷酸還原酶(HDR)。

圖4 推測的半楓荷葉和根萜類生物合成途徑AACT,acetyl-CoA acetyltransferase.HMGS,hydroxymethylglutaryl-CoA synthase.HMGR,hydroxymethylglutaryl-CoA reductase.MK,meva-lonate kinase.MVD,mevalonate diphosphate decarboxylase.DXS,1-deoxy-D-xylulose 5-phosphate synthase.DXR,1-deoxy-D-xylulose 5-phosphate reduc-toisomerase.CMS,4-diphosphocytidyl-2-C-methyl-D-erythritol synthase.CMK,4-diphosphocytidyl-2-C-methyl-D-eryth-ritol kinase.MCS,(E)-4-hydroxy-3-methyl-but-2-enyl diphosphate reductase.HDS,(E)-4-hydroxy-3-methylbut-2-enyl-diphos-phate synthase.HDR,(E)-4-hydroxy-3-methyl-but-2-enyl diphosphate reductase.Fig.4 Main terpenoid biosynthesis pathway compared with leaves and roots of S.cathayensisAACT.乙酰輔酶A ?；D移酶;HMGS.羥甲基戊二酰輔酶A 合成酶;HMGR.羥甲基戊二酰輔酶A 還原酶;MK.MVA 激酶;MVD.甲羥戊酸焦磷酸脫羧酶;DXS.脫氧木酮糖-5-磷酸合酶;DXR.脫氧木酮糖磷酸鹽還原異構酶;CMS.4-二磷酸胞苷-2-C-甲基-D-赤蘚醇合成酶;CMK.4-二磷酸胞苷-2-C-甲基赤蘚醇激酶;MCS.2-甲基赤蘚糖-2,4-環二磷酸合酶;HDS.1-羥基-2-甲基-2-(E)-丁烯基-4-二磷酸合酶;HDR.1-羥基-2-甲基-2-(E)-丁烯基-4-二磷酸還原酶。

2.5.2 苯丙烷類和黃酮類生物合成通路差異表達基因分析

參考半楓荷化學成分的相關文獻,根據本研究KEGG 途徑分析結果,對半楓荷中主要黃酮類化合物的生物合成途徑做出推測(圖5)。

苯丙烷類生物合成途徑注釋到81 個DEGs。與黃酮類化合物合成有關的通路有3個,包括黃酮類生物合成途徑(ko00941)、異黃酮類生物合成(ko00943),黃酮和黃酮醇生物合成(ko00944)注釋的30個DEGs。

從苯丙烷類生物合成和黃酮類生物合成途徑中共篩選到15個差異表達的關鍵酶,分別為苯丙氨酸解氨酶(PAL)、β-葡萄糖苷酶、4-香豆酸輔酶A 連接酶(4CL)、查爾酮合酶(CHS)、二氫黃酮醇-4-還原酶(DFR)、黃酮醇合酶(FLS)、反式-肉桂酸4-單加氧酶(C4H)、莽草酸O-羥基肉桂酰轉移酶(HCT)、5-O-(4-香豆?；?-D-喹酸酯-3′-單加氧酶(C3′H)、類黃酮3′,5′-羥化酶、肉桂醇脫氫酶(CAD)、咖啡酸3-O-甲基轉移酶(COMT)、過氧化物酶、咖啡酰輔酶A-O-甲基轉移酶(CCoAOMT)、阿魏酸-5-羥化酶。其中DFR、CHS均為上調表達,其余既有上調也有下調。

3 討論

研究通過轉錄組初步分析了半楓荷葉片和根的基因表達特征,共有39 822條Unigenes被注釋,占總數的53.84%,表明半楓荷中有大量基因信息,還有待進一步挖掘和分析。已有研究利用半楓荷葉片篩選SSR 位點,并用于半楓荷遺傳多樣性及遺傳結構的研究[43,33]中。該研究獲得的10 830 個SSRs中,二堿基重復SSRs豐度最高,占比54.68%。大多數植物的SSR 位點主要以二、三核苷酸重復類型為主,所檢測SSR 標記可為后續半楓荷遺傳圖譜的構建、分子育種等提供理論基礎。

目前半楓荷轉錄組數據都是葉片測序,根中的主要化學成分有三萜類、生物堿、黃酮類、苯丙烷類、鞣花酸衍生物、苯甲酸類、類固醇、脂肪酸、脂肪鏈烴及甾體類[9,11-12],葉片中的主要化學成分有多酚類、硬脂酸、鞣酸、β-谷甾醇等[6]。已有研究針對半楓荷葉的轉錄組數據發現有苯丙素類、黃酮類生物合成通路,多個萜類生物合成途徑關鍵基因[29-30],說明其次生代謝生物合成途徑的復雜性?；诎霔骱扇~片和根的比較轉錄組分析研究發現半楓荷不同組織部位的差異基因表達存在顯著差異,通過KEGG 數據庫比對發現半楓荷葉片與根的差異表達基因集中在萜類、胺基酸、苯丙烷類、黃酮類、脂肪酸等生物合成相關的次生代謝通路。在半楓荷葉和根的比較中,萜類化合物中的單萜、倍半萜和三萜的生物合成途徑中,下調的關鍵基因較多,可能導致半楓荷根的三萜類等有效化合物較少,在萜類化合物合成的MVA 途徑中關鍵酶HMGS、HMGR 上調和下調表達均有,而MVD 為上調表達,MEP 途徑中關鍵酶CMS、MCS、HDS和HDR都為下調表達,DXS上調和下調表達皆有,由此推測半楓荷根與葉相比,CMS、MCS、HDS和HDR可能負向調控下游產物的產生,而HMGS、HMGR和DXS這些關鍵酶既可能正向調控,也可能負向參與下游產物的生成。HMGR 和DXS分別是MVA 途徑和MEP 途徑的關鍵限速酶[45-46],這都可能影響半楓荷不同組織萜類化合物的積累模式。

半楓荷的根與葉相比,參與苯丙烷類和黃酮類化合物合成途徑的關鍵酶CHS、DFR 均為上調表達,CHS是色原酮等黃酮類物質合成的限速酶[47],DFR 是花青素合成的關鍵酶,主要可以修飾花色,其過表達會導致植物花色加深[48],在2年生和3年生苦參中,CHS和DFR 在根中的表達均高于莖中,導致根中有效成分更高[49],而CHS下調表達則會導致防風中色原酮合成速率較慢[34],由此推測半楓荷根中的黃酮類有效成分會高于葉片。其余關鍵酶則上調和下調表達均有,尤其是HCT、β-葡萄糖苷酶和過氧化物酶注釋的差異表達基因較多,這些都跟木質素、纖維素等化合物的生成相關。HCT 是木質素生物合成路徑中的關鍵酶,在苯丙烷3-羥基化上游和下游的苯丙氨酸途徑中具有雙重調節作用[50]。β-葡萄糖苷酶是纖維素酶系的重要成員,參與生物體的糖代謝[51]。過氧化物酶則可以影響木質素的合成、種類和積累[52]。因此,在本研究中半楓荷葉和根中有效成分的累積和合成、種類都存在差異,可能與合成途徑中的差異基因上、下調表達有關。不同部位次生代謝產物的合成和積累是開發新藥源的基礎,值得深入探討,后續工作可聯合代謝組學數據針對半楓荷的不同組織部位的萜類、黃酮類、多酚類物質的合成和積累展開工作。