胁迫中华蜜蜂幼虫肠道的球囊菌及其体外培养的高表达基因分析

陈大福, 王鸿权, 李汶东, 熊翠玲, 郑燕珍, 付中民, 徐细建, 黄枳腱, 郭 睿(福建农林大学蜂学学院，福建 福州 350002)

胁迫中华蜜蜂幼虫肠道的球囊菌及其体外培养的高表达基因分析

陈大福, 王鸿权, 李汶东, 熊翠玲, 郑燕珍, 付中民, 徐细建, 黄枳腱, 郭 睿
(福建农林大学蜂学学院，福建 福州 350002)

利用RNA-seq技术对胁迫中华蜜蜂(简称中蜂)6日龄幼虫肠道的球囊菌及其体外培养进行深度测序，根据基因的FPKM值筛选得到高表达基因(HEGs)，进而对HEGs进行GO及KEGG代谢通路(pathway)富集分析.Illumina测序数据经质控和过滤得到100814558条有效读段(clean reads)，平均Q30均在93.81%以上.GO富集分析结果显示处理组(AacT)的HEGs富集于39个GO term，基因富集数最多的是细胞(1872 unigene)、细胞组件(1872 unigene)和细胞进程(1748 unigenes)；对照组(AaCK)的HEGs富集于37个GO term，基因富集数最多的是细胞(785 unigenes)，其次是细胞组件 (785 unigenes)和代谢进程(776 unigenes).KEGG pathway富集分析显示，AacT的HEGs富集在119个代谢通路上，基因富集数最多的是核糖体 (176 unigenes)、碳代谢(147 unigenes)及氨基酸的生物合成(134 unigenes)；AaCK的HEGs富集在110个代谢通路上，基因富集数最多的是核糖体(179 unigenes)、氨基酸的生物合成(70 unigenes)以及碳代谢(62 unigenes).深入分析发现，对于富集在MAPK信号通路上的高表达基因数量，胁迫中蜂幼虫肠道的球囊菌远多于体外培养的球囊菌，说明病原的该通路在胁迫后期被显著激活.

中华蜜蜂; RNA-seq; 球囊菌; 高表达基因

Abstract: In this study,Ascosphaeraapisspores (AacT), from the gut of 6-day-oldApisceranaceranalarvae, andA.c.ceranawere sequenced utilizing RNA-seq technology. PureA.apisspore was used as the control (AaCK). A total of 100 814 558 clean reads with a mean Q30 of 93.81% were obtained after quality control and filtration. GO enrichment analysis suggested that the highly expressed genes (HEGs) in AacT were involved in 39 terms, which was mostly enriched in cell (1872 unigenes), followed by cellular component (1872 unigenes) and cellular process (1748 unigenes). And the HEGs in AaCK were engaged in 37 terms, with the largest group in cell (785 unigenes), followed by cellular component (785 unigenes) and metabolic process (776 unigenes). KEGG enrichment analysis showed that the HEGs in AacT were enriched in 119 pathways, dominating in ribosme (176 unigenes), followed by carbon metabolism (147 unigenes) and biosynthesis of amino acids (134 unigenes). While the HEGs in AaCK were involved in 119 pathways, and the largest groups were ribosome (179 unigenes), biosynthesis of amino acids (70 unigenes) and carbon metabolism (62 unigenes).

Keywords:Apisceranacerana; RNA-seq;Ascosphaeraapis; highly expressed genes

白垩病是蜜蜂幼虫的致死性真菌病[1]，该病病原为蜜蜂球囊菌(Ascosphaeraapis).白垩病能导致蜜蜂群势急剧下降，严重影响蜂产品的产量[2].据报道，该病可造成蜂蜜产量下降5%-37%[3].我国养蜂生产的常用蜂种是意大利蜜蜂(Apismelliferaligustica，简称意蜂)与中华蜜蜂(Apisceranacerana，简称中蜂).前者极易被球囊菌侵染而罹患白垩病，但后者几乎不受白垩病的影响.通常认为中蜂较意蜂具有更强的群体防御能力[4]，如清理行为可将被球囊菌侵染的幼虫及时清理出蜂群[5].

