中华蜜蜂幼虫肠道的高表达基因分析

解彦玲，王鸿权，张 璐，侯志贤，刀 晨，江亮亮，熊翠玲，郑燕珍，徐细建，黄枳腱，郭 睿，陈大福

(福建农林大学 蜂学学院，福建 福州 350002)

中华蜜蜂幼虫肠道的高表达基因分析

解彦玲，王鸿权，张 璐，侯志贤，刀 晨，江亮亮，熊翠玲，郑燕珍，徐细建，黄枳腱，郭 睿*，陈大福

(福建农林大学 蜂学学院，福建 福州 350002)

中华蜜蜂(中蜂)是中国特有的本土蜜蜂资源，也是养蜂生产中的常用蜂种之一，具有重要的生态和经济价值。本研究利用RNA-seq技术对中蜂5和6日龄幼虫肠道进行深度测序，得到的有效读段(clean reads)数分别为29290630与26636038，Q20分别为98.45%和98.38%。GO富集分析结果显示，Ac5中的高表达基因(HEGs)富集在44个GO term，其中基因富集数最多的为代谢进程(753 unigenes)，其次是催化活性(679 unigenes)与结合(657 unigenes)；Ac6中的HEGs富集在44个GO term，其中基因富集数最多的为代谢进程(792 unigenes)，其次是结合(760 unigenes)和催化活性(744 unigenes)。KEGG代谢通路(pathway)富集分析结果显示，Ac5的HEGs富集在129个pathway，其中基因富集数最多的是核糖体(103 unigenes)、内质网蛋白质加工(68 unigenes)和碳代谢(67 unigenes)；Ac6的HEGs富集在129个代谢通路，其中基因富集数最多的是核糖体(103 unigenes)、氧化磷酸化(88 unigenes)及内质网蛋白质加工(68 unigenes)。进一步对HEGs进行Venn分析，结果显示，二者共有的HEGs为2 015个，特有的HEGs分别为219和491个。研究结果不仅可提供中蜂幼虫肠道发育过程中的基因表达谱数据，也为分子水平上深入研究中蜂幼虫肠道发育提供重要信息。

RNA-seq；中华蜜蜂；幼虫肠道；高表达基因

东方蜜蜂(Apiscerana)主要集中分布于亚洲的中国、日本、泰国、菲律宾和伊朗等国家，以及俄罗斯远东地区[1]。中华蜜蜂(Apisceranacerana，简称中蜂)特指分布在我国的东方蜜蜂，其距今已有3 000多年的驯养历史[2]。中蜂是我国重要的本土蜜蜂资源，也是养蜂生产中的常用蜂种之一，在长期的进化过程中已高度适应我国的气候、蜜源等生态环境，具有诸多优点，例如善于利用零星蜜源、饲料消耗量少、对蜂螨和球囊菌等具有较强的抗性等[3]。

昆虫肠道具有消化、吸收和免疫防御等重要功能。有关蜜蜂肠道的研究主要集中在肠道菌群，如Gilliam[4]曾通过传统培养方法对西方蜜蜂(Apismellifera)肠道内的微生物进行分离鉴定，发现其肠道内有复杂而多样的共生微生物群；Jeyaprakash等[5]利用基于16S rRNA序列分析的分子生物学技术研究了非洲两种西方蜜蜂的肠道菌，发现其肠道内只有少数几类特定的细菌；Li等[6]在东方蜜蜂中检测到5种类型的细菌，其中Pasteurellaceae(Gamma-1)、Lactobacillus (Firm)、Neisseriaceae(Beta)和Bifidobacterium(Bifido)共占东方蜜蜂细菌总量的90%。然而，蜜蜂及其幼虫肠道发育的相关研究非常滞后，肠道过程中的基因表达规律尚不明确。前期研究中，我们已对中蜂幼虫肠道参考转录组进行了组装和注释(未发表数据)，可为其转录组学研究提供可靠的参考信息。本研究利用RNA-seq技术对中蜂5日龄和6日龄幼虫肠道进行深度测序，通过GO(gene ontology)和KEGG(Kyoto encyclopedia of genes and genomes)代谢通路富集分析对幼虫肠道高表达基因(HEGs)进行转录组学研究，研究结果可在转录组水平提供中蜂幼虫肠道发育过程中的基因表达谱数据，也为在分子水平深入研究肠道发育提供重要信息。

1 材料与方法

1.1 生物材料

本研究中使用的中蜂幼虫取自福建农林大学蜂学学院教学蜂场。

1.2 主要试剂及仪器

