西方蜜蜂工蜂蛹响应低温胁迫的差异表达基因分析

刘一名, 徐新建, 周姝婧, 姚 丹, 李 寒, 朱晨煜, 李 想,赫昱畅, 周冰峰, 朱翔杰

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

温度是影响昆虫生长发育的重要环境因子之一，蜜蜂是社会性昆虫，其蜂子(卵、幼虫、蛹的统称)发育具有狭温性且对发育温度敏感。蜜蜂具有调节巢温的能力，蜂巢有蜂子的情况下，蜂群的中心温度会稳定在35℃(Kronenberg and Heller, 1982; Stabentheineretal., 2010; 周莉等, 2015)。研究证明意大利蜜蜂Apismelliferaligustica工蜂封盖子的最适发育温区为33～35℃，当发育温度偏离封盖子发育的适温区，会导致蜜蜂发育畸形或死亡(朱翔杰, 2006; 李月, 2008)。封盖期温度的微小变化会影响成蜂的外部形态(朱翔杰, 2006)、大脑发育(Grohetal., 2004)、记忆能力(Tautzetal., 2003; Jonesetal., 2005)、行为(Becheretal., 2009)和抵抗外界环境压力的能力(Floresetal., 1996; Medrzyckietal., 2010)。

低温胁迫会提高蜜蜂体内的抗冻剂和抗氧化酶系统从而提高个体抗寒性。越冬期蜜蜂工蜂体内甘油、海藻糖、甘露醇、山梨醇等抗冻剂物质大量积累(徐凯等, 2018; Qinetal., 2019)。越冬期中华蜜蜂Apisceranacerana的总抗氧化能力以及超氧化物歧化酶和过氧化物酶活性增强(夏振宇等, 2019)。低温胁迫也会影响昆虫发育期昆虫脂肪、碳水化合物的代谢通路，进而影响昆虫的生长发育以及免疫通路(李洁琼等, 2013; 库尔班·吐松等, 2016; Xuetal., 2017; 姚丹等, 2020)。转录组分析揭示中华蜜蜂低温胁迫后糖和氨基酸的生物合成、氨基酸合成与代谢通路、钙离子活性通路相关基因差异表达，热休克蛋白、锌指蛋白和丝氨酸/苏氨酸蛋白激酶等基因可能是中华蜜蜂成蜂越冬期间抗寒能力的关键基因(Xuetal., 2017)。目前，蜂子抗寒能力和响应低温胁迫的分子机制还不清楚。

本团队前期研究了西方蜜蜂Apismellifera不同日龄封盖子在20℃条件下对低温敏感性，封盖后3 d预蛹是对低温最敏感的发育阶段；封盖后9 d处于蛹中期，对低温最不敏感；封盖后6 d蛹对低温敏感性介于封盖后3和9 d蛹之间(Wangetal., 2016)。封盖后3 d预蛹期封盖子20℃低温胁迫后发现810个基因显著上调和183个基因显著下调，其中胰岛素样肽基因ILP和CREB，结合蛋白基因CBP均上调，可能通过调控FoxO信号通路调控低温抑制3日龄封盖子的化蛹(姚丹等, 2020)。在前期研究基础上，为了探索西方蜜蜂工蜂不同蛹期封盖子对低温胁迫的响应机制，本研究在20℃低温进行短时4 h胁迫不同阶段西方蜜蜂工蜂蛹，构建转录组测序文库，根据差异表达基因(differentially expressed genes, DEGs)的功能分类和代谢通路分析，在转录组水平上探索西方蜜蜂工蜂不同蛹期封盖子对低温胁迫的响应机制。

1 材料与方法

1.1 样本材料

蜜蜂样本取自4-6月春季增长阶段的西方蜜蜂蜂群，准备蜂王产卵力强的产卵群3群和群势12足框的哺育群1群。3个产卵蜂群作为3个生物学重复群。

将空脾置于产卵蜂群进行限王产卵，以获取卵龄一致的蜜蜂卵，再将卵脾放入哺育群哺育。移入哺育群5 d后幼虫即将封盖，每2 h检查一次封盖情况，获取2 h以内封盖的封盖子，割取封盖巢脾。先在恒温恒湿箱中[施都凯仪器设备(上海)有限公司，CTHI-250B]分别培育3, 6和9 d，培育条件35±0.2℃， RH 75%(Wangetal., 2016)。对照组在该条件下继续培养4 h，记为CK3, CK6和CK9；低温处理组则移入温度20±0.2℃和RH 75%的恒温恒湿箱中处理4 h，记为T3, T6和T9。所取封盖子样本立即用液氮冷冻，放入-80℃保存备用。每个日龄的对照组和处理组分别从3个蜂群作为3个生物学重复测序，每群取10只封盖子，10只封盖子混合后作为一个测序样本，共测序18个样本。

