中华蜜蜂幼虫肠道响应球囊菌早期胁迫的转录组学

陈大福，郭睿，熊翠玲，梁勤，郑燕珍，徐细建，张曌楠，黄枳腱，张璐，王鸿权，解彦玲，童新宇



中华蜜蜂幼虫肠道响应球囊菌早期胁迫的转录组学

陈大福，郭睿，熊翠玲，梁勤，郑燕珍，徐细建，张曌楠，黄枳腱，张璐，王鸿权，解彦玲，童新宇

（福建农林大学蜂学学院，福州 350002）

【目的】白垩病是困扰养蜂生产的顽疾。目前，尚无利用二代测序技术研究中华蜜蜂（,简称中蜂）幼虫白垩病的报道。本研究利用RNA-seq技术对健康（AcCK）及球囊菌()胁迫的中蜂4日龄幼虫肠道（AcT）进行深度测序，在转录组水平研究中蜂幼虫在球囊菌胁迫早期的胁迫应答。【方法】通过Illumina HiSeq 2500平台对AcCK和AcT进行双端（PE125）测序，首先对测序数据进行质控和评估，利用edgeR软件进行差异表达基因（DEG）分析，进而对DEGs进行GO富集分析及KEGG代谢通路富集分析，最后，利用实时荧光定量PCR（qRT-PCR）验证测序数据的可靠性。【结果】 AcCK和AcT的转录组测序共得到188 457 338条原始读段（raw reads）, 经过滤得到182 088 448条有效读段（clean reads），两端Q20与Q30均在97.96%和94.97%以上，说明测序数据质量良好；主成分分析（PCA）结果显示第一与第二主成分可分别解释基因表达总体差异的75.8%和10.7%；DEG分析结果显示上调基因与下调基因的数量分别为344和239个；GO富集分析结果显示DEGs共富集在36个GO分类（term）上，其中基因富集数最多的细胞（106 unigenes）、细胞组件（106 unigenes）和代谢进程（104 unigenes）；KEGG代谢通路富集分析结果显示上调与下调基因分别富集在72个和45个代谢通路上，其中上调基因富集数最多的是核糖体（72 unigenes）、碳代谢（16 unigenes）和糖酵解（14 unigenes），而下调基因富集数最多的是碳代谢（9 unigenes）、二羧酸代谢（8 unigenes）和氨基酸生物合成（7 unigenes），进一步分析表明中蜂幼虫肠道的部分细胞免疫对球囊菌的胁迫产生应答，而体液免疫不产生应答，宿主的代谢相关基因受到球囊菌的显著抑制。【结论】揭示了中蜂幼虫在球囊菌入侵早期的胁迫应答，为深入解析中蜂幼虫的胁迫应答机制提供了重要信息，也为在分子水平研究关键应答基因打下了基础。

中华蜜蜂；幼虫肠道；球囊菌；转录组；RNA-seq

0 引言

【研究意义】蜜蜂作为一种社会学模式昆虫，对发育学、神经生物学、行为学和病原-宿主互作等研究有重要价值[1-5]。蜜蜂也是最重要的授粉昆虫，在农业生产、粮食安全和生态维持中也发挥着巨大作用[6]。白垩病是球囊菌（）特异性侵染蜜蜂幼虫的致死性真菌病，长期困扰养蜂生产，每年给养蜂业造成巨大损失。近年来，随着蜂产品全球贸易的快速发展，白垩病呈逐年上升趋势[7]。有关白垩病的研究主要集中于西方蜜蜂（）幼虫，但关于西方蜜蜂幼虫响应球囊菌胁迫的应答研究报道极少，目前尚无球囊菌侵染东方蜜蜂（）幼虫过程中的应答机制研究，也无防治白垩病的有效方法。【前人研究进展】国内养蜂生产中使用的主要蜂种是意大利蜜蜂（，简称意蜂）和中华蜜蜂（，简称中蜂），二者分别属于西方蜜蜂和东方蜜蜂。意蜂极易暴发白垩病，但中蜂极少受到该病的影响。近20年来，国内外学者对白垩病开展了一系列研究，主要集中在病原分类鉴定[8-9]、形态学[10-11]、病理学[12-13]、流行病学[14]、侵染过程[15-16]、蜜蜂防御[17-18]以及疾病防治[19]等方面。笔者课题组也在球囊菌的病理和检测方面开展了较为系统的研究[20-24]。2006年，基因组信息的公布[25]为其分子生物学研究奠定了重要基础。Evans等[26]在蜜蜂细胞内鉴定出了大多数NF-B信号通路中的基因，如的2个同系物，但都不与直系同源，表明是仅存在于果蝇的特化分支而不存在于其他昆虫；Aronstein等[27]利用cDNA-AFLP对健康及球囊菌感染西方蜜蜂幼虫进行了测序比较，发现差异表达基因（DEGs）参与了能量代谢和蛋白转运，其中的类几丁质编码基因很可能参与了宿主对球囊菌的抵抗；Cornman等[28]利用RNA-seq技术对幼虫芽孢杆菌（）感染后12 h和72 h的蜜蜂幼虫进行测序，鉴定出75个显著上调和6个显著下调基因，其中若干编码抗菌肽基因和2个编码围食膜基质基因显著上调；Julie等[29]通过RNA-seq研究了东方蜜蜂微孢子虫（）和杀虫剂单独或共同饲喂意蜂后中肠的转录组变化，发现和杀虫剂并无协同效应，长期暴露于杀虫剂环境抑制了意蜂的免疫基因的表达。2015年，韩国研究人员测序并公布了的基因组[30]，但当时并未公布基因位置和功能注释信息。为深入开展中蜂的转录组学研究，笔者课题组前期已组装并注释了中蜂幼虫肠道的参考转录组[31]。【本研究切入点】养蜂生产中，偶尔可见中蜂幼虫罹患白垩病的现象，前期研究发现中蜂幼虫可在实验室条件下被球囊菌侵染。白垩病研究基本都集中于球囊菌侵染西方蜜蜂幼虫，有关球囊菌侵染中蜂幼虫的相关研究几乎没有。本研究利用RNA-seq技术对健康及球囊菌胁迫的中蜂幼虫肠道进行深度测序，进而在转录组水平研究宿主球囊菌入侵早期的胁迫应答。【拟解决的关键问题】以健康及球囊菌胁迫的中蜂4日龄幼虫肠道作为研究对象，通过Illumina测序技术研究中蜂幼虫肠道在球囊菌入侵早期的胁迫应答，为深入解析中蜂幼虫响应球囊菌胁迫的应答机制及关键应答基因的功能研究提供重要信息。