DNaseⅠ和Oligotex mRNA Kits Midi试剂盒购自德国Qiagen公司，Dynal M280磁珠购自Invitrogen公司，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 中蜂幼虫的人工饲养

参照王倩等[7]的方法配制中蜂幼虫饲料，其中将D-果糖和D-葡萄糖换为新鲜蜂蜜。预实验中中蜂幼虫7日龄成活率达到70%以上。选择20群群势较强的蜂群(外观无白垩病症状)，对上述蜂群进行PCR检测，结果阴性的蜂群作为本研究的实验蜂群。将2日龄幼虫移至无菌的24孔细胞培养板(每孔对应1只幼虫，孔内加有35 ℃预温的幼虫饲料)，在35 ℃，相对湿度70%条件下培养至6日龄。每隔24 h更换饲料。

1.4 测序样品准备

分别剖取中蜂5日龄和6日龄的幼虫肠道，为尽量减少肠道RNA的降解，将从解剖取样到样品放入液氮速冻的时间控制在30 s以内，每剖取一组样品，液氮速冻并迅速放入-80 ℃超低温冰箱保存备用。

1.5 cDNA合成与RNA-seq

利用RNAiso Reagent试剂盒抽提上述2个幼虫肠道的总RNA，然后用RNase-free DNaseⅠ去除基因组DNA残留。RNA的质量通过琼脂糖凝胶电泳和NanoDrop ND-1000(NanoDrop， Wilmington， DE， USA)进行检测。cDNA文库构建参照张曌楠等[8]的建库方法。

上述幼虫肠道样品委托广州基迪奥生物科技有限公司进行双端测序，测序平台为Illumina Hiseq 2500。测序数据已上传NCBI SRA数据库，SRA号为：SRA456721。

1.6 基因表达丰度计算与高表达基因分析

Illumina PE125测序下机的原始读段(raw reads)，经去除含接头(adapter)的、含N比例大于10%的和低质量的reads得到有效读段(clean reads)，进而利用Tophat 软件将有效读段(clean reads)映射(mapping)至中蜂幼虫肠道参考转录组(未发表数据)，表达量的计算使用FPKM(fragments per kilobase of transcript per million mapped reads)法，其计算公式为：FPKM=(106C×103)/NL。

