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

李 寒, 徐新建,2, 周姝婧,2, 周冰峰,2, 朱晨煜, 许宏智, 姚 丹,刘一名, 王 青, 李 想, 王茗琦, 朱翔杰,2,*

(1. 福建农林大学动物科学学院(蜂学学院), 福州 350002; 2. 福建农林大学蜜蜂研究所, 福州 350002;3. 广西中医药大学, 广西高发传染病中西医结合转化医学重点实验室, 南宁 530000; 4. 泉州师范学院, 福建泉州 362000)

温度作为最重要的生态因子之一，对蜜蜂的生长发育、蜂群活动、蜜蜂个体行为、抗病抗逆力均有影响(Irwinetal., 2000; Stabentheineretal., 2021)。蜜蜂个体发育对温度敏感，与多数昆虫相比，更易受温度影响。蜜蜂属于全变态昆虫，其卵、虫、蛹统称为蜂子。蜂子的发育温区非常狭窄，蜜蜂卵的发育温区为31～38℃，封盖子的发育温区为29～38℃，最适发育温度为35℃ (朱翔杰等, 2006; Medrzyckietal., 2010; 周冰峰等, 2011)。蜜蜂卵期低温会导致发育历期延长、死亡率升高(陈琳等, 2016)。蜜蜂封盖子是老熟幼虫和蛹的总称，发育历期12 d，该期低温处理会出现封盖子死亡率升高、羽化蜜蜂寿命缩短和封盖幼虫倒置等现象(Wangetal., 2016)，翅脉形态异常(Zhuetal., 2018)，采集与学习记忆能力受到影响(Tautzetal., 2003; Jonesetal., 2005)，且低温使与记忆相关基因谷氨酸受体A型基因(Glu-RA)及N-甲基-D-天冬氨酸受体I型基因(NmdarI)显著下调表达(Yuanetal., 2016)。

前期有人关注了意大利蜜蜂Apismelliferaligustica和中华蜜蜂Apisceranacerana越冬不同时期工蜂可能具有的与抗寒基因表达量有关的低温抵御机制(Xuetal., 2017; Qinetal., 2019)，处于化蛹发育时期的蜂子进行低温胁迫的差异表达基因(differentially expressed genes, DEGs)主要与发育相关激素调控有关(姚丹等, 2020)，并未发现与越冬工蜂相同的抗寒相关的DEGs，说明发育阶段的蜜蜂与蜜蜂成虫对低温的响应差异较大，蜜蜂发育阶段忍受低温的机制还未知。研究发现处于蛹中期的8 d封盖子(复眼颜色呈深棕色, 关节处轻微着色)在受到96 h长时间低温胁迫后，恢复正常发育温度，并未出现封盖子死亡现象，该时期成为研究蜜蜂发育阶段如何忍受低温一个很好的时间窗口。同时研究发现这些经受住低温的蜜蜂羽化后，寿命却缩短50%(Wangetal., 2016)，说明低温胁迫对其发育仍有不利影响。因此，本研究对西方蜜蜂工蜂8 d封盖子低温处理96 h的样本进行转录组测序(RNA-seq)，探索蜜蜂蛹中期在基因表达层面上可能存在的低温应对机制，以及低温对其发育仍存在哪些不利影响。研究温度对蜜蜂发育的影响，有助于理解低温胁迫对狭温性昆虫的影响以及狭温性昆虫适应低温环境的机制。养蜂生产中，开箱检查、子脾调整、蜂群合并以及蜂王浆生产需移虫操作，均可能出现低温影响蜂子发育的状况。探索低温对蜂子发育的影响，有利于指导改进蜜蜂饲养管理技术(Aizen and Harder, 2009)。

1 材料与方法

1.1 供试蜜蜂