1 材料与方法

试验于2015年12月至2016年8月在福建农林大学蜂学学院蜜蜂保护学实验室进行。

1.1 材料

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

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

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

1.2 方法

1.2.1 球囊菌活化、孢子纯化及计数 将4℃保存的球囊菌培养皿置于已紫外灭菌的超净台，用酒精灯上灼烧片刻的镊子夹取少量菌丝在平板中央2 cm左右圆心区域内划线操作，盖上培养皿，置于37℃生化箱恒温培养，3 d后观察球囊菌生长情况。接种8—10 d以后，待培养皿上黑色孢子较多时，刮取孢子至干净的EP管中，用研磨棒充分研磨，按照Jensen等[32]的方法离心纯化球囊菌孢子，梯度稀释后用血球计数板进行计数。

1.2.2 中蜂幼虫的人工饲养及肠道样品准备 按照已报道的饲料配方[33]配制中蜂幼虫饲料并进行改良，预试验结果显示中蜂幼虫的7日龄成活率可达70%以上。从学院教学蜂场群势较强的健康中蜂蜂群（无白垩病症状且PCR鉴定为阴性）移取2日龄幼虫至无菌的24孔细胞培养板（每孔对应1只幼虫，孔内加有35℃预温的幼虫饲料），35℃，相对湿度70%条件下饲养。每隔24 h更换饲料。处理组3日龄幼虫饲喂含球囊菌孢子的人工饲料（孢子终浓度1×107个/mL），对照组3日龄幼虫饲喂正常饲料。饲喂孢子后24 h，分别剖取处理组（AcT）和对照组（AcCK）4日龄幼虫肠道组织，每剖取一只幼虫肠道，迅速将肠道移至RNA Free的EP管，再投入液氮速冻，待一组肠道样品（7只幼虫肠道）集齐后，迅速转移保存于-80℃。本试验进行3次生物学重复。AcCK与AcT的3个生物学重复分别为AcCK-1、AcCK-2、AcCK-3和AcT-1、AcT-2、AcT-3。因幼虫肠道是球囊菌的寄生场所，而本文旨在对中蜂幼虫在球囊菌胁迫早期的应答进行研究，故选取4日龄幼虫肠道作为测序对象。

