金线莲转录组测序及其黄酮类合成相关基因分析

邹福贤，许 文，黄泽豪，张 勋，陈抒云，林 羽，徐 伟

(福建中医药大学药学院，福州 350122)

金线莲为兰科植物花叶开唇兰(Anoectochilusroxburghii(Wall.) Lindl.)干燥全草，后更名为金线兰,具有清热凉血、祛风利湿之功效，用于治疗肾炎、支气管炎、膀胱炎、糖尿病、风湿性关节炎等[1]，为珍贵药用植物，主要分布于我国福建、台湾、云南等省份以及孟加拉国、老挝、泰国、越南等东南亚国家[2]。金线莲中主要含有多糖、黄酮、有机酸、挥发性化合物等[3]，目前多以黄酮类(芦丁、槲皮素、异鼠李素、山柰酚等)成分为指标性成分对金线莲进行质量评价[4-5]。现代药理研究表明其黄酮类成分具有抗糖尿病活性，能够降血糖[6]，并能减轻由糖尿病引起的肾损伤[7]；保护肝脏[8]等作用，因此金线莲黄酮类成分被认为是其重要药效物质基础，但是由于缺乏金线莲的基因组信息，其黄酮类成分生物合成途径尚未深入研究。

随着第2代高通量测序技术的发展，药用植物基因组学研究工作也不断深入。高通量转录组测序技术能够在没有参考基因情况下，对某一物种的转录组进行全面分析，确定基因序列和转录本，为基因组图谱尚未完成或数据信息匮乏的物种基因组和转录组提供技术支持[9]。本研究以仿野生林下种植的不同时期金线莲为样品，首次采用高通量测序技术对金线莲进行转录组分析，以期获得大量的黄酮类成分合成的基因信息，从基因表达水平阐述金线莲药用活性成分的生物合成途径，为金线莲的规范化种植以及质量评价提供参考。

1 材 料

1.1 样品采集

金线莲样品采集自福建省永泰县金线莲林下种植基地(N25°57′23.17″，E118°50′31.28″，748.625 m)，分别选取组织培养品(以下简写为AR1)、种植1个月(以下简写为AR2)和种植4个月(以下简写为AR3)的金线莲新鲜叶片作为研究材料，基原经福建中医药大学黄泽豪教授和范世明高级实验师鉴定为Anoectochilusroxburghii。样品标本存储于福建中医药大学药学院标本室，标本号分别为：35070020180915100LY，35070020180915101LY，35070020180915102LY，新鲜样品清洗后液氮速冻并保存-80 ℃下备用。

1.2 试 剂

植物总RNA快速提取试剂盒(北京百泰克生物技术有限公司)；QPCR试剂ABScript Ⅱ cDNA First-Strand Synthesis Kit和2X SYBR Green Fast qPCR Mix with High ROX(爱博泰克生物科技有限公司)；引物采用软件Primer Premier 6.0设计并由上海生工生物工程股份有限公司合成；乙腈(色谱纯，德国Merck公司)，其他试剂均为分析纯。

1.3 仪 器

Prism 7900HT型荧光定量PCR仪(美国ABI公司)；Shimadzu LC-20AT高效液相色谱仪(日本岛津公司);AR224CN型分析天平(奥豪斯仪器有限公司)；Milli-Q超纯水仪(美国Millipore公司)；TC-C18色谱柱(4.6 mm×250 mm，5 μm，美国Agilent公司)。

2 方 法

2.1 测序数据与质量控制

每个生长时期的金线莲随机选取10株的叶片委托广州基迪奥基因公司负责测序和分析工作。得到的原始序列(raw reads)中因含有带接头的，重复的，测序质量很低的 reads，这些reads会影响组装和后续分析，因此对其进行数据过滤。数据过滤的步骤如下：①去除含adaptor的reads，②去除N的比例大于10%的reads，③去除低质量reads(质量值Q≤20的碱基数占整个read的40%以上)，④获得高质量clean reads。

2.2 De novo组装

本实验利用Trinity[10]组装软件进行金线莲转录组De novo组装。Trinity首先将具有一定长度overlap的reads连成更长的片段，这些通过reads overlap 关系得到的不含N的组装片段作为组装出来的Unigene。