选取各肠道样品FPKM值大于15的高HEGs，利用在线分析工具OmicsShare(http://www.omicshare.com/tools/index.php/Home/Index/index.html)进行GO分类及KEGG pathway富集分析。

2 结果与分析

2.1 RNA-seq测序数据概述

对2个中蜂幼虫肠道样品进行深度测序，Ac5和Ac6分别得到的reads数为29893270与27182602，去除低质量reads后的clean reads数分别为29290630与26636038，Q20分别为98.45%和98.38%(表1)，说明本研究中的RNA-seq数据质量良好，可用于进一步分析。Ac5与Ac6的clean reads mapping到中蜂幼虫肠道参考转录组的比例分别为88.98%和89.87%。Ac5和Ac6样品表达基因数分别为25 987和28 004，分别占参考基因总数(43 556 unigenes)的59.66%和64.29%。

2.2 Ac5与Ac6中高表达基因的GO分类

一般认为FPKM值在0.10～3.75的是低表达基因，在3.75～15.00为中度表达基因，大于15.00的为高度表达基因(HEGs)[9]，根据FPKM值选取所有HEGs进行GO富集分析，结果显示，Ac5中的HEGs富集在44个GO term，其中基因富集数最多的为代谢进程(753 unigenes)，其次是催化活性(679 unigenes)、结合(657 unigenes)、细胞进程(591 unigenes)和单组织进程(528 unigenes)；Ac6中的HEGs富集在44个GO term，其中基因富集最多的为代谢进程(792 unigenes)，其次是结合(760 unigenes)、催化活性(744 unigenes)、细胞进程(656 unigenes)和单组织进程(582 unigenes)(图1)。

2.3 Ac5与Ac6中高表达基因的KEGG pathway富集分析

KEGG pathway 富集分析结果显示，Ac5的HEGs富集在129个代谢通路(pathway)，其中基因富集数最多的是核糖体(103 unigenes)，其次是内质网蛋白质加工(68 unigenes)和碳代谢(67 unigenes)；Ac6的HEGs富集在129个代谢通路，其中基因富集数最多的是核糖体(103 unigenes)，其次是氧化磷酸化(88 unigenes)和内质网蛋白质加工(68 unigenes)(图2)。上述结果说明，蛋白等基础物质代谢和能量代谢在中蜂幼虫肠道发育过程中非常活跃，肠道发育对物质和能量的需求旺盛。

表1RNA-seq数据信息统计

Table1Overview of RNA-seq data in this study

样本Sample原始数据Rawreads有效读段数Cleanreads(%)接头Adapter(%)低质量读段Lowquality(%)99%碱基正确率(Q20)(%)Ac52989327029290630(97.98%)271218(0.91%)328722(1.1%)3604436630(98.45%)Ac62718260226636038(97.99%)224524(0.83%)319664(1.18%)3275561178(98.38%)

A， Ac5样品HEGs的GO分类；B， Ac6样品HEGs的GO分类；1， 行为；2， 生物黏附；3， 生物调控；4， 细胞成分组织或生物合成；5， 细胞进程；6， 解毒作用 ；7， 发育进程；8， 生长；9， 免疫系统进程；10， 定位；11， 运动；12， 代谢进程；13， 多组织进程；14， 多细胞组织进程；15， 生殖；16， 生殖进程；17， 应激；18， 信号；19， 单组织进程；20， 细胞；21， 细胞连接；22， 细胞零件；23， 细胞外基质；24， 胞外区；25， 胞外区零件；26， 大分子复合物；27， 细胞膜；28，细胞膜零件；29， 细胞膜内腔；30， 细胞器；31， 细胞器零件；32， 病毒；33， 病毒零件；34， 抗氧化活性；35， 结合；36， 催化活性；37， 电子转运活性；38， 分子功能调节因子；39， 分子转换器活性；40， 核酸结合转录因子活性；41， 信号转导因子；42， 结构分子活性；43， 转录因子活性，蛋白结合；44， 转运子活性A， GO enrichment analysis of HEGs in Ac5; B， GO enrichment analysis of HEGs in Ac6; 1， Behavior; 2， Biological adhesion; 3， Biological regulation; 4， Cellular component organization or biogenesis; 5， Cellular process; 6， Detoxification; 7， Developmental process; 8， Growth; 9， Immune system process; 10， Localization; 11， Locomotion; 12， Metabolic process; 13， Multi-organism process; 14， Multicellular organismal process; 15， Reproduction; 16， Reproductive process; 17， Response to stimulus; 18， Signaling; 19， Single-organism process; 20， Cell; 21， Cell junction; 22， Cell part; 23， Extracellular matrix; 24， Extracellular region; 25， Extracellular region part; 26， Macromolecular complex; 27， Membrane; 28， Membrane part; 29， Membrane-enclosed lumen; 30， Organelle; 31， Organelle part; 32， Virion; 33， Virion part; 34， Antioxidant activity; 35， Binding; 36， Catalytic activity; 37， Electron carrier activity; 38， Molecular function regulator; 39， Molecular transducer activity; 40， Nucleic acid binding transcription factor activity; 41， Signal transducer activity; 42， Structural molecule activity; 43， Transcription factor activity， protein blinding; 44， Transporter activity图1 高表达基因的GO富集分析Fig.1 GO enrichment analysis of HEGs

2.4 Ac5与Ac6中高表达基因的Venn分析

对Ac5和Ac6中蜂幼虫肠道样品的Venn分析结果显示(图3)，二者的共有HEGs为2 015个，特有HEGs分别为219和491个，说明绝大多数HEGs在中蜂5和6日龄幼虫肠道均有高表达，推测这些共有HEGs对中蜂幼虫肠道的正常发育十分重要，而特有HEGs在幼虫肠道的不同发育时期发挥特殊作用。

图2 高表达基因的KEGG pathway富集分析Fig.2 KEGG pathway enrichment analysis of HEGs

图3 Ac5和Ac6中高表达基因的Venn分析Fig.3 Venn analysis of HEGs in Ac5 and Ac6

3 讨论

中蜂是我国特有的本土蜂种，也是我国生态系统不可或缺的组成部分，具有十分重要的生态价值[10]。此外，中蜂作为养蜂生产中的常用蜂种，具有较高的经济价值。目前，有关中蜂及其幼虫肠道发育的分子研究鲜有报道。近年来，高通量测序技术发展迅猛，因其通量大和准确度高而广泛应用于动物、植物和微生物的研究领域，在蜜蜂的相关研究中也取得一些重要进展[11-13]。本研究利用RNA-seq技术对中蜂5日龄和6日龄肠道进行深度测序，发现二者的表达基因数分别占参考基因总数的59.66%和64.29%，而前期研究显示4日龄中蜂幼虫肠道表达基因数为36 602，占84.03%(未发表数据)，说明中蜂幼虫肠道在不同发育时期所表达的基因数量差别较大，发育早期表达的基因数明显多于发育后期。本研究还发现，Ac5和Ac6中均有大量HEGs富集在核糖体和内质网蛋白加工，说明中蜂幼虫肠道的发育过程对蛋白质的需求大增，因而蛋白质代谢旺盛。还发现大量HEGs富集在氧化磷酸化，说明中蜂幼虫肠道发育过程中伴随物质代谢的增强，能量代谢也相应提高。

角质层和围食膜是昆虫免疫防御的第一道防线[1]，当病原微生物突破这道防线后，将遭遇昆虫体内细胞免疫与体液免疫的抵抗，包括内吞作用、黑化作用、吞噬作用、蛋白的酶促水解及抗菌肽等[14]。值得注意的是，对于Ac5和Ac6，相当多的HEGs富集在细胞免疫通路，如内吞作用、溶酶体和吞噬体，且随着发育时间的延长，富集在上述细胞免疫通路的HEGs数量呈增加趋势，表明在中蜂幼虫肠道发育后期，细胞免疫通路较为活跃，可能赋予中蜂幼虫较强的免疫防御能力。

西方蜜蜂的基因组早在2016年就已公布[15]，为其分子研究奠定了重要基础。前期研究中，我们已对健康及球囊菌胁迫的意大利蜜蜂(Apismelliferaligustica，简称意蜂)幼虫肠道进行了HEGs分析(待发表)。较之意蜂，中蜂的分子研究进展缓慢，其中参考基因组的缺失是最大的制约因素。直到2015年，Park等[16]通过二代测序技术对雄蜂进行测序，组装得到中蜂参考基因组，然而作者当时并未同时公布基因位置及注释信息，导致该参考基因组无法被有效利用。肠道对于蜜蜂及其幼虫而言具有十分重要的作用，它是重要的消化吸收场所和免疫防御器官。我们在前期研究中denovo组装并注释了中蜂幼虫肠道的参考转录组，可为深入研究中蜂幼虫肠道的发育机理、幼虫响应病原(如球囊菌、东方蜜蜂微孢子虫和囊状幼虫病毒)胁迫的应答机制提供重要的参考信息。

此前，中蜂幼虫肠道发育的组学研究未见报道，有关中蜂幼虫肠道发育机理的研究相对滞后。本研究在前期研究基础上，对中蜂5日龄和6日龄幼虫肠道进行HEGs分析，研究结果不仅可提供中蜂幼虫肠道发育过程中的基因表达谱数据，也为在分子水平深入研究中蜂幼虫肠道发育提供重要信息。本研究是根据基因表达量FPKM值筛选得到HEGs，进而预测其功能，若要在分子水平更加全面、深入地研究中蜂幼虫肠道的发育，则需要对不同日龄的幼虫肠道进行差异表达基因分析，进一步建立肠道发育过程中的分子调控网络、挖掘关键调控基因。这将是我们下一步的工作重点。

[1] 龚一飞. 蜜蜂分类与进化[M]. 福建: 福建科学技术出版社， 2000.

[2] 杨冠煌. 中华蜜蜂[M]. 北京: 中国农业科技出版社， 2001.

[3] 张复兴. 现代养蜂生产[M]. 北京: 中国农业大学出版社， 1998.

[4] GILLIAM M. Identification and roles of non-pathogenic microflora associated with honey bees[J].FEMSMicrobiologyLetters， 1997， 155(1): 1-10.

[5] JEYAPRAKASH A， HOY M A， ALLSOPP MH. Bacterial diversity in worker adults ofApismelliferacapensis andApismelliferascutellata (Insecta: Hymenoptera) assessed using 16S rRNA sequences[J].JournalofInvertebratePathology， 2003， 84(2): 96-103.

[6] LI JL， QIN HR， WU J， et al.The prevalence of parasites and pathogens in Asian honeybeesApisceranain China[J].PLoSOne， 2012， 7(11): e47955.

[7] 王倩， 孙亮先， 肖培新， 等. 室内人工培育中华蜜蜂幼虫技术研究[J]. 山东农业科学， 2009(11): 113-116. WANG Q， SUN L X， XIAO P X， et al. Study on technology for indoor artificial feeding ofApisceranaceranalarvae[J].ShandongAgriculturalSciences， 2009(11): 113-116. (in Chinese with English abstract)

[8] 张曌楠， 熊翠玲， 徐细建， 等. 蜜蜂球囊菌的参考转录组denovo组装及SSR分子标记开发[J]. 昆虫学报， 2017，60(1): 34-44. ZHANG Z N， XIONG C L， XU X J， et al. De novo assembly of a reference transcriptome and development of SSR markers forAscosphaeraapis[J].ActaEntomologicaSinica， 2017， 60(1): 34-44. (in Chinese with English abstract)

[9] MORTAZAVI A， WILLIAMS BA， MCCUE K， et al. Mapping and quantifying mammalian transcriptomes by RNA-Seq[J].NatureMethods， 2008， 5(7): 621-628.

[10] 丁桂玲. 中华蜜蜂(Apiscerana)群体多样性的研究[D]. 北京: 中国农业科学院， 2006. DING G L. Study on the population diversity ofApiscerana[D]. Beijing: Chinese Academy of Agricultural Sciences， 2006. (in Chinese with English abstract)

[11] CORNMAN R S， LOPEZ D， EVANS J D. Transcriptional response of honey bee larvae infected with the bacterial pathogenPaenibacilluslarvae[J].PLoSOne， 2013， 8(6): e65424.

[12] JULIE A， BARBARA MISME-AUCOUTURIER， BERNARD V， et al. Transcriptome analyses of the honeybee response toNosemaceranaeand insecticides[J].PloSOne， 2014， 9(3): e91686.

[13] CRISTINO A S， BARCHUUK A R， FREITAS F C， et al. Neuroligin-associated microRNA-932 targets actin and regulates memory in the honeybee[J].NatureCommunications， 2014， 5: 5529.

[15] HONEY BEE GENOME SEQUENCING CONSORTIUM. Insights into social insects form the genome of the honeybeeApismellifera[J].Nature， 2006， 443(7114): 931-949.

[16] PARK D， JUNE J W， CHOI B S， et al. Uncovering the novel characteristics of Asian honey bee，Apiscerana， by whole genome sequencing[J].BMCGenomics， 2015， 16: 1.

(责任编辑张 韵)

AnalysisofhighlyexpressedgenesoflarvalgutofApisceranacerana

XIE Yanling， WANG Hongquan， ZHANG Lu， HOU Zhixian， DAO Chen， JIANG Liangliang， XIONG Cuiling， ZHENG Yanzhen， XU Xijian， HUANG Zhijian， GUO Rui*， CHEN Dafu

(CollegeofBeeScience，FujianAgricultureandForestryUniversity，Fuzhou350002，China)

Apisceranaceranais not only a specific honeybee species resource in China， but also a frequently-used honeybee species in apiculture. In the present study， the 5-(Ac5) and 6-day-old (Ac6) larval guts ofA.c.ceranawere sequenced using RNA-seq. After filtration， 29290630 and 26636038 clean reads with a Q20 of 98.45% and 98.38% were obtained in Ac5 and Ac6. GO enrichment analysis results showed that the highly expressed genes (HEGs) in Ac5 were enriched in 44 GO terms， among them metabolic process (753 unigenes)， catalytic activity (679 unigenes) and binding (657 unigenes) were the largest groups; the HEGs in Ac6 were enriched in 44 GO terms， and the mostly enriched one was metabolic process (792 unigenes)， followed by binding (760 unigenes) and catalytic activity (744 unigenes). KEGG enrichment analysis results suggested that the HEGs of Ac5 were enriched in 129 pathways， among them ribosome (103 unigenes) was the largest group， followed by protein processing in endoplasmic reticulum (68 unigenes) and carbon metabolism (67 unigenes); the HEGs of Ac6 were enriched in 129 pathways， the mostly enriched ones were ribosome (103 unigenes)， oxidative phosphorylation (88 unigenes) as well as protein processing in endoplasmic reticulum (68 unigenes). Our findings can not only offer the gene expression profiles during the developmental process of the larval gut ofA.c.ceranaat transcriptome level， but also provide key information for further investigation of the larval gut’s development.

RNA-seq;Apisceranacerana; larval gut; highly expressed genes

S89

：A

：1004-1524(2017)09-1575-06

解彦玲，王鸿权，张璐，等. 中华蜜蜂幼虫肠道的高表达基因分析[J].浙江农业学报，2017，29(9)： 1575-1580.

10.3969/j.issn.1004-1524.2017.09.22

2017-03-20

现代农业产业技术体系建设专项基金(CARS-45-KXJ7)；福建农林大学科技发展基金(KF2015123)；福建省大学生创新创业项目(201610389053)

解彦玲(1996—)，女，浙江湖州人，主要从事蜂学研究。E-mail: 1103430452@qq.com

*通信作者，郭睿，E-mail: fafu_ruiguo@126.com