1.2.3 cDNA文库构建及Illumina测序 利用RNAiso Reagent试剂盒抽提处理组和对照组中蜂幼虫肠道组织的总RNA，然后用RNase-free DNaseI去除基因组DNA残留。RNA的质量通过琼脂糖凝胶电泳和NanoDrop ND-1000（NanoDrop，Wilmington，DE，USA）进行检测。利用Oligotex mRNA Kits Midi试剂盒说明书，纯化各样品总RNA中的mRNA。以10 μg mRNA作为模板，GsuI-oligo dT作为反转录引物，用1000 U Superscript II reverse transcriptase在42℃下孵育1 h合成第1链cDNA；随后利用高碘酸钠氧化mRNA的5′端帽子结构，并连接生物素；通过Dynal M280磁珠筛选连接了生物素的mRNA/cDNA，并通过碱裂解释放第1链cDNA；然后通过DNA ligase在第1链cDNA的5′末端加上接头，利用Ex Taq polymerase 合成第2链cDNA。最后，通过I酶切去除polyA和5′端接头。利用Ampure beads对上述cDNA进行纯化，cDNA文库通过TruSeqTMDNA Sample Prep Kit-Set A进行构建和TruSeq PE Cluster Kit进行扩增。上述6个幼虫肠道样品委托广州基迪奥生物科技有限公司进行测序，测序平台为Illumina HiSeq2500。本研究测得的转录组数据已上传美国国家生物技术信息中心（NCBI）SRA数据库，SRA号：SRA456721。

1.2.4 数据分析 对于下机数据，利用Perl脚本去除含有adaptor、未知核苷酸比例大于5%和低质量reads，获得有效读段（clean reads）。利用R软件（version 2.16.2）进行测序饱和度分析。使用短reads比对工具bowtie[34]将clean reads比对（mapping）到核糖体数据库（最多允许5个错配），去除mapping上核糖体的reads，将保留下来的数据用于转录组的组装及分析，进而利用SOAP aligner/soap2软件[35]将未mapping上核糖体的reads mapping到笔者课题组组装的中蜂幼虫肠道参考转录组[31]。

利用FPKM（Fragments Per Kilobase of Transcript Per Million Mapped Reads）法计算基因表达量。利用R软件（version 2.16.2）计算各样品之间的相关性系数。利用edgeR软件[36]进行DEG分析。DEG的筛选标准为FDR≤0.05且|log2fold change≥1。将DEGs向GO数据库（http://www.geneontology.org/）的各条目（term）映射，并计算每个term的基因数，从而得到具有某个GO功能的基因列表及基因数目统计，然后应用超几何检验，找出与整个基因组背景相比，在DEGs中显著富集的GO term。KEGG代谢通路显著性富集分析以KEGG pathway为单位，应用超几何检验，找出与整个基因组背景相比在差异表达基因中显著性富集的代谢通路。

1.2.5 实时荧光定量PCR（qRT-PCR）验证 为了验证RNA-seq数据，随机选取5个DEGs进行qRT-PCR。qRT-PCR反应按照SYBR Green Dye试剂盒（Vazyme公司，中国）操作说明书进行，每个反应进行3次重复。20 μL的反应体系中含SYBR Green Dye 10 μL，10.0 μmol·L-1正、反向引物各1 μL，cDNA模板DNA 1 μL，DEPC水8 μL。qRT-PCR反应在ABI 7500荧光定量PCR仪（ABI公司，美国）上进行，反应条件：95℃预变性1 min，95℃变性15 s，60℃延伸30 s，共40个循环，最后72℃延伸45 s。所选基因的相对表达量采用2-∆∆Ct法[37]计算。

2 结果

2.1 Illumina测序数据质控与评估

中蜂幼虫肠道样品的Illumina PE125测序共测得188 457 338条reads，经过滤得到182 088 448条clean reads，各样品clean reads数均在95.4%以上，两端Q20均在97.96%以上，两端Q30均在94.97%以上（表1）。饱和度分析结果显示随着测序量的增多，检测到的基因数也随之上升、增长速度趋于平缓，说明检测到的基因数趋于饱和（附图1）。上述结果表明本研究中的测序数据质量良好，可用于进一步分析。

Pearson相关性分析结果显示AcCK与AcT的组内各生物学重复之间的相关性均接近1，说明样本的重复性较高（附图2）。进一步对AcCK与AcT的各样品进行主成分分析（PCA），结果显示第一主成分（PC1）与第二主成分（PC2）共能解释样品基因表达总体差异的86.5%，AcCK与AcT的组内各生物学重复聚类良好（图1），表明AcCK与AcT的总体基因表达模式差别显著。

RNA-seq数据比对中蜂幼虫肠道参考转录组情况统计显示，各样品clean reads比对上参考转录组的比例均在88.91%以上（表2），AcCK与AcT的表达基因数目为36 604和37 108，分别占中蜂幼虫肠道参考转录组的84.05%和85.21 %。

图1 各中蜂幼虫肠道样本的PCA分析

表1 RNA-seq数据统计

表2 RNA-seq数据比对参考基因组情况统计

2.2 差异表达基因（DEG）分析

利用edgeR进行DEG分析，结果显示在球囊菌胁迫后，4日龄幼虫肠道共有583个DEGs，其中上调基因与下调基因的数量分别为344和239个（图2），说明有相当一部分基因受球囊菌胁迫而被激活，也有部分基因受到球囊菌的抑制。

2.3 差异表达基因的GO富集分析

