基于转录组测序的细风轮花青素合成途径及关键酶基因分析

赵历强 单春苗 张声祥 施圆圆 马克龙 吴家文

（1. 安徽中医药大学研究生院，合肥 230012；2. 安徽中医药大学科研实验中心、新安医学教育部重点实验室，合肥 230038；3. 安徽中医药大学中西医结合学院，合肥 230012；4. 安徽道地中药材品质提升协同创新中心，合肥 230012；5. 安徽省中医药科学院，合肥 230012）

细风轮（Clinopodium gracile（Benth.）Matsum）为唇形科（Labiatae）风轮菜属（Clinopodium）植物，又名塔花、瘦风轮、野薄荷等，主要分布于广西、福建、广东、贵州、四川等省份。细风轮的主要化学成分为三萜及其皂苷、黄酮、挥发油、苯丙素类、甾体等［1］。细风轮具有清热、解毒、活血的作用，临床上主要用于治疗白喉、咽喉肿痛、乳腺炎和过敏性皮炎等［2］疾病。

近年来已经成功从细风轮中分离出来黄酮、类黄酮和皂苷等化合物［3］。类黄酮化合物在人体内可表现出明显的生物活性，例如清除自由基、抗菌、抗溃疡、消炎、抗高脂血症、降血压及抗癌、抗病毒［4～7］等。类黄酮化合物中的花青素（Anthocyanidin）是一种天然的抗氧化剂，在蓝莓［8］、枸杞［9］和葡萄［10］等植物中含量较高，其在人体中能够起到抗氧化［11］、降血脂［12］、消除炎症［13］以及抗肿瘤［14］等作用。田野等［15］发现了紫嫣茶中的花青素具有抗结肠癌和乳腺癌的作用；刘丽敏等［16］报道了越橘中花青素能起到维护眼球健康、降低近视发生率的作用。

花青素是苯丙氨酸通过一系列的酶促反应后进一步经过酰基化、甲基化、糖基化等修饰而生成，其储存在植物的液泡中［17］。DFR（二氢黄酮醇还原酶）是参与花青素生物合成的限速酶［18］，可以催化二氢黄酮醇生成花色素［19］，大量研究表明DFR 在花青素的生物合成中起到了重要的调控作用［20～22］。

高通量测序技术可以在没有参考基因组的情况下对中草药转录组进行全面分析，确定其转录本 和 功 能 基 因 序 列［23～25］。目 前 黑 果 枸 杞［26］、麻黄［27］、紫背天癸［28］等中药已经完成了其转录组数据库的构建，并对花青素生物合成途径进行了研究，但目前尚未有细风轮花青素生物合成途径研究的报道。本研究以细风轮花、叶、茎、根4个组织为样品，首次采用高通量测序技术对细风轮进行转录组分析，从基因水平分析细风轮药用活性成分花青素的生物合成途径，为将来利用基因工程或代谢工程技术提高花青素的产量、进一步开发利用花青素奠定了实验基础。

1 材料与方法

1.1 实验材料

本实验所用细风轮植株采集于安徽中医药大学的药物园，由安徽中医药大学药用植物分类学老师杨青山鉴定。新鲜的细风轮植株洗净后，分离根、茎、叶、花各器官，用灭菌蒸馏水多次冲洗后滤纸擦干，分别收集放入50 mL 的离心管中，并迅速置于液氮中冻存［29］。

1.2 细风轮总RNA的提取

将所需器具高温灭菌后，把细风轮的组织分别放入研钵中，边加液氮边研磨，粉末充分混匀后放入离心管中，离心后取上清液，用RNA 提取试剂盒（北京天根生化科技有限公司）对RNA 进行纯化，用Agilent 2100 生物分析仪检测其浓度及完整性后利用BGISEQ-500 平台进行转录组测序。

1.3 数据过滤及组装

将转录组测序所得到的原始读序（raw reads）进行数据过滤，去除低质量、未知碱基N 含量大于5%和包含接头的读序（reads）后得到干净读序（clean reads）。使用Trinity 软件对干净读序进行重头组装，然后使用Tgicl 将组装的转录本进行聚类去冗余，最后得到Unigene。

1.4 Unigene的分类及注释

