伞裙追寄蝇滞育关联基因的转录组学分析

张 博, 韩海斌, 徐林波, 高书晶, 甘 霖, 岳方正, 刘爱萍,*

(1.中国农业科学院草原研究所, 呼和浩特 100010; 2.阿拉善盟科技信息研究所, 内蒙古阿拉善 750399;3.国家林业和草原局森林和草原病虫害防治总站, 沈阳 110034)

伞裙追寄蝇Exoristacivilis是一种牧草优势寄生性天敌，隶属于双翅目(Diptera)寄蝇科(Tachinidae)，能够对粘虫Mythimnaseparata、草地螟Loxostegesticticalis、草地贪夜蛾Spodopterafrugiperda等多种牧草害虫起到有效的防控作用(王建梅等, 2013, 2015)。昆虫滞育是一种生理过程，其特征在于行为修饰、代谢抑制以及增强的胁迫耐受性(King and Macrae, 2015)。近年来，关于伞裙追寄蝇滞育类型、滞育温光周期的诱导以及滞育关联蛋白的研究取得了一定的成果(韩海斌等, 2018; 姜珊, 2018)，为揭示其滞育的分子机制奠定了基础。

新一代测序技术的出现，极大促进了昆虫转录组学研究(张棋麟和袁明龙, 2013)。目前利用高通量技术对双翅目昆虫的滞育研究也取得了一定的进展，利用高通量技术发现了在柑橘大实蝇Bactroceraminax滞育调控的相关研究中筛选出12个显著上调和6个显著下调的滞育关联基因以及重要的滞育关联蛋白hsp23, hsp70和hsp90(Luetal., 2014; Wangetal., 2019)。张倩等(2019)利用iTRAQ技术筛选淡色库蚊Culexpipienspallens越冬后期和滞育期的差异表达蛋白，共鉴定出244个差异表达蛋白，包括126个上调蛋白和118个下调蛋白，这些差异表达蛋白的功能与脂质代谢、细胞骨架重塑、碳水化合物代谢等有关。

目前，调控伞裙追寄蝇滞育的分子调控机制仍然不完全明确。为此，本研究对伞裙追寄蝇正常发育及滞育的蛹进行转录组测序，对筛选出的滞育关联基因进行KEGG通路富集分析，了解滞育关联基因主要参与的代谢途径和信号转导机制，为进一步研究伞裙追寄蝇的滞育调控基因以及分子机理奠定基础，并为通过遗传信息研发滞育剂、滞育解除剂，从而实现人工调控伞裙追寄蝇等天敌昆虫的寿命，实现在田间精准释放，发挥生物防治的最佳效果提供科学依据。

1 材料与方法

1.1 供试昆虫及滞育诱导

伞裙追寄蝇蛹为中国农业科学院草原研究所实验室所饲养。将正常发育的非滞育蛹在25℃±1℃室温条件下培养1 d，经用液氮处理后放入-80℃冰箱保存备用；将正常发育的非滞育蛹置于人工气候箱内，在温度11±1℃、相对湿度70%±10%、光周期8L∶16D的条件下，经过40 d诱导可获得滞育态蛹(姜珊, 2018)，40 d后选取重量与正常发育非滞育蛹一致的存活滞育蛹，经用液氮处理后放入-80℃冰箱保存备用。然后，将备用的伞裙追寄蝇正常发育非滞育蛹和滞育蛹进行RNA样品检测，主要方法为：(1)琼脂糖凝胶电泳分析RNA降解程度以及是否有污染；(2)NanoDrop检测RNA的纯度(OD260/OD280比值)；(3)Qubit对RNA浓度进行精确定量；(4)Agilent 2100精确检测RNA的完整性。

1.2 转录组测序和组装