DEGs的GO分类结果显示，这些DEGs分为3类：生物学进程、细胞组分和分子功能，分布于36个GO term上，分别为代谢进程、生长、发育进程、生殖、生殖进程、运动、细胞成分组织或生物起源、免疫系统进程、多组织进程、定位、单组织进程、细胞进程、应激反应、信号、生物调控、结构分子活性、电子载体活性、抗氧化活性、脒基-核苷酸交换因子活性、蛋白结合转录因子活性、催化活性、分子功能调节器、核酸结合转录因子活性、转运器活性、结合、分子传感器活性、细胞、细胞组件、大分子复合物、细胞器组件、细胞器、突触、细胞膜内腔、细胞膜和细胞膜组件（图3）。GO富集分析结果显示，基因富集数最多的GO term为细胞（106 unigenes），其次为细胞组件（106 unigenes）和代谢进程（104 unigenes）。

2.4 差异表达基因的KEGG代谢通路富集分析

分别对上调基因和下调基因进行KEGG代谢通路富集分析，结果表明上调基因富集在72个代谢通路上，其中基因富集数最多的是核糖体（ribosome）（72 unigenes），其次为碳代谢（carbon metabolism）（16 unigenes）和糖酵解（glycolysis/gluconeogenesis）（14 unigenes）（图4-A），而下调基因富集在45个代谢通路上，其中基因富集数最多的是碳代谢（9 unigenes），其次为二羧酸代谢（glyoxylate and dicarboxylate metabolism）（8 unigenes）和氨基酸生物合成（biosynthesis of amino acids）（7 unigenes）（图4-B）。值得关注的是，上调基因中分别有3、2、1和1个基因富集在内吞作用（endocytosis）、吞噬体（phagosome）、MAPK信号通路（MAPK signaling pathway）和泛素介导的蛋白水解（ubiquitin mediated proteolysis），说明意蜂幼虫的免疫相关通路被球囊菌胁迫所激活；多达32个下调基因富集在新陈代谢相关通路，如碳代谢（9 unigenes）、氧化磷酸化（oxidative phosphorylation）（5 unigenes）、精氨酸生物合成（arginine biosynthesis）（3 unigenes）、氮代谢（nitrogen metabolism）（2 unigenes）和丙酮酸代谢（pyruvate metabolism）（2 unigenes），说明对中蜂幼虫的新陈代谢相关通路受到球囊菌的抑制，球囊菌入侵对宿主的代谢系统产生广泛的影响。

图2 球囊菌胁迫中蜂幼虫肠道的DEGs

1：生物调控 Biological regulation；2：细胞成分组织或生物合成 Cellular component organization or biosynthesis；3：细胞进程 Cellular process；4：发育进程 Developmental process；5：生长 Growth；6：免疫系统进程 Immune system process；7：定位 Localization；8：运动 Locomotion；9：代谢进程 Metabolic process；10：多组织进程 Multi-organism process；11：多细胞生物进程 Multicellular organismal process；12：生殖 Reproduction；13：生殖进程 Reproductive process；14：应激反应 Response to stimulus；15：信号Signaling；16：单组织进程 Single-organism process；17：细胞 Cell；18：细胞组件 Cell part；19：大分子复合物 Macromolecular complex；20：细胞膜 Membrane；21：细胞膜组件 Membrane part；22：细胞膜内腔 Membrane-enclosed lumen；23：细胞器Organelle；24：细胞器组件 Organelle part；25：突触 Synapse；26：抗氧化活性 Antioxidant activity；27：结合 Binding；28：催化活性 Catalytic activity；29：电子载体活性 Electron carrier activity；30：脒基核苷酸交换因子活性 Guanyl-nucleotide exchange factor activity；31：分子功能调节器 Molecular function regulator；32：分子转换器活性 Molecular transducer activity；33：核酸结合转录因子活性 Nucleic acid binding transcription factor activity；34：蛋白结合转录因子活性 Protein binding transcription factor activity；35：结构分子活性 Structural molecule activity；36：转运器活性 Transporter activity

A：上调基因Up-regulated genes；B：下调基因Down-regulated genes

进而对上述免疫相关基因和代谢相关基因的表达模式进行分析，结果显示球囊菌胁迫后，中蜂幼虫肠道的免疫相关基因表达量呈不同程度的上调（图5-A），而代谢相关基因表达量不同程度地下调（图5-B），进一步表明中蜂幼虫肠道的免疫相关基因受球囊菌胁迫而激活、代谢相关基因受到球囊菌的抑制。

2.5 RNA-seq数据的qRT-PCR验证

利用qRT-PCR检测随机选取的5个DEGs（4个上调基因，1个下调基因），结果显示这些DEGs的表达水平变化趋势与RNA-seq数据中的基因表达水平变化趋势一致（图6），证明了本研究转录组数据的可靠性。