1.2 高通量测序及测序数据质控

样本文库构建和高通量测序委托广州基迪奥生物科技有限公司使用Illumina HiSeqTM4000平台进行。下机数据使用fastp进行过滤，去除含adapter的reads、含N比例大于10%的reads、全部都是A碱基的reads、质量值Q≤20的碱基数占整条read 50%以上的reads后得到clean reads(Chenetal., 2018)，对比到西方蜜蜂的基因组Amel_HAv3.1(NCBI Assembly: GCF_003254395.2)。所得原始数据已上传NCBI的SRA数据库(获取号：SRR15258477-SRR15258487, SRR15258489-SRR15258493, SRR15258497-SRR15258498)。

1.3 差异表达基因(DEGs)筛选及功能分析

使用DESeq2软件对生物学重复的处理组与对照组间进行差异分析(Loveetal., 2014)，DEGs筛选条件为|log2FC|≥1[fold change(FC)=处理组FPKM/对照组FPKM]，FDR≤0.05。分别对T3vsCK3, T6vsCK6和T9vsCK9这3组数据进行DEGs分析，得到3个发育阶段低温处理后的DEGs。分别对每个发育阶段的DEGs做Venn分析，进一步进行GO功能分类和KEGG通路分析。

1.4 DEGs的RT-qPCR验证

为验证RNA-seq数据的可靠性，从各比较组T3vsCK3, T6vsCK6和T9vsCK9分别选取8, 6和 5个DEGs(表1)进行RT-qPCR验证，利用NCBI Primer Blast设计引物(表1)。将1.2节测序所用的RNA反转录得到cDNA。以cDNA作为模板进行RT-qPCR，反应体系和程序按全式金生物科技有限公司的PerfectStartTMGreen qPCR SuperMix说明书进行。PCR反应体系(10 μL): H2O 3.4 μL, 2×PerfectStartTMSuperMix 5 μL, 上下游引物(10 μmol/L)各0.2 μL, cDNA模板1 μL, Dye Ⅱ 0.2 μL。反应条件: 95℃预变性30 s; 95℃变性5 s, 60℃退火和延伸30 s, 共40个循环；熔解曲线默认QuantStudio 5系统程序(95℃ 5 s, 60℃ 30 s, 95℃ 1 s)。以Actin作为内参基因，利用2-ΔΔCt法计算目的基因相对表达量(Livak and Schmittgen, 2001)，每个目的基因表达量检测进行3次生物学重复和3次技术重复，取平均值。各基因的相对表达量在与测序数据对比时，先经过log2(平均值)的转换。

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

2 结果

2.1 测序数据质控

通过测序，CK3, CK6, CK9, T3, T6和 T9分别得到原始读段平均为90 772 537, 89 198 197, 83 691 883, 95 418 175, 84 371 909和93 973 924，过滤后得到高质量有效读段平均分别为90 656 963, 89 031 343, 83 586 031, 95 297 370, 84 218 866和93 813 699。所有样品两端的平均Q20和Q30分别为97.28%和92.43%。表明本研究测序数据质量良好，可用于进一步分析。

2.2 DEGs的筛选

与对照组(35℃下)相比，受到低温胁迫后，封盖后3 d预蛹,以及封盖后6和9 d蛹分别有220, 50和26个DEGs。对低温更敏感的封盖后3 d预蛹在受低温胁迫后DEGs数较多；低温影响相对较小的封盖后6和9 d蛹期DEGs数量也相对较少。通过差异表达基因筛选及Venn分析，发现不同日龄的封盖子存在3个共同响应低温的基因，分别为上调表达的3-磷酸肌醇依赖性蛋白激酶1基因(PDK1)(GenBank登录号: LOC410732)和下调表达的Toll样受体基因(Tollo)(GenBank登录号: LOC410231)、H+/Cl-交换转运蛋白7基因(CLCN7)(GenBank登录号: LOC413069)；封盖后3 d预蛹和封盖后6 d蛹有5个共有基因，为肌联样蛋白基因(Titin)(GenBank登录号: LOC107965630)，mab蛋白21基因(mab21)(GenBank登录号: LOC550677)，胞质铁硫蛋白组装蛋白Ciao1基因(Ciao1)(GenBank登录号: LOC551415)，RNA聚合酶Ⅱ转录亚基的介体28基因(MED28)(GenBank登录号: LOC726666)和囊泡抑制性氨基酸转运体基因(VIAAT)(GenBank登录号: LOC409089)；封盖后6 d蛹和封盖后9 d蛹有1个共有基因，为双功能精氨酸脱甲基酶和赖氨酰羟化酶基因(PSR, GenBank登录号： LOC411069)；封盖后3 d预蛹以及封盖后6和9 d蛹分别有211, 41和21个特有响应低温胁迫的基因(图1)。

