低温胁迫下意大利蜜蜂预蛹的差异表达基因趋势分析

姚 丹, 周姝婧, 周冰峰, 徐新建, 朱翔杰

(福建农林大学动物科学学院(蜂学学院), 福州 350002)

蜜蜂不仅为人类提供了丰富多样的蜂产品，而且是农业中最重要的授粉昆虫，在维护生态平衡发挥重要作用(Allen-Wardelletal., 1998)。近年来，全球气候变化和病虫害危害导致蜜蜂的种群数量大幅下降(Milbrathetal., 2015)，蜜蜂蜂群的减少引起了人们对蜜蜂健康的极大关注。

环境温度是影响昆虫分布、形态、生理、行为等的主要生态因素之一(Burnett, 1949; Ponsonby and Copland, 1996)。 独居性昆虫可通过休眠或滞育等方式耐受低温等不良环境，生态幅一般较宽(Tauber and Tauber, 1976)。而蜜蜂的卵、虫、蛹(统称为蜂子)是在蜂群维持较稳定的巢温下发育的，因此形成了狭温性特点，发育温区为29～38℃，最适发育温度为35℃(Tautzetal., 2003; Grohetal., 2004; 朱翔杰等, 2006)。偏离最适温度，羽化后成蜂的外部形态(Medrzyckietal., 2010)、翅脉发育(周冰峰等, 2011; Zhuetal., 2018)、神经突触数量(Grohetal., 2004)、寿命(Wangetal., 2016)和行为(Jonesetal., 2005; Becheretal., 2009)等均受到影响。蜜蜂发育对温度的敏感性，使得蜂子成为研究温度对个体发育影响的理想材料，有助于认识温度与生命的关系。

低温对昆虫的影响主要表现在糖代谢、脂类代谢、氨基酸代谢、嘌呤代谢和硫胺素代谢等代谢通路上富集大量差异表达基因(differentially expressed genes, DEGs)(库尔班·吐松等, 2016; Tangetal., 2017; Wangetal., 2017; Xuetal., 2017; Qinetal., 2019)。在0℃下20日龄中华蜜蜂Apisceranacerana成蜂转录组测序检测到与环境进程相关的Ca2+信号通路、FoxO信号通路、Hippo信号通路等，同时热激蛋白、锌指蛋白和丝氨酸/苏氨酸蛋白激酶等差异表达(Xuetal., 2017)。

相对其他日龄的蜜蜂封盖子，3日龄封盖子处于幼虫-蛹变态的过渡期，表现出对低温胁迫最敏感。前期研究发现，不同日龄封盖子经20℃低温胁迫后3日龄封盖子最早出现大量封盖子死亡，羽化蜜蜂的寿命也最短(Wangetal., 2016)。到目前为止，尚无研究对蜜蜂封盖子发育期间低温敏感的机制给出解释。为了深入研究低温对蜜蜂化蛹影响的分子机制，本研究采用转录组学方法对低温处理不同时间的意大利蜜蜂Apismelliferaligustica3日龄封盖子的共有差异表达基因进行趋势分析，从分子水平探讨蜜蜂化蛹时受低温影响的关键基因和关键通路。本研究对丰富昆虫生态学和昆虫发育生物学有重要意义。

1 材料与方法

1.1 供试蜜蜂

样本取自春季增长阶段的意大利蜜蜂A.m.ligustica蜂群，通过限制蜂王产卵空间，获得卵龄一致的蜜蜂卵。3 d后将卵脾放入同一个哺育群哺育。移入哺育群5 d后这些幼虫即将封盖，获取4 h内同时封盖的封盖子(封盖巢房内的幼虫和蛹)，割取样本巢脾，放入恒温恒湿箱(35±0.2℃, RH 75%)[CTHI-250B, 施都凯仪器设备(上海)有限公司]中培育至3日龄(对照组，CK)。再将一部分3日龄封盖子放入温度20±0.2℃、RH 75%(Wangetal., 2016)的恒温恒湿箱中分别处理18 h和36 h(分别标记为T18和T36)，作为低温处理组。所取的封盖子样本立即用液氮冷冻，放入-80℃保存备用。每个样本用30只封盖子(3个蜂群，10头/群)混合建库和高通量测序。

