伞裙追寄蝇滞育关联基因和滞育关联代谢物的联合分析

李媛媛，德力格尔，张 博，林 晨，石 凯*

(1.内蒙古民族大学农学院，内蒙古通辽 028043；2.中国农业科学院草原研究所，呼和浩特 010010；3.内蒙古师范大学生命科学与技术学院，呼和浩特 010022)

伞裙追寄蝇Exoristacivilis隶属于双翅目寄蝇科，在我国北京、内蒙古、河北等地区都有分布(赵洪有等,1985)，是一种牧草优势寄生性天敌，能寄生于多种鳞翅目害虫。但是该天敌产品的存放时间较短，并且运输也会对天敌产品造成一定的伤害。为了减少对天敌产品的损害，延长产品的货架期，人工诱导天敌滞育是解决以上问题的重要方式，能够帮助昆虫度过不利环境，维持个体生存，保证种群繁衍(徐卫华,2008)。滞育是一种应对不利环境出现之前的一种遗传性适应，滞育的发生受外界和内部环境因子的调控，固定在某一虫态发育停滞的现象，一旦昆虫进入滞育状态，即使给以适生条件亦不复苏(Denlinger and Armbruster,2014)。滞育为昆虫提供了巨大的适应性优势，允许昆虫在不适合的季节性环境中继续生存，使其生命周期与适合生长、发育和繁殖的周期同步，有助于发挥生物防治在害虫管控中的作用(Handetal.,2016)。滞育过程包括诱导、准备、启动、维持、终止以及滞育后发育这6个阶段(Kostal,2006)。关于伞裙追寄蝇滞育调控，目前开展的研究主要集中于：滞育诱导的温光周期、低温贮藏条件、主要物质变化含量以及滞育关联蛋白等方面。韩海斌等对伞裙追寄蝇滞育诱导的温光周期、低温贮藏条件进行了相关研究，明确了伞裙追寄蝇的滞育虫态是蛹，温度11±1℃，光周期L ∶D=8 ∶16(韩海斌等,2018)。孙程鹏等对伞裙追寄蝇滞育期间的主要物质变化做了研究，比较分析了滞育与非滞育虫态的糖原、海藻糖、蛋白质、甘油含量变化，认为糖原是滞育期间的主要抗寒物质(孙程鹏等,2017)。

转录组学是从RNA水平研究基因表达的情况，是研究细胞功能的重要手段(Griffiths and Wang,2009)；代谢组学作为其延伸，是通过小分子代谢物的变化得出的规律来揭示其参与生命活动的机制。代谢组学为转录组学提供了更加表观显著的数据支持(Dunnetal.,2011)。对两者进行联合分析，是因为细胞内的生命活动由基因、蛋白质、以及小分子代谢产物来共同承担，大分子的功能性变化最终会体现于代谢层面，基因表达在功能水平上的微小变化会在代谢物上得到放大，从而使检测更容易(Tengetal.,2009)。除此之外，许多基因的非功能性变化不会在代谢物上反映出来，看不到上游信息向下游信息传递过程中的功能性变化。

目前对家蚕Bombyxmori、七星瓢虫Coccinellaseptempunctata、大豆食心虫Leguminivoraglycinivorella、白纹伊蚊Aedesalbopictus、葱蝇Deliaantigua等昆虫均开展了滞育相关的转录组学研究(Qietal.,2015;蒋涛,2017;任爽,2018;胡甜,2019;李新畅,2020;陈艳荣等,2021)，而利用代谢组学开展的滞育研究相对较少，主要包括柑橘大实蝇Bactroceraminax、丽蝇蛹集金小蜂Nasoniavitripennis(李玉艳,2015;王攀,2020)，但是尚缺乏结合多种组学对昆虫进行滞育的研究。本文通过转录组和代谢组对伞裙追寄蝇滞育蛹和非滞育蛹的差异表达基因和差异代谢物进行联合分析，为明确其滞育的分子机制奠定理论基础。

1 材料与方法

