基于转录组测序的襄麦冬块根甾体皂苷合成关联基因表达差异分析

隋 哲,程 旭,王海燕,余海忠

(1.湖北文理学院食品科学技术学院·化学工程学院,湖北 襄阳 441053;2.襄阳市农业科学院,湖北 襄阳 441057)

【研究意义】甾体皂苷是百合科山麦冬属植物襄麦冬(Liriopespicatavar.proliferaY. T. Ma)块根中的主要活性成分之一,具有抗肿瘤等多种药理活性,在临床医学、食品营养学等研究中日益受到关注[1-4]。【前人研究进展】目前,甾体皂苷的生产以中药材(如襄麦冬等)为原料,通过分离纯化工艺进行制备[5]。提高襄麦冬的产量,对实现甾体皂苷的稳定供应具有重要意义。然而,由于麦冬块根的生产周期长、皂苷含量低、大田的连作障碍等因素制约了襄麦冬产业的可持续发展[6]。研究表明,植物皂苷类成分的含量及组成,不仅随植物遗传背景、组织类型、生长年龄、生理状况、环境因子的改变而发生显著变化,还取决于皂苷合成酶及其相关基因的表达水平[7-10]。进一步研究显示,随着生育期的延长,襄麦冬块根中甾体皂苷的含量会逐步升高,在快速膨大期和采收期达到峰值[5, 11-13]。【本研究切入点】据此可认为,甾体皂苷合成酶基因在不同的襄麦冬块根生长发育阶段的表达具有差异性。若结合现有研究方法如高通量测序等,则能从分子水平对甾体皂苷的形成过程进行解析,挖掘甾体皂苷合成相关酶的基因及相关基因的表达特点,以更深入全面了解甾体皂苷的形成途径及影响因素[14-17]。【拟解决的关键问题】以襄麦冬基因组为研究对象,挖掘麦冬的甾体皂苷生物合成调控酶的基因,分析其表达情况。采集始见期和膨大期的襄麦冬块根样品,通过高通量测序比较分析2个转录组中的表达差异基因,并从分子水平上探究襄麦冬甾体皂苷的生物合成机理,从而为提高甾体皂苷的产量提供理论依据。

1.1 试验材料

实验材料采自汉江中游湖北襄阳欧庙襄麦冬良好农业规范(Good agricultural practices,GAP)种植基地,分别在襄麦冬块根始见期、块根膨大期选择生长良好的植株,剪取块根样品(图1)加入液氮速冻后,置于-80 ℃保存。

A:始见期样品;B:膨大期样品。A:Samples of initial onset stages(XMDSJQ);B:Samples of expansion stages(XMDPDQ).图1 襄麦冬块根始见期和膨大期的样品Fig.1 Samples from the root tubers of Liriope spicata var. prolifera at the initial onset stage or expansion stage

1.2 主要仪器和试剂

实时荧光定量PCR仪(Bio-Rad,美国);Agilent 2100生物分析仪(Agilent,美国)。主要试剂:Plant RNA Purification Reagent (Invitrogen,美国)、TruseqTM RNA sample prep Kit(Illumina,美国)、TBS380 Picogreen (Invitrogen,美国)、Certified Low Range Ultra Agarose(Bio-Rad,美国)、cBot Truseq PE Cluster Kit v3-cBot-HS(Illumina,美国)、Hiseq2000 Truseq SBS Kit v3-HS (200-cycles)(Illumina,美国)。

1.3 实验方法

1.3.1 总RNA的提取与Illumina测序 按照Invitrogen公司的Plant RNA Purification Reagent kit提取襄麦冬块根总RNA。对总RNA进行定量并电泳分析完整性后,置于-80 ℃保存备用。委托上海美吉生物医药科技有限公司,采用HiSeqTM2000平台对构建的2个cDNA文库进行Illumina测序。

