中华蜜蜂工蜂中肠响应东方蜜蜂微孢子虫胁迫的高表达基因分析

付中民， 周丁丁， 陈华枝 ， 耿四海 ， 陈大福 ， 郑燕珍 ， 熊翠玲 ， 徐国钧 ， 张 曦 ， 郭 睿

(1. 福建农林大学动物科学学院蜂学学院，福州 350002； 2. 枣庄职业学院，枣庄 277100)

1 引 言

蜜蜂是自然界中最重要的授粉昆虫，能显著提高作物、果蔬的产量和品质[1]，具有无可比拟的经济和生态价值.蜜蜂也是一种高度进化的社会性昆虫，在社会行为学和种群遗传学等研究领域应用广泛[2].中华蜜蜂(Apisceranacerana，简称中蜂)是东方蜜蜂(Apiscerana)的指名亚种，也是适应我国本土自然环境的特有蜂种，具有善于利用零星蜜源、采集力强及抗病性强等优势[3].

微孢子虫是专性胞内寄生的真菌，宿主范围包括哺乳动物、鱼类和昆虫等脊椎动物和无脊椎动物[4].微孢子虫对家蚕和蜜蜂等经济昆虫的影响尤为严重[5].东方蜜蜂微孢子虫(Nosemaceranae)特异性侵染成年蜜蜂的中肠上皮细胞，可导致蜜蜂免疫力下降，工蜂离巢飞行提前，寿命缩减，以及产蜜量下降，严重危害蜜蜂健康[6-7].随着东方蜜蜂微孢子虫基因组(大小：7.86 Mb)的公布和东方蜜蜂基因组(大小：228.332 Mb)完善，二者互作的组学研究不断深入[8-10].较多的研究结果表明，N.ceranae是导致蜂群崩溃综合症的关键因素之一[11-12].前人对N.ceranae与西方蜜蜂(Apismellifera)的相互作用进行了较多的组学研究，取得了一些重要进展.Aufauvre等[13]利用二代测序技术对N.ceranae、杀虫剂胁迫及二者共同胁迫1 d和7 d的西方蜜蜂的中肠样品进行转录组测序，分析发现7 d时宿主的中肠免疫调节系统受到了强烈抑制，表皮防御功能以及海藻糖代谢通路受到了显著影响.Badaoui等[14]利用RNA-seq技术研究了N.ceranae感染5、10和15 d的西方蜜蜂的基因表达谱，发现差异表达基因可归入3种表达模式，进一步分析发现N.ceranae感染对宿主的氨基酸代谢产生干扰及抗菌肽基因的表达产生抑制，表明N.ceranae感染抑制了西方蜜蜂的免疫系统.但有关N.ceranae-东方蜜蜂互作的组学研究较为滞后，相关信息匮乏.

笔者所在课题组前期对N.ceranae胁迫意大利蜜蜂(Apismelliferaligustica，简称意蜂)工蜂过程中宿主和病原的高表达基因(highly expressed gene, HEG)进行全面分析，揭示了二者的HEG的表达谱及潜在作用[15].本研究利用Illumina测序技术对正常及N.ceranae胁迫的中蜂7 d和10 d工蜂中肠进行测序，通过生物信息学方法对宿主的HEG进行深入分析，以期在转录组水平探究中蜂对N.ceranae的胁迫应答.研究结果可揭示HEG在中蜂胁迫应答中的作用，为探明背后的分子机制提供有益的信息和线索.

2 材料与方法

2.1 实验材料

中蜂工蜂取自福建农林大学动物科学学院蜂学学院教学蜂场；N.ceranae孢子由福建农林大学动物科学学院蜂学学院蜜蜂保护学实验室保存并纯化[16].

2.2 方 法

2.2.1 孢子纯化、接种感染及人工饲养N.ceranae孢子的粗提和纯化方法如下：(1)从N.ceranae感染严重的意蜂蜂群巢门口抓取外勤蜂若干只，拉取中肠放入研钵加水充分研磨，差速离心后弃上清液，纯水悬浮沉淀，继续多次差速离心以去除杂质，将得到的孢子粗提液，按照1∶10比例加入无菌水后分装到1.5 mL灭菌的EP管内，4 ℃，3000 r/min离心5 min，吸去上清液后将5管浑浊液合为1管，4 ℃，1 200 r/min离心3 min；(2)分别配制25%、50%、75%、100%的percoll纯化液.各浓度的percoll液取200 μL按照密度从高到低的浓度依次分装到干净的EP管中，最后加入离心后的浑浊液，4 ℃，14 000 r/min离心30 min；(3)将无菌的注射器小心吸取孢子，针头插入乳白色的孢子带中，注意不要吸到孢子带上下方的液体.

