基于RNA-seq技术分析意大利蜜蜂幼虫肠道响应球囊菌胁迫的高表达基因

熊翠玲， 黄枳腱， 梁 勤， 徐细建， 张曌楠， 骆 群， 刘 敏， 陈大福， 郭 睿

(1.福建农林大学蜂学学院，福建福州 350002； 2.江西省养蜂研究所，江西南昌 330200)

蜜蜂白垩病是由蜜蜂球囊菌侵染蜜蜂幼虫而导致的致死性真菌病[1]。该病主要分布在欧洲、北美和中国等地区，早在1913年，Maassen等就报道了白垩病的发生，我国大陆1961年也曾报道发生类似蜜蜂白垩病的病害[2]。蜜蜂白垩病主要是使老熟幼虫或封盖幼虫死亡，雄蜂幼虫最易感染，因为雄蜂幼虫多在巢脾的边缘，易受低温的影响。白垩病会导致大批幼虫死亡，造成群势下降，严重病群失去产蜂蜜和蜂王浆的能力，严重影响蜂群的生产力。球囊菌属于球囊菌门球囊菌属，迄今已鉴定出球囊菌属的22个种，它们以致病性真菌或腐生菌存在于蜜蜂体内；在蜂群里，患病幼虫的尸体以及被污染的饲料与巢脾是疾病传播的主要来源[3]。国内养蜂生产中使用的主要蜂种是意大利蜜蜂和中华蜜蜂，意蜂极易受到球囊菌的侵染而暴发白垩病，而中蜂几乎不受影响。2006年，Evans等公布了西方蜜蜂的基因组，为在分子水平研究西方蜜蜂奠定了基础[4]。对于白垩病的前期研究主要是集中在生化[5]、病理[6]及检测[7]等方面，分子水平的研究相对较少。目前，高通量测序技术应用于白垩病的研究报道极少，蜜蜂幼虫响应球囊菌胁迫的应答研究仍不深入。

以RNA seq为代表的高通量测序技术发展迅猛，现已广泛应用于昆虫包括蜜蜂的相关研究[8-11]，然而，尚无利用二代测序技术研究白垩病过程蜜蜂幼虫应答的相关研究，本研究利用RNA-seq技术对球囊菌胁迫和正常意蜂幼虫肠道进行深度测序，对表达水平最高的100个基因进行GO分类和KEGG代谢通路富集分析，并比较分析球囊菌侵染幼虫肠道与正常幼虫肠道之间高表达基因的差异。

1 材料与方法

1.1 生物材料

本研究中使用的意大利蜜蜂幼虫取自福建农林大学蜂学学院教学蜂场，蜜蜂球囊菌菌株由福建农林大学蜂学学院蜜蜂保护实验室保存并活化。

1.2 主要试验试剂及仪器

RNase-free水购自中国上海生工生物公司，DNaseⅠ和Oligotex mRNA Kits Midi试剂盒购自德国Qiagen公司，Dynal M280磁珠购自Invitrogen公司，高碘酸钠购自美国Sigma公司，DNA ligase购自美国Thermo公司，RNA Reagent抽提试剂盒、ExTaqpolymerase及Superscript Ⅱ reverse transcriptase均购自日本TaKaRa公司。纯化cDNA的Ampure beads为美国Agencourt产品，cDNA文库构建试剂盒TruSeqTMDNA Sample Prep Kit -Set A为美国Illumina公司产品。

恒温恒湿气候箱购自中国宁波江南仪器厂，超纯水仪购自中国四川沃特尔水处理设备有限公司，高速冷冻离心机购自德国Eppendorf公司，倒置显微镜为中国上海光学仪器五厂产品，超净工作台为中国苏州安泰空气技术有限公司产品，PCR仪为美国Bio Rad公司产品，凝胶成像系统为中国上海培清科技有限公司产品，超低温冰箱为中科美菱低温科技股份有限公司产品。

1.3 试验方法