1.3.2 生物信息学分析 使用Trinity软件(http://trinityrnaseq.github.io/)对质检处理的序列进行从头测序与拼接。使用RSEM(http://deweylab.biostat.wisc.edu/rsem)和edge R软件(http://www. bioconductor.org/packages/release/bioc/html/edgeR.html)对差异表达的转录本进行统计和分析。使用Trinity软件对拼接结果进行开放阅读框(Open Reading Frame,ORF)的预测,使用Blast软件完成序列功能注释。使用blast2go软件(http://www.blast2go.com/b2ghome)进行基因本体(Gene Ontology,GO)的功能分类。使用blastx 2.2.24+软件比对数据库STRING 9.0 (http://string-db.org/),完成基因产物的分类(Clusters of Orthologous Groups,COG)。使用blastx/blastp 2.2.24+软件进行代谢通路的分析。使用goatools软件(https://github.com/tanghaibao/goatools)对差异基因进行GO功能显著性富集分析。使用KOBAS软件(http://kobas.cbi.pku.edu.cn/home.do)对各组基因差异的表达进行代谢通路(Kyoto Encyclopedia of Genes and Genomes,KEGG)的显著性富集分析。

1.3.3 基因验证和分析 根据样本转录组数据的分析结果,从筛选获得的差异表达基因中选取8个关键基因,采用反转录荧光定量PCR技术(Real time quantitative PCR,RT-qPCR)验证Illumina测序结果,并对襄麦冬块根膨大期的差异表达基因进行表达谱分析,包括糖基转移酶保守结构域(Glycosyl transferase family 8 protein,GTF8,comp11269_c0_seq4)、环阿齐醇合成酶(Cycloartenol synthase,CS,comp1761_c0_seq5)、角鲨烯合成酶(Squalene synthase,SS,comp13813_c0_seq1)、单萜烯合成酶(Monoterpene synthases,MS,comp15786_c0 _seq2)、F26G(comp309_c0_seq3)、甾醇14α-去甲基化酶(Obtusifoliol 14alpha-demethylase,OF14D,comp444_c0_seq4)、甾酮C27单加氧酶2(Sterol 4-alpha-methyl-oxidase 2,SMO2,comp939_c0_seq8)和环丙基固醇异构酶(Cycloprop Yl sterol isomerase1,CPI1,comp4924_c0_seq9)。各个基因在转录组的表达情况见表1。利用Primer Express 3.0设计定量PCR特异性引物,并将翻译延长因子(Translation elongation factor,TEF)和(Translation β-tubulin,TUB)作为内参基因,引物如表2所示。

表1 襄麦冬转录组中8个关键基因的测序结果Table 1 The sequencing results of 8 key isogenes from the transcriptome of Liriope spicata var. prolifera

表2 RT-qPCR分析的引物Table 2 The primers designed for RT-qPCR

构建25 μL RT-qPCR反应体系:SYBR®Premix ExTaqTM II,12.5 μL;上游引物(10 μmol/L)1 μL;下游引物(10 μmol/L)1 μL;cDNA 2 μL;ddH2O 8.5 μL。RT-qPCR反应程序:95 ℃ 预变性2 min;95 ℃ 变性 30 s,58 ℃ 退火 30 s,72 ℃ 延伸 30 s,45 个循环;95 ℃ 30 s,55 ℃ 30 s,80次循环;55 ℃降温至0.5 ℃ 10 s。

2.1 测序数据的去杂及拼接分析

完成襄麦冬块根始见期和膨大期2个样品XMDPDQ、XMDSJQ的Iiilumia测序,共得到12.9 GB的原始数据。去除测序质量值低于25的reads、接头污染的reads、含N较多的reads以及经过以上步骤后长度小于25 nt的reads,分别获得高质量序列59 405 731和71 478 348条。对上述质控后的高质量序列进行从头测序并拼接,共得到56 067个genes和118 265个isogenes。isogene长度在351～15 626 bp,平均长度787 bp,其中以401～600 bp长度区间内分布的Contig数量最多,占36.17%,拼接结果如表3所示,isogene大小分布趋势如图2所示。

表3 襄麦冬样品的Illumina 测序拼接结果Table 3 Summary of the sequence assembly after Illumina sequencing of Liriope spicata var. prolifera

图2 襄麦冬样品的isogene大小分布趋势Fig.2 Size distribution of isogene of Liriope spicata var. prolifera

2.2 表达量统计及表达差异分析

首先,利用RSEM软件分别统计与isogenes和genes相匹配的reads数量,以此数量来计算isogenes和genes的表达量。随后,利用edgeR软件对统计的reads数量进行均一化,据此对基因的差异表达情况进行分析,确定P-value、FDR值均小于0.05的基因为差异表达基因。

从图3可知,在襄麦冬块根迅速膨大生育期XMDPDQ样品中,共有993条差异表达的genes和5562条差异表达的isogenes,上调表达的genes有139条,上调表达的isogenes有2157条。

A:在样本XMDPDQ与XMDSJQ差异表达的gene;B:在样本XMDPDQ与XMDSJQ中差异表达的isogene。红色点表示在两个样本中都表达的差异基因(P<0.05,FDR<0.05);黑色点表示在两个样本中都表达的非差异基因;黄色部分表示仅在一个样本中表达的基因。A:Differentially expressed genes from two stage samples;B:Differentially expressed isogenes from two stage samples.The red pots were differentially expressed genes between samples (P<0.05, FDR<0.05); The black pots were the genes expressed in all the samples; The yellow part indicated the genes only expressed in one sample.图3 样本XMDPDQ与XMDSJQ差异表达基因Fig.3 Visualization of differentially expressed genes between XMDPDQ and XMDSJQ

2.3 拼接序列的基因预测与注释

使用Trinity软件对拼接序列进行ORF预测,其中,预测出ORF的序列有62 507条,未预测出ORF的序列有55 758条。预测核酸cds序列的长度最短为150 bp,最长为15 357 bp,具体长度分布如图4所示。

图4 预测的cds长度分布Fig.4 Size distribution of predicted cds

同时,将上述预测出ORF的蛋白序列使用Blastp、未预测出ORF的DNA序列使用Blastx分别与NCBI NR、String、KEGG Genes数据库进行比对,完成所有拼接的注释,比对参数设置期望值E-value为1e-5。能注释到三大数据库中的拼接序列数量如表4所示。预测出ORF的蛋白序列有14 665条可注释到Genes库;有45 043条可注释到NR库,相似性介于35%～100%,其中335条达100%;有10 007条可注释到String库,相似性介于75%～100%,其中111条达100%。未预测出ORF的DNA序列有3802条可注释到Genes库;有13 515条可注释到NR库,相似性介于28%～100%,其中315条达100%;有1207条可注释到String库,相似性介于75%～100%,其中42条达100%。

表4 拼接序列与NR, String和KEGG数据库中公共数据的对比Table 4 BLAST result of isogenes against the data from three databases (NR, String and KEGG)

2.4 GO功能分类与富集显著性分析

比对GO数据库,完成42 388条拼接序列的GO注释。其中,归属于生物学过程的25个功能组基因中,细胞过程类群注释的基因最多(图5),有24 409条,占比57.59%;在细胞组分19个功能组中,细胞类群注释的基因最多,有28 131条,占比66.37%;在分子功能16个功能组中,催化活性类群注释的基因最多,有22 028条,占比51.97%。说明,本研究所建立并分析的转录组功能性基因,分别涉及襄麦冬生长及皂苷代谢的核心环节,因此可作为GO富集显著性分析的对象。

图5 GO功能的分类Fig.5 GO function classification of isogene

GO富集显著性分析发现,在2个转录组中存在5563个差异表达基因,分布在357个功能类群中,其中差异基因分布数量最多的功能类群为分子功能(2657个差异基因)、生物过程(2409个)、细胞组成(2342个)、细胞部分(1954个)、催化活性(1873个)、胞内组分(1802个)、细胞过程(1713个)、代谢过程(1684个)、结合(1669个)、胞内细胞器(1389个)。差异基因分布数量最少的功能类群为类胡萝卜素双加氧酶活性、9-顺式环氧类胡萝卜素双加氧酶活性以及表皮结构组成,均只有4个基因。说明,始见期与膨大期襄麦冬在生命活动能力上差异显著,暗示了皂苷合成与代谢水平间的关系:即膨大期的细胞功能更活跃,皂苷合成水平更高。

2.5 差异基因的COG/KOG注释

比对COG数据库来预测差异基因的可能功能,共有10 077条序列可注释到23个COG类群(图6-A)。基因数目最多的是一般功能预测基因,预测出2129功能基因条,占比21.13%,为最大类群,其后3位依次是转录(9.94%)、复制、重组及修复(9.21%)以及信号转导机制(8.47%),而细胞动力仅有0.12%,为最小类群。同时,比对KOG数据库发现,共有8187条序列可注释到25个KOG类群(图6-B)。基因数最多的是一般功能预测基因,预测基因1388条,占比16.95%,为最大类群。其后3位依次是转录后修饰与蛋白伴侣(11.10%)、信号转导机制(9.16%)、转录(6.88%),而细胞动力仅有0.06%,为最小类群。

A.COG分类比例的统计;B.KOG分类比例的统计。A.COG function classification of the two-staged samples; B.KOG function classification of the two-staged samples.图6 同源蛋白簇的分析Fig.6 Histogram presentation of COG classification and KOG classification

在信息存储和加工(Information storage and processing)类别中,归类到核酸复制、修复及转录过程(Replication, repair and transcription)的差异基因最多;在细胞过程与信号转导(Cellular processing and singaling)类别中,归类到蛋白修饰与折叠的差异基因最多;在细胞代谢(Metabolism)类别中,归类到碳水化合物运输与代谢、氨基酸运输与代谢的差异基因最多。说明,襄麦冬在不同生长期的生命活动差异显著,并由此影响麦冬皂苷的合成。

2.6 KEGG功能分类与富集显著性分析

采用KOBAS软件,对各组差异表达基因进行KEGG代谢通路显著性富集分析,发现共有26 626条拼接序列可映射到不同的KEGG途径中。KEGG富集显著性分析(Qvalue≤0.05)发现,在襄麦冬块根始见期、膨大期,共有1981个差异表达基因映射到227条KEGG途径,其中细胞周期通路(60个差异基因)、酵母细胞周期(52个差异基因)和淀粉及蔗糖代谢(51个差异基因)为3个主要的途径。

综上,上述转录组基因注释、分类及显著性分析结果均表明,襄麦冬在不同生长期的生命代谢活动差异显著,基础代谢差异对皂苷合成的相关代谢产生影响。

2.7 与甾体皂苷生物合成相关的基因

通过基因功能分类及富集显著性分析,初步梳理出43个可能参与甾体皂苷生物合成的差异表达基因,主要分布在倍半萜和三萜生物合成(Sesquiterpenoid and triterpenoid biosynthesis)(图7)、类固醇生物合成(Steroid biosynthesis)、油菜类固醇生物合成(Brassinosteroid biosynthesis)、淀粉和蔗糖代谢(Starch and sucrose metabolism)等KEGG途径中。

图7 参与甾体皂苷生物合成的基因在相关KEGG途径中的分布预测Fig.7 Predicted distribution of genes involved in Steroidal saponin biosynthesis in the related KEGG pathways

在上述初步筛选的43个基因中,有8个基因的表达在不同生长阶段的襄麦冬块根中具有显著差异。在襄麦冬块根膨大期,上调表达的有角鲨烯合成酶基因(comp13813_c0_seq1、comp13813_c0_seq2)、环阿屯醇合成酶基因(comp1761_c0_seq1、comp1761_c0_seq2、comp1761_c0_seq7、comp1761_c0_seq4、comp1761_c0_seq5)等,下调表达的有SMO1基因(comp669_c0_seq1)、SMO2基因(comp939_c0_seq2、comp939_c0_seq1)、呋甾皂苷26-O-β-葡萄糖苷酶基因(Furostanol glycoside 26-O-beta-glucosidase,F26G,comp309_c0_seq1、comp309_c0_seq2、comp309_c0_seq3、comp5904_c0_seq6)等。

2.8 RT-qPCR验证

通过RT-qPCR验证选择的8个基因,发现其在襄麦冬的不同时期表现出不同的表达量,且结果与转录组分析结果相符(图8)。

柱状图表示基因的相对表达量,误差线表示每次实验的差异。*表示在以MM为对照的0.05水平相对显著。其中,A～D表示在膨大期表达上调的基因:A表示Glycosyl transferase family 8 protein基因在两阶段襄麦冬块根中的表达差异;B表示Squalene synthase基因在两阶段襄麦冬块根中的表达差异;C表示Cycloartenol synthase基因在两阶段襄麦冬块根中的表达差异;D表示Monoterpene synthases基因在两阶段襄麦冬块根中的表达差异;E～H表示在块根膨大期下调的基因:E表示F26G基因在两阶段襄麦冬块根中的表达差异;F表示1,4-beta-D-xylan synthase基因在两阶段襄麦冬块根中的表达差异;G表示SMO2基因在两阶段襄麦冬块根中的表达差异;H表示CPI1基因在两阶段襄麦冬块根中的表达差异。The up-regulated genes extracted from swelled samples were shown in A-D, which were GTF8, SS, CS and MS, respectively;The down-regulated genes extracted from swelled samples were shown in E-H, which were gene F26G, OF14D, SMO2 and CPI1, respectively.图8 RT-qPCR验证XMDSJQ和XMDPDQ的基因表达结果Fig.8 The results of RT-qPCR between XMDSJQ and XMDPDQ

本研究通过对始见期襄麦冬块根和膨大期襄麦冬块根转录组的比较分析,发现在襄麦冬的不同生育期,参与调控代谢甾体皂苷生成的蛋白和酶基因存在差异性表达,初步梳理出43个可能直接参与甾体皂苷生物合成的差异表达基因。通过GO功能分类和COG/KOG功能标注,发现除一般性的功能基因外,差异基因主要参与核酸及蛋白质合成过程及氨基酸和碳水化合物运输与代谢等相关基因,一方面说明不同生长阶段的襄麦冬块根发育过程中蛋白质和核酸的合成情况有较大差异,另一方面说明它们对碳水化合物、氨基酸等营养物质的代谢与运输能力也有较大差异。即不同生长阶段的襄麦冬在合成代谢与分解代谢两方面的基因表达差异明显,由此推测不同阶段襄麦冬的代谢状况差异明显。

在KEGG富集分析中,对26 626条拼接序列进行了富集差异性分析,结果显示主要有3条代谢通路在分析中都得到了显著性富集,大量的细胞周期相关基因、酵母细胞周期相关基因和能量代谢相关基因参与襄麦冬块根组织的生长发育,调节襄麦冬块根组织细胞的增殖以及能量的供应,是襄麦冬组织中甾体皂苷表达的先决条件。在此基础上,针对甾体皂苷生物合成相关基因的分析结果显示,共有43个差异基因可能参与甾体皂苷的生物合成,主要映射到倍半萜和三萜生物合成、类固醇生物合成等多个次生代谢产物合成通路注释基因,说明不同生长阶段的襄麦冬块根组织次级代谢水平存在显著差异。

在初步筛选的43个基因中,发现8个差异性表达显著的基因,转录组分析结果与RT-qPCR的验证结果相符。相较于始见期样本,在膨大期样本中有4个基因表达上调。①糖基转移酶(Glycosyl transferase family 8 protein,GTF8):可催化环化三萜类化合物的糖基化,是皂苷生物合成的关键酶;②角鲨烯合成酶(Squalene synthase,SS):角鲨烯合成的关键酶,催化异戊二烯途径中的碳流向三萜生物合成,而且角鲨烯合成酶的活性与角鲨烯等三萜化合物的含量相关,而2,3-氧化角鲨烯是甾体皂苷与萜类化合物的分支和前体物质;③环阿屯醇合成酶(Cycloartenol synthase,CS):催化2,3-氧化角鲨烯成环形成甾体皂苷元的初始碳架-环阿屯醇,随后在酶的连续催化下形成胆固醇;④单帖合成酶(Monoterpene synthases,MS):单萜生物合成的关键酶。根据现有研究,上述4种酶均在甾体皂苷及皂苷元的生物合成中发挥重要作用[18-23],其表达上调,表明膨大期的襄麦冬中进行着更旺盛的甾体皂苷合成代谢。而在下调的基因中,包含消化呋甾皂苷C26葡萄糖水解酶的基因(Furostanol glyscoside 26-O-β-glucosidase,F26G)、裂解细胞壁并促进甾体皂苷释放的木聚糖酶的基因(1,4-beta-D-xylan synthase,OF14D)等[24-25]。除了挖掘并验证合成的甾体皂苷关键酶基因外,上调基因与下调基因的功能对比还表明,襄麦冬块根内皂苷的生物合成与其生长状态紧密关联,甾体皂苷主要合成于膨大期的襄麦冬块根内。

初步揭示与襄麦冬甾体皂苷生物合成有潜在关联的基因和代谢通路,为今后进一步研究襄麦冬甾体皂苷的主要关联基因及代谢通路的功能以阐明其合成机理奠定了基础。

猜你喜欢甾体块根差异基因非节肢类海洋无脊椎动物中蜕皮甾体及其功能研究海洋通报(2020年2期)2020-09-04基于RNA 测序研究人参二醇对大鼠心血管内皮细胞基因表达的影响 （正文见第26 页）心电与循环(2020年1期)2020-02-27海洋微生物来源的甾体化合物及其生物活性研究进展天然产物研究与开发(2018年7期)2018-08-21紫玉簪活性甾体皂苷的制备工艺研究天然产物研究与开发(2018年5期)2018-06-13德钦乌头块根化学成分的研究中成药(2018年2期)2018-05-09块根块茎类植物细胞悬浮培养技术与应用现代园艺(2017年19期)2018-01-19改良CLSI M38-A2应用于皮肤癣菌对甾体皂苷敏感性的测定中成药(2017年4期)2017-05-17紫檀芪处理对酿酒酵母基因组表达变化的影响江苏农业科学(2017年5期)2017-04-15木薯块根拔起的最大应力数值模拟及试验农业工程学报(2016年6期)2016-05-17药用植物珠子参新鲜块根DNA提取方法研究云南中医学院学报(2014年4期)2014-07-31