意大利蜜蜂工蜂中肠响应东方蜜蜂微孢子虫胁迫的可变剪接基因分析

周倪红, 王海朋, 周丁丁, 付中民, 祝智威, 范元婵, 张 曦, 熊翠玲, 郑燕珍, 陈大福, 郭 睿

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

蜜蜂作为自然界中最重要的授粉昆虫，能够大幅提高农产品的产量，并且在维持生态系统平衡方面具有不可替代的作用[1].近年来，由于病原菌、寄生虫、病毒、农药、杀虫剂及自然环境等因素的影响，蜜蜂健康受到严重威胁，数量正在急剧下降[2-4].东方蜜蜂微孢子虫(Nosemaceranae)是一种广泛存在的蜜蜂真菌病原，可通过蜜蜂的采集和交哺等行为感染蜜蜂，主要寄生于宿主的中肠上皮细胞，可引起蜜蜂发生微孢子虫病[5].感染N.ceranae的蜜蜂伴有寿命缩短、王浆腺发育不全、采集能力明显降低以及方向认知能力减弱等现象[6-7]，可导致染病蜂群的抗病性和抗逆性大幅减弱、蜂群群势严重下降[8].此外，N.ceranae被认为是蜂群崩溃综合征的主要诱因之一[9].

可变剪接是指mRNA前体通过不同的剪接方式形成不同的剪接异构体，经过调控元件和剪接因子的相互作用，形成不同组合，最终生成不同的蛋白质，是真核生物转录水平调控的重要机制，丰富了转录组和蛋白质结构与功能的多样性[10-11].内含子保留、可变3′剪接位点、可变5′剪接位点和外显子缺失是4种常见的可变剪接类型.人类基因组包含约25 000个基因，但可变剪接产生的转录本能翻译出多达90 000种蛋白质[12].高等真核生物中存在丰富的基因可变剪接事件，动物中约20%的多外显子基因中有新的剪接位点，并且根据RNA-seq数据和EST数据估计95%的外显子基因都能发生可变剪接[13].线虫中的已注释基因都存在可变剪接[14].拟南芥中发生可变剪接的基因占比达到61%[15].此外，可变剪接在果蝇的生长发育、神经分布和性别决定等方面发挥着重要作用[16].目前，可变剪接相关研究多集中于哺乳动物及少数模式昆虫如果蝇和家蚕[17-19].然而，对于意大利蜜蜂(Apismelliferaligustica，简称意蜂)响应N.ceranae的胁迫应答过程中发生的可变剪接事件、可变剪接基因(alternatively spliced gene, ASG)应答及作用，仍缺乏相关研究报道.笔者所在课题组前期已在mRNA组学水平探究了中华蜜蜂(Apisceranacerana，简称中蜂)幼虫肠道响应球囊菌胁迫的ASG应答[20]；此外，利用Illumina测序技术对正常及N.ceranae胁迫的意蜂工蜂中肠进行了全转录组测序，基于高质量的测序数据和生物信息学分析方法对宿主中肠的微小RNA(miRNA)、长链非编码RNA(long non-coding RNA, lncRNA)和环状RNA(circular RNA, circRNA)差异表达谱及调控网络以及宿主响应胁迫的高表达基因表达谱进行了全面解析[21-24].

在前期研究基础上，本研究基于已获得的高质量全转录组数据，利用生物信息学分析方法进一步对意蜂工蜂中肠响应N.ceranae胁迫过程中的ASG数量、种类、表达谱及潜在作用进行探究.研究结果可在组学水平揭示ASG在意蜂工蜂响应N.ceranae胁迫过程中的作用，为在分子水平深入研究ASG的功能、宿主的胁迫应答机制提供依据.

1 材料与方法

1.1 意蜂工蜂中肠的全转录组数据来源