1.1 供试昆虫

伞裙追寄蝇是在中国农业科学院草原研究所沙尔沁基地采集获得的，以20%蜂蜜水喂养；寄主粘虫Mythimnaseparata为中国农业科学院草原所实验室饲养，以人工饲料喂食，其幼虫主要在长为20 cm，高为10 cm的塑料养虫盒中饲养，待粘虫长至5龄，转入长为35 cm，宽为30 cm，高为10 cm亚克力养虫盒饲养，同时放入交配好的寄蝇。每盒粘虫与交配后的寄蝇的释放比例为50 ∶1，在温度25±1℃，相对湿度70%±1%，光周期16 L ∶8 D的环境条件下进行饲养，在亚克力盒中放入适量人工饲料，在盒中悬挂20%蜂蜜水的棉球。待寄蝇成功寄生后，将获得寄生的末龄粘虫放入装有灭菌土的养虫圆盒中饲养，每日喷清水，保持灭菌土的水分。

1.2 样品收集

根据实验室前期的研究基础，获得了伞裙追寄蝇滞育诱导的温度11±1℃、光周期8 L ∶16 D、滞育的敏感虫态为蛹。寄蝇成功寄生粘虫后，在温度25±1℃，光周期16 L ∶8 D，相对湿度 70%±10%环境条件下饲养，直至寄蝇3龄幼虫出现。将一定量的寄蝇3龄幼虫放入设置条件：温度11±1℃，光周期8 L ∶16 D，相对湿度 70%±10%、光照强度为8 800 Lx人工气候箱中持续诱导40 d获得滞育蛹，将滞育蛹经过液氮处理，放入-80℃冰箱保存；将另一部分寄蝇3龄幼虫放入：温度25±1℃，光周期16 L ∶8 D，相对湿度70%±10%条件继续饲养，直至伞裙追寄蝇蛹出现，收集非滞育蛹作为对照组，将非滞育蛹经过液氮处理，放入-80℃冰箱保存。

1.3 转录组学测序及差异表达基因的筛选

伞裙追寄蝇非滞育和滞育蛹的转录组测序工作由北京诺禾致源科技股份有限公司完成。采用3个生物学重复，每个生物学重复以2头蛹为样本。本研究为无参考基因组的转录组分析，测序完成后，过滤raw reads，去掉低质量的reads后，将获得的clean reads采用Trinity(Grabherretal.,2011)拼接成转录本，利用Corset程序对转录本进行层次聚类，获得最长的cluster序列进行后续的分析(Davidson and Oshlack,2014)。

采用DESeq2对输入的read count数据进行分析，设置筛选域值为：(1)padj<0.05；(2)|log2 FoldChange|>1；对大量的基因进行独立的统计假设检验，得到的Pvalue进行校正，P-adjusted为校正后的Pvalue，最后获得转录组测序中的差异基因，P-adjusted越小，表示基因表达差异越显著。

1.4 代谢组学检测及差异代谢物的筛选

伞裙追寄蝇非滞育和滞育蛹的代谢组测序采用6个生物学重复，每个生物学重复以2头蛹为样本。基于液质联用技术进行非靶向代谢组学研究，基于质谱检测得到的Raw文件，先将原始文件导入CD软件中，进行谱图处理并进行数据库搜索，然后对数据进行质控保证数据结果的准确度、可靠性，代谢物分子式预测使用mzCloud数据库进行比对。代谢物主要进行多元统计分析包括主成分分析(Principal component analysis，PCA)、偏最小二乘法分析(PLS-DA)等揭示不同比对组别代谢物组成的差异。利用层次聚类和代谢物-代谢物相关性分析揭示了代谢物和样本之间的关系。最后，通过代谢通路等功能分析发现代谢物相关的生物学意义。

采用PLS-DA模型第一主成分的变量投影重要度(Variable Importance in theProjection，VIP)值，并结合T-test的P值来寻找差异性表达代谢物，设置阈值为VIP>1.0，差异倍数FC>2.0或FC<0.5且Pvalue<0.05。