利用NCBI Blastx 软件将Unigene 分别注释到NT（非冗余核酸序列数据库）、KOG（NCBI 直系同源家族蛋白数据库）、NR（NCBI非冗余蛋白序列数据库）、SwissProt（高质量非冗余的蛋白数据库）、KEGG（京都基因和基因组百科全书）以及GO（基因本位论数据库）和Pfam（多序列比对和隐马尔科夫模型构建的蛋白家族数据库）。

1.5 结构特征分析

利用软件Translate tool 确定细风轮DFR 的开放读码框（ORF）序列，利用MEGA 5 软件和CLUSTALX 1.83 软件对筛选出的DFR 的氨基酸序列进行比对，寻找其保守的活性位点并构建其系统 发 育 进 化 树，最 后 采 用Swiss-Model（https：//swissmodel.expasy.org/）模拟出DFR的三级结构，再利用Pymol软件描绘出其三级结构［30］。

2 结果和分析

2.1 细风轮总代谢和次生代谢概述

本研究利用BGISEQ-500 平台一共测了43.81 Gb数据，经过从头组装后得到128 856个Unigene，利用Mapman 软件对细风轮FPKM（Fragments per Kilobase Million）值大于1 的Unigene 在总代谢和次生代谢途径中进行富集，发现在总代谢中基因主要富集在脂质、核苷酸代谢途径（见图1）；而在次生代谢中基因则主要富集在黄酮、木质素木酚素以及苯丙素类代谢途径（见图2）。

图2 细风轮中FPKM＞1的Unigene在次生代谢中的富集及表达水平Fig.2 Enrichment and expression level of unigenes（FPKM＞1）from C.gracilein secondary metabolism

图1 细风轮中FPKM＞1的Unigene在总代谢中的富集及表达水平Fig.1 Enrichment and expression level of unigenes（FPKM＞1）from C.gracile in total metabolism

2.2 细风轮中次生代谢途径概述

通过KEGG 数据库分析，我们对细风轮可能参与各种代谢途径的Unigene 进行了统计，结果表明2 312 个Unigene 参与了次生代谢产物的生物合成，包括14 个代谢途径（见表1），其中苯丙素生物合成途径（ko00940）含1 406 个Unigene，类黄酮生物合成途径（ko00941）含455 个Unigene，吲哚生物碱生物合成途径（ko00901）含187个Unigene，花青素的生物合成途径（ko00942）含72 个Unigene。

表1 参与细风轮次生代谢途径的Unigene数目Table 1 Number of unigenes involved in secondary metabolism pathways in C.gracile

表2 花青素合成关键酶Table 2 Key enzymes for the anthocyanin

2.3 鉴定与花青素生物合成相关的Unigene

通过分析细风轮转录组数据我们发现了6 个与花青素生物合成相关的酶（见表2），分别由不同数目的Unigene 编码。二氢黄酮醇还原酶（bifunctional dihydroflavonol 4-reductase，DFR）由7个unigene 编码、黄烷酮3 羟化酶（flavanone-3-dioxygenas，F3H）由7个Unigene编码、黄烷酮合成酶（chalcone synthase，CHS）由5 个Unigene 编码、苯 丙 氨 酸 解 氨 酶（phenylalanine ammonia-lyase，PAL）由7 个Unigene编码、4-香豆酰辅酶A 连接酶（4-coumarate-CoA ligase，4CL）由10 个Unigene 编码、查耳酮异构酶（chalcone isomerase，CHI）由4个Unigene编码。

图3 细风轮生物合成代谢的途径分类Figure 3 Pathway classifications for biosynthesis metabolism in C.gracile

PAL 是花青素合成起始阶段的关键酶，催化苯丙氨酸生成反式肉桂酸，4CL催化4-香豆酰生成香豆酰CoA，CHS催化4-香豆酰CoA生成柚皮素查尔酮，CHI 催化柚皮素查尔酮生成柚皮素，F3H 则催化柚皮素生成二氢山萘酚，DFR 催化二氢山萘酚生成白天竺葵苷元。苯丙氨酸经过这一系列酶促反应最终生成了3种类型花青素即天竺葵色素、花翠素和矢车菊素。在细分轮花青素的生物合成途径中，DFR 基因在叶和花中的表达量明显高于根和茎，PAL 和F3H 基因在根和茎中的表达量明显高于叶和花，而CHI 基因主要在花中高表达（见图4）。