A：中蜂幼虫肠道的免疫相关基因A. c. cerana larval gut’s immune-related DEGs；B：中蜂幼虫肠道的代谢相关基因A. c. cerana larval gut’s metabolism-related DEGs

A：类围食膜因子3-C前体表皮蛋白编码基因cuticular protein analogous to peritrophins 3-C precursor encoding gene；B：表皮蛋白1前体编码基因cuticular protein 1 precursor encoding gene；C：类表皮蛋白同种型X2编码基因cuticle protein-like isoform X2 encoding gene；D：类表皮蛋白38编码基因cuticle protein 38-like encoding gene；E：G蛋白耦联受体Mth-like 3编码基因probable G-protein coupled receptor Mth-like 3 encoding gene

3 讨论

白垩病对蜂群的危害严重，每年给养蜂业造成巨大损失。目前，白垩病的组学研究十分有限。前期研究中，笔者课题组组装了中蜂幼虫肠道的参考转录组，并对其进行了功能及代谢通路注释[31]。在此基础上，本研究利用RNA-seq技术对健康及球囊菌胁迫的中蜂4日龄幼虫肠道进行深度测序，进而研究其响应球囊菌胁迫的应答。肠道是昆虫的重要器官，在消化、发育和免疫等诸多方面中发挥重要作用。本研究针对中蜂幼虫肠道进行测序，其转录组变化能更为精确地反映宿主响应病原胁迫的应答。在AcCK和AcT中分别检测到已知基因36 604和37 108个，约占中蜂幼虫肠道转录组基因总数的84.05%和85.21%，推测部分未检测到的基因可能在中蜂4日龄幼虫肠道不表达。

养蜂生产中，意蜂幼虫极易被球囊菌侵染而暴发白垩病，而中蜂幼虫极少受其影响。前期研究发现在球囊菌胁迫后的意蜂幼虫肠道中共检测到24个DEGs，其中上调基因与下调基因分别为4和20个（未发表数据）。本研究发现球囊菌胁迫后的中蜂幼虫肠道有583个基因差异表达（上调基因344个，下调基因239个），远多于球囊菌胁迫后的意蜂幼虫肠道，这些DEGs或许与中蜂的球囊菌抗性密切相关。此前，大多数研究都集中于蜜蜂的群体防御，通常认为清理行为在蜜蜂抵御球囊菌入侵过程中至关重要。前期试验结果表明，实验室条件下对中蜂幼虫和意蜂幼虫接种球囊菌，两种幼虫皆能被侵染而发生白垩病，但前者的发病率远低于后者，说明在中蜂幼虫与意蜂幼虫在个体水平上对球囊菌也存在较大的抗性差异。角质层和围食膜是昆虫免疫防御的第一道防线[38]，当病原微生物突破这道防线后，将遭遇昆虫体内细胞免疫与体液免疫的抵抗，包括内吞作用、黑化作用、吞噬作用、蛋白的酶促水解及抗菌肽等[39-40]。内吞作用和吞噬作用在蜜蜂抵御真菌侵染的过程中发挥着重要作用[40]。本研究发现分别有3个和2个上调基因富集在内吞作用和吞噬体，说明这两个代谢通路受球囊菌胁迫而激活。泛素介导的蛋白水解对清除衰老和损伤细胞意义重大，同时也在调节细胞进程的诸多方面发挥着关键作用，包括细胞周期、分化、发育和免疫等[41]。本研究发现1个上调基因富集在泛素介导的蛋白水解，说明该通路同样被球囊菌所激活。上述结果表明中蜂幼虫的部分细胞免疫对球囊菌的胁迫产生应答，推测它们在球囊菌入侵的早期发挥重要作用。昆虫的体液免疫也会因病原微生物的入侵而被激活，从而合成分泌抗菌肽杀灭病原，但在本研究中，未发现有DEGs富集在体液免疫如NF-B和Jak-STAT信号通路，也未发现蜜蜂4种抗菌肽（abaecin、apidaecin、defensin和hymenoptaecin）编码基因的差异表达，表明中蜂幼虫的体液免疫在球囊菌胁迫早期尚未被球囊菌激活。碳代谢对生物体内氨基酸的合成与转化至关重要。本研究发现在球囊菌胁迫的中蜂4日龄幼虫肠道中，富集在碳代谢上的DEGs包括16个上调基因和9个下调基因，推测多数碳代谢相关基因上调表达可能与宿主的细胞免疫水平增强有关，而部分数碳代谢相关基因表达水平下调反映出球囊菌-中蜂幼虫间互作的复杂性。本研究还发现多达72个上调基因富集在核糖体，说明在球囊菌胁迫早期，宿主的蛋白合成活跃，推测宿主通过提升蛋白合成满足抵御球囊菌入侵及自我修复的需要。