测序样本取自西方蜜蜂蜂群，取样时间为春季2021年4-6月蜂群增长阶段，通过隔王栅限王产卵，获得卵龄一致的蜜蜂卵。将卵脾放入哺育群哺育，以获得发育良好且一致的封盖子。移入哺育群5 d后大幼虫即将封盖，取4 h内封盖的封盖子，割取封盖巢脾，在恒温恒湿箱35±0.2℃(RH 75%)[CTHI-250B, 施都凯仪器设备(上海)有限公司, 精度±0.1℃]中培育8 d作为对照组，记作CK。再将8 d封盖子(蛹中期)放入低温20℃(±0.2℃, RH 75%)(Wangetal., 2016)进行96 h处理为处理组，记作T。低温处理后的封盖子样本立即用液氮冷冻，放入-80℃保存备用。每个蜂群采集样本10头，3群作为3个生物学重复进行转录组测序分析。

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

委托广州基迪奥生物科技有限公司完成构建cDNA文库以及高通量测序，测序平台为Illumina HiSeqTM。对下机数据的原始读段(raw reads)进行过滤，去除含接头的读段(reads)、含N比例大于10%的读段和低质量的读段，最终得到高质量读段(clean reads)(Chenetal., 2018)，然后利用Tophat2(2.1.1)软件将读段比对到西方蜜蜂的参考基因组上Amel_HAv3.1(NCBI Assembly: GCF_003254395.2)，并利用Cufflinks组装转录本。转录组原始数据已上传NCBI的SRA数据库(获取号: SRR15710549-SRR15710554)。

1.3 DEGs的筛选及其GO及KEGG富集分析

使用DESeq2软件对生物学重复的处理组与对照组间进行DEGs分析(Loveetal., 2014)，DEGs筛选条件为|log2FC|≥1 (FC, fold change, 指基因表达水平的变化, FC=T/CK, 其中T为处理组的基因表达水平, CK为对照组的基因表达水平)，错误发现率(false discovery rate, FDR)≤0.05。对CKvsT的DEGs进行GO功能和KEGG通路注释及富集分析(Omicshare Tools, http:∥www.omicshare. com/tools/home/index/index.html)。FDR≤0.01为差异极显著基因。

1.4 DEGs的RT-qPCR验证

为验证RNA-Seq数据的可靠性，从DEGs中随机挑选了5个DEGs,即二氢嘧啶酶(dihydropyrimidinase, DPYS)基因、丝氨酸/苏氨酸蛋白激酶(UNC-51-like kinase 3, ULK3)基因、过氧化物酶体酰基辅酶A氧化酶3(peroxisomal acyl-coenzyme A oxidase 3, ACOX3)基因、叉头框蛋白(forkhead box protein O, FoxO)基因和脂肪酰基辅酶A还原酶CG5065(putative fatty acyl-CoA reductase CG5065, FACRCG5065)基因，进行RT-qPCR验证并利用NCBI Primer Blast设计引物(表1)。从3个蜂群再次取样，取样方法同1.1节。每个样本设置3个蜂群的生物学重复和3次技术重复。利用全式金生物科技有限公司的Transzol Up Plus RNA Kit试剂盒提取样本总RNA，利用反转录试剂盒合成cDNA第1链(北京全式金生物技术有限公司)。以上述的cDNA作为模板进行qPCR。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℃预变性30 s; 95℃变性5 s, 60℃退火和延伸30 s, 共40 个循环；熔解曲线默认系统程序95℃ 5 s; 60℃ 30 s; 95℃ 1 s(QuantStudio5, 美国)。以Actin作为内参基因。

表1 引物信息

1.5 数据分析

RT-qPCR实验基因的表达水平按照2-ΔΔCt法进行计算(Livak and Schmittgen, 2001)，每个目的基因进行3次生物学重复和3次技术重复，取平均值，并计算数据标准差；应用SPSS V21.0软件对RT-qPCR实验处理组与对照组之间平均相对表达量的差异显著性进行非配对T检验。RT-qPCR数据结果在与RNA-Seq测序数据对比时，先经过log2(平均值)的转换。

2 结果

2.1 测序结果质量评估

通过对CK和T测序，分别得到原始读段平均为51 929 304.67和51 198 348.67，过滤后得到高质量有效读段平均分别为50 122 846和50 803 814。统计两端的Q20和Q30平均分别在98.20%和94.85%以上，表明本研究测序数据质量良好，可用于进一步分析。

2.2 低温胁迫后的DEGs

8 d封盖子20℃低温胁迫96 h后，CKvsT间DEGs总数为1 101个，上调基因总数为330个，下调基因总数771个，下调基因数大于上调基因数。