1.2 高通量测序及测序数据质量评估

利用Trizol法提取1.1节样品全虫总RNA，用带有Oligo(dT)的磁珠富集mRNA后，加入fragmentation buffer使其片段成为短片段的mRNA为模板。用六碱基随机引物(random hexamers)合成cDNA第1链，并加入缓冲液、dNTPs、RNase H和DNA polymerase I合成cDNA第2链，经过QiaQuick PCR试剂盒纯化并加EB缓冲液洗脱经末端修复、加碱基A，加测序接头，再经琼脂糖凝胶电泳回收目的大小片段，并进行PCR扩增，从而完成整个文库制备工作，构建好cDNA文库，高通量测序委托广州基迪奥生物科技有限公司完成，测序平台为Illumina HiSeqTM。对下机数据的raw reads进行过滤，去除含adapter的reads、含N比例大于10%的reads和低质量的reads，最终得到clean reads，然后利用Tophat2(2.1.1)软件将reads比对到西方蜜蜂Apismellifera的参考基因组(ncbi_GCF_000002195.4)上，并利用Cufflinks组装转录本。 转录组原始数据已上传NCBI的SRA数据库(获取号: SRR10912937-SRR10912945)。

1.3 差异表达基因趋势分析

使用DESeq2软件(V1.20.0)对有生物学重复的处理组与对照组间(T18vsCK, T36vsCK)基因表达量进行差异分析(Loveetal., 2014)， DEGs筛选条件为FDR(错误发现率)≤0.05, |log2Fold change| (Fold change, FC, 指基因表达水平的变化倍数, FC=T/CK，其中T为处理组的基因表达水平， CK为对照组的基因表达水平)≥1。利用STEM(Short Time-series Expression Miner)，选择参数{-pro 20 -ratio 1.0[log2(2)=1, log2(1.5)=0.5849625, log2(1.2)=0.2630344]}(Ernst and Bar-Joseph, 2006),对2个比较组间共有的差异表达基因进行趋势分析，在软件中设置P≤0.05为差异显著的表达模式。

1.4 差异表达基因GO及KEGG pathway富集分析

对各个趋势中的基因进行GO(Gene Ontology)和KEGG(Kyoto Encyclopedia of Genes and Genomes)代谢通路进行功能注释(Omicshare tools, http:∥www.omicshare. com/tools/home/index/index.html)，并通过假设检验计算得到P值。得到的P值通过FDR校正之后，以P值≤0.05为阈值，满足此条件的GO条目(term)和通路(pathway)定义为在该趋势中显著富集的GO term和Pathway。

1.5 RT-qPCR基因验证

随机挑选5个趋势显著的基因进行RT-qPCR验证,β-actin为内参基因(引物见表1)。利用RNA抽提试剂盒提取样本总RNA，利用反转录试剂盒合成cDNA第1链(北京全式金生物技术有限公司)。以上述的cDNA作为模板进行qPCR，反应体系(20 μL): H2O 6.6 μL, 2×TransStart SuperMix 10 μL, 上下游引物(5 μmol/L)各1 μL, cDNA模板1 μL,Dye II 0.4 μL。qPCR反应条件: 95℃预变性4 min; 95℃变性15 s, 60℃退火和延伸45 s，共40个循环；熔解曲线默认系统程序(QuantStudio5, 美国)。样本取自1.1节，对照组和处理组的封盖子均来自3个蜂群，每个蜂群取1头进行混合提取总RNA。每个样本设置3次目的基因的重复扩增。

表1 RT-qPCR引物信息Table 1 Primers for RT-qPCR

1.6 数据分析

目的基因的RT-qPCR表达水平按照内参基因的2-ΔΔCt法进行计算，即ΔCt=Ct目的基因-Ct内参基因(Livak and Schmittgen, 2001)，各基因的相对表达量经log2(X)转换后与高通量测序中的表达水平进行比较。应用SPSS V21.0软件对每个基因处理组与对照组之间平均相对表达量的差异显著性进行非配对T检验。