中蜂幼虫在球囊菌入侵早期并不表现出明显的白垩病症状，但宿主和病原之间发生着复杂的互作，本研究发现有多达76个DEGs富集在新陈代谢（32个新陈代谢相关通路），本研究的结果表明此时中蜂幼虫的代谢系统已受到球囊菌一定程度地干扰和破坏，推测病原通过抑制宿主的新陈代谢减少营养物质的消耗、降低宿主的免疫应答水平，从而促进病原自身的侵染。在球囊菌入侵早期，这种干扰和破坏作用的累积，可能会导致入侵晚期白垩病的发生。此前的白垩病研究集中于意蜂且聚焦在病程晚期，忽略了病程早期的研究。本研究利用RNA-seq技术开展中蜂幼虫肠道响应球囊菌胁迫的转录组学研究，发现在球囊菌胁迫的早期阶段，中蜂幼虫肠道也伴随着复杂的应答。目前尚无一种杀真菌剂被批准应用于养蜂生产[1]，养蜂生产实践中主要通过选育抗病品系、改善养蜂管理和保持清洁卫生来防治白垩病[41]，但是效果并不理想，亟需有效的防治策略。未来的研究方向是利用RNA-seq技术对球囊菌胁迫后的不同日龄中蜂幼虫肠道进行深度测序，全局研究中蜂幼虫在转录组水平的应答，并应用趋势分析或WGCNA分析筛选出关键应答基因，进而利用RNAi等技术手段验证其功能。

4 结论

利用RNA-seq技术对健康及球囊菌胁迫的中蜂4日龄幼虫肠道进行深度测序，通过对DEGs的GO和KEGG富集分析，发现在球囊菌胁迫早期，中蜂幼虫的部分细胞免疫被球囊菌激活，而体液免疫未被激活，宿主的新陈代谢系统受到球囊菌的显著抑制。研究结果为深入解析中蜂幼虫响应球囊菌胁迫的应答机制提供了重要信息，也为白垩病的有效防治打下了基础。

References

[1] GALIZIA C G, EISENHARDT D, GIURFA M, MENZEL R.. Dordrecht Netherlands. New York: Springer, 2012.

[2] BEGNA D, Han B, FENG M, FANG Y, LI J. Differential expression of nuclear proteomes between honeybee (L.) queen and worker larvae: a deep insight into caste pathway decisions., 2012, 11(2): 1317-1329.

[3] ZAYED A, ROBINSON G. Understanding the relationship between brain gene expression and social behavior: lessons from the honey bee., 2012, 46(6): 591-615.

[4] FORET S, KUCHARSKI R, PELLEGRINI M, FENG S H, JACOBSEN S E, ROBINSON G E, MALESZKA R. DNA methylation dynamics, metabolic fluxes, gene splicing, and alternative phenotypes in honey bees., 2012, 109(13): 4968-4973.

[5] KURZE C, ROUTTU J, MORITZ R F. Parasite resistance and tolerance in honeybees at the individual and social level., 2016, 119(4): 290-297.

[6] Committee on the Status of Pollinators in North America.. National Academies Press, 2007.

[7] AIZEN M A, GARIBALDI L A, CUNNINGHAM S A, KLEIN A M. How much does agriculture depend on pollinators? Lessons from long-term trends in crop production., 2009, 103(9): 1579-1588.

[8] SPILTOIR C F. Life cycle of., 1955, 42(6): 501-518.

[9] LUMBSCH H T, Huhndorf S M. Outline of Ascomycota., 2007, 13: 1-58.

[10] BISSETT J. Contribution toward a monograph of the genus., 1988, 66(12): 2541-2560.

[11] CHORBIŃSKI P. Enzymatic activity of strains of., 2003, 59(11): 1019-1022.

[12] BAMFORD S, HEATH L A F. The effects of temperature and pH on the germination of spores of the chalkbrood fungus,., 1989, 28(1): 36-40.

[13] WINSTON M L. The biology of the honey bee//. Harvard University Press, Cambridge, USA, 1991.

[14] FLORES J M, SPIVAK M, GUTIÉRREZ I. Spores ofcontained in wax foundation can infect honeybee brood., 2005, 108(1/2): 141-144.

[15] BAILEY L.. Academic Press, London, UK, 1991.

[16] THEANTANA T, CHANTAWANNAKUL P. Protease and beta-N- acetylglucosaminidase of honey bee chalkbrood pathogen., 2008, 47(1): 68-76.

[17] EVANS J D, SPIVAK M. Socialized medicine: individual and communal disease barriers in honey bees.,2010, 103(Suppl. 1): S62-S72.

[18] TANJI T, HU X, WEBER A N, IP Y. Toll and IMD pathways synergistically activate an innate immune response in., 2007, 27(12): 4578-4588.

[19] HORNITZKY M. Literature review of chalkbrood. A report for the RIRDC. Publication No. 01/150, Kingston, ACT, AU. 2001.

