海菜花转录组测序及生物信息学分析

李健玲，秦 波，黄 欣，蒋日红，孙 苗，梁圣华，黄耀恒，韦广绥

（1.广西壮族自治区林业科学研究院，广西特色经济林培育与利用重点实验室，广西 南宁530002；2.北京林业大学，国家花卉工程技术研究中心，北京 100083；3.广西国有高峰林场，广西 南宁 530025）

海菜花（Ottelia acuminata）属水鳖科水车前属多年生沉水植物，为中国特有种，国家二级保护植物。该植物对水质要求较高，是一种环境指示性植物，主要分布区域为广东、海南、广西、四川、贵州、云南等地[1]。海菜花是一种典型的高钾低钠型蔬菜，含有钙、铁、蛋白质、抗坏血酸、酚类等丰富的营养物质，具有较高的食用价值；其中酚类物质是天然的抗氧化剂，对DNA 损伤有良好的修复作用。另一方面，海菜花还能富集重金属元素铅，具有一定的生态修复价值[2-3]。目前，关于海菜花的研究多集中在栽培[4-5]、进化和遗传学分析[6-7]、化学成分分析[3]等方面，其分子水平的研究仅有叶绿体基因组的报道[8]。

分子生物学技术在植物领域尤其是农作物方面应用的较早且深入，但是在水生植物及药用植物上的研究报道仍然有限。高通量转录组作为初步掌握代谢通路和生物合成基因信息的一个手段，在许多物种中已经得到了广泛的应用。为了弥补海菜花转录组上的空白，研究对海菜花叶片进行转录组测序，挖掘其重要基因和调控信息，以期为海菜花功能基因和遗传多样性分析及分子育种奠定基础。

1 材料与方法

1.1 试验材料

采集新鲜的海菜花叶片，并立即投入液氮中保存。

1.2 海菜花cDNA 文库构建及测序

使用Trizol 法对海菜花的RNA 进行提取。利用Illumina TruSeqTM RNA sample prep Kit（Illumina，美国）方法构建海菜花的RNA 文库，首先使用带有Oligo（dT）的磁珠对带有polyA 尾巴的mRNA 进行富集，并用超声波将mRNA 进行片段化打断。加入随机寡核苷酸为引物，在M-MuLV 逆转录酶体系下反转合成第一条链的cDNA，之后用RNaseH 清除反应中的RNA，在DNA polymerase I 体系下，以dNTPs 为原料进行cDNA 第二条链的合成。对合成后的双链cDNA 进行纯化、末端修复、加A 尾以及连接测序接头，利用AMPure XP beads 对处理过的双链cDNA 进行筛选，挑出200 bp 左右的片段进行扩增，之后对得到的PCR 产物进行纯化，最终产物即为文库。

利用琼脂糖凝胶电泳、NanoPhotometer spectrophotometer、Qubit2.0 Fluorometer、Agilent 2100 bioanalyzer 对RNA 的质量进行检验，样品合格后进行测序。在测序的flow cell 中加入不同标记的4 种dNTP 以及DNA 聚合酶和接头引物，当测序链在延伸互补的时候，每个dNTP 会发出相应的荧光，从而被测序仪捕捉获得序列信息。

海菜花的RNA 提取、文库构建及测序均由上海凌恩生物科技有限公司完成。

1.3 转录组组装及分析

使用Trimmomatic 对原始测序数据进行质控，Trinity 进行转录组拼接，并使用Salmon 进行拼接数据的比对，之后使用Diamond 对测序数据进行NR数据库、KEGG 数据库、Swiss-Port 数据库的注释，利用eggNOG-mapper 进行eggNOG 数据库和GO 数据库的注释。

2 结果与分析

2.1 海菜花转录组测序及组装

对海菜花的叶片进行转录组测序，共检测到26 658 026 条reads，共计3 998 703 900 个碱基。对原始数据进行质控进一步去除reads 中的接头序列、质量较低的碱基以及短序列后，得到clean data 的reads 有 26 497 174 条，共计3 961 376 102 个碱基，其中Q20%高达98.90%，Q30%高达96.19%，GC含量占比为50.57%，表明海菜花转录组数据质量较高，可以开展后续的数据分析。