2 结果

2.1 测序结果质量

通过对CK, T18和T36测序，分别得到原始读段平均为56 759 342, 45 977 301和54 196 127，过滤后得到高质量有效读段平均分别为55 578 624, 45 032 615和53 121 097。统计两端的Q20和Q30分别平均在98.20%和94.85%以上。表明本研究测序数据质量良好，可用于进一步分析。

2.2 差异表达基因趋势

通过差异表达基因(DEGs)分析，T18vsCK有1 394个DEGs，T36vsCK有2 722个DEGs，说明随着低温处理的时间增长，受影响的基因数明显增加。其中两比较组共有DEGs 1 062个，T18vsCK特有DEGs有332个，T36vsCK特有DEGs有1 660个，上调趋势的DEGs显著多于下调趋势的DEGs。

对两比较组共有的1 062个DEGs进行趋势分析，在8个基因表达模式中有3个为显著基因表达模式(P<0.05)，包括2个上调表达模式，Profile 6模式为基因表达量先增加后保持不变，有539个基因；Profile 7模式为基因表达量随低温胁迫时间延长一直上调，有271个基因；1个下调表达模式Profile 1，为随低温胁迫时间延长基因表达量先下调表达，之后保持不变，有183个基因(图1)。多数基因归入上述3个显著基因表达模式。

图1 低温胁迫不同时间3日龄意大利蜜蜂封盖子差异表达基因表达趋势Fig. 1 Trend analysis of differentially expressed genes (DEGs) in the 3 d-old post-capped brood of Apis mellifera ligustica subjected to low temperature stress for different timeA: DEGs的趋势分析结果，趋势类型方框内数字为P值The result of trend analysis of DEGs, the number in the box denoting the P-value; B: Profile 6包含DEGs的表达模式Expression profile of DEGs in Profile 6; C: Profile 7包含DEGs的表达模式Expression profile of DEGs in Profile 7; D: Profile 1包含DEGs的表达模式Expression profile of DEGs in Profile 1. CK: 对照组(未经低温胁迫的3日龄封盖子)Control group (3 d-old post-capped brood unexposed to low temperature); T18: 20℃低温处理18 h的3日龄封盖子 3 d-old post-capped brood exposed to low temperature (20℃) for 18 h； T36: 20℃低温处理36 h的3日龄封盖子 3 d-old post-capped brood exposed to low temperature (20℃) for 36 h. 下图同The same below.

2.3 差异表达基因的GO分类

GO富集分析结果显示，3日龄意蜂封盖子在受到低温胁迫后DEGs归类到生物学进程的基因数最多。Profile 6趋势上的DEGs分别富集在30个GO条目上，Profile 7趋势上的DEGs富集在28个GO条目上。这两种趋势从生物学进程分析，主要富集在细胞进程、代谢进程、单有机体进程；从分子功能分析，主要富集在结合、催化活性、分子传感器活性；从细胞组件分析，主要富集在细胞膜组分、细胞膜、细胞。除此之外Profile 6趋势中的DEGs还富集在核酸结合转录因子活性、生物学进程负调控、发育进程、突触和突触组分上，而Profile 7趋势中还富集在行为(1个)、转录因子活性、蛋白质结合(1个)、细胞连接等条目上(图2)。

Profile 1为随低温胁迫时间延长基因表达量先下调表达，之后保持不变，该趋势上富集183个DEGs，这些基因富集在25个GO条目上，排在前10的GO条目分别为代谢进程(28个)、单有机体进程(26个)、结合(22个)、细胞进程(21个)、催化活性(18个)、细胞膜组分(12个)、细胞膜(12个)、细胞(12个)、细胞组分(12个)、应激(8个)条目上(图2)。说明低温胁迫主要影响3日龄封盖子的代谢进程和细胞结构等。同时机体核酸结合转录因子活性、发育相关的基因不会随着低温胁迫时间延长表达量改变，与行为相关的基因表达量随着低温胁迫时间延长而升高。