1.5 转录组和代谢组联合分析

1.5.1差异基因和差异代谢物的相关性分析

基于皮尔森相关系数对伞裙追寄蝇显著差异基因、差异代谢物进行相关性分析，以此来衡量差异基因与代谢物之间的关联程度。当相关系数<0时，为负相关；当相关系数>0时，为正相关。展示Top 50的差异代谢物(按Pvalue从小到大排序)和Top 100的差异基因(按Pvalue 从小到大排序)；若差异代谢物数目<50或差异基因数目<100，则展示所有的差异代谢物或所有的差异基因。

1.5.2差异基因与差异代谢物的KEGG数据库映射

将所有差异表达基因与差异表达代谢物同时向KEGG映射，获得两者共同的pathway信息，明确差异表达基因与差异代谢物共同参与的主要生化途径和信号转导途径。

2 结果与分析

2.1 测序与组装

滞育组中的原始序列为176 530 572，非滞育组中的原始序列为169 609 048；过滤后，获得滞育组171 785 530 clean reads，获得非滞育组165 045 546 clean reads。经过拼接，共获得137 725转录本，其平均长度为1 367 bp，然后转录本经过进一步组装，共获得58 050 unigenes(如图1)。

图1 拼接转录本与基因序列长度分布图Fig.1 Transcripts with unigenes sequence length distribution

2.2 差异表达基因的筛选

通过筛选条件padj<0.05且|log2FoldChange|>1，确定了伞裙追寄蝇滞育蛹和非滞育蛹的差异表达基因(differentially expressed genes，DEGs)个数为7 513个，其中4 355个上调表达，3 158个下调表达(图2)。而经过KEGG数据库注释成功的基因共有2 642个，映射到223条通路中，被分成代谢、遗传信息加工、环境信息处理、细胞过程和有机体系统这5个类别。

图2 伞裙追寄蝇滞育和非滞育蛹的差异表达基因火山图Fig.2 Volcanoplot of differentially expressed gene between D and ND pupaes of Exorista civilis

2.3 样品检测及差异代谢物筛选

经检测，伞裙追寄蝇非滞育(non-diapause,ND)与滞育(diapause,D)蛹的总样品，无论在正离子模式下还是在负离子模式下，在PCA图上，合格的质量控制(quality control,QC)样本，用于测定进样前仪器状态及平衡色谱-质谱系统，并用于评价整个实验过程中系统稳定性。差异较小，呈现聚集分布，说明对伞裙追寄蝇非滞育、滞育蛹样品的分析方法稳定性较好，数据质量较高(图3-A,B)。

图3 伞裙追寄蝇非滞育与滞育蛹总样品的主成分分析Fig.3 Principal component analysis of total samples of non-diapause and diapause pupae of Exorista civilis 注：A，负离子模式；B，正离子模式。图中横坐标PC1和纵坐标PC2分别表示排名第一和第二的主成分的得分，散点颜色表示样本的实验分组，置信椭圆为95%。Note:A,Negative model;B,Positive model.In the figure,the abscissa PC1 and the ordinate PC2 represented the scores of the first and second principal components respectively,and the color of the scatter indicated the experimental groups of the samples,with a confidence ellipse of 95%.

在两种模式下共检测获得差异代谢物501个，在负离子模式下，筛选出差异代谢物169个，其中99个上调，70个下调。在正离子模式下，筛选出332个差异代谢物，其中268个上调，64个下调。

2.4 差异基因与差异代谢物表达相关性分析