伞裙追寄蝇非滞育和滞育蛹的转录组测序工作由北京诺禾致源科技股份有限公司完成。测序进行3个生物学重复，每个生物学重复以4头蛹为样本，两种蛹共构建6个文库。样品检测合格后，用带有Oligo(dT)的磁珠富集真核生物mRNA。随后加入fragmentation buffer将mRNA打断成短片段，以mRNA为模板，用六碱基随机引物(random hexamers)合成cDNA第1链，然后加入缓冲液，再加纯化的双链cDNA。纯化的双链cDNA先进行末端修复、加A尾并连接测序接头，再进行片段大小选择。最后进行PCR扩增，得到最终的文库。文库构建完成后，先使用Qubit2.0进行初步定量，稀释文库至1.5 ng/μL，随后使用Agilent 2100对文库的insert size进行检测，insert size符合预期后，使用Q-PCR方法对文库的有效浓度进行准确定量(文库有效浓度>2 nmol/L)，以保证文库质量。库检合格后，把不同文库按照有效浓度及目标下机数据量的需求pooling后进行高通量测序。测序完成后，高通量测序(Illumina HiSeqTM2000)得到的原始图像数据文件经碱基识别(base calling)分析转化为原始测序序列，然后进行测序错误率分布检查，获得高质量的序列(clean reads)后，采用Trinity软件对clean reads进行拼接，Corset参差聚类后得到最长的序列进行后续法分析。

1.3 生物信息学分析

为获得全面的基因功能信息，对伞裙追寄蝇进行了七大数据库的基因功能注释，包括：NCBI非冗余蛋白序列数据库(NCBI non-redundant protein sequences, Nr)，NCBI非冗余核苷酸序列数据库(NCBI non-redundant nucleotide sequences, Nt)，蛋白家族(protein family, Pfam)，直系同源基因簇(cluster of orthologous groups of proteins, KOG/COG), Swiss-Prot数据库(a manually annotated and reviewed protein sequence database), 京都基因与基因组百科全书(KEGG ortholog database, KO)，基因本体论(gene ontolooy, GO)。拼接得到的伞裙追寄蝇转录组测序结果与UniProt 或者Nr数据库进行Blastx比对(将核苷酸序列翻译为蛋白，再进行比对)，筛选条件E-value<1e-5，得到注释基因。通过RSEM软件得到每个样品比对到每个基因上的readcount数目，并对其进行FPKM(expected number of fragments per kilobase of transcript sequence per millions base pairs sequenced)转换，分析基因的表达水平；对转录组测序的非滞育蛹和滞育蛹的样品进行了样品间的相关性检查，采用R语言进行Pearson相关系数的计算；将正常非滞育蛹和滞育蛹的转录组测序结果进行比较，采用校正后的P值(p-adjusted, padj)与差异表达倍数(fold change, FC)方法进行差异表达基因挑选，padj是对P值进行了多重假设检验，能更好地控制假阳性率。对筛选出的差异表达基因进行聚类分析，使用聚类软件包pheatmap，针对差异表达基因的并集，以基因的相对表达水平值进行聚类。采用相应的距离算法计算每个基因之间的距离，然后通过反复迭代，计算基因之间的相对距离，最后根据基因的相对距离远近分成不同的subcluster，从而实现聚类。对伞裙追寄蝇非滞育蛹和滞育蛹差异表达基因的筛选条件为: padj<0.05且fold change≥32 或 padj<0.05 且fold change≤1/32。与非滞育蛹的相比，滞育蛹基因差异表达倍数在32倍以上的滞育关联基因，并进行重点分析。

1.4 滞育关联基因 KEGG通路分析

将滞育关联基因序列制成FASTA格式文件，应用KAAS(KEGG Automatic Annotation Server.http:∥www.genome.jp/kegg/)在线pathway比对分析工具进行KEGG通路富集分析。

2 结果

2.1 伞裙追寄蝇非滞育和滞育蛹转录组测序结果

本研究对伞裙追寄蝇非滞育和滞育蛹进行转录组测序，高质量样品数据超过7.42 Gb，每个样品的Q20均大于96.08%，Q30均大于90.01%(表1)，通过转录组数据结果发现其准确度较高，可用于后续分析。

表1 伞裙追寄蝇非滞育和滞育蛹转录组测序数据质量评估统计