[20] 梁勤, 陈大福, 王建鼎. 营养生态条件对蜜蜂球囊菌生长及产孢的影响. 中国生态农业学报,2001, 9(4): 31-34.

LIANG Q, CHEN D F, WANG J D. Effects on the mycelia growth and spore-forming ofunder ecological condition of nutrients.,2001, 9(4): 31-34. (in Chinese)

[21] 李江红, 郑志阳, 陈大福, 梁勤. 影响蜜蜂球囊菌侵染蜜蜂幼虫的因素及侵染过程观察. 昆虫学报, 2012, 55(7): 790-797.

LI J H, ZHENG Z Y, CHEN D F, LIANG Q. Factors influencinginfection on honeybee larvae and observation on the infection process., 2012, 55(7): 790-797. (in Chinese)

[22] 熊翠玲, 陈大福, 付中民, 马晓云, 梁勤. 蜜蜂球囊菌检测分子标记的灵敏性测定. 中国蜂业, 2010, 61(11): 16-18.

XIONG C L, CHEN D F, FU Z M, MA X Y, LIANG Q. Determination of the sensitivity of DNA molecular marker for detectings., 2010, 61(11): 16-18. (in Chinese)

[23] 席伟军, 李江红, 陈大福, 梁勤. 环介导等温扩增(LAMP)技术检测蜜蜂球囊菌. 中国农业科学, 2016, 49(4): 765-774.

XI W J, LI J H, CHEN D F, LIANG Q. Diagnosis of theby the loop-mediated isothermal amplification., 2016, 49(4): 765-774. (in Chinese)

[24] 郑志阳, 李江红, 梁勤, 陈大福. 蜜蜂球囊菌分泌多种胞外酶侵染蜜蜂幼虫. 福建农林大学学报 (自然科学版), 2011, 40(3): 280-284.

ZHENG Z Y, LI J H, LIANG Q, CHEN D F.secretes multiple extracellular enzymes to infect honeybee larvae., 2011, 40(3): 280-284. (in Chinese)

[25] Honeybee Genome Sequencing Consortium. Insights into social insects from the genome of the honeybee., 2006, 443(7114): 931-949.

[26] EVANS J D, ARONSTEIN K A, CHEN Y P, HETRU C, IMLER J L, JIANG H, KANOST M, THOMPSON G J, ZOU Z, HULTMARK D. Immune pathways and defense mechanisms in honey bees,., 2006, 15(5): 645-656.

[27] ARONSTEIN K A, MURRAY K D, SALDIVAR E. Transcriptional responses in honey bee larvae infected with chalkbrood fungus., 2010, 11: 391.

[28] CORNMAN R S, LOPEZ D, EVANS J D. Transcriptional response of honey bee larvae infected with the bacterial pathogen., 2013, 8(6): e65424.

[29] JULIE A, BARBARA M A, BERNARD V, CATHERINE T, FRÉDÉRIC D, NICOLAS B. Transcriptome analyses of the honeybee response toand insecticides., 2014, 9(3): e91686.

[30] PARK D, JUNG J W, CHOI B S, JAYAKODI M, LEE J, LIM J, YU Y, CHOI Y S, LEE M L, PARK Y, CHOI I Y, YANG T J, EDWARDS O R, NAH G, KWON H W. Uncovering the novel characteristics of Asian honey bee,, by whole genome sequencing., 2015, 16: 1.

[31] 徐细建, 郭睿, 骆群, 熊翠玲, 梁勤, 张串联, 郑燕珍, 张曌楠, 黄枳腱, 张璐, 李汶东, 陈大福. 中华蜜蜂幼虫肠道参考转录组的组装及SSR分子标记鉴定. 中国农业科学, 2017, 50(6): 1157-1166.

Xu X J, Guo R, Luo Q, Xiong C L, Liang Q, Zhang C L, Zheng Y Z, Zhang Z N, Huang Z J, Zhang L, Li W D, Chen D F.transcriptome assembly forlarval gut and identification of SSR molecular markers., 2017, 50(6): 1157-1166. (in Chinese)

[32] JENSEN A B, ARONSTEIN K, FLORES J M, VOJVODIC S, PALACIO M A, SPIVAK M. Standard methods for fungal brood disease research., 2013, 52(1): 79-88.

[33] 王倩, 孙亮先, 肖培新, 刘锋, 康明江, 胥保华. 室内人工培育中华蜜蜂幼虫技术研究. 山东农业科学, 2009(11): 113-116.

WANG Q, SUN L X, XIAO P X, LIU F, KANG M J, XU B H. Study on technology for indoor artificial feeding oflarvae., 2009(11): 113-116. (in Chinese)

[34] LANGMEAD B, TRAPNELL C, POP M, SALZBERG S L. Ultrafast and memory-efficient alignment of short DNA sequences to the human genome., 2009, 10(3): R25.