1.3.1 意蜂幼虫饲养 中蜂幼虫的饲料配方参照王倩等的方法[12]，其中将D-果糖和D-葡萄糖换为新鲜蜂蜜。预试验中空白对照组中蜂幼虫7日龄成活率达到90%以上，说明本实验室建立的中华蜜蜂幼虫人工饲养模型可为本研究提供符合要求的中蜂幼虫。从福建农林大学蜂学学院蜂场选择20群健康蜂群(无白垩病症状且蜜蜂球囊菌特异性引物常规PCR检测为阴性)。用移虫针挑取2日龄幼虫，放入无菌的24孔细胞培养板(每孔对应1头幼虫，孔内加有35 ℃预温的幼虫饲料)，将24孔板放入恒温恒湿培养箱，35 ℃、70% RH条件下培养。每隔24 h吸去旧饲料，加入新饲料。3日龄时，处理组饲喂含有球囊菌孢子的人工饲料(终浓度为1×107孢子/mL)，空白对照组饲喂正常人工饲料。适宜条件下饲喂幼虫至7日龄，处理组和对照组均设置3个生物学重复。

1.3.2 测序样品准备 肠道组织是球囊菌孢子萌发的主要场所，也是蜜蜂抵御病原入侵的重要免疫器官，故选择幼虫肠道作为测序材料。分别剖取处理组和对照组4、5、6日龄幼虫肠道，为尽量减少肠道RNA的降解，将从解剖取样到样品放入液氮速冻的时间控制在30 s以内，每剖取1组样品，立即液氮速冻并迅速放入-80 ℃超低温冰箱保存，通过Illumina Hiseq 2500平台对处理组和对照组的12个幼虫肠道样品进行深度测序，其中处理组样品包括4日龄幼虫肠道(AmT1-1、AmT1-2、AmT1-3)、5日龄幼虫肠道(AmT2-1、AmT2-2、AmT2-3)和6日龄幼虫肠道(AmT3-1、AmT3-2、AmT3-3)，对照组为4日龄幼虫肠道(AmCK-1、AmCK-2、AmCK-3)。

1.3.3 RNA抽提和cDNA合成 液氮冷却研磨，按照RNA抽提试剂盒的说明书进行RNA抽提，最后将RNA溶解于200 μL RNase-free水中。取1 μL溶解的RNA，用Evolution 600分光光度计(Thermo Scientific公司，美国)定量确定浓度。

1.3.4 mRNA纯化与cDNA合成 抽提的总RNA首先利用10 U DNase Ⅰ在37 ℃条件下消化1 h，然后利用Oligotex mRNA Kits Midi进行mRNA纯化，纯化后的mRNA用分光光度计进行定量。以10 μg mRNA作为模板，GsuI-oligo dT作为反转录引物，用1 000 U SuperscriptⅡ逆转录酶在42 ℃下孵育1 h合成第1链cDNA；随后利用高碘酸钠氧化mRNA的5′端帽子结构，并连接生物素；通过Dynal M280磁珠筛选连接了生物素的mRNA/cDNA，并通过碱裂解释放第1链cDNA；然后通过DNA连接酶在第1链cDNA的5′末端加上接头，利用ExTaq聚合酶合成第2链cDNA。最后通过GsuⅠ酶切去除polyA和5′端接头。

1.3.5 cDNA测序 将3、4、5、6日龄意蜂幼虫中肠组织的cDNA利用705型Fisher Scientific超声波破碎仪打断至 300～500 bp，利用Ampure beads(Illumina公司，美国)进行纯化。随后纯化的cDNA利用TruSeqTMDNA Sample Prep Kit -Set A(Illumina公司，美国)构建文库，并利用TruSeq PE Cluster Kit(Illumina公司，美国)进行扩增。最后Illumina Hiseq 2500测序仪进行测序反应。测序数据已上传NCBI SRA数据库，SRA Num：SRA456722。

1.3.6 基因表达丰度统计 通过Illumina PE125测序，获得大量原始读段，去除含接头、含N比例大于10%的和低质量的读段，得到清晰读段，进而利用Tophat软件将清晰读段比对至意蜂的参考基因组(Amel_4.5)，表达量的计算使用FPKM法[10]，其计算公式为：FPKM=(106C×103)/(NL)。