图1 西方蜜蜂工蜂蛹响应低温胁迫的差异表达基因(DEGs)韦恩图Fig.1 Venn map of differentially expressed genes (DEGs)in Apis mellifera worker pupae in responseto low temperature stressCK3, CK6和CK9分别为35℃培养的封盖后3 d工蜂预蛹以及封盖后6和9 d工蜂蛹对照组；T3, T6和T9分别为20℃培养的封盖后3 d工蜂预蛹以及封盖后6和9 d工蜂蛹低温处理组。CK3, CK6 and CK9 were the 3-d-old post-capped prepupae of workers and 6- and 9-d-old post-capped pupae of workers, respectively reared at 35℃ as the control groups, while T3, T6 and T9 were the 3-d-old post-capped prepupae of workers and 6- and 9-d-old post-capped pupae of workers, respectively exposed to 20℃ as the cold treatment groups.下同The same below.

2.3 DEGs的GO功能分类

封盖后3 d预蛹响应低温胁迫的DEGs的GO功能分类涉及生物学进程(41个)、分子功能(34个)和细胞组分(18个)，富集的二级GO条目多于封盖后6和9 d蛹的，富集基因数较多的为生物学进程大类的代谢进程(metabolic process)、细胞进程(cellular process)、单有机体进程(single-organism process)、生物学进程调控(regulation of biological process)、生物学调控(biological regulation)；分子功能大类的催化活性(catalytic activity)、结合(binding)；细胞组分大类的细胞(cell)、细胞部分(cell part)、细胞器(organelle)；封盖后6 d蛹的DEGs富集数分布为生物学进程(8个)、分子功能(8个)和细胞组分(6个)；封盖后9 d蛹的基因富集数分布为生物学进程(9个)、分子功能(9个)和细胞组分(2个)(图2)。

图2 西方蜜蜂工蜂蛹响应低温胁迫的DEGs的GO功能分类Fig.2 GO functional classification of DEGs in Apis mellifera worker pupae in response to low temperature stress

封盖后3 d预蛹的DEGs分别显著富集生物学进程和分子功能大类的13和9条GO条目，多数与机体和细胞对激素和刺激的反应、转录活性以及酶活性调节和催化活性调节相关(图3)；封盖后6 d蛹的DEGs显著富集7条GO条目，生物学进程大类的4条显著富集条目均与酚类代谢相关，分子功能大类的3条显著富集条目均与氧化还原酶活性相关(图4)，下调DEG酪氨酸羟化酶基因TyHyd均富集在这7条条目上。封盖后9 d蛹的DEGs显著富集的条目只有分子功能大类的结合。

图3 西方蜜蜂工蜂封盖后3 d预蛹响应低温胁迫的DEGs显著富集的GO条目Fig.3 Significantly enriched GO terms of DEGs in 3-d-old post-capped prepupae of Apis mellifera workersin response to low temperature stress

图4 西方蜜蜂工蜂封盖后6 d蛹响应低温胁迫的DEGs显著富集的GO条目Fig.4 Significantly enriched GO terms of DEGs of 6-d-old post-capped pupae of Apis mellifera workersin response to low temperature stress

2.4 DEGs的KEGG通路富集

KEGG B级通路富集显示，受到低温胁迫后封盖后3 d预蛹以及封盖后6和9 d蛹的DEGs均有在代谢上的整体概述图(global and overview maps)和氨基酸代谢(amino acid metabolism)；环境信息加工的信号转导(signal transduction)；细胞进程的运输和分解代谢(transport and catabolism)有富集(图5)。3个不同阶段封盖子的共有DEG 3-磷酸肌醇依赖性蛋白激酶1基因PDK1同时富集在自噬-动物、mTOR信号通路和FoxO信号通路(表2)。

