基于转录组分析椰心叶甲啮小蜂复壮和保种种群差异表达基因

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

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

重大危险性外来有害生物椰心叶甲Brontispa longissima原产于印度尼西亚和巴布亚新几内亚，主要危害植物心叶组织，造成心叶大量死亡，严重影响棕榈植物的景观度和植株的健康生长[1]。2002年首次在海南省海口市椰子树上发现椰心叶甲[2]，随后短短几年便传播至海南全省[3]。近年来，椰心叶甲在海南省各地槟榔园频繁发生，扩散蔓延快[4－5]，对椰心叶甲防治研究仍需加强。椰心叶甲啮小峰Tetrastichus brontispae是椰心叶甲的蛹期寄生蜂[6]，国内外学者已对椰心叶甲啮小蜂的成虫习性、寄主适应性、生物生态学特性，以及室内大量繁殖和野外释放等方面进行了相关研究[6－10]。饲养并释放椰心叶甲啮小蜂可作为防治椰心叶甲的有效手段。

自然界中的个体在种群进化的过程中，为了适应生存环境的变化会产生一系列的变异，当这种变异朝着各项生理指标下降的方向发展，就称之为退化[11]。室内饲养的昆虫退化往往表现在繁殖能力下降、种群整齐度下降、寿命缩短和抗病能力下降等方面。昆虫退化现象是普遍存在的，目前关于黑腹果蝇Drosophila melanogaster[12]、烟蚜茧蜂Aphidius gifuensis[13]和管氏肿腿蜂Scleroderma guani[14]种群的退化和复壮已经进行了相关研究，都表明长时间的近亲繁殖会导致昆虫种群的退化，通过回交和外源种群杂交的方法可以有效复壮。为探究椰心叶甲啮小蜂种群退化的分子机制，笔者基于椰心叶甲啮小蜂2个种群的转录组数据，对差异基因进行分析。

1 材料与方法

1.1 供试蜂种 选用中国热带农业科学院椰子研究所天敌工厂种群(每年复壮1次)和中国热带农业科学院环境与植物保护研究所天敌工厂种群(自2004年引进后隔离繁殖，未经复壮)，各设置3次重复，每重复20头，均为1～2日龄的椰心叶甲啮小蜂(不分雌雄)。

1.2 转录组数据 用于对比分析的转录组数据已提交至NCBI SRA数据库中，Genbank登录号:PRJNA678031。原始数据质量分析详见刘华伟等[15]。测序的质量较好，可用于后续差异基因的分析。

