基于转录组的瑶药半枫荷根和叶的差异分析

周 云,韦玉丹,黄日涛

(广西医科大学 药学院,南宁 530021)

半枫荷始载于《岭南采药录》[1],以半枫荷同名同用途的植物有5科7属14种[2]。《广西壮族自治区瑶药材质量标准》记载的半枫荷(或半荷风、扁荷崩)是半枫荷(SemiliquidambarcathayensisChang)的干燥地上部分,隶属于蕈树科(Altingiaceae)半枫荷属(Semiliquidambar)[3-4]或金缕梅科(Hamamelidaceae)枫香树亚科(Subfam.Liquidambaroideae)半枫荷属[5],因富含黄酮类、萜类等活性成分而具有祛风除湿、活血舒筋等疗效,且治疗效果优于其他来源植物[6-9]。半枫荷的根、茎、树皮、枝叶甚至花蜜都具有祛风除湿、活血消肿的独特疗效[2,10]。现有的半枫荷化学成分研究表明其根和叶主要差异物质为黄酮类、萜类和多糖类等活性成分[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。基于拼接所得Unigenes序列信息采用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 项目。纵轴为该分类内基因个数(右)及其占被注释上基因总数的百分比(左)。纵轴(右)黑色数字注释到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 酰基转移酶;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]。因此,在本研究中半枫荷叶和根中有效成分的累积和合成、种类都存在差异,可能与合成途径中的差异基因上、下调表达有关。不同部位次生代谢产物的合成和积累是开发新药源的基础,值得深入探讨,后续工作可联合代谢组学数据针对半枫荷的不同组织部位的萜类、黄酮类、多酚类物质的合成和积累展开工作。