参照本实验室前期已建立的方法[17]进行中蜂工蜂的接种感染和人工饲养：(1)配制含孢子的50%(w/v)糖水，孢子浓度为2×108/mL；(2)选取群势较强的成熟封盖子脾迅速提至实验室，放入34±0.5 ℃，(60%～70%)RH培养箱，将刚出房的工蜂放入灭菌的塑料盒(塑料盒四周打孔以通风，30只/盒)，将装有50%(w/v)无菌糖水的饲喂器插入每个盒子上方；(3)将工蜂出房当天记为0 d，34±0.5 ℃培养24 h；将上述中蜂工蜂单头固定后饥饿处理2 h，处理组工蜂单只饲喂含孢子的糖水5 μL，对照组工蜂单只饲喂不含孢子的糖水5 μL；(4)24 h后均饲喂不含孢子的糖水，每日更换糖水，及时清理未食尽糖水的工蜂及死去的工蜂，食尽糖水的工蜂用于后续实验，用50%(w/v)的糖水饲喂工蜂，适宜条件饲喂至10 d.处理组和对照组均包含3个生物学重复，每组重复包含9只工蜂中肠.适宜条件饲喂至7 d时，处理组中的9只工蜂中肠，分别放入3个EP管，该混合样品记为Ac7T(Ac7T-1、Ac7T-2、Ac7T-3)；对照组的9只工蜂中肠，分别放入3个EP管，该混合样品记为Ac7CK(Ac7CK-1、Ac7CK-2、Ac7CK-3)；每次取样在20 s内迅速放于液氮后，迅速转移到-80℃超低温冰箱保存备用.10 d的处理和对照组中肠样品也按照上述方法进行取样和保存，分别记为Ac10T(Ac10T-1、Ac10T-2、Ac10T-3)和Ac10CK(Ac10CK-1、Ac10CK-2、Ac10CK-3).

2.2.2 RNA提取、cDNA建库及RNA-seq 首先利用RNAiso Reagent试剂盒分别抽提上述12个中肠组织样品的总RNA，然后用RNase-Free DNase I去除残留的基因组DNA.接着利用Oligotex mRNA Kits Midi对抽提的总RNA进行RNA纯化，通过分光光度计对纯化RNA进行定量检测.RNA的质量通过琼脂糖凝胶电泳和NanoDrop ND-1000( Nano Drop公司，美国)进行检测.参照陈大福等[18]的方法进行cDNA文库构建.委托广州基迪奥生物技术有限公司，采用Illumina Hiseq 4000测序平台进行双端测序.

2.2.3 测序数据质控、HEG筛选及生物信息学分析 对于测序得到的原始读段(raw reads)，去除含接头(adapter)、含未知碱基(N)的比例大于10%、低质量的读段，用Tophat软件将过滤得到的有效读段(clean reads)映射核糖体数据库，进而将未比对上的clean reads映射东方蜜蜂基因组(assembly ACSNU-2.0)，将比对上的clean reads用于后续的数据分析.

