椰心叶甲啮小蜂转录组分析及基因功能注释

刘华伟，李朝绪，李 芬，吕朝军，吴少英*，覃伟权*

（1.海南大学热带作物学院，海口 570228；2.中国热带农业科学院椰子研究所，海口 571399；3.海南大学植物保护学院，海口 570228；4.海南省院士团队创新中心，海口 571339；5.海南省热带油料作物生物学重点实验室，海口 571399）

椰心叶甲Brontispalongissima(Gestro)属鞘翅目Coleoptera铁甲科Hispidae害虫，其成虫和幼虫均以棕榈科植物未展开的心叶为食，并且几乎所有的棕榈科植物都可受其为害[1,2]。自2002年6月首次在海南省发现椰心叶甲为害，至2006年仅四年的时间，已蔓延至海南全省[3]。椰心叶甲啮小蜂TetrastichusbrontispaeFerrière是椰心叶甲的优势天敌寄生蜂。目前国内外已经对该蜂的生物生态学、室内繁殖和野外释放等方面进行了相关研究[4-8]，但由于缺乏遗传信息，椰心叶甲啮小蜂的遗传多样性、基因功能等的研究还较为滞后。二代测序技术的转录组测序技术（RNA-seq），具有低成本、快速和高准确性等优势，能够在缺乏基因信息的条件下获得物种的代谢和生长规律，并揭示其基因与生物学特性内在关联，同时还可获得物种大多数的转录产物[9]。本研究通过RNA-seq技术对两个种群椰心叶甲啮小蜂进行转录组测序，并进行了功能基因注释和分析，旨在为椰心叶甲啮小蜂的功能基因挖掘等分子生物学研究中提供理论基础。

1 材料与方法

1.1 材料

为后续基于转录组数据分析连续自交是否会引起椰心叶甲啮小蜂的种群退化，退化程度以及退化对其影响。本试验中所采用的两个种群椰心叶甲啮小蜂分别选用中国热带农业科学院椰子研究所天敌工厂提供的经过复壮（每年一次）的种群和中国热带农业科学院环境与植物保护研究所天敌工厂提供的自引进后隔离繁殖未经复壮蜂种。两个种群各设置3次试验重复，每个重复由20头椰心叶甲啮小蜂（雌雄蜂混合的1～2日龄成蜂）混合而成。

1.2 总RNA提取，文库构建和转录组测序

本试验将收集的两个种群椰心叶甲啮小蜂经液氮快速冷冻后，用干冰将样品送至上海派森诺生物科技股份有限公司进行总RNA提取、文库构建和转录组测序。

1.3 序列组装和转录本拼接

RNA-seq测序完成后，统计总reads的数量和长度、Q30值、模糊碱基（N）所占比例以及Q20和Q30所占比例等。经数据过滤后获得 clean reads，继而统计 clean reads的数量、总长度和占比。通过 Trinity Software（http:///trinityrnaseq.Github.io/）软件对clean reads进行拼装获得转录本（Transcript），选取每个基因最长的转录本作为代表序列（Unigene）[9,10]，并进一步统计Transcript和Unigene的总长度、序列总数、最大长度、平均长度、N50、N50%、N90、N90%以及GC含量等指标。

1.4 Unigene功能注释

根据基因的相似性，通过 BLAST比对工具（参数设置：E≤1e-5）将椰心叶甲啮小蜂转录组获得的Unigene与 NCBI non-redundant protein sequences（NR）数据库、Gene Ontology（GO）数据库、Kyoto Encyclopedia of Genes and Genome（KEGG）数据库和 evolutionary genealogy of genes: Non-supervised Orthologous Groups（eggNOG）数据库比对，进行功能注释。

2 结果与分析

2.1 椰心叶甲啮小蜂转录组测序数据及转录本拼接质量情况

转录组数据集已保存在NCBI SRA数据库中，编号：PRJNA678031。对每个样品的原始数据（Raw data）进行统计（表1），两个种群的6个样品分别得到48514092、43095558、44266534、49405442、47396830、47326810条reads，且每个样本的碱基总量均在6.46 Gb以上；碱基百分比Q20均大于97.47%、碱基百分比Q30均大于93.44%，含量相近；模糊碱基占比在0.001485%～0.001512%的较低水平，通过数据过滤后统计发现（表2），Clean Reads占reads总数的93.94%～94.37%。对Transcript和Unigene序列进行统计（表3），共得到78930条Transcript和29535条Unigene，Unigene的长度分布统计图见图1，总长度为51330466 bp，平均长度为1737.95 bp，N50值为3547 bp。