2.3 DEGs的GO功能和KEGG通路注释及富集分析

对获得DEGs的GO功能注释和富集分析结果表明，DEGs主要富集在生物学过程(biological process)(518个)，分子功能(molecular function)(307个)和细胞组分(cellular components)(276个)3大类的27个二级GO条目上，其中基因富集数最多的是代谢过程(metabolic process)(142个)和细胞进程(cellular process)(142个)，其次是结合(binding)(131个)、催化活性(catalytic activity)(120个)和单有机体过程(single organism process)(118个)。上调基因富集最多的为结合(45个)，其次是细胞进程(43个)和代谢过程(41个)；下调基因富集最多的为代谢过程(101个)和细胞进程(99个)，其次是单有机体过程(92个)和结合(86个)。在不同的GO条目上下调基因数均多于上调基因数(图1)。

图1 低温胁迫西方蜜蜂工蜂8 d封盖子DEGs的GO分类

KEGG通路注释和富集分析结果显示，DEGs富集的通路分为细胞进程(30个)、环境信息加工(environmental information processing)(30个)、遗传信息加工(genetic information processing)(49个)、代谢(metabolism)(456个)、有机体系统(organismal system)(9个)5大类(图2)；DEGs富集在113个通路上，富集数排在前10的全部富集于代谢模块，其中富集数最多通路的是代谢通路(signaling pathway)(99个)，其次是次生代谢物的生物合成(biosynthesis of secondary metabolites)(36个)、抗生素的生物合成(antibiotic biosynthesis)(26个)、 微生物在不同环境中的代谢(metabolism of microorganism in different environments)(24个)、碳代谢(carbon metabolism)(15个)、嘌呤代谢(purine metabolism)(13个)、氨基糖和核苷酸糖代谢(amino sugar and nucleotide sugar metabolism)(12个)、脂肪代谢(lipid metabolism)(11个)、氨基酸的生物合成(biosynthesis of amino acids)(11个)。与抗氧化系统相关的过氧化物酶体通路差异极显著基因数有10个，均下调表达，其中控制短链脱氢酶家族成员4(DHRS4)合成的相关基因有两个(LOC100577607和LOC412304)；控制脂酰辅酶A还原酶合成相关的基因有3个(LOC100576414,LOC100578329和LOC725187)。而与谷胱甘肽相关的D-谷氨酰胺和D-谷氨酸代谢通路、谷胱甘肽代谢通路上的谷氨酸脱氢酶基因(LOC409253)和谷胱甘肽S-转移酶基因(LOC552283)下调表达。与发育密切相关的昆虫激素代谢通路中，控制合成保幼激素环解酶基因(JHEH)下调表达，mTOR信号通路控制合成丝氨酸/苏氨酸蛋白激酶基因(ULK3)下调表达，且FoxO信号通路叉头框蛋白基因(FoxO)及胰岛素受体基因(INSR)上调表达(表2)。

图2 低温胁迫西方蜜蜂工蜂8 d封盖子DEGs的KEGG 通路分析

2.4 DEGs的RT-qPCR验证

为验证测序数据的可靠性，选取5个DEGs的表达量进行RT-qPCR定量检测，结果与对应的转录组测序结果中这些差异基因的趋势一致。RNA测序结果的基因表达水平变化趋势一致，证实RNA测序结果可信(图3)。

图3 低温胁迫西方蜜蜂工蜂8 d封盖子DEGs的RNA-Seq测序数据的RT-qPCR验证

3 讨论