图2 低温胁迫不同时间3日龄意大利蜜蜂封盖子差异表达基因GO富集分析Fig. 2 GO enrichment analysis of differentially expressed genes in the 3 d-old post-capped brood of Apis mellifera ligustica subjected to low temperature stress for different time

2.4 差异表达基因的KEGG pathway富集分析

通过KEGG pathway富集分析结果显示，Profile 6趋势上的DEGs富集在70条信号通路上，其中有48条与代谢相关的信号通路。基因数富集最多的为代谢途径(22个)、神经活性受体(11个)、次生代谢物的生物合成(11个)、不同环境中的微生物代谢(8个)、抗生素的生物合成(6个)、细胞外基质-受体相互作用(5个)、氨基酸的生物合成(5个)、拼接体(5个)、胞吞作用(5个)等；找到基因ILP(insulin-like peptide receptor gene, ID411297)和基因FoxO(forkhead box protein O gene, ID 727091)同时富集在FoxO信号通路和寿命调节信号通路上(表2)。

表2 Profile 6中DEGs基因数富集前12的KEGG通路Table 2 Top 12 KEGG pathways enriched by DEGs of Profile 6

Profile 7趋势中的271个DEGs富集在51条信号通路上，其中与代谢相关的信号通路有31条。基因数富集最多的为代谢进程(16个)、次生代谢物的生物合成(7个)、嘌呤代谢(4个)、赖氨酸降解(3个)、Jak-STAT信号通路(2个)、Notch信号通路(2个)、淀粉和蔗糖代谢(2个)、TGF-β信号通路(2个)、基础转录因子(2个)、Hippo信号通路(2个)；其中基因CBP(CREB-binding protein gene, ID726332)同时富集在Jak-STAT, Notch, TGF-β, FoxO和Wnt 5条信号通路上(表3)。

表3 Profile 7中DEGs基因数富集前12的KEGG通路Table 3 Top 12 KEGG pathways enriched by DEGs in Profile 7

Profile 1模式中183个DEGs富集在43条信号通路上。DEGs主要富集在代谢进程模块中，包括新陈代谢总览中的代谢途径(13个)、不同环境中微生物代谢(3个)、次生代谢物的生物合成(3个)；核苷酸代谢中有嘌呤代谢(3个)和嘧啶代谢(3个)；氨基酸代谢中有酪氨酸代谢(2个)、丙氨酸、天冬氨酸和谷氨酸代谢(2个)；有1个基因CYP450306a1 (phm, cytochrome P450 306a1 gene, ID08398)富集在昆虫激素的生物合成信号通路上，还发现基因WNT1(protein Wnt-1 gene, ID为413502)同时富集在Hippo, Wnt和hedgehog(Hh)3条信号通路上(表4)。

表4 Profile 1中DEGs基因数富集前13的KEGG通路Table 4 Top 13 KEGG pathways enriched by DEGs in Profile 1

2.5 差异表达基因的RT-qPCR验证

利用RT-qPCR技术检测随机挑选的5个DEGs (kibra,FoxO,CBP,Myo20和arm)的表达模式。结果显示这些DEGs的表达水平与RNA-Seq测序数据的基因表达水平变化趋势一致(图3)，说明本次转录组的测序数据真实可靠。

图3 低温胁迫不同时间3日龄意大利蜜蜂封盖子差异表达基因RNA-Seq测序数据的RT-qPCR验证Fig. 3 Verification of RNA-Seq sequencing data of differentially expressed genes in the 3 d-old post-capped brood of Apis mellifera ligustica subjected to low temperature stress for different time by RT-qPCRA: T18 vs CK; B: T36 vs CK. 图中数据为平均值±标准差，柱上双星号代表同一基因处理组与对照组之间相对表达量差异极显著(P≤0.01, 非配对T检验)。Data in theFigure are mean±SD, and the double asterisk above bars represents extremely significant difference between the treatment group and the control group of the same gene in the relative expression level (P≤0.01, unpaired T-test).

3 讨论