图1 椰心叶甲啮小蜂转录组Unigene组装长度分布统计Fig.1 Statistics of Unigene length distribution of T.brontispae transcriptome

表1 椰心叶甲啮小蜂转录组测序数据及质量情况统计Table 1 Transcriptome sequencing data and quality statistics of T.brontispae

表2 椰心叶甲啮小蜂转录组高质量Reads与碱基统计Table 2 High quality Reads and base statistics of T.brontispae transcriptome

表3 椰心叶甲啮小蜂转录组Transcript和Unigene序列统计表Table 3 Statistics of transcript and Unigene sequences of T.brontispae transcriptome

以上数据表明样本的测序质量，文库构建质量和拼接的完整性都较好，可以用于后续的分析。

2.2 Unigene功能注释情况

通过Unigene和四个公共数据库的比对并进行功能注释（表4）。共有13796条Unigene被注释，其在4个数据库中的注释如下：13401条基因（45.37%）在NR数据库中获得了注释，3834条基因（12.98%）在GO数据库中获得了注释，12707条基因（43.02%）在eggNOG数据库中获得了注释，5999条基因（20.31%）在KEGG数据库中获得了注释。

表4 椰心叶甲啮小蜂转录组注释结果汇总Table 4 Summary of annotated results of T.brontispae transcriptome

2.3 椰心叶甲啮小蜂转录组Unigene的NR功能分类

通过 Blastx将椰心叶甲啮小蜂转录组中的Unigene与NR数据库进行比对，注释到13401条Unigene。统计注释结果并绘制物种分布图（图2），结果显示椰心叶甲啮小蜂与蝇茧蜂Diachasmaalloeum、佛罗里达弓背蚁Camponotusfloridanus、黑褐毛蚁Lasiusniger、赤眼蜂Trichogrammapretiosum、多胚跳小蜂Copidosomafloridanum、榕小蜂Ceratosolensolmsimarchali、丽蝇蛹集金小蜂Nasoniavitripennis都有同源序列分布，其中与丽蝇蛹集金小蜂相似序列最多，占50.34%，与榕小蜂、佛罗里达跳小蜂和赤眼蜂分别有13.3%、7.42%和5.95%的相似序列，与黑蚁、佛罗里达弓背蚁和蝇茧蜂中的同源序列较少（0.98%～1.57%）。

图2 椰心叶甲啮小蜂转录组NR注释结果统计图Fig.2 NR annotated species distribution of T.brontispae transcriptome

2.4 椰心叶甲啮小蜂转录组Unigene的GO分类注释

对椰心叶甲啮小蜂转录组Unigene进行GO功能分类，共有3834条（12.98%）Unigene获得注释。将注释得到的Unigene划分为三大类（生物过程、细胞成分和分子功能）67个分支（图3），统计注释到每一类的基因数量，发现在生物过程类中细胞过程（1527条Unigene）的占比最大（39.83%）；在细胞成分大类中膜（1324条Unigene）的占比最大（34.53%）；在分子功能大类中结合（1733条Unigene）的占比最大（45.20%）。

图3 椰心叶甲啮小蜂转录组Unigene的GO功能分类Fig.3 GO analysis of T.brontispae transcriptome Unigene

2.5 椰心叶甲啮小蜂转录组Unigene的eggNOG功能注释

将椰心叶甲啮小蜂转录组的Unigene与eggNOG数据库进行比对，注释到12707条Unigene，根据功能可将其划分为26个功能区域（图4），统计注释到各类功能的基因数量，发现无特征基因占比最多（63.06%），其次是新陈代谢类（33.49%），细胞过程和信号传递，信息存储与处理类分别占 20.82%和19.67%。

图4 椰心叶甲啮小蜂转录组Unigene的eggNOG功能分类Fig.4 eggNOG analysis of T.brontispae transcriptome Unigene

2.6 椰心叶甲啮小蜂转录组Unigene的KEGG代谢途径分类分析