2.3 Unigene基本功能注释

将Unigene通过NCBI中Blastx比对到蛋白质数据库Nr、Swiss-Prot、KEGG和KOG(E-value<0.000 01)，得到序列相似性最高的蛋白被用于进行金线莲Unigene的功能注释。剩余部分，则采用ESTScan[11]软件预测其编码区。利用Blast2GO[12]软件将在Nr数据库中获得注释到的Unigene进行GO注释；并运用WEGO[13]软件对所有Unigene进行GO功能分类统计，确定金线莲相关功能基因分布特征。物种分布统计是利用Blastx将组装出来的Unigene序列与Nr数据库进行比对后，取每个Unigene在Nr库中比对结果最好(E值最低)的那一条序列为对应同源序列确定同源序列所属物种，统计比对到各个物种的同源序列数量。

2.4 基因表达量计算

用RPKM法(Reads Per Kb per Million reads)计算基因表达水平，见公式(1)：

RPKM=(1 000 000×C)/(N×L/1 000)

(1)

设RPKM为Unigene A的表达量，则C为比对到Unigene A的reads数，N为比对到所有Unigene的总reads数，L为Unigene A的碱基数。RPKM法能消除基因长度和测序量差异对计算基因表达的影响，计算得到的基因表达量可直接用于比较不同样品间的基因表达差异。

2.5 差异基因筛选

利用经过FDR(false discovery rate)[14]校正后的P与 log2FC来筛选差异基因，筛选条件为：校正后的P<0.05且|log2FC|>1。其中log2FC为样品1与样品2 RPKM差异倍数的对数值，以2为底数计算得到的值，见公式(2)。

log2FC=log2(RPKMsample_2/RPKMsample_1)

(2)

2.6 通路显著性富集分析

以KEGG通路为单位，应用超几何检验，找出与整个基因组背景相比，在差异表达基因中显著性富集的通路。

2.7 Q-PCR扩增

Q-PCR引物利用Primer Premier 6.0软件进行设计。Q-PCR操作方法按照试剂盒操作说明，基因相对表达量采用2-ΔΔCt法进行计算。

2.8 含量测定

采用HPLC的方法对金线莲中6种主要黄酮类成分进行含量测定[15]。

3 结 果

3.1 金线莲RNA-seq测序及质量评价

通过IlluminaHiseqTM2000高通量测序，3个样品的原始reads数在50 629 778～54 654 486之间，每个reads含有150个碱基，GC含量在48.46%～49.63%之间。经过数据过滤之后，3个样品平均剩余98.09%的高质量clean reads，且碱基质量值Q30(碱基错误率小于0.01%)均在90%以上，说明测序质量可靠(见表1)。

Table 1 Comparison of data before and after filtration in transcriptome of Anoectochilus roxburghii (AR)

3.2 金线莲转录组组装与分析

利用Trinity软件对上述的高质量reads根据序列之间的重叠部分进行混合组装。组装结果质量评估可从N50数值和测序深度(Sequencing Depth)来评估(一般认为N50长度超过800 bp就可认为组装序列完整性较好)；本实验所组装得到的序列N50均超过1 400 bp，测序深度均在100×以上(一般认为测序深度应达到10～15×以上[16])，测序覆盖度和测序错误率控制均得以保证，组装序列完整性较好(见表2)，可以用于后续分析和信息挖掘工作。测序深度为测序得到的碱基总量(bp)与基因组大小(Genome)的比值。

Table 2 Data assembly for unigenes and transcript in transcriptome of AR

3.3 基因功能注释及分类

将3个样品组装得到的基因进行合并去除重复，共得到51 370条基因，并将这些基因序列分别与Nr、GO、Swiss-Prot、KO和KOG数据库进行比对。经过比对，发现在Nr数据库中被注释到的基因最多，有22 970条；最少的被注释到的数据库为KEGG数据库，其中有5 555条基因在5个数据库中均被注释到(图1)。

Figure 1 Statistics of annotation results in Nr,GO,Swiss-Prot,KEGG and KOG databases