真社会性昆虫蜜蜂蜂子具有典型的狭温性发育特点，蜂群维持稳定的合适巢温对蜂子的健康发育具有十分重要意义。3日龄封盖子是蜜蜂进行化蛹起始阶段(预蛹期)，其低温敏感性强。本研究采用趋势分析对两个梯度时间的低温胁迫预蛹转录本表达模式进行聚类分析，得到了梯度低温胁迫条件下的DEGs及归类，筛选出具有生物学意义的基因集和蜜蜂响应低温胁迫的关键基因。在此基础上，我们重点讨论了FoxO信号通路中的FoxO上调表达在低温胁迫后蜜蜂预蛹无法蜕皮化蛹中的可能机制。

低温胁迫可能导致蜕皮激素(20-hydroxyecdysone, 20E)的合成减少。20E是甾醇类化合物，是在昆虫前胸腺(prothoracic glands, PG)中以胆固醇或植物固醇类物质为底物，经一系列酶，特别是碳22-羟化酶、碳2-羟化酶、碳25-羟化酶的催化作用合成的(Gilbertetetal., 2003; Gilbertet, 2004)。在对下调模式Profile 1中的DEGs进行KEGG富集分析显示，在激素合成信号通路上存在一个下调基因CYP450306a1 (phm)，它通过调节碳25-羟化酶的活性来参与20E合成，phm表达量下调将会导致碳25-羟化酶活性降低，进而影响合成20E(Niwaetal., 2004; Warrenetal., 2004)。因此推测低温会影响20E合成，导致其滴度降低。

低温胁迫后，3日龄封盖子无法化蛹，可能是由于ILP调控FoxO磷酸化留在细胞质内，无法行使其功能而导致无法蜕皮。昆虫化蛹时ILP调控FoxO完成蜕皮的过程如下：在昆虫中存在脊椎动物的胰岛素类似物ILP，它的重要靶标蛋白为FoxO，在ILP上调时，FoxO发生磷酸化，无法进入细胞核而留在细胞质内，促进细胞生长，延迟化蛹(Panetal., 2018)；当ILP下调时，FoxO去磷酸化入核与20E和EcR-Usp聚合体激活后，诱导Br-C,E75A和E75B等初级及DHR3,DHR4和ftz-f1等次级应答基因表达，来激活羧肽酶-A(CPA)表达完成蜕皮发生(Gonzyetal., 2002; Zhou and Riddiford, 2002; Wu and Brown, 2006; Minakuchietal., 2008; Jinderetal., 2013)。在研究中上调模式Profile 6的差异表达基因的富集分析结果显示，在蜕皮相关的FoxO信号通路(3个)、与寿命相关的寿命调节信号通路(2个)同时富集到了ILP和FoxO基因。ILP的表达量上调会导致FoxO磷酸化后无法入核留在细胞质内，无法被20E和EcR-Usp聚合体激活，导致20E的初级应答基因Br-C无法接受20E信号，蜕皮信号传递受阻。另外，20E具有拮抗胰岛素的类似物ILP的作用，同时与ILP共同调控FoxO信号通路，从前一段论述低温胁迫可能导致20E滴度降低，也正好与该结果相符。因此推测在3日龄封盖子在受到低温胁迫后，FoxO发生磷酸化留在细胞质内不能入核，是阻碍蜜蜂化蛹的一个关键基因。化蛹是一个极其复杂的过程，蜕皮只是其中的一个环节，除了蜕皮受到影响之外，低温胁迫带来的影响应该还有很多，仍然需要继续关注。

我们前期研究结果显示3日龄意蜂封盖子在受到低温胁迫后寿命显著缩短(Wangetal., 2016)，FoxO在果蝇的寿命调节起到重要作用(Tataretal., 2004)。转录组数据显示FoxO基因在受到低温胁迫时表达量发生改变，这也可能是受到低温胁迫后，蜜蜂寿命缩短的原因之一。

致谢福建农林大学科研助手刘也齐参与样本获取工作，硕士研究生刘一名、冯睿蓉和李寒参与部分文献研究工作，本科生冯海成、蒋心怡、赫昱畅、常斐然参与部分实验，在此表示感谢！