本研究发现，低温胁迫可导致调控发育系统发生异常，可能是蜜蜂蛹发育停滞或减缓的原因，封盖子以此降低代谢，降低能量消耗，保存实力应对低温带来的不良影响；同时低温胁迫可能影响谷胱甘肽代谢、过氧化物酶体等影响抗氧化水平，蜜蜂蛹累积氧化损伤，从而羽化蜜蜂受到影响；蜜蜂蛹中期低温胁迫，蜜蜂蛹发育停滞或减缓以应对低温带来的不良影响，可能与保幼激素(JH)降解受到影响，蜕皮激素(20E)合成减少，FoxO磷酸化，细胞自噬被抑制等系列变化有关。蜜蜂的生长发育由20E和JH及胰岛素等多种激素共同调控的(Pursleysetal., 2000; Wangetal., 2013)，激素的滴度变化共同诱导着蜜蜂的发育过程(刘川冬等, 2017)。研究结果显示差异表达基因保幼激素环解酶基因(JHEH)下调表达(表2)，另一方面，FoxO基因上调表达(表2)，FoxO与JHEH启动子序列结合，可抑制JHEH表达(Zengetal., 2017)，这两个基因的差异表达均可导致与JH降解相关的JHEH的合成减少。JHEH的作用是降解JH(Zhangetal., 2005; Tsubotaetal., 2010)，调控JH滴度降低。因此推测低温胁迫导致JH的降解受到影响，JH滴度与正常发育个体相比相对较高，有助于蜜蜂发育维持现状。我们还发现低温处理后参与20E合成的基因P450302a1(Dib)(Chávezetal., 2000)下调表达(表2)，同时，JH滴度相对较高也会抑制20E合成(Gaoetal., 2022)，这说明低温胁迫封盖子会使20E合成减少。低温处理后胰岛素受体基因(INSR)上调表达(表2)，INSR具有磷酸化FoxO蛋白的功能，FoxO蛋白如果发生磷酸化就无法入核，而留在细胞质被降解。昆虫在发育过程中，FoxO蛋白具有启动蜕皮激素下游初级应答基因Br-C传递蜕皮信号和发挥皮层融离作用(Heeetal., 2012; Süren-Castilloetal., 2012; Caietal., 2016)，是蜕皮激素发挥功能的重要环节(Gonzyetal., 2002)。说明低温抑制了蜕皮激素功能，从而发育缓慢或停止发育。这与本团队对低温处理蜜蜂化蛹期封盖子转录组进行分析发现胰岛素类似物基因(ILP)的表达量上调会导致FoxO磷酸化后无法激活20E下游基因(姚丹等, 2020)相似，说明在不同发育阶段低温胁迫均会造成FoxO磷酸化，但其上游的胰岛素信号通路中差异基因不同。此外，还检测到mTOR信号通路丝氨酸/苏氨酸蛋白激酶基因(ULK3)下调表达。mTOR信号通路是生长调节的重要环节，mTOR通路主要通过生长因子，胰岛素与胰岛素样生长因子(IGF)和氨基酸和葡萄糖营养来调节细胞生长，参与细胞自噬(Sarkaretal., 2007; Yang and Guan, 2007)。而此酶在细胞自噬过程中发挥重要作用，介导自噬的起始(Russelletal., 2013)，ULK3基因表达受到低温影响后，细胞自噬也将受到不同程度的影响。因此，低温胁迫影响蛹中期的细胞自噬而使蜜蜂发育减缓。

表2 低温胁迫西方蜜蜂工蜂8 d封盖子的DEGs与发育及抗氧化相关的KEGG通路

蜜蜂发育阶段通过激素调控停止或减缓发育，而非抗寒小分子物质积累或抗冻蛋白来度过低温环境。很多昆虫处于低温环境中时，积累小分子物质以抵御低温环境(Baustetal., 1984; 欧阳芳和戈峰, 2014)。昆虫利用小分子物质积累及抗冻蛋白(AFPs)使冰点及过冷却点降低以抵御低温，利用热激蛋白(HSPs)恢复肽链功能。沙葱萤叶甲Galerucadaurica虫体内抗寒系统为小分子碳水化合物-氨基酸类-脂肪(李浩等, 2014)；二化螟Chilosupperssalis滞育幼虫以蛋白质及小分子糖类积累来抵御低温(张拥军, 2007)。AFPs通过与冰晶不可逆结合降低昆虫过冷却点，云杉蚜虫等昆虫利用AFPs越冬(周艳霞等, 2006)。HSPs可促进低温影响的细胞内肽链重新折叠恢复功能，叶蝉在-7.5～-5℃下hsp20及hsp70表达量上升以抵御低温(Huangetal., 2007)。发育过程中的昆虫多通过滞育或休眠的方式抵御低温，这一过程中主要涉及滞育激素(DH)和JH来调控发育。家蚕Bombyxmori卵期通过DH的作用将体内海藻糖转化为葡萄糖，通过小分子糖类物质积累抵御低温；熊蜂蜂王主要通过JH滴度调控卵巢活性使发育停滞(Fueusawaetal., 1999; Hartfelder, 2000)。本研究未发现HSP, AFP和DH基因差异表达，且正常工蜂不参与生殖，这说明处于发育阶段的蜜蜂工蜂具有狭温性，与其他昆虫具有不同的抵御低温机制，这种机制可能是狭温性昆虫所特有的。