将椰心叶甲啮小蜂转录组的Unigene与KEGG数据库进行对比，共有5999条Unigene获得注释，涉及的代谢通路可归为五个大类（新陈代谢、遗传信息处理、环境信息处理、细胞进程和有机系统）35个子类（图 5）。统计注释到各类通路的基因数量，发现在新陈代谢一类中获得注释最多的是碳水化合物代谢（299条Unigene）；遗传信息处理一类中翻译获得注释最多（361条Unigene）；环境信息处理一类中获得注释最多的是信号转导（735条Unigene），细胞进程和有机系统两类中获得注释最多的分别为运输和分解代谢（382条Unigene）和内分泌系统（388条Unigene）。

图5 椰心叶甲啮小蜂转录组KEGG注释统计图Fig.5 KEGG annotation statistical chart of T.brontispae transcriptome

3 讨论

第二代高通量测序技术（Next-Generation Sequencing，NGS）因测序时间短、成本低、高准确性和所获得数据量大等优点，被广泛应用于非模式生物分子生物学研究中[9,10]。椰心叶甲啮小蜂是椰心叶甲的优势寄生蜂，但是由于其遗传信息数据的缺乏，对其分子生物学的研究仍然较少。本研究旨在通过该技术分别对经过复壮和未经复壮两个种群的椰心叶甲啮小蜂进行转录组测序和拼装，揭示椰心叶甲啮小蜂整体基因表达特征，并为后续基于转录组数据分析复壮对椰心叶甲啮小蜂的影响，种群遗传结构等研究提供数据支撑。

一般来说，Q30在80%以上，N50值不小于800 bp（N50值越大表示长片段越多）就可以认为测序质量可靠，序列组装的完整性较好[11-13]。本研究通过无参考基因组分析椰心叶甲啮小蜂转录组的特异性，序列拼接后共获得29535条Unigene，Q30和N50值分别为93.44%和3547 bp，可以认为本研究中所得到的测序结果质量和序列拼装的完整性都较好，可以满足后续分析的基本要求。

对组装后所获得的29535条Unigene在NR、GO、KEGG和eggNOG数据库进行基因功能注释，共有13796条Unigenes被注释，仍有15739条Unigenes未被注释。这一结果在许多生物的转录组测序中都有存在，如植物的彩色马蹄莲[10,14]，动物中的云锦杜鹃[15]和昆虫中的丽蝇蛹集金小蜂[16]都有出现，可能由于椰心叶甲啮小蜂缺乏基因组方面研究的基础资料，使得部分Unigene在数据库中无法得以注释[11]，也可能是由于部分Unigene片段太短或椰心叶甲啮小蜂中存在新的功能基因而导致的[17]。

与Nr公共数据库比对发现，有13401条Unigene获得注释（45.37%），其中与丽蝇蛹集金小蜂相似的序列最多（50.34%）。在椰心叶甲啮小蜂的转录组数据中有19.36%的Unigene属于其他序列，这在很多物种的转录组测序中也有出现[18-20]，可能属于椰心叶甲啮小蜂自身特有的与大多数物种不同的序列，也可能这些基因属于非编码RNA或不是功能基因。此外，还有54.63%的Unigene未获得相关注释信息，此现象在其他昆虫中也存在，如二化螟盘绒茧蜂，菜粉蝶等[21,22]，其原因可能是Unigene片段过短，基因信息的暂时缺乏，或者椰心叶甲啮小蜂中存在新的功能基因[23]。在G0数据库中有3834条Unigene获得注释，仅占总Unigene的12.98%，相对注释率较低，这进一步说明了椰心叶甲啮小蜂基因信息的缺乏[24]，也可能与G0数据库信息不够完善有关[25]。eggNOG数据库中注释到的Unigene也较多（43.02%），除无特征序列外，注释到与新陈代谢有关的Unigene较多（24.44%），说明椰心叶甲啮小蜂自身具有较强的代谢能力。在 KEGG数据库中注释到 5999条 Unigene，同样与新陈代谢有关的代谢通路最多，有 1842条 Unigene（30.7%），与eggNOG分析结果相一致。这些注释结果为在分子层面研究椰心叶甲啮小蜂提供了充足的依据。

综上，本研究利用转录组测序技术，对椰心叶甲啮小蜂复壮和未复壮两个种群进行了转录组测序，并利用四大数据库对获取的Unigene进行了功能注释、分类和代谢途径预测等，为下一步深入研究椰心叶甲啮小蜂两个种群的差异表达基因和其他分子生物学研究提供了参考。