封盖后3 d预蛹响应低温胁迫的DEGs所富集的信号通路最多，较多地富集在代谢上的碳水化合物代谢(carbohydrate metabolism)、脂质代谢(lipid metabolism)、糖的生物合成与代谢(glycan biosynthesis and metabolism)等；细胞进程上的细胞生长与死亡(cell growth and death)；且只有封盖后3 d预蛹在有机体系统上的通路有富集(图5)。封盖后3 d预蛹特有的DEG胰岛素受体底物1-B基因IRS1-B和雷帕霉素复合物靶点基因LST8富集在自噬-动物和mTOR信号通路(表2)。

图5 西方蜜蜂工蜂封盖后3 d预蛹以及封盖后6和9 d蛹响应低温胁迫的DEGs的KEGG通路富集Fig.5 KEGG pathway enrichment of DEGs in 3-d-old post-capped prepupae and 6- and 9-d-old post-capped pupaeof Apis mellifera workers in response to low temperature stressA: 代谢Metabolism; B: 遗传信息加工Genetic information processing; C: 环境信息加工Environmental information processing; D: 细胞进程Cellular process; E: 有机体系统Organismal system.

2.5 DEGs的RT-qPCR验证

RT-qPCR检测结果与RNA测序结果的基因表达水平变化趋势一致，证实RNA测序结果可信(图6)。

3 讨论

在社会性进化过程中，蜂子自身适应环境温度的能力较弱，多个研究已经证实蜜蜂封盖子发育期受低温或高温胁迫后影响蜜蜂个体发育健康(朱翔杰, 2006; 李月, 2008; Wangetal., 2016)。但低温胁迫影响蜜蜂个体发育的机制仍未知，因此本研究通过对封盖后3 d预蛹以及封盖后6和9 d蛹不同发育阶段的蜂子响应低温胁迫后的比较转录组学分析，寻找低温胁迫影响蜂子发育的关键基因和关键通路，为进一步探索其分子机制奠定基础。

通过比较蜜蜂蛹不同发育阶段的差异基因，发现了这3个阶段共有的差异基因PDK1，根据该基因所在的3条通路分析，找到了低温胁迫条件下，FoxO信号通路是影响自噬、保幼激素和蜕皮激素调控网络的核心。进一步结合封盖后3 d预蛹所特有的差异基因分析，找到了更多的参与自噬、保幼激素和蜕皮激素调控的差异基因(表2)，这不仅进一步证明了低温胁迫会影响以FoxO信号通路为核心，涉及到自噬和激素调控网络的存在，也很好地解释了封盖后3 d预蛹相对其他发育阶段对低温最敏感的事实。同时，在封盖后6 d的蛹所特有基因中，检测到与昆虫表皮着色和免疫相关的DEGTyHyd，这与前期研究发现该时期蛹低温处理后出现较高比例的黑色个体的现象(朱翔杰, 2006; Wangetal., 2016)一致，为进一步探索该现象的分子机制提供了思路。而封盖后9 d蛹特有的少数基因中，未检测到与预期目标相关的关键基因和关键通路，可能因为该阶段对低温最不敏感，在同样长的低温胁迫强度下，受到的影响较小。

封盖后3 d预蛹以及封盖后6和9 d蛹响应低温胁迫后，DEGPDK1均显著上调表达(图6)，可能调控蜜蜂蛹在低温胁迫下的细胞自噬和蜕皮受到抑制。PDK1是蛋白激酶B(PKB)的上游激酶，PDK1具有两个保守的结构域，即PH结构域和氨基端激酶结构域，能够通过激酶结构域磷酸化一系列激酶保守区域而激活这些激酶，来调节细胞代谢、生长、扩散、生存和抗凋亡等诸多生理过程(Biondi, 2004)。KEGG结果显示，PDK1同时富集在mTOR信号通路、FoxO信号通路和自噬-动物通路(表2)。PDK1是调节细胞自噬的关键因子之一，PDK1会增强PI3K信号通路活性并激活mTOR，使自噬过程受到抑制，细胞正常的代谢水平受到影响(McManusetal., 2016)。PDK1还可通过影响FoxO信号通路抑制蜕皮。PDK1磷酸化Akt使其激活，磷酸化的Akt作用于下游的FoxO， 使 FoxO 磷酸化而不能进入细胞核，滞留在细胞质内的FoxO不能与其靶基因结合，失去转录功能，无法完成蜕皮(Cantrell, 2001; Sarbassovetal., 2005; Caietal., 2016; Linetal., 2018)。PDK1上调表达可引起的Akt磷酸化和FoxO磷酸化， 结果使 FoxO 无法入核启动下游基因的转录。这与我们前期发现的胰岛素样肽基因ILP和CREB结合蛋白基因CBP抑制FoxO入核启动下游基因的转录结果(姚丹等, 2020)一致。这进一步证实了FoxO信号通路介导低温胁迫对蜜蜂蛹发育的影响。