前期已利用RNA-seq技术对正常及N.ceranae胁迫的意蜂工蜂中肠进行了全转录组测序，获得了高质量的测序数据(包括mRNA组学数据).样品制备及测序过程简述如下：(1)将成熟封盖子脾提至实验室，收集刚出房的意蜂工蜂(0 d)，放入四面打孔的塑料盒，置于恒温恒湿培养箱培养，培养条件为(34±0.5) ℃，60%～70% RH；(2)待工蜂经过24 h的适应后，利用已制备的N.ceranae纯净孢子(1 mL含106个)对处理组刚出房的意蜂工蜂(0 d)进行饲喂接种(每只5 μL)；对照组刚出房的意蜂工蜂每只饲喂不含孢子的糖水5 μL；该试验进行3次生物学重复；(3)在实验室条件下对上述处理组和对照组工蜂进行笼养，培养条件为(34±0.5) ℃，60%～70% RH，分别于接种后7和10 d拉取对照组和处理组的工蜂中肠，各组样品分别为：Am7CK(Am7CK-1、Am7CK-2、Am7CK-3)，Am7T(Am7T-1、Am7T-2、Am7T-3)，Am10CK(Am10CK-1、Am10CK-2、Am10CK-3)，Am10T(Am10T-1、Am10T-2、Am10T-3).上述12个中肠样品经液氮速冻后委托广州基迪奥生物科技有限公司的进行全转录组测序，测序平台为Illumina HiseqTM4000.相关的原始数据已上传至NCBI SRA数据库，SRA号：SRA456721.

1.2 测序数据质控

对于下机的原始读段(raw reads)，利用Per脚本去除含有adaptor、未知核苷酸(N)比例大于5%和低质量测序数据，得到的高质量有效读段(clean reads)[22]用于本研究中可变剪接基因的相关分析.利用公式FPKM=total exon fragment/[mapped fragment (millions)×exon length (KB)]计算所有表达基因的FPKM值.

1.3 ASG类型分析

参照郭睿等[20]的方法，利用短reads比对工具bowtie分别将各测序样品的clean reads映射到核糖体数据库，去除mapping上核糖体的reads后，将剩余的clean reads进一步比对到西方蜜蜂的参考基因组(Amel_4.5)，比对结果包含了剪接位点信息.

将TopHat2比对结果中的所有剪接位点信息按不低于5个测序数据过滤掉大部分偶然事件造成的噪音，再与参考的已知剪接位点进行比较(允许1 bp的误差)，找出已知的剪接位点，并对剩下的新剪接位点进行可变剪接事件的分类统计.

1.4 ASG的相关生物信息学分析