样品间基因表达水平相关性是检验实验可靠性和样本选择是否合理的重要指标。样品间基因表达水平的相关性反映了样品间的相似程度，即不同处理或组织的样品在基因表达水平方面的相似度。相关系数越接近1，表明样品之间基因表达模式的相似度越高，样品间的差异表达基因也越少(张宇镭等, 2005)。本研究中生物学重复间的样品的相关系数大于生物学重复外的样品的相关系数(图1)，表明滞育蛹样品与非滞育蛹样品基因表达模式存在明显差异。

图1 伞裙追寄蝇非滞育和滞育蛹转录组样品间基因表达水平Pearson相关系数热图

差异表达基因聚类分析用于判断不同实验条件下差异基因表达量的聚类模式(杨春梅等, 2003)。图2红色表示高表达，蓝色表示低表达。本研究检测结果中蓝色和红色部分能够明显区分，说明样品整体的测序质量好。

图2 伞裙追寄蝇非滞育和滞育蛹转录组中差异表达基因聚类图

2.2 滞育关联基因 KEGG通路分析

伞裙追寄蝇非滞育蛹和滞育蛹转录组测序在KEGG直系同源系统(KEGG orthdogy, KO)中成功注释了6 643个基因，KO系统将各个KEGG注释系统联系在一起，可完成转录组的功能注释。根据筛选条件：padj<0.05且fold change≥2或padj<0.05且fold change≤0.5，获得在伞裙追寄蝇非滞育蛹和滞育蛹转录组间差异2倍以上的差异表达基因共有2 648个，滞育蛹中上调表达基因1 518个，下调表达基因1 130个。对这些差异表达基因进行KEGG分类(图3)，所属KEGG代谢通路分为5类，包括细胞过程、环境信息处理、遗传信息处理、代谢(metabolism)和有机系统，这些基因主要集中在信号转导、内分泌系统、细胞群落和碳水化合物代谢等途径。与伞裙追寄蝇非滞育蛹相比，滞育蛹转录组中差异在32倍以上的滞育关联基因为454个，可分为代谢、环境信息处理、遗传信息处理、细胞过程和有机系统五大类，每类中滞育关联基因参与的通路占比前3的如表2所示。其中上调基因为406个，下调基因为48个。通过KEGG数据库分析，发现454个滞育关联基因共参与到141个通路，其中参与焦点粘连通路的滞育关联基因最多，为13个；其次是参与碳代谢通路的，为12个。

表2 伞裙追寄蝇滞育蛹转录组中滞育关联基因KEGG通路分析

图3 伞裙追寄蝇滞育蛹转录组中与非滞育蛹相比差异表达基因的KEGG功能分类

2.2.1代谢:根据KO分类结果，滞育关联基因共参与代谢分类中的碳水化合物代谢(carbohydrate metabolism)、能量代谢(energy metabolism)、脂代谢(lipid metabolism)、核苷酸代谢(nucleotide metabolism)、氨基酸代谢(amino acid metabolism)等11个途径，涉及47个通路。

分析结果显示，在糖酵解/糖异生(glycolysis/gluconeogenesis)通路中，丙酮酸激酶(pyruvate kinase, PK)基因表达上调101.60倍，在柠檬酸循环(citrate cycle, TCA cycle)中，丙酮酸脱氢酶的E1组分(pyruvate dehydrogenase E1 component, PDHA)α亚基基因表达上调151.49倍，PK是糖酵解过程中发生第二次底物磷酸化的关键酶，催化反应生成丙酮酸。在有氧条件下，丙酮酸脱氢酶复合体(pyruvate dehydrogenase complex, PDHC)催化由糖酵解产生的丙酮酸进行氧化脱羧反应，促进乙酰辅酶A的生成，进入柠檬酸循环。在脂代谢中，脂酰辅酶A还原酶(fatty acyl-CoA reductase, FAR)基因表达上调，上调倍数最高可达146.49倍。在能量代谢中，参与氧化磷酸化的滞育关联基因有9个，氧化磷酸化是一个电子传递过程，对线粒体呼吸电子链的传递十分重要。细胞色素C氧化酶是线粒体中电子传递链的末端酶，催化电子从细胞色素C转移到氧。细胞色素氧化酶的5个亚基基因表达上调，其中细胞色素C氧化酶亚基4(cytochrome c oxidase subunit 4, Cox4)是1个关键亚基。在昆虫激素的生物合成(insect hormone biosynthesis)通路中，保幼激素酯酶(juvenile-hormone esterase, JHE)基因表达下调，下调了82.28倍。