图4 花青素生物合成途径及关键酶在组织中表达水平分析聚类热图中列F、L、R、S分别表示花、叶、根、茎；行表示Unigene；颜色表示表达水平Fig.4 Anthocyanin biosynthesis pathway and Cluster heat map of the expression of key enzymes in tissuesF，L，R and S represent flower，leaf，root，and stem samples，respectively.Red and green represent high and low expression levels，respectively

2.4 细风轮花青素生物合成中基因差异表达分析

将细风轮叶和花组织进行比较，参与花青素生物合成的Unigene 有41 个差异表达，其中18 个在叶中上调、23个在叶中下调；将细风轮叶和根组织进行比较，参与花青素生物合成的Unigene有49个差异表达，其中22 个在叶中上调，27 个在叶中下调；将细风轮茎和花组织进行比较，参与花青素生物合成的Unigene 有41 个差异表达，其中24 个在茎中上调，17个在茎中下调（见图5）。

2.5 细风轮花青素生物合成相关Unigene 的GO分析

GO 分类（见图6）显示有69 个参与花青素生物合成的Unigene 被成功注释，并将其分为21 类，其中有9 个Unigene 参与“代谢过程”、9 个Unigene与“催化活性”相关、8 个Unigene 具有“氧化还原酶”活性、8个Unigene具“结合金属离子”特性。

2.6 参与花青素生物合成的苯丙氨酸解氨酶结构特征

图5 细分轮4个组织中参与花青素生物合成的Unigene差异基因表达A. 在细分轮4种组织中表达的Unigenes聚类热图；B. 细分轮两两组织比较基因表达上调和下调的数目Fig.5 Differential gene expression of anthocyanin in four tissuesA.Heatmap of unigenes of C.gracile involved in anthocyanin biosynthesis；B.The number of unigenes up-regulated or down-regulated in the specified samples are shown

图6 细分轮花青素生物合成相关的Unigene GO富集Fig.6 GO enrichment of unigenes involved in anthocyanin biosynthesis in C.gracile

二氢黄酮醇还原酶（DFR）是类黄酮生物合成途径的关键酶，催化二羟基黄酮醇转变为花色素，序列比对分析表明细风轮中编码DFR 酶的序列同源性较高，其中CL3369-1 与CL3369-2 这两条氨基酸序列一致性达87.86%，它们通过Swiss-Model 构建的空间结构也非常相似。我们选择与模板蛋白同源性最高的DFR（CL-3369-1）描绘其二级结构和空间结构。DFR 二级结构主要由α-螺旋和β-折叠构成，共含有15 个α-螺旋，11 个β-折叠，α-螺旋包裹着β-折叠，另外还有4 个TT 结构和2 个η 结构（见图7）。

DFR 的空间结构模型为同源二聚体，呈现为致密的球状结构，DFR 酶的NAD+特异性结合位点位于第一个α-螺旋与第一个β-折叠之间的loop 区域V9-Y29（VTGGSGYIASFLISHLLERGY）（见图7）。

2.7 DFR同源关系及系统进化树分析

通过NCBI BLAST 比对不同物种间DFR 氨基酸序列（见表3），结果表明细风轮与一串红、芝麻、风铃木等植物的DFR 蛋白序列一致性均达到65%以上。DFR 系统进化树显示，在已知数据库中细风轮与一串红（Salvia splendens）遗传距离最短，亲缘关系最近，这与它们同属于唇形科植物的是一致的（见图8）。

图7 细风轮DFR的二级结构和三级结构示意图A～B.细风轮DFR的三级结构（B紫色区域为NAD+结合位点）；C.细轮DFR的二级结构（黑色下划线显示NAD+结合位点）Fig.7 The secondary structure and tertiary structure of DFR in C.gracileA.The tertiary structure of DFR in C.gracile with helices colored in red，strands colored in yellow；B.The purple region shows the binding site for NAD+；C.The secondary structure of DFR in C.gracile（The black underline shows the NAD+binding site）

图8 细风轮DFR系统发育进化树0.05表示遗传距离；黑色三角突出显示细风轮Fig.8 Phylogenetic analysis of DFR in C.gracile0.05 shows genetic distance；The black triangle highlights the C.gracile