利用Omics Share在线工具集合(http://www.omicshare.com/tools/Home/Soft/venn)进行ASG的Venn分析、GO分类及KEGG代谢通路富集分析，参数采用默认参数.

2 结果与分析

2.1 意蜂工蜂中肠ASG类型

在Am7CK、Am10CK、 Am7T和Am10T中检测到可变剪接事件数分别为73 383、85 618、79 813和74 693次，对应的基因数分别为9 586、10 063、9 829和9 622个，各样品中平均每个ASG发生的可变剪接事件数分别为7.66、8.50、8.12和7.76次.各样品的可变剪接事件中外显子跨越数量最多，分别达到4 411、5 247、4 993和4 629个，所占可变剪接事件总数比例分别为6.01%、6.13%、6.29%和6.25%(图1).

2.2 意蜂工蜂中肠ASG的Venn分析

Venn分析结果显示，Am7CK、Am10CK、Am7T和Am10T的共有ASG数为9 056个，占ASG总数的85.7%，4个中肠样品的特有ASG数分别为100、153、275和95个(图2A).此外，Am7CK与Am7T的共有ASG数为9 302个，特有的ASG数分别为284和527个；Am10CK与Am10T的共有ASG数为9 412个，特有的ASG数分别为651和210个；Am7CK与Am10CK的共有ASG数为9 380个，特有的ASG数分别为206和683个(图2B-D).

2.3 意蜂工蜂中肠共有和特有ASG的GO分类

GO分类结果显示Am7CK、Am10CK、Am7T和Am10T的共有ASG涉及分子功能、细胞组分和生物进程相关的50个GO条目(图3)，富集基因数最多的为结合(2 439个)、代谢进程(2 361个)以及细胞进程(2 313个)；Am7CK特有的ASG主要富集在结合(47个)、细胞进程(43个)和代谢进程(41个)等26个功能条目，Am7T特有的ASG主要富集在结合(71个)、细胞进程(67个)和单组织进程(60个)等30个功能条目，Am10CK特有的ASG主要富集在结合(93个)、细胞进程(74个)和代谢进程(68个)等33个功能条目，Am10T特有的ASG主要富集在结合(22个)、催化活性(20个)和细胞进程(18个)等个24功能条目.

2.4 意蜂工蜂中肠共有和特有ASG的KEGG代谢通路富集分析

进一步对Am7CK、Am10CK、Am7T和Am10T的共有ASG进行KEGG分析，结果表明，有2 259个ASG富集在325个代谢通路(图4)，富集基因数前10位的为内质网蛋白加工(120个)、核糖体(118个)、RNA转运(116个)、内吞作用(114个)、嘌呤代谢(112个)、癌症路径(109个)、剪接体(108个)、亨廷顿氏舞蹈症(102个)、泛素介导的蛋白水解作用(91个)、氧化磷酸化(90个)；还有部分ASG富集在PI3K-Akt(74个)、Hippo(55个)、Hedgehog(18个)等与生长发育相关的信号通路，Jak-STAT(20个)、Toll-like受体(15个)和NF-kappa B(13个)等与免疫相关的信号通路.

Am7CK的特有ASG富集在甘油磷脂代谢(3个)、甘氨酸、丝氨酸和苏氨酸代谢(2个)、半胱氨酸和甲硫氨酸代谢(2个)等31条代谢通路；Am7T的特有ASG富集在神经活性配体—受体相互作用(5个)、糖酵解/糖异生(4个)、缬氨酸、亮氨酸和异亮氨酸降解(3个)等104条代谢通路；Am10CK的特有ASG富集在RAP1信号通路(5个)、亨廷顿氏舞蹈症(5个)、胃酸分泌(4个)等126条代谢通路；Am10T的特有ASG富集在神经活性配体—受体相互作用(1个)、糖苷脂类生物合成—乳酸和新乳糖系列(1)、谷氨酸、亮氨酸和异亮氨酸降解(1个)等22个代谢通路.上述4个样品的共有ASG在糖酵解/糖异生信号通路的富集详情如图5所示.

3 讨论与结论

意蜂是我国养蜂生产中使用的主要蜂种，具有重要的经济价值，此外，意蜂作为主要的授粉昆虫，对农业生产意义重大.N.ceranae的原始寄主是东方蜜蜂，后跨种感染西方蜜蜂，现已成为危害蜜蜂健康和养蜂生产的真菌性病原.较多的研究结果表明[26-29]，可变剪接在动植物的生长发育、胁迫响应中起到非常重要的作用.近20年来，人们围绕西方蜜蜂-N.ceranae互作进行了较多的生物化学及分子生物学研究，随着高通量测序技术的快速发展和应用，相关的组学研究日益增多.然而，对于蜜蜂在N.ceranae胁迫过程的ASG应答表现，仍缺乏相关研究报道，ASG在胁迫响应中的作用尚不清楚.为深入探究意蜂-N.ceranae互作，本课题组前期结合链特异性建库的RNA-seq技术对正常及N.ceranae胁迫的意蜂7和10 d工蜂中肠进行全转录组测序，基于高质量的组学数据全面解析了宿主和病原的高表达基因表达谱[22]，意蜂工蜂中肠发育过程的差异lncRNA表达谱及调控网络、差异circRNA表达谱及调控网络[23-24].基于已获得的高质量mRNA组学数据，本研究通过生物信息学方法对意蜂工蜂中肠响应N.ceranae胁迫的ASG应答进行深入分析，在Am7CK、Am10CK、 Am7T和Am10T样品中检测到可变剪接事件数分别为73 383、85 618、79 813和74 693次，对应的基因数分别为9 586、10 063、9 829和9 622个，各样品中平均每个ASG上发生的可变剪接事件数分别为7.66、8.50、8.12和7.76次；其中，各样品总事件中外显子跨越事件类型均占比最大，分别为4 411(6.01%)、5 247(6.13%)、4 993(6.29%)和4 629(6.25%)个.本课题组前期对中华蜜蜂(Apisceranacerana)幼虫肠道响应蜜蜂球囊菌(Ascosphaeraapis)胁迫的ASG应答进行了探究，鉴定出发生于9 124个基因的57 327次可变剪接事件，并发现以基因间区(17.68%)、可变3′端剪切(15.32%)外显子跨越(14.12%)和可变5′端剪切(12.81%)类型为主[20].本研究的结果与之存在差异，表明在不同蜂种对不同真菌病原的胁迫应答中，不同数量的ASG和不同的可变剪接类型发挥着不同的作用.Am7CK与Am10CK中的可变剪接事件数分别为73 383和85 618次，表明可变剪接事件随意蜂工蜂中肠发育时间延长而逐渐活跃；Am7T中的可变剪接事件数较Am7CK增加8.76%，但Am10T中的可变剪接事件数较Am10CK下降12.76%，说明意蜂工蜂中肠的响应N.ceranae胁迫的前期和后期，ASG通过上调和下调可变剪接事件的数量参与宿主的胁迫应答，但具体的机制仍需进一步深入研究.此外，还发现9 065个ASG为Am7CK、Am10CK、Am7T和Am10T所共有；Am7T和Am10T的特有ASG数分别为153和95个；推测上述共有ASG在意蜂工蜂中肠的生长发育、新陈代谢和胁迫应答中具有基础性的作用，特有ASG在宿主胁迫应答的不同阶段发挥特殊功能.

本研究中，Am7CK、Am10CK、Am7T和Am10T的共有ASG涉及50个功能条目，以结合(2 439个)、代谢进程(2 361个)、细胞进程(2 313个)为主，说明共有ASG在意蜂工蜂中肠的生长发育和新陈代谢等生命活动中发挥广泛作用.此外，Am7T和Am10T中分别有22和2个特有ASG富集在应激反应，表明ASG参与宿主对N.ceranae的免疫防御，但参与程度随胁迫时间的延长而下降.进一步对上述共有ASG进行代谢通路富集分析，发现它们富集在多达325条信号通路，其中包含112条物质和能量代谢相关通路，例如三羧酸循环(34个)和半乳糖代谢(20个)等14条碳水化合物代谢相关通路，脂肪酸降解(24个)和甘油酯代谢(34个)等16条脂质代谢相关通路，嘌呤代谢(112个)和嘧啶代谢(75个)等2条核苷酸代谢相关通路，以及6条能量代谢相关通路(150个)；此外，还有多达117个ASG富集在免疫系统相关通路.上述结果再次表明意蜂工蜂响应N.ceranae胁迫的应答过程中，共有ASG广泛参与宿主的物质和能量代谢以及免疫防御.

昆虫体内存在两种天然免疫防御系统，即体液免疫和细胞免疫.体液免疫主要包括Toll、IMD和JAK/STAT 3条信号通路，通过信号转导及免疫途径调控免疫相关基因的表达，诱导产生抗菌肽和其他效应分子[30-31].细胞免疫由血细胞介导，主要完成对病原物的包裹、吞噬和集结等[32-33].吞噬作用的发生依赖一些可识别病原物的受体如PRRs等，病原物被昆虫表面受体识别后，受体被激活，诱导胞间信号级联反应，完成对病原物的吞噬[34-35].集结作用与包裹作用类似，通过血细胞粘附、聚集在细菌表面，进而形成黑化或非黑化结节，以此来对病原物产生应答[36-37].Am7T的特有ASG富集在104条代谢通路，其中包括内吞作用(1)、溶酶体(1)和黑化作用(1)等3条细胞免疫通路；Am10T的特有基因仅富集在吞噬作用这1条细胞免疫通路；上述结果表明意蜂工蜂中肠响应N.ceranae胁迫的不同时期，特有ASG通过参与不同的细胞免疫通路对病原产生免疫应答，胁迫前期的应答以内吞作用、溶酶体和黑化作用为主，吞噬作用在胁迫后期的应答中占据主要位置.

糖酵解作为生物体最重要的分解代谢途径之一，是很多哺乳动物组织和细胞(红细胞、脑等)唯一代谢能量来源.研究表明[38]，昆虫消化吸收的部分葡萄糖在脂肪体中转化为海藻糖或作为糖原储存，机体需要供能时，海藻糖再转化为葡萄糖进入糖酵解途径；在柞蚕蛹发育过程中，己糖激酶基因的表达对柞蚕蛹解除滞育有重要调节作用，在蛹发育中期己糖激酶活力及葡萄糖含量之间呈现出互补关系，即己糖激酶的活性增强，血淋巴中葡萄糖的含量减少；被N.ceranae感染后的西方蜜蜂饥饿水平上升，能量消耗增多，蜜蜂死亡率升高[39].本研究表明，Am7T中有4个特有ASG富集在糖酵解/糖异生通路，表明宿主在响应N.ceranae胁迫的前期，通过ASG调节糖酵解/糖异生通路以适应N.ceranae造成的能量胁迫.

综上所述，本研究通过生物信息方法对意蜂工蜂中肠响应N.ceranae胁迫的ASG应答及ASG的潜在作用进行了细致深入的分析，研究结果不仅提供了宿主ASG的数量和种类信息，而且揭示了共有和特有ASG在宿主胁迫响应中的重要作用，为ASG的分子克隆和功能研究提供了必要的数据资源，也为深入认识意蜂对N.ceranae的胁迫应答提供了新见解.下一步将通过更高级的分析方法如基因权重共表达分析、融合mRNA组学数据和非编码RNA组学数据进一步挖掘ASG的相关信息，并筛选出关键ASG，用于克隆其全长序列.