2.2.2有机系统：根据KO注释结果，有机系统中的滞育关联基因参与免疫系统、内分泌系统、循环系统、排泄系统、神经系统、感觉系统以及消化系统等7个系统，其中有58个滞育关联基因参与内分泌系统，涉及内分泌系统的15个通路。

钙调蛋白(calmodulin, CALM)基因上调表达主要涉及内分泌系统，钙调蛋白本身无任何酶活性，可以通过结合Ca2+形成Ca2+-CaM复合体，调控细胞代谢过程以及下游的靶蛋白、靶酶。Hsp70的一种：heat shock 70kDa protein 1/8(HSPA1s)在内分泌系统的雌激素信号通路(estrogen signaling pathway)、免疫系统的抗原处理和展示(antigen processing and presentation)通路中既有上调也有下调，上调倍数最高可达98.53倍，下调倍数为2.56倍。

2.2.3细胞过程：根据KO注释结果，参与细胞过程的73个滞育关联基因主要涉及焦点粘连(focal adhesion)、肌动蛋白细胞骨架的调节(regulation of actin cytoskeleton)、紧密连接(tight junction)、缝隙连接(adherens junction)、吞噬体(phagosome)、溶酶体(lysosome)、细胞凋亡(apoptosis)、卵母细胞减数分裂(oocyte meiosis)、间隙连接(gap junction)、过氧化物酶体(peroxisome)、p53信号通路(p53 signaling pathway)、细胞周期(cell cycle)、胞吞作用(endocytosis)这13个通路。

在吞噬体通路中肌动蛋白基因表达占比最高，共有9个基因被注释为肌动蛋白β/γ1(actin beta/gamma 1, ACTB_G1)，肌动蛋白是真核生物中最丰富的蛋白质之一，对维持细胞形态具有重要的作用。在溶酶体中组织蛋白酶表达占比最高，共有6个基因分别被注释为组织蛋白酶G(cathepsin G, CTSG)、组织蛋白酶S(cathepsin S, CTSS)、组织蛋白酶K(cathepsin K, CTSK)。

在细胞凋亡通路中，细胞色素C(cytochrome c, CYC)基因表达上调，CYC介导呼吸链中的电子传递，影响ATP的产生，其释放可诱导细胞凋亡。在过氧化氢物酶体通路中，过氧化氢酶(catalase, CAT)基因表达上调，上调123.17倍。

2.2.4环境信息处理：根据滞育关联基因的KO比对结果，环境信息处理共涉及到信号转导(signal transduction)、信号分子及其相互作用(signaling molecules and interaction)、膜转运(membrane transport)3个途径。60个滞育关联基因共参与环磷酸腺苷(cAMP signaling pathway)，rap1信号通路(rap1 signaling pathway)，MAPK信号通路(MAPK signaling pathway)，钙信号通路(calcium signaling pathway)，cGMP-PKG信号通路(cGMP-PKG signaling pathway)，hippo信号通路(hippo signaling pathway)等在内的19个通路。在ECM-受体相互作用通路中，肌腱蛋白(tenascin, TN)基因表达下调，下调倍数最高可达51.98倍。

2.2.5遗传信息处理：根据KO注释结果，在遗传信息处理第二层次的分类中涉及复制与修复(replication and repair)，折叠、分类和降解(folding, sorting and degradation)，翻译(translation)和转录(transcription)4个途径，第3层次分类中涉及内质网蛋白加工(protein processing in endoplasmic reticulum)、RNA运输(RNA transport)、泛素介导的蛋白质水解(ubiquitin-mediated proteolysis)等7个通路。泛素蛋白连接酶E3 XIAP连锁凋亡抑制蛋白(E3 ubiquitin-protein ligase XIAP, XIAP)基因，ring结构域蛋白1(ring-box protein 1, RBX1)基因共同参与了泛素介导的蛋白质水解这一通路，RBX1基因表达上调，上调293.29倍；XIAP基因表达下调，下调42.81倍。