[35] HURGOBIN B. Short read alignment using SOAP2//. New York, NY , United States: Humana Press, 2016: 241-252.

[36] ROBINSON M D, MCCARTHY D J, SMYTH G K. EdgeR: a bioconductor package for differential expression analysis of digital gene expression data., 2010, 26(1): 139-140.

[37] LIVAK K J, SCHMITTGEN T D. Analysis of relative gene expression data using real-time quantitative PCR and the 2-∆∆Ctmethod., 2001, 25: 402-408.

[38] ORIHEL T C. The peritrophic membrane: its role as a barrier to infection of the arthropod host//Maramorosch K, Shope R E.. Academic Press, New York, USA, 1975: 67-73.

[39] GLIŃSKI Z, JAROSZ J. Infection and immunity in the honey bee., 2001, 36: 12-24.

[40] GLINSKI Z, BUCZEK K. Response of theto fungal infections., 2003, 38: 183-189.

[41] McBride W H, Iwamoto K S, Syljuasen R, Pervan M, Pajonk F. The role of the ubiquitin/proteasome system in cellular responses to radiation., 2003, 22: 5755-5773.

（责任编辑 岳梅）

A: AcCK-1; B: AcCK-2; C: AcCK-3; D: AcT-1; E: AcT-2; F: AcT-3

附图1 各中蜂幼虫肠道样品的测序饱和度

Fig. S1 Sequencing saturation of everylarval gut sample

附图2 各中蜂幼虫肠道样品不同生物学重复间的相关性

Fig. S2 Pearson correlation between every two biological repeats within eachlarval gut sample

Transcriptome oflarval gut Underthe Stressof

CHEN DaFu, GUO Rui, XIONG CuiLing, LIANG qin, ZHENG YanZhen, XU XiJian, ZHANG ZhaoNan, HUANG ZhiJian, ZHANG Lu, WANG HongQuan, XIE YanLing, TONG XinYu

(College of Bee Science, Fujian Agriculture and Forestry University, Fuzhou 350002)

【Objective】So far, no study on application of next-generation sequencing technology for the research of chalkbrood disease was reported. In the present research, RNA-seq technology was utilized to deep sequence of normal and-treated 4th instarlarval gut to mine larvae’s responses tochallenge.【Method】AcCK (un-treated group) and AcT (-treated group) were sequenced on Illumina HiSeq 2500 platform. After evaluation and filtration of raw data from RNA-seq, differentially expressed gene (DEG) analysis was performed using edgeR software, further, Gene Ontology (GO) and KEGG pathway enrichment analyses were carried out, and finally, real-time quantitative PCR (qRT-PCR) was conducted to validate the RNA-seq data.【Result】 In total, RNA-seq produced 188 457 338 raw reads, and after filtration, 182 088 448 clean reads were obtained, Q20 and Q30 of each sample were above 97.96% and 94.97%, respectively, indicating that RNA-seq data in this research were with high quality. principle component analysis (PCA) was performed on all genes level and the result showed PC1 and PC2 were able to account for 75.8% and 10.7% of the expressed genes’ overall differences, respectively. DEG analysis result displayed that there were 344 up-regulated genes and 239 down-regulated genes in AcCK VS AcT. GO enrichment analysis result showed that the DEGs were enriched in 36 GO terms, among them, the mostly enriched ones were cell (106 unigenes), cell part (106 unigenes) and metabolic process (104 unigenes). KEGG pathway enrichment analysis result suggested that up- and down-regulated genes were enriched in 72 and 45 pathways, respectively, and the mostly enriched pathways for up-regulated genes were ribosome (72 unigenes), carbon metabolism (16 unigenes) and glycolysis/gluconeogenesis (14 unigenes), while the mostly enriched pathways for down-regulated genes were carbon metabolism (9 unigenes), glyoxylate and dicarboxylate metabolism (8 unigenes) and amino acids biosynthesis (7 unigenes). further analysis demonstrated that the immune-related genes inlarval gut were activated, while the metabolism-related genes were greatly inhibited.【Conclusion】The findings of the study not only uncovered thelarval gut’s responses toduring the early stage of invasion, but also provided key information for clarifying the mechanism underlying the host’s responses to, thus laying a foundation for functional investigation of key responding genes.

; larval gut;; transcriptome; RNA-seq

2016-12-26；接受日期：2017-03-06

国家自然科学基金（30800806）、国家现代农业产业技术体系（蜜蜂）建设专项（CARS-45-KXJ7）、福建农林大学科技发展资金（KF2015123）

陈大福，Tel：0591-83726835；E-mail：dfchen826@fafu.edu.cn。通信作者郭睿，Tel：0591-87640197；E-mail：ruiguo@fafu.edu.cn