低温引起蛹中期工蜂抗氧化酶表达量下降，这可能会导致氧化还原平衡被破坏，内部氧化水平升高，从而积累氧化损伤。在本研究中蜜蜂蛹中期经过长时间低温胁迫后，进行差异表达基因的KEGG分析，结果显示谷胱甘肽代谢通路谷胱甘肽脱氢酶基因(GDH)、谷胱甘肽硫转移酶基因(GST)以及过氧化物酶体通路短链脱氢酶家族成员4基因(DHRS4)均下调表达(表2)，这些基因均与抗氧化相关。谷胱甘肽代谢相关酶通过使过氧化物失活抑制细胞内活性氧含量，从而维持细胞内氧化还原平衡。GDH位于动植物以及微生物体内的线粒体中，对于细胞内氧化还原平衡起着至关重要的作用(Goldin and Frieden, 1971; Xuetal., 2016)。而GSTs参与调控昆虫的抗逆功能，通过对中华蜜蜂GST基因进行基因沉默后，发现蜜蜂存活率下降，体内氧化水平升高，蜜蜂受到了氧化损伤(Enayatietal., 2005; Mengetal., 2014)。由此推测蜜蜂蛹中期受到低温胁迫会使谷胱甘肽代谢受到影响，引起氧化损伤。此外，辅酶依赖性视黄醛脱氧还原酶(NRDR)作为一种羰基还原酶参与细胞内多种物质的氧化还原反应(张静, 2014)，NRDR由DHRS4基因翻译而成(Perssonetal., 1995)，存在于吞噬细胞质膜，能生成用于清除病原微生物的活性氧(Steck and Grassl, 2014)。低温胁迫后封盖子体内控制NRDR合成的基因下调表达时，抗氧化酶作用受到影响，细胞可能积累氧化损伤。

低温状况下昆虫需要抗氧化酶系统相关基因高表达以抵御低温带来的氧化损伤。昆虫体内的抗氧化酶主要包括超氧化物歧化酶(SOD)、过氧化氢酶(CAT)和过氧化物酶(POD) (李周直等, 1994)。本研究结果显示发育阶段的蜜蜂抗氧化酶系统与成蜂及其他昆虫不同，主要为谷胱甘肽脱氢酶(GDH)、谷胱甘肽还原酶(GR)和谷胱甘肽S-转移酶(GSTs)等。荒漠昆虫小胸鳖甲Microderapunctipennis在低温条件下抗氧化酶基因SOD高表达(肖荣, 2015)；西藏飞蝗在低温条件下体壁与消化道的抗氧化酶SOD, POD和CAT活性增强，以维持虫体氧化水平(吴蕾, 2010)；小金蝠蛾Hepialusxiaojinensis幼虫受低温影响4 h后，体内丙二醛(MDA)含量增加，其GSTs活性增强(张青等, 2016)。从这些研究结果看出，一方面对低温耐受能力强的昆虫抗氧化酶种类与蜜蜂蛹中期不同，另一方面这些昆虫通过抗氧化酶基因上调表达提高抗氧化能力，而有研究表明越冬期中华蜜蜂体内磷酸海藻糖合成酶(TPS)及SOD基因表达水平随低温胁迫时间延长先升高后降低(夏振宇等, 2019)，结合本研究中抗氧化酶表达量下降，96 h的低温处理强度超出蜜蜂抗氧化能力，推测其他昆虫延长低温处理时间也会和蜜蜂一样出现抗氧化酶表达量下降的现象，但这个时间可能需要很长。