图4中聚类枝越短表示差异表达基因或差异代谢物的相似性越高。蓝色表示负相关，红色表示正相关。通过分析发现基因与代谢物是正相关的是：在负离子模式下，注释为肌动蛋白β/γ1，肌球蛋白重链，丙烯基硫酸酶B，原肌球蛋白1等80个差异基因与注释为L-组氨酸等22个差异代谢物呈正相关，注释为酪氨酸酶，乙酰胆碱酯酶等20个差异基因和注释为2-羟基戊二酸、6-胍基己酸、N-乙酰-DL-谷氨酸、高肌肽、肌肽等28个差异代谢物呈正相关。注释为肌动蛋白β/γ1，肌球蛋白重链，丙烯基硫酸酶B，原肌球蛋白1等80个差异基因和注释为2-羟基戊二酸、6-胍基己酸、N-乙酰-DL-谷氨酸、高肌肽、肌肽等28个差异代谢物呈负相关，注释为酪氨酸酶，乙酰胆碱酯酶20个差异基因和注释为2-羟基戊二酸、6-胍基己酸、N-乙酰-DL-谷氨酸、高肌肽、肌肽等28个差异代谢物呈负相关。

图4 伞裙追寄蝇滞育和非滞育蛹差异基因与差异代谢物表达相关性分析热图Fig.4 Heat map of correlation analysis between DEMs and DEGs of Exorista civilis注：A，负离子模式；B，正离子模式。纵向表示差异表达基因的聚类，横向表示差异代谢物的聚类。Note:A,Negative ion mode;B,Positive ion mode.Vertical axis indicated clustering of differentially expressed genes and horizontal axis indicated clustering of differential metabolites.

在正离子模式下，肌球蛋白重链、丙烯基硫酸酶B、肌动蛋白β/γ1，钙调蛋白等80个差异基因和注释为L-苯丙氨酸、5-甲基胞嘧啶、D-(+)-脯氨酸、1-甲基-L-组氨酸等39个差异代谢物呈正相关，酪氨酸酶，乙酰胆碱酯酶等20个差异基因和DL-二棕榈酰磷脂酰胆碱等11个差异代谢物呈正相关，酪氨酸酶，乙酰胆碱酯酶等20个差异基因和L-苯丙氨酸、5-甲基胞嘧啶、D-(+)-脯氨酸、1-甲基-L-组氨酸等39个差异代谢物呈负相关，80个差异基因和DL-二棕榈酰磷脂酰胆碱等11个差异代谢物呈负相关。

2.5 差异基因与差异代谢物pathway分析

在负离子模式下，差异表达基因与差异代谢物在氨基酸代谢、消化系统、免疫系统、神经系统方面关联性较高，差异表达基因和差异表达代谢物共富集到31条通路，其中差异基因和代谢物富集个数在2种以上的通路为16条，富集的差异代谢物主要为：UDP-N-乙酰葡糖胺、D-葡糖胺-6-磷酸盐、肌肽、L-组氨酸、组胺、花生四烯酸、顺式乌头酸、柠檬酸、高肌肽、4-氧代脯氨酸、磷酸精氨酸、L-脯氨酸、L-苯丙氨酸、赤藓糖醇、UDP-N-乙酰葡萄糖胺、D-氨基葡萄糖6-磷酸。其中，柠檬酸循环中，富集到28个差异表达基因，2种差异代谢物；在氨基糖和核苷酸糖代谢(Amino sugar and nucleotide sugar metabolism)中，富集到17个差异表达基因，2种差异代谢物；在精氨酸和脯氨酸代谢(Arginine and proline metabolism)中，富集到21个差异表达基因，5种差异代谢物。

在正离子模式下，差异表达基因与差异代谢物在消化系统、信号转导途径中关联性较高。DEGs和DEMs共富集到39条通路，其中差异基因和代谢物富集个数在2种以上的通路为15条。富集的差异代谢物主要为：L-苯丙氨酸、钙三醇、L-组氨酸、L-酪氨酸、乙酰胆碱、腺苷-5′-单磷酸盐、1-甲基组氨酸、苯甲酸。在两种离子模式下，差异代谢物为氨基酸所占的比重最多。其中，在叉头转录因子通路中，富集了11差异表达基因，1种差异代谢物；在cAMP信号通路中，富集到31个差异表达基因，1种差异代谢物(如图5)。