1.3 差异基因表达分析及功能富集分析 皮尔逊相关系数(pearson′s correlation coefficient，r2)作为样品间相关性的评估指标[16]，r2越接近1，说明两个样品之间的相关性越强。通过聚类分析判断差异表达基因在不同实验条件下的表达模式，在样品间表达量相关性高的基因被归为一类，使用R语言Pheatmap软件包对所有比较组的差异基因的并集和样品进行双向聚类分析，根据同一基因在不同样品中的表达水平和同一样品中不同基因的表达模式进行聚类，采用Euclidean方法计算距离，采用层次聚类最长距离法(Complete Linkage)进行聚类。通过DESeq(http://www.bioconductor.org/packages/release/bioc/html/DESeq.html)软件对复壮种群(WC)和保种种群(DZ)进行差异表达分析，筛选差异表达基因条件为:表达差异倍数log2FoldChange＞1，显著性 P-value＜0.05。 使用topGO进行 GO富集分析，显著富集的标准为P-value＜0.05，确定差异基因行使的主要生物学功能。将差异基因与KEGG数据库比对，显著富集的标准为P-value＜0.05，获得差异基因显著性富集的通路注释信息。

1.4 RT-qPCR验证 选择与椰心叶甲啮小蜂各种行为相关的 7个基因(编号为 DN935，DN4191，DN4534，DN26397，DN14402，DN535，DN1980)进行荧光定量PCR，以验证转录组数据的可靠性。

2 结果与分析

2.1 椰心叶甲啮小蜂差异基因表达量分析 通过对复壮种群3个重复(WCA、WCB和WCC)与保种种群3个重复(DZA、DZB和DZC)间的差异基因表达量进行分析，各样本间的相关系数r2均大于0.95，表示种群内各样本间及复壮与保种种群之间的相关性较强(表1)。通过聚类分析发现，在复壮种群和保种种群的6个样本中，高表达基因和低表达基因的分布在不同种群之间出现分离，在同一种群的不同生物学重复中聚集在一起。表明所用样本在组内的基因表达模式相似，生物学重复性较好，且样本分组较合理(图1)。

表1 样品基因表达量相关性Tab.1 Correlation of sequenced samples based on gene expression level

图1 差异表达基因表达模式聚类图Fig.1 Expression pattern clustering of DEGs

相较于复壮种群，保种种群(DZ/WC)共有553个差异表达基因，其中上调基因206个、下调基因347个。两组样本间基因差异表达的分布情况见图2。其中，有3个基因的注释结果为蜂王浆主蛋白(major royal jelly proteins，MRJPs)，3 个基因的注释结果为细胞色素P450(cytochrome P450，CYP450)，它们在保种种群中下调表达；此外，还发现几种嗅觉相关蛋白。

图2 基因差异表达火山图分析Fig.2 Volcano plot of DEGs

2.2 差异表达基因的GO富集分析和KEGG通路分析 差异表达基因共有14个基因显著富集到分子功能(molecular function，MF)中，并且全部下调，富集到单氧酶活动(monooxygenase activity)的差异基因数量最多，占供试基因总数的28.57%。有25个基因显著富集到细胞组分(cellular component，CC)中，其中10个基因上调，15个基因下调，富集到细胞膜(cell membrane)中的差异基因数量最多，占供试基因总数的68%。有25个基因显著富集到生物过程(biological process，BP)中，其中13个基因上调，12个基因下调，此分类中各项显著富集到的差异基因都较少，均为1～2个。在显著富集的差异基因中，包括几种气味受体和细胞色素 P450(图3)。

图3 差异表达基因GO分类Fig.3 GO classification of DEGs

差异表达基因共有4条显著富集的代谢通路，分别为Hh信号通路(Hedgehog signaling pathway)、Hh信号通路－fly(Hedgehog signaling pathwayfly)、胰腺分泌(pancreatic secretion)和神经活性配体受体相互作用(neuroactive ligand-receptor interaction)。因4条通路涉及的基因有重合，所以共计10个差异基因(表2)。

表2 差异表达基因的显著富集代谢通路Tab.2 KEGG analysis of DEGs

2.3 RT-qPCR验证转录组数据 在椰心叶甲啮小蜂保种和复壮2个种群中，7个基因的表达变化趋势与转录组的测序结果是一致的，可以认为转录组的数据是可靠准确的(图4)。

图4 各基因转录组及RT-qPCR统计图Fig.4 Transcriptome and RT-qPCR of each gene

3 结论与讨论

有研究证明胰蛋白酶(TRINITY＿DN334＿c0＿g1，TRINITY＿DN7050＿c0＿g1)与自身免疫防御有关[17]，还发现离子型谷氨酸受体(TRINITY＿DN29200＿c0＿g1)属于嗅觉相关蛋白中的离子型受体(IRs)，参与化学感受基因感觉[18]。在差异表达基因的KEGG通路富集分析中发现有2条通路中的基因全部下调，表示保种种群的椰心叶甲啮小蜂的自身免疫防御和嗅觉系统可能出现退化。

蜂王浆主蛋白是蜂王浆中的水溶性蛋白，与蜂王浆的生物学功能和医疗保健作用可能有密切的关系，比如蜂王浆主蛋白参与雌蜂的级型分化，为幼虫发育提供可利用的氮元素并可能在蜜蜂行为调控中起作用。另外还发现蜂王浆主蛋白在其他动物中具有重要的医疗保健作用，主要表现为高效的抗菌活性、较强的免疫活性和促进特殊细胞增殖等[19]。在椰心叶甲啮小蜂转录组数据中，共检测出25个蜂王浆主蛋白基因，其中3个基因差异表达且全部下调，可以说明保种种群椰心叶甲啮小蜂的抗菌活性、免疫活性等出现退化。

细胞色素 P450(cytochrome P450，CYP450)是一个超家族酶系，几乎存在于昆虫的所有组织中，参与外源化合物的代谢、蜕皮激素和保幼激素的合成与降解等过程，还可以调节昆虫的生长和发育[20]。在椰心叶甲啮小蜂转录组中发现124个P450基因，其中有3个基因差异表达，全部下调。在对飞蝗CYP450的研究中发现，CYP303A1可能与蜕皮激素代谢调控和影响新表皮增厚有关[21]，进而搜索椰心叶甲啮小蜂转录组数据时发现一个幼虫表皮蛋白也同样下调，表明保种种群的椰心叶甲啮小蜂的幼虫发育过程可能受阻。但在果蝇的研究中表示，CYP303A1与果蝇刚毛结构和机械感觉有关[22]，在茶小绿叶蝉 Empoasca flavescens[23]中还发现CYP303A1参与药物代谢解毒过程，出现这种情况的原因可能是物种的亲缘关系较远，同一个基因在不同物种间可能存在差异。

嗅觉系统是帮助昆虫识别各种化学气味分子，引导昆虫进行觅食、交配、选择产卵地点和躲避天敌等生命活动的重要系统[24－25]。研究表明，有多种嗅觉相关蛋白参与了嗅觉识别过程，主要有气味结合蛋白(odorant binding proteins，OBPs)、气味降解酶(odor degrading enzymes，ODEs)、嗅觉受体(olfactory receptors，OLRs)、气味受体(odorant receptors，ORs)、化学感受蛋白(chemosensory proteins，CSPs)、感觉神经元膜蛋白(sensory neuron membrane proteins，SNMPs)及离子型受体(ionotropic receptors，IRs)等[26]。在椰心叶甲啮小蜂转录组数据中，共发现153个 ORs基因，394个OBPs基因，18个SNMPs基因和16个IRs基因，并未发现ODEs、OLRs、CSPs相关基因。在这些基因中，共有7个基因在两个种群中差异表达，包括4个ORs基因上调，2个OBPs基因和1个IRs基因下调。一般认为OBPs和ORs都是在嗅觉识别进程中起作用，气味结合蛋白结合和转运气味分子，气味受体完成信号的激活和转导，但在本研究中ORs(上调)和OBPs(下调)的表达模式相反。在对果蝇的研究中发现，OR82a基因会抑制果蝇的求偶行为[27]，但其他ORs的功能仍不明确，需要继续研究。

通过对椰心叶甲啮小蜂2个种群的转录组数据进行差异基因分析，发现保种种群在一些方面已经出现了明显的退化，这也表明在室内大量繁殖过程中，合理的繁殖技术和适时的复壮是非常有必要的，后续应针对退化规律、复壮时间等展开研究，以期为更新现有的繁蜂技术提供数据支撑。