3 讨论

人工诱导伞裙追寄蝇蛹进入滞育，可以有效延长其货架期。本研究以伞裙追寄蝇非滞育和滞育蛹转录组为研究对象，通过对比分析，筛选出454个滞育关联基因。这些基因主要集中在信号转导、内分泌系统、碳水化合物代谢等途径。

在伞裙追寄蝇滞育期间，丙酮酸激酶基因、丙酮酸脱氢酶的E1组分α亚基表达上调，有利于中间产物——乙酰辅酶A生成，然后进入TCA循环(Naitoetal., 1998)。细胞色素氧化酶相关基因的表达上调，利于呼吸电子链的传递，为昆虫的滞育期提供能量。由此可以推测糖酵解、柠檬酸循环、氧化磷酸化在伞裙追寄蝇滞育期间的能量供给并不是维持在一个低水平。在脂代谢中，FAR催化脂酰辅酶A还原成脂肪醇，这是生成脂肪醇的主要途径。关于FAR的研究主要集中在鳞翅目昆虫信息素，在家蚕Bombyxmori和巢蛾属Yponomeuta的研究中表明一种FAR催化性腺中多种性信息素，而欧洲玉米螟Ostrinianubilalis两个亚种的研究则表明2个FAR基因生成不同的性信息素(Motoetal., 2003; Wahlenetal., 2009; Lassanceetal., 2010; 杨璞等, 2012)。JHE由脂肪体细胞和上皮细胞产生，是主要的保幼激素(juvenile hormone, JH)特异性降解酶之一，裂解JH酯键产生甲醇和JH酸，在调节JH代谢中发挥重要作用(李胜等, 2004)。JHE基因表达下调，由此可以推测滞育期间伞裙追寄蝇体内保幼激素降解受到阻碍，JH含量可能不是保持在一个低水平。

肌腱蛋白属于大细胞外基质糖蛋白家族，在神经系统发育过程中以特定模式表达，并在成熟后下调。大量的体外研究表明，肌腱蛋白在神经系统中的作用与前体细胞迁移，轴突生长和引导有关(Joester and Faissner, 2001)。可以推测,在滞育期间伞裙追寄蝇蛹的肌腱蛋白可能处于发育成熟阶段。

过氧化氢酶基因在过氧化氢物酶体通路发挥着保护酶的作用，使伞裙追寄蝇在滞育期间抵抗逆境胁迫，避免受到损伤(于德玲和王昌留, 2016)。在菜蛾盘绒茧蜂Cotesiavestalis滞育期间，发现低浓度的过氧化氢和高活性的CAT有利于增强该蜂的抗逆能力(郝仲萍等, 2012; 胡甘雨等, 2014)。

细胞凋亡是一种细胞程序性死亡过程，由许多过程相互关联，共同调控。在伞裙追寄蝇滞育期间，XIAP基因、RBX1基因共同参与泛素介导的蛋白质水解通路，RBX1是泛素连接酶E3的一个催化亚基，可以控制许多短寿命调节蛋白的更新(Jiaetal., 2011; 王步勇等, 2013)。在秀丽隐杆线虫Caenorhabditiselegans的胚胎发育，生殖细胞丰度以及成年活力或寿命具有重要的调控作用(Moore and Boyd, 2004)。XIAP是具有抑制细胞凋亡的蛋白，RBX1与XIAP发挥了拮抗作用。RBX1具有泛素连接酶E3活性，可以引起XIAP自身泛素化而被降解，同时也能够引起底物发生泛素化，抑制细胞凋亡的线粒体途径，发挥抗凋亡的作用(李中海等, 2003)。

本研究通过对伞裙追寄蝇非滞育蛹与滞育蛹的转录组差异表达基因KEGG富集通路的分析比较，揭示了其滞育的关键基因，为阐明其滞育机理奠定了一定的理论基础。关键滞育关联基因在伞裙追寄蝇滞育过程中的差异性表达及作用，还有待于进一步验证与探究。