图5 伞裙追寄蝇差异基因与差异代谢物KEGG通路关联分析Fig.5 Association analysis of differential genes and differential metabolite KEGG pathway in Exorista civilis注：A，负离子模式；B，正离子模式。横坐标为该通路中富集到的差异代谢物或差异基因与该通路中注释到的代谢物或基因个数的比值(Ratio)，纵坐标为代谢组-转录组共同富集到的KEGG通路。数量，通路中富集的代谢物或基因的个数。Note:A,Negative ion mode;B,Positive ion mode.The horizontal coordinate was the ratio of differential metabolites or differential genes enriched in the pathway to the number of metabolites or genes annotated in the pathway,and the vertical coordinate was the metabolome-transcriptome co-enriched to the KEGG pathway.Count,the number of metabolites or genes enriched in the pathway.

通过对伞裙追寄蝇滞育关联基因和滞育关联代谢物的联合分析，发现在正离子模式下，cAMP通路包括的31个DEGs，其中22个DEGs上调表达，主要有蛋白激酶A、钙调蛋白；9个DEGs下调表达，主要有RAC1。DEMs主要为腺苷-5′-单磷酸，乙酰胆碱。差异代谢物腺苷-5′-单磷酸(Adenosine 5′-monophosphate,AMP)还参与调控胰岛素信号途径、FOXO信号通路、糖酵解，进而影响到钙调蛋白、过氧化氢酶(Catalase,CAT)、超氧化物歧化酶(superoxide dismutase,SOD)、磷酸烯醇式丙酮酸(phosphoenolpyruvate carboxykinase,PEPCK)的表达(如图6)。在负离子模式下，发现柠檬酸循环中，获得KEGG成功注释的差异表达基因为28个，上调的差异表达基因为27个，下调表达的差异基因为1个，即PEPCK，下调倍数为9.43倍，在代谢组学中发现滞育关联代谢物为柠檬酸，下调表达2.97倍。

图6 伞裙追寄蝇滞育关联基因和滞育关联代谢物富集通图Fig.6 Enrichment pathways for diapause-associated genes and diapause-associated metabolites of Exorista civilis注：红方框表示差异基因上调，绿方框表示差异基因下调，红色椭圆表示差异代谢物上调。Note:Red boxes indicated differential gene up-regulated,green boxes indicated differential gene down-regulated,and red ovals indicated differential metabolite up-regulated.

在氨酰-tRNA的生物合成中，苯丙氨酸在正、负离子两种模式下均表达上调，上调倍数最高约为2.80倍，苯丙氨酰-tRNA合成酶α链(phenylalanyl-tRNA synthetase alpha chain,pheS)表达上调，上调倍数为2.37倍；另外在正离子模式下，酪氨酸表达上调，上调倍数2.44倍。

但是在两种离子模式下，未发现同时在转录组和代谢组中显著富集的通路，在转录组和代谢组中任一存在显著性富集的通路见表1。在负离子模式下，显著性富集的通路主要有柠檬酸循环、矿物质吸收、癌症的中枢碳代谢、精氨酸和脯氨酸代谢、蛋白质的消化和吸收、ABC转运子、组氨酸代谢、精氨酸和脯氨酸代谢，在正离子模式下，主要有氨基乙酰-tRNA的生物合成、组氨酸代谢、癌症的中枢碳代谢。

3 结论与讨论