近20年来，国内外学者对球囊菌开展了大量研究，主要涉及培养方法[6-7]、形态学[8-9]、流行病学[10-11]、增殖方式[12-13]、病理学[14-15]、免疫防御[16-17]以及疾病防治[18-20]等方面.本课题组也在球囊菌的生理、病理和检测方面开展了较多研究[21].基因组信息的缺失严重制约球囊菌的分子研究.直到2006年，Qin et al[22]测序公布了球囊菌ARSEF7406菌株的基因组序列信息，为在分子水平上研究球囊菌奠定了基础，但作者当时并没有公布球囊菌的基因位置和功能注释信息，导致近10年来球囊菌的分子研究进展缓慢.Cornman等利用Roche 454焦磷酸测序技术对来自培养基的球囊菌菌丝和来自西方蜜蜂(Apismellifera)幼虫感染组织的球囊菌菌丝进行了转录组测序，功能分析表明球囊菌的差异表达基因(DEGs)参与了球囊菌交配类型、细胞内信号转导以及应激反应[23].

球囊菌侵染的相关研究多集中在西方蜜蜂幼虫，有关球囊菌侵染东方蜜蜂幼虫(Apiscerana)的研究报道极其有限.前期研究中，我们从中蜂白垩状幼虫尸体中分离培养出真菌病原，经形态学和分子生物学鉴定该病原即为球囊菌(未发表数据).目前，有关球囊菌对中蜂幼虫的致病机理的研究未见报道.本研究利用RNA-seq技术对体外培养及胁迫中蜂6日龄幼虫肠道的球囊菌进行转录组学分析，筛选得到对照组和处理组球囊菌的高表达基因(differentially expressed genes, HEGs)，进一步通过GO(Gene Ontology)及KEGG(Kyoto Encyclopedia of Genes and Genomes)代谢通路(pathway)富集分析对HEGs进行深入挖掘，研究结果可为球囊菌致病的分子机理提供参考信息，也为白垩病的治疗奠定初步基础.

1 材料与方法

1.1 生物材料

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

1.2 主要实验试剂及仪器

DNaseI和Oligotex mRNA Kits Midi试剂盒购自德国Qiagen公司，Dynal M280磁珠购自Invitrogen公司，DNA ligase购自美国Thermo公司，RNA Reagent抽提试剂盒、Ex Taq polymerase及Superscript Ⅱ reverse transcriptase均购自日本TaKaRa公司.纯化cDNA的Ampure beads为美国Agencourt产品，cDNA文库构建试剂盒TruSeqTMDNA Sample Prep Kit-Set A为美国Illumina公司产品.其它试剂为国产分析纯.

恒温恒湿气候箱购自中国宁波江南仪器厂，倒置显微镜为中国上海光学仪器五厂产品，超净工作台为中国苏州安泰空气技术有限公司产品，PCR仪为美国Bio Rad公司产品，凝胶成像系统为中国上海培清科技有限公司产品，超低温冰箱为中科美菱低温科技股份有限公司产品.

1.3 试验方法

1.3.1 中蜂幼虫的人工饲养及球囊菌孢子纯化 中蜂幼虫的饲料参照王倩等的配方[24]配制并进行改良.预试验结果显示中蜂幼虫7日龄成活率达到70%以上.选择20群健康蜂群(无白垩病症状且PCR检测为阴性)作为供试蜜蜂.将2日龄幼虫移入预置饲料的24孔细胞培养板(1只·孔-1)， 35 ℃，70% RH条件下培养.每隔24 h更换饲料.3日龄幼虫饲喂含球囊菌孢子的人工饲料(孢子终浓度为1×107孢子·mL-1) ，适宜条件饲养幼虫至6日龄.本试验设3个生物学重复.

1.3.2 测序样品准备及Illumina测序 按照Jesen et al[25]的方法纯化球囊菌孢子(AaCK)，纯化孢子经液氮速冻后转移至-80 ℃超低温冰箱保存备用.同时，剖取6日幼虫肠道(AacT)，立即投入液氮速冻，再转移至-80 ℃保存备用.利用RNAiso Reagent试剂盒抽提上述3个球囊菌孢子样品(AaCK-1、AaCK-2、AaCK-3)与3个幼虫肠道样品(AacT-1、AacT-2、AacT-3)的总RNA，然后用RNase-free DNaseI去除基因组DNA残留.RNA的质量通过琼脂糖凝胶电泳和NanoDrop ND-1000(NanoDrop, Wilmington, DE, USA)进行检测.cDNA文库的构建参照张曌楠等[26]的方法进行.

委托广州基迪奥生物科技有限公司对上述6个样品进行深度测序，测序平台为Illumina Hiseq 2500.测序数据已上传NCBI SRA数据库，SRA号：SRA454366.