表2 西方蜜蜂工蜂封盖后3 d预蛹响应低温胁迫的DEGs富集数大于2的KEGG通路Table 2 Number of KEGG pathways with the enrichment number of DEGs more than two in 3-d-old post-cappedprepupae of Apis mellifera workers in response to low temperature stress

图6 西方蜜蜂工蜂蛹响应低温胁迫的DEGs RNA-Seq测序结果的RT-qPCR验证Fig.6 Validation of the RNA-Seq sequencing results of DEGs in Apis mellifera worker pupaein response to low temperature stress by RT-qPCR图中RT-qPCR数据为平均值±标准差；柱上星号代表处理组与对照组之间基因表达量差异显著(*P<0.05; **P<0.01; ***P<0.001)(非配对T检验)。RT-qPCR data in the figure are means±SD.Sterisks above bars represent significant difference in the gene expression level between the treatment group and the control group (*P<0.05; **P<0.01; ***P<0.001)(unpaired T-test).

封盖后的预蛹和不同阶段的蛹对低温的耐受性不同，前期的研究结果显示，封盖后3 d预蛹因为处于化蛹阶段，昆虫内部发生剧烈组织重组，所以比封盖后6 d和9 d蛹受低温影响更大(Wangetal., 2016)。从本研究的结果推测，封盖后3 d预蛹与后面的蛹期相比，细胞自噬和蜕皮受到抑制程度更大，因此导致处于化蛹阶段蜜蜂对低温更敏感。封盖后3 d预蛹一方面通过海藻糖转运蛋白基因下调表达，自噬受到抑制程度增强。除了主要通过mTOR途径调节自噬外，海藻糖也是自噬的增强因子(Sarkaretal., 2007)。结果显示在预蛹时期受到低温胁迫后，海藻糖转运蛋白1基因Tret1和海藻糖转运蛋白1-2基因Tret1-2显著下调表达。海藻糖进出细胞发挥其特定的功能，需要转运蛋白的帮助才能完成(於卫东等, 2020)，因此这两个基因下调表达可能会导致海藻糖的功能发挥受到抑制，从而使自噬受到限制。另一方面，与PDK1共同属于胰岛素信号通路的基因IRS1-B在低温胁迫后上调表达，胰岛素信号通路位于PI3K信号通路上游，胰岛素可以促进PDK1的表达(Panetal., 2018)。可能通过IRS1-B上调表达，促进PDK1的表达，蜕皮抑制程度增加。除此之外，蜜蜂预蛹期蜕皮抑制程度增加，还可能是因为低温导致激素调控异常，从而导致蜜蜂致死性蜕皮缺陷。蜜蜂预蛹期要经历从幼虫态到蛹态的变态发育，而昆虫的变态发育主要受到保幼激素和蜕皮激素调控(Thummel, 2001)。封盖后3 d的预蛹在受到低温胁迫后，检测到3个与激素调控相关的DEGs:Kr-h1,Ftz-F1和Eth。在受到低温胁迫后，检测结果显示基因Kr-h1的表达量显著升高，Kr-h1是位于保幼激素途径下游的关键因子，具有抑制20E应答基因表达的作用，阻止20E诱导变态发生(Minakuchietal., 2011; Belles, 2020)，因此由蜕皮激素诱导的变态可能被抑制。同时，还检测到在低温胁迫后，由蜕皮激素诱导的转录因子核激素受体Ftz-F1基因Ftz-F1和蜕皮激素诱导激素基因Eth均下调表达，这2个基因均作用于20E下游的，具有激活蜕皮行为的功能(Parketal., 2002; Dubrovskyetal., 2011)，Eth基因缺失还会导致昆虫在化蛹时死亡(Lenaertsetal., 2017)，说明蜕皮激素诱导变态的功能受到抑制，可能导致蜜蜂无法蜕掉幼虫皮而化蛹失败。

封盖后6 d蛹在受到低温胁迫时，检测到TyHyd的表达量显著降低，可能导致酪氨酸羟化酶活性降低。有研究表明酪氨酸羟化酶是全变态昆虫生长发育所必需的酶，且与昆虫角质层着色与免疫都有关系(Qiaoetal., 2016; Zhangetal., 2019; Liuetal., 2020)。低温胁迫使TH的表达受阻，可能会对封盖子的表皮着色和免疫反应造成影响。