物种分布统计表明，金线莲与油棕Elaeisguineensis和海枣Phoenixdactylifera具有的同源序列最多，分别占到总序列的24.27%和22.39%，其他的主要为小果野蕉Musaacuminata(9.89%)、菠萝Ananascomosus(4.83%)、木豆Cajanuscajan(3.29%)、水芙蓉Nelumbonucifera(2.53%)、粳稻群体OryzasativaJaponicaGroup(2.21%)、可可Theobromacacao(2.04%)、葡萄Vitisvinifera(1.51%)、野胡萝卜Daucuscarota(1.21%)，所匹配到的物种均和金线莲一样同为单子叶植物。经GO数据库注释到的基因根据其功能可以分为3大类：生物过程(biological process)、细胞组分(cellular component)和分子功能(molecular function)；其中生物过程中细胞进程(cellular process)和代谢进程(metabolic process)最高，细胞组分中细胞(cell)和细胞部分(cell part)最高，分子功能中绑定(binding)和催化活性(catalytic activity)最高(图2)。另外将金线莲转录组数据与台湾银线莲(Anoectochilusformosanus)转录组数据进行对比，各分类中占比最高的两个基因功能与台湾银线莲的趋势一致[17]。

3.4 差异基因表达量比较

通过“2.5”项下筛选条件发现共有19 065个基因具有显著差异，占总基因数37.11%。3个样品之间两两比较，发现AR1-VS-AR3之间的差异基因最多，有12 931个；而AR2-VS-AR3两组之间的差异基因最少，为5 946个；其中764个基因在任意两组之间均有显著性差异(图3)。

3.5 差异基因的通路富集分析

在生物体内，不同基因相互协调行使其生物学功能，而基于通路的分析有助于更进一步了解基因的生物学功能。在显著差异基因中共有1643条基因被注释到KEGG数据库中，AR1-VS-AR2、AR1-VS-AR3和AR2-VS-AR3中被注释的基因数分别为1 279、991和378条。通过统计各组之间具有显著性差异(P<0.05)的通路，3组之间共有27条具有显著性富集的通路，其中黄酮类生物合成(Flavonoid biosynthesis)、苯丙素类生物合成(Phenylpropanoid biosynthesis)、α-亚麻酸代谢(alpha-Linolenic acid metabolism)、氰氨基酸代谢(Cyanoamino acid metabolism)4条通路在3组两两之间的比较均具有显著性富集(表3)，其中黄酮类生物合成差异基因改变最为明显。

黄酮类生物合成通路(flavonoid biosynthesis)直接与次生代谢产物黄酮类成分具有直接关系，该代谢通路中共有22个基因具有显著性差异，涉及到6种酶，分别为反式肉桂酸4-单加氧酶(trans-cinnamate 4-monooxygenase，C4H)、咖啡酰辅酶AO-甲基转移酶(caffeoyl-CoAO-methyltransferase，CCOMT)、查尔酮合成酶(chalcone synthase，CHS)、黄酮醇合成酶(flavonol synthase，FLS)、莽草酸O-羟基肉桂酰转移酶(shikimateO-hydroxycinnamoyltransferase，HST)、黄酮类3′，5′-羟化酶(flavonoid 3′,5′-hydroxylase，F3′5′H)，各类酶及其对应基因的总FPKM值叠加图见图4。随着生长时间的增加，种植后的样品黄酮类生物合成相关酶的表达量基本高于组织培养样品，其中C4H、CCOMT、HST和F3′5′H的表达量呈现上调趋势，总体表现出种植后的表达量高于组织培养的趋势。

Figure 2 GO function classification of unigenes in transcriptome of AR

Figure 3 Number of up-regulated and down-regulated in differentially expressed genes between groups

Table 3 Significant enrichment pathways in differentially expressed genes between groups and corresponding P value

Figure 4 FPKM results of six unigenes involved in flavonoid biosynthesis of AR

3.6 Q-PCR分析