在双翅目昆虫体内，蜕皮激素的合成不以cAMP为第二信使，而是依赖于钙离子，钙离子和钙调蛋白结合，在家蚕、烟草天蛾Manducasexta的研究中发现，cAMP是活化PTTH的第二信使(Hua and Koolman,1995)，导致依赖于cAMP的蛋白激酶活化，最终导致蜕皮激素的合成和分泌的增加(王荫长等,2004)。在正离子模式下，在伞裙追寄蝇cAMP信号途径中，cAMP转化为差异代谢物AMP，AMP上调表达，其作为中间产物参与胰岛素信号途径、FoXO信号通路、糖酵解途经的代谢过程。本研究中AMP通过调控FoXO，FoXO作为糖酵解的上游基因，对PEPCK产生影响，PEPCK活性主要反映在转录组水平上，其表达水平直接对糖异生产生影响，在调控机体代谢中发挥着重要的作用(Hanson and Reshef,1997)。研究表明胰岛素信号途径主要通过促进糖原或海藻糖的合成，抑制糖异生的作用(彭竹清和郝友进,2019)。在伞裙追寄蝇的转录组测序中，确实发现PEPCK表达下调，作为糖异生重要限速酶，抑制了糖异生的第一步，说明糖异生途径受到抑制。但与之不同的是肉蝇Sarcophagaperegrina和果蝇Drosophilamelanogaster在滞育期间，PEPCK基因的表达显著增强；滞育解除后，其表达量明显下降(Hahn and Denlinger,2011)。在苹果实蝇Rhagoletispomonella滞育后期，代谢由抑制状态恢复的24 h内，PEPCK表达量下降(Raglandetal.,2011)。并且与白纹伊蚊、麻蝇Sarcophagacrassipalpis滞育期间的PEPCK表达趋势也不同，其PEPCK在滞育期间表达量较为丰富(Raglandetal.,2010;Poelchauetal.,2011)。然而伞裙追寄蝇PEPCK表达结果与丽蝇蛹集金小蜂在滞育期间的表达趋势相同(高飞等,2020)。胰岛素作为重要的内分泌激素，其通路还能够参与果蝇等一些昆虫的滞育诱导调控。Tatar等(2001)干扰果蝇胰岛素信号后，发现果蝇繁殖停止、体内能量储备开始增加。在尖音库蚊Culexpipiens、家蚕滞育的研究中发现，胰岛素信号通路不独立对昆虫滞育调控产生作用，而是作为中介接受上游通路信号，或是将信号传递给下游内分泌通路(梁瀚清等,2014)。这些结果与本研究胰岛素信号通路参与的滞育调控方式具有一致性。

在负离子模式下，柠檬酸循环中，除了PEPCK表达下调外，其他差异表达基因均表达上调，在糖酵解中，己糖激酶表达上调，促进6-磷酸葡萄糖进入糖酵解，通过上调表达的丙酮酸脱氢酶促进丙酮酸转化为乙酰辅酶A，从而进入TCA循环。在该通路中，差异代谢物顺式乌头酸、柠檬酸均为表达下调，这个下调表达趋势与差异表达基因PEPCK一致，虽然不能明确其下调原因，但是可以推测PEPCK对糖异生起着抑制的作用，而糖酵解途径是被激活的。

组成蛋白质分子的20种氨基酸在合成蛋白质之前，均必须活化以获取能量。活化反应是在专门的氨酰-tRNA合成酶催化下进行的。在氨酰-tRNA的生物合成中，苯丙氨酸表达上调，在西藏飞蝗Locustamigratoriatibetensis抗寒物质的研究中发现，低温驯化和短光照条件下，其苯丙氨酸也出现累积的情况(朱昱翰,2014)，这与伞裙追寄蝇滞育期间苯丙氨酸的代谢组结果具有一致性。另外，黑色素源于酪氨酸，苯丙氨酸是合成酪氨酸的前体物质，伞裙追寄蝇滞育期间酪氨酸表达上调，表明伞裙追寄蝇滞育期间与苯丙氨酸相关的蛋白质合成在活化过程中是较为活跃的。在家蚕幼虫斑纹黑色素合成的主要来源是苯丙氨酸羟化酶催化苯丙氨酸生成的酪氨酸(陈萍,2011)，推测苯丙氨酸以及相关基因的上调表达，可能对伞裙追寄蝇滞育期间蛹表皮的黑色素沉积有着重要的影响。

然而目前对于伞裙追寄蝇滞育过程，即诱导、准备、启动、维持、终止时期，起到关键作用的特异性表达基因仍不清楚，滞育中某一过程的滞育关联基因有待于进一步研究。