利用FPKM(Fragment Per Kilobase of per Million mapped reads)算法计算基因表达量，公式为FPKM =Total exon fragment/[Mapped fragment (millions) × exon length (KB)]，按照FPKM >15的标准从所有表达基因中筛选出HEG[19].利用在线工具平台OmicsShare(http://www.omicshare.com/tools/index.php/Home/Index/index.html)对筛选出的共有及特有HEG进行GO分类及KEGG pathway富集分析等相关生物信息学分析.

表1 RT-PCR的引物信息

2.2.4 HEG的RT-PCR验证 为验证本研究RNA-seq的数据的可靠性，从4组中肠样品的共有HEG中随机选取8个HEG进行RT-PCR验证.根据上述8个基因的序列，利用Primer Premier 5设计特异性引物，委托福州擎科生物技术有限公司进行引物合成，相关引物序列信息见表1.利用试剂盒AxyPrepTMMultisource Total RNA Miniprep(AXYGEN公司，中国)抽提上述4组中肠样品的总RNA.以上述总RNA作为模板，利用试剂盒HiScript 1st Strand cDNA Synthsis Kit(Vazyme公司，中国)进行反转录，反应步骤按照说明书进行.PCR反应在T100热循环仪(Bio-Rad公司，美国)上进行，反应体系包括：Mixture 10 μL，ddH2O 7 μL，上下游引物各1 μL，cDNA模板1 μL.反应条件：94 ℃预变性5 min；94 ℃变性50 s，60 ℃退火30 s，72 ℃延伸1 min，共34个循环；最后72 ℃延伸10 min，4 ℃保存.PCR产物经1.5%的琼脂糖凝胶电泳检测，核酸凝胶成像仪(上海培清，中国)观察拍照.

3 结果与分析

3.1 高通量测序数据概述

本文中蜂工蜂中肠样品共测得1 809 736 786条raw reads，经严格质控得到1 562 162 742条clean reads，将未比对上核糖体数据库的clean reads映射到东方蜜蜂基因组，Ac7CK、Ac7T、Ac10CK及Ac10T各组的平均比对率为75.78%、55.01%、78.13%和44.19%；Q30均达93.34%及以上(表2).此外，Ac7CK、Ac7T、Ac10CK及Ac10T的Pearson相关系数都在0.872 8及以上(图1).上述结果说明测序样品的重复性较好，测序数据的质量可满足进一步分析需要.

表2 RNA-seq数据信息统计

图1 中蜂工蜂中肠各组样品内不同生物学重复间的Pearson相关性

Fig.1 Pearson correlations among different biological repeats within eachA.c.ceranaworker’s midgut sample group

3.2 HEG的Venn分析

根据FPKM >15的标准，分别从Ac7CK、Ac7T、Ac10CK及Ac10T中筛选出3 162、3 311、3 305和2 463个HEG.Venn分析结果显示四组样品的共有HEG为2 074个，特有HEG分别为89、283、156和78个(图2).

图2 中蜂工蜂中肠各组样品HEG的Venn分析

Fig.2 Venn investigation of HEGs in everyA.c.ceranaworker midgut sample group

3.3 共有及特有HEG的GO分类

GO富集分析结果显示，共有HEG主要涉及生物学进程、细胞组分和分子功能三类，分布于45个GO条目.生物进程中，富集基因数最多的是代谢进程(409)、细胞进程(380)和单组织进程(314)(图3A)；细胞组分中，富集基因数最多的是细胞(254)、细胞组件(254)和细胞膜(153)(图3B)；分子功能中，富集基因数最多的是结合(385)、催化活性(329)和转运器活性(60)(图3C).

Ac7CK特有HEG涉及37个GO条目，基因富集数最多的是细胞进程(204)、代谢进程(196)、结合(191)、催化活性(181)和单组织进程(165)；Ac7T特有HEG涉及35个GO条目，基因富集数最多的是细胞进程(219)、代谢进程(210)、结合(193)、催化活性(189)和单组织进程(183)；Ac10CK特有HEG涉及37个GO条目，基因富集数最多的是细胞进程(230)、代谢进程(223)、结合(217)、催化活性(203)和单组织进程(183)；Ac10T特有HEG涉及28个GO条目，基因富集数最多的是细胞进程(71)、结合(63)、单组织进程(59)、代谢进程(58)和催化活性(51).

图3 共有HEG的GO分类Fig.3 GO categorization of the shared HEGs

3.4 共有及特有HEG的KEGG代谢通路富集分析

KEGG富集分析结果显示，共有HEG富集在39条通路，其中基因富集数最多的前3位分别是核糖体(90)、氧化磷酸化(63)和蛋白酶体(25)(图4).

Ac7CK的特有HEG富集在39条通路，富集基因数最多的前4位分别是内质网蛋白加工(33)、氧化磷酸化(27)、核糖体(25)和内吞作用(25)；Ac7T的特有HEG富集在39条通路，富集基因数最多的前4位分别是核糖体(32)、氧化磷酸化(29)、内质网蛋白加工(26)和剪接体(23)；Ac10CK的特有HEG富集在39条通路，富集基因数最多的前4位分别是内质网蛋白加工(38)、氧化磷酸化(32)、核糖体(27)、和碳代谢(26)；Ac10T的特有HEG富集在37条通路，富集基因数最多的前4位分别是内质网蛋白加工(8)、PI3K-Akt信号通路(7)、Wnt信号通路(6)和TGF-β信号通路(5).

图4 共有HEG的KEGG代谢通路富集分析

圆圈大小代表富集在某一通路的基因数多少，越大代表基因数越多；圆圈颜色代表某一通路的富集基因的显著性高低，越红代表显著性越高

Fig.4 KEGG pathway enrichment analysis of the shared HEGs

The size of the circle indicates the number of enriched genes in a certain pathway, the bigger the more; the color of the circle indicates the significance of enriched genes in a certain pathway, and the redder the color, the higher the more significant

图5 8个共有HEG的RT-PCR验证

A:类易化海藻糖转运体tret1编码基因; B: 类胰蛋白酶-1编码基因; C: 类磷脂酶A1编码基因; D: 精氨酸激酶亚基X1编码基因; E: 类α胰蛋白酶-3编码基因; F: 40S核糖体蛋白S15编码基因; G: 类锌羧肽酶编码基因; H: 类生命必需蛋白lethal(2)编码基因

Fig.5 RT-PCR verification of eight shared HEGs

A:facilitated trehalose transporter Tret1-like encoded gene; B: trypsin-1-like encoded gene; C: phospholipase A1-like encoded gene; D: arginine kinase isoform X1 encoded gene; E: trypsin alpha-3-like encoded gene; F: 40S ribosomal protein S15 encoded gene; G: zinc carboxypeptidase-like encoded gene; H: protein lethal(2) essential for life-like encoded gene

从4个样品共有的HEG中随机挑选8个进行RT-PCR验证，电泳结果显示8对引物均能扩增出符合预期的目的条带(图5)，初步证实了本研究预测的HEG真实表达.

4 讨 论

中蜂是我国养蜂生产的主要蜂种之一，经长期进化已高度适应我国的生态环境.中蜂是N.ceranae的原始寄主，二者经过长期的协同进化已相互适应，N.ceranae经跨种侵染已扩散至世界各地的西方蜜蜂蜂群.目前，关于N.ceranae与蜜蜂互作方面的研究主要集中在西方蜜蜂[20-22].前人也对中蜂与N.ceranae的互作进行了一些研究，但组学研究较为滞后，相关信息颇为匮乏.Chaimanee等[21]曾对N.ceranae侵染3、 6和12 d的西方蜜蜂(Apismellifera)工蜂进行研究，通过检测4种抗菌肽基因(defensin、abaecin、apidaecin和hymenoptaecin)、eater基因和卵黄原蛋白编码基因(vitellogenin)的表达水平，发现N.ceranae侵染对宿主的体液免疫造成了较强抑制.Sinpoo等[22]对西方蜜蜂和东方蜜蜂进行了蜜蜂微孢子虫(Nosemaapis)和N.ceranae交叉感染试验，发现二种蜜蜂的抗菌肽编码基因和细胞免疫相关基因均显著上调表达，但东方蜜蜂的上皮细胞表现出更强的体液免疫水平且含有更少的孢子数，作者推测东方蜜蜂对N.apis和N.ceranae可能具有更强的抵抗力.本研究利用RNA-seq技术对正常及N.ceranae感染的中蜂工蜂中肠进行测序，从正常的7和10 d工蜂中肠分别筛选出3 162及3 305个HEG，从N.ceranae胁迫的7和10 d工蜂中肠分别筛选出3 311及2 463个HEG.Ac7CK、Ac7T、Ac10CK和Ac10T的共有HEG数为2 074个，推测这些共有HEG对中蜂工蜂中肠的生长和发育中具有重要作用.此外，Ac7T和Ac10T的特有HEG数分别为283及78个，说明随着胁迫时间的延长，宿主应答的特有HEG数呈下降趋势，推测这些HEG在宿主响应N.ceranae胁迫的不同阶段发挥特殊作用.本研究中，Ac7CK、Ac7T、Ac10CK及Ac10T各组的平均比对率为75.78%、55.01%、78.13%和44.19%，处理组的平均比对率明显低于对照组，究其原因，处理组中肠样品包含大量的病原孢子，其测序数据为宿主和病原的混合数据，因病原的clean reads无法比对上东方蜜蜂的参考基因组，势必影响宿主clean reads的比对率；而对照组中肠样品的测序数据中绝大部分都为宿主本身的数据，因而clean reads的比对率较高.进一步分析发现，处理组和对照组测序得到的raw reads，以及质控的clean reads，基本处于同一数量级，因此对后续生物信息学分析的影响很小.

本研究中，Ac7T和Ac10T的特有HEG中分别有219和71个注释到细胞进程，210和58个注释到代谢进程，189和51个注释到催化活性，显示伴随N.ceranae的入侵，中蜂工蜂中肠的新陈代谢和细胞生命活动处于不同幅度的激活状态，但激活程度随着病原增殖而呈下降趋势.病原微生物入侵昆虫时会遭到宿主细胞免疫和体液免疫的联合抵抗，包括内吞作用、吞噬作用、酶促水解作用、黑化作用以及抗菌肽的合成与释放等[15].本研究发现，Ac7T和Ac10T的特有HEG中分别有95和68个富集在应激反应，1和1个富集在免疫系统进程，表明中蜂工蜂对N.ceranae产生了免疫应答且应答水平随胁迫时间的延长有所下降.

本研究中，Ac7T和Ac10T的210和58个特有HEG可注释到77和39条物质代谢通路，包含碳水化合物代谢、氨基酸代谢、辅助因子与维生素代谢、核苷酸代谢、多糖生物合成与代谢等，进一步说明中蜂工蜂中肠的物质代谢因N.ceranae胁迫而激活，但激活程度随胁迫时间延长而降低.前期研究中，笔者所在课题组对胁迫中蜂7和10日龄工蜂中肠的N.ceranae的共有HEG进行了代谢通路分析，发现与物质代谢相关的信号通路多达43条，包含碳代谢(15)、糖酵解/糖异生(11)、氧化磷酸化(8)和磷酸戊糖途径(5)等，表明N.ceranae在侵染过程中的物质和能量代谢较为旺盛(数据未发表).结合本研究中的分析结果，推测N.ceranae在其增殖过程中一方面增强自身物质和能量代谢以满足核酸和蛋白合成需要，另一方面通过某种操纵策略提高宿主的代谢水平，促使其摄入更多的物质和能量满足中蜂和N.ceranae的需要，从而利于病原的增殖.Li等[23]对新羽化的成年意蜂工蜂接种感染N.ceranae，通过检测宿主体内的脂类含量发现宿主的脂质含量随感染时间的延长而显著降低，作者推测脂质可能是用于维持宿主生长和病原增殖的首选营养物质.本研究发现，Ac7T和Ac10T中分别有31和8个特有HEG富集到12和6条脂质代谢通路，包括脂肪酸降解、磷脂酰甘油代谢、鞘脂代谢等，表明宿主通过提高脂质的合成来满足自身和N.ceranae的需要，与前人的研究结果一致.前人研究发现在N.ceranae感染的西方蜜蜂对低浓度的蔗糖溶液可产生伸吻反应，取食更多的蔗糖溶液，氧化磷酸化富集基因的数量增加，表明N.ceranae对西方蜜蜂造成了能量胁迫[24-26].本研究中，Ac7T和Ac10T中分别有29和1个特有HEG富集在氧化磷酸化，表明N.ceranae对宿主造成能量胁迫，与西方蜜蜂的相关研究结果一致.

Toll样受体是一类富含亮氨酸的重复序列的模式识别受体蛋白质家族，广泛分布于免疫细胞及某些体细胞的表面，作为宿主防御反应中的重要组成部分以及先天性免疫和获得性免疫的连接点[27].MAPK属于丝氨酸-苏氨酸蛋白激酶，MAPK信号通路是一种包括MAPK、MAPK激酶激酶、MAPK激酶的保守三级激酶模式，共同调节着细胞的生长、分化、增殖及对环境的应激反应、免疫应答、炎症反应等多种细胞生理及病理的过程[28].NF-κB信号通路通过调节某些抗死亡基因和分子的表达，介导调节巨噬细胞中的炎症反应及免疫应答[29].本研究中，Ac7T中分别有3、6和1个特有HEG富集在Toll样受体信号通路、MAPK信号通路和NF-κB信号通路；Ac10T中有1个特有HEG富集在NF-κB信号通路，但未发现有HEG富集在其余两条体液免疫通路.上述结果表明Toll样受体等3条信号通路在中蜂工蜂响应N.ceranae胁迫的前期表现活跃，可能发挥重要的免疫防御作用，随着胁迫时间的延长，宿主的NF-κB信号通路在免疫防御中的角色渐趋关键.细胞因子介导的Jak-STAT信号通路参与细胞的增殖、分化、凋亡以及免疫调节等诸多生物学过程.本研究中，Ac10T中有2个特有HEG富集在Jak-STAT信号通路，表明该信号通路被N.ceranae激活.笔者团队在前期研究中发现对于N.ceranae胁迫10 d的意蜂工蜂中肠，没有HEG富集Jak-STAT和NF-κB通路[15]，表明中蜂和意蜂响应N.ceranae胁迫的免疫应答存在差异，值得进一步深入探讨.N.ceranae感染可抑制西方蜜蜂中肠内稳态和细胞凋亡相关基因的表达水平[15].本研究发现，Ac7T中有2个特有HEG富集在细胞凋亡，Ac10T中无HEG富集在此通路，推测中蜂工蜂中肠细胞可通过提高相关基因的表达促进细胞凋亡，以限制N.ceranae的增殖.

综上所述，本研究对正常及N.ceranae胁迫的中蜂工蜂中肠HEG进行了细致深入的生物信息学分析，提供了宿主的HEG表达谱，揭示了HEG在宿主响应N.ceranae胁迫中的潜在作用.研究结果为在分子水平深入解析中蜂响应N.ceranae胁迫的应答机制、宿主-病原互作机制提供了有益的信息和线索.