引物设计结果见表4，采用Zhang等[18]确定的ACT2作为金线莲内参基因。对上述每一类酶的变化差异表达基因进行Q-PCR验证，以AR1为参照样品，AR2和AR3的各基因相对表达量见图5。Q-PCR结果表明验证的所有基因中，AR3的相对表达量均显著性高于AR1，AR2的相对表达量大多数与AR1无显著性差异，但综合整个AR1，AR2和AR3的相对表达量来说，基本呈现出随着种植时间的增加，基因相对表达量也相应增加的趋势，与RNA-seq的数据基本符合。

3.7 金线莲黄酮类成分含量测定

选择相同批次样品，对其主要的黄酮类成分：芦丁(rutin)、异槲皮苷(isoquercitrin)、水仙苷(narcissin)、槲皮素(quercetin)、山柰酚(kaempferol)和异鼠李素(isorhamnrtin)进行含量测定，HPLC图谱见图6，含量结果见图7。金线莲黄酮含量随着种植时间的增加而增加，种植时期的含量极显著高于组织培养的，种植4个月的金线莲黄酮类成分含量显著性高于种植1个月的金线莲黄酮类成分含量。

Table 4 Primers for Q-PCR to verify six genes involved in flavonoid biosynthesis pathway in AR

Figure 5 Q-PCR verification for six genes involved in flavonoid biosynthesis pathway in AR

Figure 6 HPLC chromatograms of mixed reference solution (A),and analytical sample solution (B)

Figure 7 Content changes of six flavonoids in AR1,AR2 and AR3

3.8 金线莲中黄酮类成分生物合成相关基因鉴定

根据KEGG数据库的注释，筛选出和黄酮类成分合成相关的基因及其对应的酶，主要包括CHS、查尔酮异构酶(Chalcone isomerase，CHI)、柚皮素3-双加氧酶(Naringenin 3-dioxygenase，F3H)、FLS、F3′H、黄酮类甲氧基转移酶(flavoneO-methyltransferase，FOMT)、黄酮醇-3-O-糖基转移酶(Flavonol 3-O-glucosyltransferase，F3GT)、黄酮醇-3-O-葡萄糖苷L-鼠李糖基转移酶(Flavonol-3-O-glucoside L-rhamnosyltransferase，FG2)，结合金线莲中已经被报道的黄酮类成分，推测其生物合成途径见图8。

Figure 8 Probable flavonoids biosynthesis pathway of AR

4 讨 论

本研究首次对金线莲采用第2代转录组测序技术进行测序。通过该技术共得到金线莲的51 370 条Unigenes，在Nr、GO、Swiss-Prot、KO和KOG 5个数据库中均有大量基因被注释，为金线莲的育种开发及有效成分生物合成研究提供基础。课题组前期开展了金线莲黄酮在根、茎、叶的分布，其中金线莲中叶的总黄酮含量占全草的总黄酮含量超过80%，而茎和根约各占不到10%，为了更好的表征金线莲黄酮的合成相关基因的变化，本研究选择其表达量最为丰富的叶部位进行转录组研究。此外，本研究还首次同时结合Q-PCR技术和化学成分含量测定对金线莲转录组测序结果进行验证。

在不同生长时期的金线莲两两比较中，有较多的基因具有显著性差异，通过KEGG通路富集分析，主要差异基因主要涉及黄酮类生物合成、苯丙素生物合成、α-亚麻酸代谢和氰氨基酸代谢4条通路。黄酮类生物合成通路中大部分关键酶的表达量基本呈现上升趋势，经过Q-PCR验证和化学成分含量检测，AR3的黄酮类成分合成相关酶的表达高于AR1，与总黄酮含量变化一致。而AR1和AR2相比，C4H酶的表达量显著升高，CCOMT、CHS、F3′5′H等酶的相关基因表达量具有上升趋势。黄酮生物合成关键酶具有遗传多样性，基因表达比较复杂，整个黄酮类成分的生物合成具有复杂性，合成途径存在许多交叉点[19-20]，化合物的含量变化可能是由于多个酶进行调控的，有待更进一步深入研究金线莲的黄酮类成分生物合成网络调控。

本研究依靠转录组学技术揭示金线莲黄酮类成分生物合成途径，能够为金线莲的黄酮类成分生物合成及调控机制的相关研究奠定基础。

致 谢：本项目在福建省科技厅重点实验室和福建省中药资源研究与开发利用重点实验室完成。