1.3.3 基因表达丰度计算与高表达基因的分析 对于测序下机的原始读段(raw reads)，去除含adapter、N比例大于10%及低质量的reads，得到有效读段(clean reads)，利用Tophat软件将clean reads比对至球囊菌的参考转录组[26]，表达量的计算使用FPKM(Fragments Per Kilobase of transcript per Million mapped reads)法，其计算公式为：FPKM=(106C×103)/NL.利用在线分析工具OmicsShare(http://www.omicshare.com/tools/index.php/Home/Index/index.html)对AaCK与AacT中FPKM值大于15的HEGs进行GO及KEGG代谢通路富集分析.

2 结果与分析

2.1 RNA-seq测序数据概述

本研究中，各样品的RNA-seq得到的clean reads数均在4398766条以上，两端Q20均在97.05%以上，两端Q30均在92.55%以上(表1).将clean reads映射(mapping)核糖体数据库和球囊菌参考转录组，各样品mapping率均在87%以上.Pearson相关性分析结果显示AaCK与AacT的组内各生物学重复之间的相关性均在0.8134以上(0.8134～0.9999)，说明样本的重复性较好.进一步对AaCK与AacT进行主成分分析(PCA)，结果显示第1主成分(PC1)与第2主成分(PC2)可分别解释样品基因表达总体方差的97.9%和1.5%，AaCK与AacT的基因表达模式差异明显且各样品的组内重复聚类良好.上述结果说明本研究中的转录组数据质量良好，可用于进一步分析.

表1 RNA-seq数据信息统计表Table 1 Summary of RNA-seq datasets

2.2 HEGs的GO富集分析

GO富集分析结果显示，AaCK与AacT的HEGs均富集于3个大类，即生物学进程、细胞组分和分子功能.其中，AaCK的HEGs富集于37个GO term，基因富集数最多的是细胞(785 unigenes)、细胞组件(785 unigenes)和细胞进程(776 unigenes)；AacT的HEGs富集于39个GO term，基因富集数最多的是细胞(1872 unigenes)、细胞组件(1872 unigenes)和细胞进程(1748 unigenes).上述结果显示，AacT的HEGs较之AaCK的多数GO term富集基因数均有较大幅度提高，说明胁迫后期的球囊菌的新陈代谢活跃、生命活动增强.

2.3 HEGs的KEGG代谢通路富集分析

KEGG代谢通路富集分析结果显示，AaCK的HEGs富集在110个代谢通路，其中基因富集数最多的是核糖体(179 unigenes)，其次是氨基酸的生物合成(70 unigenes)和碳代谢(62 unigenes)(图1A).AacT的HEGs富集在119个代谢通路，其中基因富集数最多的是核糖体(176 unigenes)，其次是碳代谢(147 unigenes)和氨基酸的生物合成(134 unigenes)(图1B).

进一步分析结果显示，AaCK与AacT中均有HEGs富集在MAPK信号通路(图2)，但AacT中富集在此通路的HEGs数量(27 unigenes)远多于AaCK(9 unigenes)，说明球囊菌的MAPK信号通路在胁迫后期被显著激活.

2.4 HEGs的Venn分析

进一步对AaCK和AacT中的HEGs进行Venn分析，结果显示 AaCK和AacT的共有HEGs为459个，二者的特有HEGs分别为1962和6185个.推测这些共有HEGs对球囊菌的生命活动至关重要，AacT特有HEGs在球囊菌胁迫后期发挥特殊作用.

A：AaCK中的HEGs；B：AacT中的HEGs.图1 HEGs的KEGG pathway富集分析Fig.1 KEGG pathway enrichment analysis of HEGs

3 讨论与结论

为对球囊菌进行深入的转录组学研究，本课题组前期已组装并注释球囊菌的参考转录组[23]，在此基础上，本研究对胁迫中蜂6日龄幼虫肠道的球囊菌及其体外培养进行测序，进而对不同状态下的球囊菌进行HEGs分析.在Cornman et al[23]的研究中，测序对象是来自培养基的球囊菌菌丝和来自蜜蜂幼虫感染组织的球囊菌菌丝，其中后者是感染幼虫肠道内萌发的菌丝突破肠壁后继续生长最后突破体表的菌丝，而此时的菌丝不再与宿主互作，菌丝的转录组变化无法真实反映球囊菌在侵染过程中的基因表达谱.本研究的测序对象是纯培养的球囊菌和球囊菌胁迫的中蜂幼虫肠道，前者尚未与宿主发生互作，而后者因处于中蜂幼虫肠道内而与宿主发生复杂的互作，因此本研究中的转录组分析可更真实地反映球囊菌在侵染过程中的基因表达谱.此外，Cornman等人的研究是对人工培养的球囊菌菌丝和幼虫感染组织的菌丝进行DEGs分析，因此，无法将他们的研究结果与本研究的结果直接对比.在球囊菌侵染中蜂幼虫的后期，菌丝已大量生长并快速穿透肠壁，技术上无法将宿主和病原完全分离.因此，为得到纯净的球囊菌转录组数据，本研究采取的策略是对胁迫肠道进行测序，得到意蜂幼虫肠道、肠道内球囊菌以及肠道内其它微生物的混合转录组数据，首先通过映射意蜂参考基因组去除宿主本身的转录组数据，再将球囊菌和肠道内微生物的混合数据映射前期组装并注释的球囊菌参考转录组，去除比对不上的数据，保留能够比对上的数据，即为球囊菌本身的转录组数据.

本研究中，AacT中有60.74%的HEGs(2098 unigenes)富集在新陈代谢, 包括碳代谢(ko01200)、脂肪酸降解(ko00071)以及甘氨酸、丝氨酸及苏氨酸代谢(ko00260)等107个物质代谢，还包括氧化磷酸化(ko00190)、氮代谢(ko00910)、硫代谢(ko00920)以及甲烷代谢(ko00680)等4个能量代谢通路.而对照组只有55.18%的HEGs(858 unigenes)富集在代谢进程中，表明球囊菌在胁迫中蜂幼虫肠道的后期，通过提高自身的物质和能量代谢满足增殖需求.真菌具有精细而保守的信号通路，可通过调节黑色素形成、配合、毒力以及形态发生等来感知外界环境并迅速做出应答[23].在真菌的信号通路中，MAPK信号通路至关重要，在真菌致病性及响应环境胁迫方面发挥关键作用[27].本研究发现，胁迫后期的球囊菌中有27个HEGs富集在MAPK信号通路，较球囊菌纯化孢子有大幅提高，表明该通路在胁迫后期的球囊菌中非常活跃，说明MAPK信号通路在球囊菌侵染宿主的过程中可能发挥关键作用.下一步将人工合成siRNA对富集在此通路上的高表达基因进行RNAi，有望为白垩病的治疗提供潜在靶标.

本研究只对胁迫中蜂6日龄幼虫肠道的球囊菌及其体外培养进行HEGs分析，这是从基因表达量的角度分析球囊菌的基因表达谱，若要全面解析球囊菌致病的分子机理，则需要从基因相对变化倍数的角度进一步挖掘球囊菌在胁迫中蜂幼虫肠道过程中的转录组信息，本课题组目前正在对胁迫中蜂幼虫肠道的球囊菌进行DEGs分析，接下来将对球囊菌胁迫其它日龄的幼虫肠道进行转录组测序，并通过趋势分析或基因权重共表达分析(WGCNA)从全局水平进行深入研究.

本研究利用RNA-seq技术对胁迫中蜂幼虫肠道的球囊菌及其体外培养进行深度测序，筛选并深入分析各样品的HEGs，研究结果可为解析球囊菌致病的分子机理提供重要信息，为关键致病基因的功能验证奠定基础，有望为白垩病治疗提供潜在分子靶标.

[1] 陈淑静,康雪冬,冯峰.蜜蜂白垩病研究初报[J].中国蜂业,1994(6):4-7.

[2] BAILEY L. Infectious Diseases of the Honeybee[M]. London: Land Books, 1963.

[3] ZAGHLOUL O A, MOURAD A K, KADY M B E, et al. Assessment of losses in honey yield due to the chalkbrood disease, with reference to the determination of its economic injury levels in Egypt[J]. Communications in Agricultural & Applied Biological Sciences, 2005,70(4):703-714.

[4] INVERNIZZI C, RIVAS F, BETTUCI L. Resistance to chalkbrood disease inApismelliferaL. (Hymenoptera: Apidae) colonies with different hygienic behavior[J]. Neotropical Entomology, 2011,40(1):28-34.

[5] PARK D, JUNG JW, CHOI BS, et al. Uncovering the novel characteristics of Asian honey bee,Apiscerana, by whole genome sequencing[J]. BMC Genomics, 2015,16(1):1.

[6] LAF H. Chalk brood pathogens: a review[J]. Bee World, 2015,16(6):865-869.

[7] ANDERSON D L, GIBSON N L. New species and isolates of spore-cyst fungi (Plectomycetes: Ascosphaerales) from Australia[J]. Systematic Botany, 1998,11(1):53-72.

[8] SPILTOIR C F. Life cycle ofAscosphaeraapis(Pericystisapis)[J]. American Journal of Botany, 1955,42(6):501-508.

[9] SKOU J P. More details in support of the class Ascosphaeromycetes[J]. Mycotaxson, 1988,51(1):191-198.

[10] FLORES J M, SPIVAK M, GUTIERREZ I. Spores ofAscosphaeraapiscontained in wax foundation can infect honeybee brood[J]. Veterinary Microbiology, 2005,108(1-2):141-144.

[11] FLORES J M, GUTIERREZ I, ESPEJO R. The role of pollen in chalkbrood diseaseinApismellifera: transmission and predisposing conditions[J]. Mycologia, 2005,97(6):1 171-1 176.

[12] PÖGGELER S. Mating-type genes for classical strain improvements of Ascomycetes[J]. Applied Microbiology and Biotechnology, 2001,56(5-6):589-601.

[13] CHORBIÉSKI P. Enzymatic activity of strains ofAscosphaeraapis[J]. Medycyna Weterynaryjna-Veterinary Medicine-Science and Practice, 2003,59(11):1 019-1 022.

[14] BAILEY L, BALL B V. Honey Bee Pathology[M]. London & San Diego: Academic Press, London, 1991.

[15] THEANTANA T, CHANTAWANNAKUL P. Protease and β-N acetylglucosaminidase of honey bee chalkbrood pathogenAscosphaeraapis[J]. Journal of Apicultural Research, 2008,47(1):68-76.

[16] STANLEY D, MILLER J, TUNAZ H. Eicosanoid actions in insect immunity[J]. Innate Immunity, 2009,1(4):282-290.

[17] EVANS J D, SPIVAK M. Socialized medicine: individual and communal disease barriers in honey bees[J]. Journal of Invertebrate Pathology, 2010,103(1):62-72.

[18] HORNITZKY M. Literature review of chalkbrood, a fungal disease of honeybees[J]. Rural Industries Research and Development Corporation, 2001,1:150.

[19] TARPY D R. Genetic diversity within honeybee colonies prevents severe infections and promotes colony growth[J]. Proceedings of the Royal Society B-Biological Sciences, 2002,270,1510:99-103.

[20] MOURAD A K, ZAGHLOUL O A, KADY E L, et al. A novel approach for the management of the chalkbrood disease infesting honeybeeApismelliferaL. (Hymenoptera: Apidae) colonies in Egypt[J]. Communication of Agricultural and Applied Biology Science, 2005,70(4):601-611.

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

[22] QIN X, EVANS J D, ARONSTEIN K A, et al. Genome sequences of the honey bee pathogensPaenibacilluslarvae, andAscosphaeraapis[J]. Insect Molecular Biology, 2006,15(5):715-718.

[23] CORNMAN R S, BENNETT A K, MURRAY K D, et al. Transcriptome analysis of the honey bee fungal pathogen,Ascosphaeraapis: implications for host pathogenesis[J]. BMC Genomics, 2012,13(7):285.

[24] 王倩,孙亮先,肖培新,等.室内人工培育中华蜜蜂幼虫技术研究[J].山东农业科学,2009,11:113-116.

[25] JENSEN A B, ARONSTEIN K, FLORES J M, et al. Standard methods for fungal brood disease research[J]. Methods, 2013,52(1):1-20.

[26] 张曌楠,熊翠玲,徐细建,等.蜜蜂球囊菌的参考转录组denovo组装及SSR分子标记开发[J].昆虫学报,2017,60(1):34-44.

[27] PEREZ P, CANSADO J. Cell integrity signaling and response to stress in fission yeast[J]. Current Protein & Peptide Science, 2010,11(8):680-692.

(责任编辑:吴显达)

AnalysisofhighlyexpressedgenesofAscosphaeraapisinfectingthegutofApisceranaceranalarvaeanditsinvitroculture

CHEN Dafu, WANG Hongquan, LI Wendong, XIONG Cuiling, ZHENG Yanzhen, FU Zhongmin, XU Xijian, HUANG Zhijian, GUO Rui
(College of Bee Science, Fujian Agriculture and Forestry University, Fuzhou, Fujian 350002, China)

S895.1

A

1671-5470(2017)05-0562-07

10.13323/j.cnki.j.fafu(nat.sci.).2017.05.014

2017-03-14

2017-05-02

现代农业产业技术体系建设专项资金(CARS-44-KXJ7)；福建省中青年教师教育科研项目(JAT17018)；福建农林大学科技发展基金(KF2015123)；国家级大学生创新创业项目(201610389016).

陈大福(1973-),男,副教授,博士.研究方向:蜜蜂保护学.Email:dfchen826@fafu.edu.cn.通讯作者郭睿(1987-),男,讲师,博士.研究方向:蜜蜂分子生物学.Email:ruiguo@fafu.edu.cn.