1.4 高表达基因的GO分类及KEGG代谢通路富集分析

根据基因的FPKM值，选取各肠道样品表达量水平最高的前100个基因，利用在线分析工具OmicsShare(http：//www.omicshare.com/tools/index.php/Home/Index/index.html)进行GO分类及KEGG代谢通路富集分析。

2 结果与分析

2.1 RNA seq测序数据概述

对12个肠道样品进行深度测序，获得的读段数都在 26 405 020 条以上，去除低质量读段后的清晰读段数达到 26 405 020 条以上，Q20达到97%以上，说明RNA-seq数据质量较好。首先将清晰读段比对核糖体数据库，最多允许5个错配，去除比对上核糖体的读段，剩余读段用于比对意蜂参考基因组(表1、表2)。各样品比对上的比例均在86%以上，说明测序数据真实可靠。并对各肠道样品进行主成分分析(PCA)(图1)，结果显示处理组和对照组样品的聚类情况良好，说明测序样品具有较好的重复性。

2.2 各肠道样品前100个高表达基因的GO分类

一般认为FPKM 值在0.10～3.75之间的是低表达基因，在>3.75～15之间为中度表达基因，大于15的为高表达基因，根据FPKM值对各样品基因进行排序，选取前100位的高表达基因进行下一步分析。

表1 RNA-seq数据过滤信息统计

注：括号内为占比，表2同。

表2 Clean reads与参考基因组的比对信息统计

将意蜂幼虫肠道表达水平最高的100个基因进行GO注释(图2)，这100个高表达基因主要集中在细胞组分、生物学进程和分子功能上发挥作用，AmCK、AmT1和AmT2的高表达基因GO分类数均为29个，AmT3的分类数有27个，在细胞组分方面，基因主要富集在细胞整体、细胞部分、胞外区、大分子复合物、隔膜、隔膜部分、细胞器和细胞器部分；在生物学进程方面，基因主要富集在细胞进程、行为、生物调节、细胞组分起源、发育进程、代谢进程、定位、多组织进程、多细胞组织进程、繁殖、繁殖进程、刺激反应、单组织进程；在分子功能方面，基因主要富集在抗氧化活性、电子载体活性、结合(Binding)、催化活性、结构分子活性。此外，富集在细胞组分中的基因比率高于在分子功能、生物学进程中的基因比率。

AmT1与AmCK比较，同是4日龄，高表达基因富集的GO分类数量相同，且分别富集在细胞组分、生物学进程和分子功能上的基因数量基本一致，说明在胁迫的早期阶段，球囊菌尚未唤起意蜂幼虫的应答，AmT1、AmT2、AmT3之间比较，高表达基因富集的GO分类数量逐渐下降，且分别富集在细胞组分、生物学进程和分子功能上的基因数量也呈逐渐下降的趋势，说明随着时间延长，球囊菌胁迫增强，意蜂幼虫的应答也逐渐增强。

同时，GO富集分析结果显示，随着胁迫时间的延长，富集在代谢过程和催化活性上的基因数量呈下降趋势，表明随着球囊菌的增殖破坏了幼虫肠道的正常代谢功能。

2.3 各肠道样品前100个高表达基因的KEGG代谢通路富集分析

对样品AmCK、AmT1、AmT2、AmT3表达水平最高的前100个基因进行KEGG在线分析，结果显示这些基因分别富集在94、93、93、83个代谢通路上(图3)。统计分析发现这些高表达基因主要集中在新陈代谢、有机体系统、遗传信息加工、环境信息加工等6个途径上。

2.3.1 新陈代谢 根据KEGG直系同源号进行分类，AmCK新陈代谢相关基因涉及脂代谢、碳水化合物代谢、其他氨基酸代谢、氨基酸代谢、能量代谢、萜类和酮类化合物的代谢、外来物质的降解和代谢等8个代谢途径。分析结果显示，细胞色素C氧化酶Ⅰ、线粒体ATP合酶脂质结合蛋白、细胞色素b、细胞色素C氧化酶Ⅱ(线粒体)4个编码蛋白参与能量代谢，酮脂酰CoA硫解酶、谷胱甘肽S-转移酶S1、类视黄醛脱氢酶参与脂类代谢过程、α-葡萄糖苷酶前体、山梨糖醇脱氢酶、果糖二磷酸醛缩酶参与碳代谢过程。精氨酸激酶同工酶、酮脂酰CoA硫解酶、类视黄醛脱氢酶参与氨基酸代谢过程，在细胞色素P450通路上的谷胱甘肽S-转移酶S1基因、谷胱甘肽S-转移酶D1型X2基因共同作用，细胞色素P450是一种含高铁血红素的蛋白，作为一种末端氧化酶，细胞色素P450主要催化机体内源和外源性物质在体内的氧化反应，参与大部分药物的氧化代谢[13-14]。可以得到多个基因可以参与多个通路。

2.3.2 遗传信息加工 遗传信息加工通路中主要涉及转录、翻译、折叠、组装和降解，其中参与翻译通路高表达基因占据一个很大的比例，其中主要有编码40S核糖体蛋白基因和60S核糖体蛋白，总共有37个基因编码这2类蛋白，翻译通路中的RNA转运的延伸因子1-α同工酶基因，该基因参与许多重要的细胞过程和疾病，包括信号传导、翻译控制、凋亡、细胞骨架组成、病毒复制等[15]。热休克蛋白在涉及转录、转移、折叠、组装和降解过程，该基因参与新生蛋白折叠、组装加工成有生物学活性的蛋白质。

2.3.3 环境信息加工 在KEGG分析中，环境信息加工过程具有2个方面，分别是信号转导和信号分子及其相互作用。40S核糖体蛋白S6、 肌动蛋白相关蛋白等参与相关信号转导表达。

2.3.4 细胞进程 细胞进程中高表达基因主要富集在转运和分解代谢、细胞运动、细胞通讯等3个方面。天冬氨酸蛋白酶、GL12416、NPC2蛋白质、肌动蛋白、热休克蛋白等参与转运和分解代谢。通路中GL12416基因产物参与吞噬体的形成以及细胞通讯过程中的间隙连接，相关基因还参与溶酶体的形成和作用。肌动蛋白参与吞噬体通路过程，是细胞的一种重要骨架蛋白。同时肌动蛋白细胞分泌、吞噬、移动、胞质流动和胞质分离等细胞通讯过程中起重要作用。

2.3.5 有机体系统 有机体系统主要包括免疫系统、循环系统、消化系统等9个系统。超氧化物歧化酶、热休克蛋白等蛋白参与有机体系统。热休克蛋白参与MPKM通路，丝裂原活化蛋白激酶(MAPK)是生物体内重要的信号转导系统之一，参与介导生长、发育、分裂、分化、死亡以及细胞间的功能同步等多种细胞过程，MAPK能被多种炎性刺激激活，并对炎症的发生、发展起重要调控作用[16-17]。

2.4 对照组与处理组间比较

意蜂幼虫肠道AmCK、AmT1、AmT2之间的差异性不大，高表达基因富集在新陈代谢、有机体系统、遗传信息加工、环境信息加工和细胞进程上的代谢通路数基本一致。AmT3相对于AmCK具有一些特异性基因，特异性表达出一些抗菌肽、防御素前体蛋白物质等，抗菌肽指昆虫体内经诱导而产生的一类具有抗菌活性的碱性多肽物质，抗菌肽对部分真菌、原虫、病毒及癌细胞等均具有强有力的杀伤作用[18-19]，说明随着球囊菌胁迫时间的增加，幼虫肠道诱导产生抗菌肽。

2.5 处理组间比较

在球囊菌胁迫前期(AmT1、AmT2)，高表达基因富集在新陈代谢、有机体系统、遗传信息加工、环境信息加工和细胞进程上的代谢通路数基本一致，而在AmT3胁迫后期，高表达基因富集在新陈代谢上的代谢通路数大幅下降，并且免疫性物质增加，再次说明球囊菌对意蜂幼虫的生长代谢产生了严重影响，防御机制开始发挥作用。

KEGG代谢通路富集分析结果(表3至表6)显示，AmCK、AmT1、AmT2和AmT3中基因富集量最多的都为核糖体通路，但是富集的基因数量呈逐渐下降的趋势，说明球囊菌破坏了宿主的蛋白合成系统，推测球囊菌能通过互作减少宿主的营养供给，从而促进自身增殖。

表3 AmCK组基因富集量前10位KEGG代谢通路统计

表4 AmT1样品基因富集量最多的前10个KEGG代谢通路

表5 AmT2样品基因富集量最多的前10个KEGG代谢通路

2.6 各肠道样品前100个高表达基因的Venn分析

对处理组和对照各幼虫肠道样品的Venn分析结果显示， 72个高表达基因为AmCK、AmT1、AmT、AmT3所共有，各样品中的特有高表达基因分别为1、2、2、13个，特有的高表达基因大多是编码抗菌肽、防御素前体、铁蛋白、抑制剂等免疫方面的基因，推测这些共有基因在球囊菌侵染的全过程都发挥着重要功能，而少数特有基因则在病程的各个阶段分别发挥作用，对这些共有与特有高表达基因的深入研究将有助于揭示意蜂幼虫的应答机制(图4)。

表6 AmT3样品基因富集量最多的前10个KEGG代谢通路

3 讨论

本研究选取意蜂幼虫肠道作为测序材料，因为肠道是球囊菌侵染的主要场所，其转录组变化能更准确地反映意蜂幼虫对球囊菌入侵的应答。为了从全局水平了解意蜂幼虫受球囊菌胁迫前后的基因表达谱，本研究利用RNA-seq技术对正常幼虫肠道和球囊菌胁迫的幼虫肠道进行深度测序，通过比对意蜂参考基因组，共发现有13 592个基因在幼虫肠道中表达，其中已知基因有11 918个(87.68%)，另有197个新基因，这些新基因将为完善意蜂基因组信息提供一定补充。

Venn分析结果显示，4个肠道样品中的共有高表达基因有72个，比如部分参与核糖体通路的基因，合成加工成供幼虫自身生长发育需要的物质，也有参与核糖体通路的U-box结构域蛋白基因，可作为分子伴侣或辅分子伴侣，还可能作为剪接体的部分因子，识别细胞内错误折叠或异常蛋白质，从而使其正确折叠或降解。另外，U-BOX蛋白质还与癌病变以及自身免疫性疾病等病理过程有关，该类蛋白质通常还在将底物转运给蛋白体的过程中起作用，调节着体内蛋白质的降解[20]。各样品中的特有高表达基因分别有1、2、2、13个，对比处理组样品和对照组样品中的特有高表达基因，将有可能是意蜂幼虫响应球囊菌胁迫的关键应答基因，基因编码的核苷酸结合蛋白β亚基1在AmT3上特异性表达，该蛋白特异性结合抗原，可以起到抗体的作用，在抵御侵染上起到一定作用。对比处理组各样品中的特有高表达基因，或许能为明确球囊菌胁迫过程不同阶段的关键应答基因提供重要线索，这些特有高表达基因值得进一步研究。

由于经费限制，本研究只设4日龄正常幼虫肠道作为对照组，会导致一些基因表达信息的缺失，因此，对4、5、6日龄正常幼虫肠道分别取样并同时进行深度测序，便于更全面地获得意蜂幼虫响应球囊菌胁迫的基因表达谱信息。

本研究中，在进行RNA-seq数据分析时，发现有部分reads不能定位于已有的意蜂参考基因组，它们有可能是尚未克隆的新基因，这些新基因也许在意蜂幼虫肠道发育中发挥重要作用。本研究利用RNA-seq技术对正常和球囊菌胁迫意蜂幼虫肠道进行深度测序，通过对各样品表达量水平最高的前100个基因进行GO分类和KEGG代谢通路富集分析，初步解析了意蜂幼虫响应球囊菌胁迫的应答，推测在6日龄之后免疫机制迅速增强，应答机制的深入阐释有赖于差异表达基因的进一步研究。