由于海菜花没有参考基因组，因此利用Trinity（2.8.6）对clean data 进行从头组装，得到组装后的unigene 有27 040 条，长度为27 761 688 bp，最长的unigene 有12 094 bp，N50 为1 340 bp，平均unigene 为1 026.69 bp，其中N50 高于平均长度说明组装效果较好。同时，将测序获得的高质量序列与拼接序列进行比对，比对率达到86.89%，也表明拼接效果较好。在得到的unigene 中，401～600 bp 的序列长度占比最大，达到7 134 条（26.38%），序列长度分布见图1。

图1 海菜花unigene 长度分布

2.2 海菜花转录组注释

将获得的27 040 条unigene 分别与不同数据库进行比对，结果如图2 所示。有18 910 条unigene分别被NR、GO、COG、KEGG、SWISS 这5 大数据库注释到，占unigene 总数的69.93%。其中，NR数据库比对到的unigene 最多，有18 874 条，占总数的69.80%；其次是SWISS 数据库，比对到14 020条unigene，占总数的51.85%；第三是COG 数据库，比对到13 531 条unigene，占总数的50.04%；KEGG 数据库比对到的unigene 数量最少，为7 552条，占总数的27.93%；在所有数据库中均比对到的unigene 有4 819 条，占总数的17.82%，推测为新基因。

图2 海菜花unigene 在各数据库的注释结果

2.3 NR 数据库注释

NR数据库的功能注释能够了解海菜花转录组中序列与哪些物种有同源分布。由图3 可知，海菜花与摇蚊（Clunio marinus）比对到的同源序列最多，为2 622 条，占所有比对到序列的13.89%；其次是油棕（Elaeis guineensis），达到2 062 条，占10.93%；海枣（Phoenix dactylifera）中比对到1 646条，占8.72%；另外，芋（Colocasia esculenta）、大叶藻（Zostera marina）、莲（Nelumbo nucifera）、凤梨（Ananas comosus）、小果野蕉（Musa acuminatasubsp.malaccensis）、沉水樟变型（Cinnamomum micranthumf.kanehirae）、石刁柏（Asparagus officinalis）中 分 别 比 对 到1 435、947、750、637、570、570、544 条，分别占比7.60%、5.02%、3.97%、3.38%、3.02%、3.02%、2.88%；而其他物种中比对到的序列数均低于500 条，总数达7 091 条，占比为37.57%。

图3 海菜花unigene 在NR 数据库中比对到的同源序列排名前10 的物种

2.4 COG 数据库注释

由图4 可知，海菜花转录组数据在COG 数据库中共比对到unigene 13 531 条，但功能未知的基因占多数，排在KOG 预测的第1 位，表明海菜花中有许多序列还未得到有效的验证和研究；其中，转录后修饰、蛋白质转换、伴侣蛋白（O：Posttranscriptional modification，protein turnover，chaperones），信号传导机制（T：Signal transduction mechanisms），翻译、核糖体结构与生物发生（J：Translation，ribosome structure and biogenesis），RNA转录和修饰（A：RNA processing and modification）分别排在KOG 预测的第2～5 位，代表了海菜花生长发育和信号转导中基本过程活动比较活跃。

图4 海菜花叶片转录组COG 功能注释

2.5 GO 注释

在GO 注释中，海菜花共注释到10 505 条unigene，分为生物过程（biological process）、细胞组分（cellular component）和分子功能（molecular function）3 个大类和43 个功能亚类（图5）。生物过程注释到的unigene 最多，亚类中细胞进程和代谢过程中unigene 数量最多，而碳利用和生物矿化所占比例最少。在分子功能中，催化活性和结合比对到的unigene 数量最多，蛋白标签、运货受体活性、养分库活性和小分子传感器活性比对到的数量较少。在细胞组分中包括了4 个亚类，最高的为细胞解刨学实体，最少的为其他器官部分。