表3 不同物种间DFR比对结果Table 3 Homologous alignment of DFR protein

3 讨论

本研究利用BGISEQ-500 平台对细风轮花、叶、茎、根四个组织进行了转录组分析。利用Mapman软件发现FPKM＞1的Unigene主要富集在总代谢中的脂质、细胞壁代谢以及次生代谢中的黄酮、苯丙素类代谢途径。KEGG 分析表明，有2 312 个Unigene参与到次生代谢通路的14条途径中，其中主要的代谢途径是苯丙素、黄酮和花青素代谢途径。GO 分析显示参与花青素生物合成的Unigene主要富集在“代谢过程”和“催化活性”两类功能。

前人研究证明PAL、CHS、CHI、DFR、F3H 和4CL 均是花青素合成路径中的关键酶，它们在控制花青素生物合成方面发挥重要作用［31～35］。虎娟等利用转录组测序技术从黑果枸杞果皮中鉴定出与花青素合成相关的关键基因CHS、CHI、DFR、F3′H（类黄酮3′羟化酶）、ANS（花青素合酶）和F3′5′H（类黄酮3′5′羟化酶），进一步对他们的核苷酸序列、蛋白质一级序列及空间结构特征进行了分析，并利用RT-PCR 技术克隆了这些酶基因，验证了它们是真实存在的［28］；蒋会兵等通过转录组测序分析紫芽茶树的PAL、CHS、ANS、和UFGT（类黄酮糖基转移酶）基因在不同种类中的差异表达，认为它们可能是紫芽茶树中花青素合成途径的关键酶基因［36］；Liu 等通过对不同花色毛茛的转录组测序鉴定出其花青素合成关键酶及其在不同花色毛茛中的表达量，发现关键酶基因中的F3H、F3′H、3GT（类黄酮-3-O 糖苷转移酶）、5GT（类黄酮-5-O 糖苷转移酶）和FMT2（类黄酮-O-甲基转移酶2基因）在颜色为红色品种中高表达，而F3′5′H和3MAT（花色苷3-O-葡萄糖苷-6-O-丙二酰基转移酶基因）在紫色品种中高表达［37］，说明不同花色毛茛中的高表达关键酶基因种类不同。

在细风轮花青素合成途径中，40个Unigene编码6 个关键酶（PAL、CHS、CHI、DFR、F3H、4CL），其中DFR 基因在叶和花的表达量明显高于根和茎，CHS 和CHI 基因在花和根中的表达量高于叶和茎，而PAL 和F3H 基因主要在根和茎中高表达，4CL 在茎中表达最低。细风轮中花青素在各组织中的含量可能是这些关键基因协同作用的结果，具体哪些基因发挥最关键的作用，有待于我们进一步研究。

二氢黄酮醇还原酶（DFR）是花青素生物合成下游途径中的第一个关键酶［20］，细风轮中DFR 酶是由7 个Unigene 编码的，它们氨基酸序列的一致性为42.08%，我们选择与蛋白质结构数据库（PDB）中已知结构的葡萄DFR 酶（3bxx.1.A［38］）同源性最高的细风轮DFR 酶序列（编码此酶的Unigene 编号为CL3369-1）构建其空间结构，发现DFR 酶二级结构主要由α-螺旋和β-折叠构成，α-螺旋包裹着β-折叠形成致密的“夹心饼干”样结构，NAD+特异性结合位点位于第一个α-螺旋与第一个β-折叠之间的loop 区域，在植物中有高度的保守性。目前，有研究表明在花青素的生物合成中DFR 既是诱导酶［39］，又是限速酶［40］，在红叶山茶、叶菜型甘薯和花生等植物中DFR 活性增加，花青素含量也增加［41～43］。细风轮中的DFR 酶很可能在调控细风轮中花青素的合成与积累中也发挥重要作用。

本研究对细风轮转录组数据进行了深入挖掘，对细风轮花青素生物合成途径进行了详尽的分析，鉴定了花青素生物合成中关键酶基因，为细风轮花青素生物合成途径关键酶基因的克隆及其结构和功能的进一步研究提供了依据，也为将来花青素生物合成途径的解析及调控机制的研究奠定了基础。