图5 海菜花叶片转录组GO 注释结果

2.6 KEGG 通路分析

为了研究基因具体的功能，对海菜花转录组中的unigene 进行KEGG 通路分析，共分为5 大类和20 个亚类（图6）。第1 大类代谢过程中比对到的unigene 数量最多，包括11 个亚类，全局和概述图谱比对到的数量最多，达4 678 条unigene，其次是碳水化合物代谢，比对上1 047 条unigene，氨基酸代谢、能量代谢、酯类代谢分别比对到587、498、486 条，表明在细胞过程中关于代谢的基因较多。第2 大类是遗传信息处理，包括5 个亚类，翻译过程中比对到的unigene 最多，达到784 条，其次是折叠、分类和降解，比对到662 条，病毒信息中比对到的数量最少，为46 条。第3 大类为环境信息处理，包括2 个亚类，分别为信号转导（比对到309条）和膜转运（比对到25 条）。第4 大类为细胞过程，包括2 个亚类，为运输和分解代谢（比对到427 条）和细胞运动（比对到84 条）。最后1 类为有机系统，仅包含1 个亚类，为环境适应，比对到200 条unigene。

图6 海菜花叶片转录组KEGG 注释分类统计

2.7 SSR 分析

在海菜花转录组中共检索到4 217 个SSR 位点，包含了6 种核苷酸重复类型，具体情况如图7 所示。单核苷酸重复数量最多，为1 518 个，占比为36.00%，A/T 重复比例最高；其次为三核苷酸重复，有1 437 个，占比为34.08%，且类型较多；双核苷酸重复为1 214 个，占比为28.79%，以AG/CT 类型比例最高；四核苷酸位点重复有34 个，占比为0.81%；六核苷酸位点重复有9 个，占比为0.21%；五核苷酸位点重复数最少，仅5 个，占比0.11%。

图7 海菜花叶片转录组SSR 分布情况

3 讨 论

海菜花在景观营造、食用和药用方面有很高的开发利用价值。目前，对水生植物转录组的研究种类较少，包括圆海链藻（Thalassiosira rotula）[9]、莲（Nelumbo nucifera）[10-11]、齿叶睡莲（Nymphaea lotus）[12]、东方泽泻（Alisma orientale）[13]、芡实（Euryale feroxSalisb.ex Konig et Sims）[14]等，通过转录组测序技术对其次生代谢物的合成、花开放规律、胁迫响应和遗传标记等进行了深入研究。海菜花在生态适应性上对环境的要求较高，且有不同的生态适应类群，次生代谢物丰富，对其进行分子层面的深入研究具有重要意义。

该研究在高通量转录组测序技术的支持下，对海菜花叶片的转录组信息进行解读，共检测到unigene 27 040 条，其中有18 910 条unigene 在数据库中比对上了，获得了海菜花中大量的基因信息。在NR 数据库得到了最多的基因注释，共18 874 个，与植物同源基因比对最多的是油棕，另外包括芋、大叶藻、莲和凤梨等，这些植物类型均属于湿生环境，表明海菜花可能在转录调控和生态适应方面与湿生植物有相似的基因调控类型。在KOG 数据库中，海菜花与已知基因比对上最多的基因涉及转录后修饰、蛋白质转换、伴侣蛋白等方面。对基因通路进行KEGG 分析，共比对到7 552 个unigene，以细胞过程比对到的unigene 最多，且与代谢途径密切相关；另外，在GO注释中，生物过程注释到的unigene最多，表明在海菜花叶片中新陈代谢活动较为旺盛。SSR分析中单核苷酸、双核苷酸和三核苷酸重复是主要类型，占比达98.87%，这些位点的发现对海菜花的分子标记开发和利用具有重要价值，能为植物鉴定、遗传育种和药理研究提供直接的证据。

通过对海菜花叶片进行高通量转录组测序，获得了大量的基因序列信息和注释信息，为海菜花的功能基因和遗传多样性分析及分子育种奠定了基础，也为海菜花的生态适应性、食用和药用价值的进一步开发利用提供了依据。