转录组分析揭示东方蜜蜂微孢子虫侵染意大利蜜蜂的分子机制

耿四海, 周丁丁, 范小雪, 蒋海宾, 祝智威, 王 杰, 范元婵,王心蕊, 熊翠玲, 郑燕珍, 付中民, 陈大福, 郭 睿

(福建农林大学动物科学学院(蜂学学院), 福州 350002)

微孢子虫是一种细胞内寄生的单细胞真菌，广泛寄生脊椎动物和无脊椎动物(Cuomoetal., 2012)。自1857年首次在患病家蚕Bombyxmori中发现家蚕微孢子虫Nosemabombycis(Nageli, 1857)以来，至今已经有1 400多种微孢子虫被鉴定出来(吴杰等, 2017)。微孢子虫的基因组紧凑，缺少典型的线粒体，取而代之是一种微粒体，因而高度依赖宿主细胞提供能量(Burrietal., 2006)。

东方蜜蜂微孢子虫Nosemaceranae是专性寄生蜜蜂中肠上皮细胞的单细胞真菌病原，可引起蜜蜂微孢子虫病(Forsgren and Fries, 2010)。Freis等(1996)首次在北京地区的东方蜜蜂Apiscerana蜂群中分离和鉴定出东方蜜蜂微孢子虫，之后该孢子虫逐渐扩散到欧洲(Higesetal., 2006)和中国台湾地区(Huangetal., 2007)的西方蜜蜂Apismellifera蜂群，现已蔓延至世界各养蜂国家和地区(Shumkovaetal., 2018)。前人的研究结果表明东方蜜蜂微孢子虫感染可导致成年蜜蜂肠道破坏、消化系统紊乱,脂质代谢加快和脂质储存减少,免疫基因和细胞凋亡被抑制,能量胁迫,采集日龄提前及采集能力下降等；东方蜜蜂微孢子虫感染严重可造成蜂群的群势降低、采蜜量减少以及越冬失败(Higesetal., 2007; Antúnezetal., 2009; Mayack and Naug, 2009; Botíasetal., 2013; Chenetal., 2013; Mikeetal., 2013; Kurzeetal., 2015)。此外，有研究表明东方蜜蜂微孢子虫与蜂群崩溃失调症之间存在密切的关联(Higesetal., 2009; Simeunovicetal., 2014)。东方蜜蜂微孢子虫在体外以休眠态的孢子存在，孢子壁由几丁质和孢壁蛋白组成，能够抵御外界的不良环境，孢子被蜜蜂摄取之后，在蜜蜂中肠的特殊理化环境的刺激下发芽，内部高度压缩盘旋的极丝弹射而出，刺穿中肠上皮细胞并将有侵染性的胞原质注入，随后进入繁殖期，利用宿主的物质与能量合成系统进行增殖，孢子数量逐渐增多，最终导致宿主细胞裂解，释放出的成熟孢子继续感染临近的上皮细胞，或者随粪便排出体外，感染新的宿主(Simeunovicetal., 2014)。

随着高通量测序技术和生物信息学算法的发展，RNA-Seq技术不仅广泛应用于人类(Violletetal., 2017)、动物(李鹏飞等, 2018)和植物(Munetal., 2017)的研究，在微生物的相关研究中也得到了较好的应用，包括绿僵菌Metarhiziumanisopliae(Wangetal., 2014)、里氏木霉Trichodermareesei(Riesetal., 2013)、家蚕微孢子虫(Liuetal., 2016a)和蜜蜂球囊菌Ascosphaeraapis(郭睿等, 2017)等。Huang等(2016)对感染东方蜜蜂微孢子虫后1-6 d的西方蜜蜂中肠进行测序，从混合数据中筛选出东方蜜蜂微孢子虫的数据，进而通过统计东方蜜蜂微孢子虫的基因表达量，发现1 122个基因可分为4种表达模式；随后，该团队对东方蜜蜂微孢子虫的Dicer进行RNAi，然后通过转录组测序和分析发现Dicer调控东方蜜蜂微孢子虫的基因表达(Huangetal., 2018)。笔者团队前期结合链特异性建库方法和二代测序技术对东方蜜蜂微孢子虫孢子进行了深度测序，在全基因组水平预测、鉴定和分析了东方蜜蜂微孢子虫的长链非编码RNA(long non-coding RNA, lncRNA)和环状RNA(circular RNA, circRNA)(Guoetal., 2018a, 2018b)。但对于东方蜜蜂微孢子虫的组学研究较之蜜蜂仍然较为滞后，东方蜜蜂微孢子虫的侵染和增殖机制还未阐明。

笔者团队前期已对东方蜜蜂微孢子虫感染7和10 d的意大利蜜蜂Apismelliferaligustica工蜂中肠进行转录组深度测序，通过连续比对西方蜜蜂和东方蜜蜂微孢子虫的基因组过滤得到东方蜜蜂微孢子虫的纯净数据，结合东方蜜蜂微孢子虫孢子的转录组数据(Guoetal., 2018a)开展了东方蜜蜂微孢子虫的高表达基因(highly expressed gene, HEG)分析，解析了东方蜜蜂微孢子虫侵染意大利蜜蜂过程的HEG表达谱及潜在作用(熊翠玲等, 2019)。本研究在前期研究基础上对东方蜜蜂微孢子虫进行深入的转录组分析，进一步通过分析病原的毒力因子和侵染相关因子，揭示东方蜜蜂微孢子虫侵染意大利蜜蜂的分子机制。研究结果可为探明东方蜜蜂微孢子虫的侵染机制提供基础，为蜜蜂微孢子虫病的治疗提供潜在靶点。

1 材料与方法

1.1 转录组数据来源

前期研究中，笔者团队利用基于链特异性建库的lncRNA-Seq技术对东方蜜蜂微孢子虫的纯化孢子进行了深度测序(NcCK-1, NcCK-2和NcCK-3为3个生物学重复)，获得了高质量的mRNA和lncRNA组学数据(Guoetal., 2018a)，其中的mRNA组学数据可作为本研究中的对照组数据；此外，还利用二代测序技术对东方蜜蜂微孢子虫感染7 d(NcT1-1, NcT1-2和NcT1-3为3个生物学重复)和10 d(NcT2-1, NcT2-2和NcT2-3为3个生物学重复)的意大利蜜蜂工蜂中肠进行了全转录组测序，获得了高质量的测序数据(熊翠玲等, 2019)，其中的mRNA组学数据包含宿主来源的数据和病原来源的数据，因此通过连续比对核糖体数据库、西方蜜蜂基因组和东方蜜蜂微孢子虫基因组对上述混合数据进行过滤，筛滤得到的病原的mRNA组学数据(即侵染意大利蜜蜂工蜂的东方蜜蜂微孢子虫的mRNA组学数据)作为本研究中的处理组数据。

上述转录组测序的原始数据已上传美国国家生物技术信息中心(NCBI)SRA数据库，Bioproject号分别为PRJNA395264(NcCK)和PRJNA406998(NcT1和NcT2)。

1.2 差异表达基因(DEG)分析

采用FPKM算法计算和归一化基因表达量。利用edgeR软件筛选NcCKvsNcT1, NcCKvsNcT2和NcT1vsNcT2比较组的DEG，标准为P≤0.05且|log2(Fold change)|≥1。利用OmicShare在线平台(www.omicshare.com)对各比较组中上调基因和下调基因进行Venn分析，采用默认参数。参照郭睿等(2017)的方法，利用BLASTall软件将DEG比对到GO和KEGG数据库，进行功能分类和代谢通路注释，采用默认参数。

1.3 毒力因子和侵染相关因子分析

东方蜜蜂微孢子虫是专性寄生成年蜜蜂中肠上皮细胞的真菌病原，其增殖所需的物质和能量完全依赖宿主细胞提供(Burrietal., 2006)。根据DEG在Nr数据库的注释信息，同时参考前人关于东方蜜蜂微孢子虫的基因组学、转录组学和分子生物学的研究报道(Cornmanetal., 2009; Paldietal., 2010; Huangetal., 2016, 2018; Pelinetal., 2016; Rodríguez-Garcíaetal., 2018)以及其他微孢子虫的相关研究报道(Corradi and Slamovits, 2011; Katinkaetal., 2011; Nakjangetal., 2013; Simeunovicetal., 2014)，筛选出几丁质酶、孢壁蛋白、蓖麻毒素B凝集素、极管蛋白等毒力因子相关DEG，以及己糖激酶、丙酮酸激酶、磷酸果糖激酶、ATP/ADP移位酶、ABC转运蛋白等侵染相关因子相关DEG，并统计上述相关DEG在各组样品中的表达量(FPKM)及在各比较组中的log2(Fold change)。

1.4 DEG的RT-qPCR验证

从NcCKvsNcT1, NcCKvsNcT2和NcT1vsNcT2比较组中分别随机选取7个DEG，其中NcCKvsNcT1中包含6个上调基因和1个下调基因，NcCKvsNcT2中包含6个上调基因和1个下调基因，NcT1vsNcT2中包含5个上调基因和2个下调基因，以actin(GenBank登录号: 36320842)作为内参基因，利用Primer Premier 5设计引物(表1)，委托福州擎科生物有限公司合成。利用RNA提取试剂盒(TaKaRa公司，大连)分别提取NcCK以及NcT1和NcT2中肠总RNA(熊翠玲等, 2019)，然后利用cDNA第1链合成试剂盒(Vazyme公司，南京)反转录成cDNA，作为模板进行qPCR，反应在QuantStudio 3(ThermoFisher公司，美国)上进行。qPCR反应按照SYBR Green Dye试剂盒(Vazyme公司，南京)操作说明书进行，本实验进行3次生物学重复，每个生物学重复包含6头意大利蜜蜂工蜂中肠。反应体系(20 μL)：正反向引物(2.5 μmol/L)各1 μL, cDNA模板1 μL, SYBR Green Dye 10 μL, DEPC水7 μL。反应程序: 95℃预变性1 min; 95℃变性15 s, 60℃延伸30 s，共40个循环。通过2-△△Ct法对上述基因的相对表达量进行计算。利用Graph Prism 7软件进行数据分析和绘图。

2 结果

2.1 东方蜜蜂微孢子虫转录组的DEG分析

差异分析结果显示，NcCKvsNcT1, NcCKvsNcT2和NcT1vsNcT2比较组分别得到1 397, 1 497和52个DEG，其中上调和下调基因分别为497和900个, 517和977个，38和13个。Venn分析结果显示，有10个基因在上述3个比较组中共同上调表达(图1: A)，它们涉及碳水化合物代谢、蓖麻毒素B凝集素、细胞表面糖基化蛋白、信号肽酶和假定蛋白(表2)；仅有1个下调表达的假定蛋白编码基因为各比较组所共有(图1: B; 表2)。

表1 本研究使用的引物Table 1 Primers used in this study

图1 各比较组中东方蜜蜂微孢子虫转录组差异表达基因(DEGs)Venn分析Fig. 1 Venn analysis of differentially expressed genes (DEGs) in Nosema ceranae transcriptome between various comparison groupsA: 上调基因Up-regulated genes; B: 下调基因Down-regulated genes. NcCK: 微孢子虫纯化的孢子Purified spores of N.ceranae; NcT1, NcT2: 分别为感染意大利蜜蜂工蜂7和10 d中肠内的东方蜜蜂微孢子虫N. ceranae in the midgut of Apis mellifera ligustica workers at 7 and 10 d post infection, respectively. 下同The same below.

表2 各比较组中东方蜜蜂微孢子虫转录组共同上调与下调基因信息统计Table 2 Summary of commonly up- and down-regulated genes in Nosema ceranae transcriptome between various comparison groups

2.2 东方蜜蜂微孢子虫转录组中DEG的GO分类及KEGG 通路富集分析

GO分类结果表明，NcCKvsNcT1中DEG涉及生物学进程、细胞组分和分子功能相关的23个功能条目，其中代谢进程(221)、催化活性(195)、结合(193)、细胞进程(193)、细胞(116)、细胞组件(116)、单组织进程(97)和细胞器(63)的条目富集基因数最多；NcCKvsNcT2中DEG也涉及23个功能条目，同样有较多的DEG富集在代谢进程(218)、催化活性(196)、细胞进程(195)、结合(193)、细胞(119)、细胞组件(119)、单组织进程(94)和细胞器(51)等；NcT1vsNcT2中DEG涉及生物学进程(4)、细胞组分(5)和分子功能(2)相关的11个功能条目，其中催化活性、细胞进程、代谢进程、单组织进程和结合的条目富集基因数最多，分别为11, 10, 10, 7和6个(图2)。

进一步对各比较组的DEG进行KEGG通路富集分析，结果显示NcCKvsNcT1中有355个DEG富集在80条通路，涉及代谢、遗传信息加工、环境信息加工、细胞进程和有机体系统等5个大类，其中富集基因数最多的是新陈代谢途径(93)、次级代谢产物合成(37)、核糖体(33)、抗生素生物合成(29)和真核生物核糖体生物合成(27)(表3)；NcCKvsNcT2中有356个DEG富集在79条通路，其中富集基因数最多的是新陈代谢途径、次级代谢产物合成、核糖体、抗生素合成和真核生物核糖体生物合成，分别为92, 39, 32, 31和26个(表3)。此外NcCKvsNcT1中有7个DEG富集在MAPK信号通路，其中有5个基因上调表达，GenBank登录号分别为9423510, 9424462, 9422987, 9423963和9422306；另有2个基因下调表达，GenBank登录号分别为9422598和9424007。 NcCKvsNcT2中同样也有7个DEG富集在MAPK信号通路，包括5个上调基因，GenBank登录号分别为9423510, 9424610, 9424462, 9422987和9422306；以及2个下调基因，GenBank登录号分别为9411598和9424007。在NcCK, NcT1和NcT2中，GenBank登录号为9424462的基因的FPKM值分别为966, 829和22，GenBank登录号为9422306的基因的FPKM值分别为209, 205和5。

图2 各比较组中东方蜜蜂微孢子虫转录组差异表达基因(DEGs)的GO分类Fig. 2 GO categories of differentially expressed genes (DEGs) in Nosema ceranae transcriptome between various comparison groups

2.3 东方蜜蜂微孢子虫的毒力因子分析

根据NcCKvsNcT1, NcCKvsNcT2和NcT1vsNcT2的DEG的Nr数据库注释情况对东方蜜蜂微孢子虫的毒力因子进行统计，分别发现5个孢壁蛋白基因、1个孢壁蛋白前体基因和1个孢壁和锚定盘复合蛋白基因，其中孢壁蛋白9基因(GenBank登录号: 9422782)在NcCKvsNcT1和NcCKvsNcT2中下调表达[log2(Fold change)值分别为-2.10和-2.06]，孢壁蛋白12基因(GenBank登录号: 9422437)也下调表达[log2(Fold change)值分别为-0.86和-0.67]；孢壁蛋白8基因(GenBank登录号: 9422311)在NcCKvsNcT1中表达量下调[log2(Fold change)值为-2.05]，但在NcCKvsNcT2中其表达量不变；在NcCKvsNcT1和NcCKvsNcT2中均上调表达的基因包括2个孢壁蛋白基因[GenBank登录号为9422653的基因的log2(Fold change)值分别为7.11和7.12; GenBank登录号为9423533的基因的log2(Fold change)值分别为6.32和6.93]、1个孢壁蛋白前体基因[GenBank登录号: 9422389, log2(Fold change)值分别为5.03和5.46]和1个孢壁和锚定盘复合蛋白基因[GenBank登录号: 9424403, log2(Fold change)值分别为5.54和4.61]。此外，另有1个内切几丁质酶基因[GenBank登录号: 9422749, log2(Fold change)值分别为2.94和3.60]、1个几丁质合酶1基因[GenBank登录号为9422288的基因的log2(Fold change)值分别为6.08和5.56]、3个极管蛋白基因[GenBank登录号为9423198的基因的log2(Fold change)值分别为5.62和6.34; GenBank登录号为9423201的基因的log2(Fold change)值分别为5.74和6.37; GenBank登录号为9422340的基因的log2(Fold change)值分别为3.94和4.39]和7个蓖麻毒素B凝集素基因[GenBank登录号为9422649的基因的log2(Fold change)值分别为4.76和5.98; GenBank登录号为9422719的基因的log2(Fold change)值分别为4.89和6.32; GenBank登录号为9422642的基因的log2(Fold change)值分别为4.77和4.98; GenBank登录号为9423906的基因的log2(Fold change)值分别为2.32和3.23; GenBank登录号为9422645的基因的log2(Fold change)值分别为3.43和4.74; GenBank登录号为9422644的基因的log2(Fold change)值分别为2.72和3.21; GenBank登录号为9422647的基因的log2(Fold change)值分别为5.62和3.38]在NcCKvsNcT1和NcCKvsNcT2中均上调表达(表4)。

表3 NcCK vs NcT1和NcCK vs NcT2中差异表达基因(DEGs)富集数前18的KEGG通路Table 3 Top 18 KEGG pathways enriched by differentially expressed genes (DEGs) in NcCK vs NcT1 and NcCK vs NcT2

2.4 东方蜜蜂微孢子虫的侵染相关因子分析

进一步对东方蜜蜂微孢子虫的其他侵染相关因子进行统计，发现在NcCKvsNcT1和NcCKvsNcT2中均上调表达的糖酵解关键酶基因包括己糖激酶基因[GenBank登录号: 9423428, log2(Fold change)值分别为3.38和2.79]、丙酮酸激酶基因[GenBank登录号: 9424331, log2(Fold change)值分别为1.81和2.69]和6-磷酸果糖激酶基因[GenBank登录号: 9424162, log2(Fold change)值分别为1.91和2.79]等(表5)。在NcCKvsNcT1和NcCKvsNcT2中均上调表达的有3个ATP/ADP移位酶的编码基因[GenBank登录号为9424586的基因的log2(Fold change) 值分别为2.31和2.34; GenBank登录号为9422880的基因的log2(Fold change)值分别为3.60和3.44; GenBank登录号为9424363的基因的log2(Fold change)值分别为0.72和0.45]。此外，有1个ATP/ADP移位酶的编码基因[GenBank登录号:9422375, log2(Fold change)值分别为-3.65和-3.90]在NcCKvsNcT1和NcCKvsNcT2中均下调表达。另发现7个DEG涉及ABC转运蛋白，其中有2个[GenBank登录号为9423041的基因的log2(Fold change)值分别为1.83和1.85; GenBank登录号为9422518的基因的log2(Fold change)值分别为0.93和0.64]在NcCKvsNcT1和NcCKvsNcT2中均上调表达；还有1个DEG(GenBank登录号: 9423588)在NcCKvsNcT1中下调表达[log2(Fold change)值为-0.79]，而在NcCKvsNcT2中表达量上调[log2(Fold change)值为0.67]；其余的4个ABC转运蛋白基因[GenBank登录号为9423050的基因的log2(Fold change)值分别为-3.65和-3.34; GenBank登录号为9423588的基因的log2(Fold change)值分别为-1.72和-0.13; GenBank登录号为9422298的基因的log2(Fold change)值分别为-1.91和-1.56; GenBank登录号为9422297的基因的log2(Fold change)值分别为-0.19和-0.35]在NcCKvsNcT1和NcCKvsNcT2中表现为表达量下调(表5)。

表4 东方蜜蜂微孢子虫的毒力因子统计Table 4 Summary of virulence factors of Nosema ceranae

表5 东方蜜蜂微孢子虫的侵染相关因子统计Table 5 Summary of infection-associated factors of Nosema ceranae

2.5 东方蜜蜂微孢子虫转录组DEG的RT-qPCR验证

为验证RNA-Seq结果，在NcCKvsNcT1, NcCKvsNcT2和NcT1vsNcT2比较组中分别随机挑取7个DEG进行RT-qPCR验证，实验结果与转录组数据的差异表达趋势一致(图3)，证明了测序数据和基因差异表达情况的可靠性。

图3 各比较组中东方蜜蜂微孢子虫转录组差异表达基因(DEGs)的RT-qPCR验证Fig. 3 RT-qPCR verification of differentially expressed genes (DEGs) in Nosema ceranae transcriptome between various comparison groupsA: NcCK vs NcT1; B: NcCK vs NcT2; C: NcT1 vs NcT2. 误差棒表示各DEG的RT-qPCR结果的方差。Error bars represent the variance of RT-qPCR results of each DEG.

3 讨论

前人研究报道东方蜜蜂微孢子虫在感染西方蜜蜂后10-20 d内都处于一个持续增殖的过程，在此过程病原的孢子数持续增多(Huang and Solter, 2013)。前期研究中，我们发现东方蜜蜂微孢子虫感染的意大利蜜蜂工蜂的死亡率在感染后1-10 d范围内不断上升，且与未感染工蜂的死亡率差异极显著(Chenetal., 2019)，表明随着宿主中肠内东方蜜蜂微孢子虫的不断增殖，病原数量持续增多，对宿主的胁迫压力持续上升。在另一项研究中，我们发现东方蜜蜂微孢子虫的孢子数在4-13 d的范围内持续增加(未发表数据)。本研究是属于完整课题的一部分，完整课题的主要研究目的是探究意大利蜜蜂与东方蜜蜂微孢子虫的相互作用机制，一方面是解析宿主的胁迫应答和免疫应答机制，另一方面是解析病原的侵染机制。鉴于东方蜜蜂微孢子虫的增殖周期尚无定论，以及处于增殖过程的病原存在孢原质、裂殖体和成熟孢子等各种状态，综合考虑后选取东方蜜蜂微孢子虫感染后7和10 d作为本研究取样、测序和数据分析的两个时间点。目前仍缺乏不依赖于蜜蜂宿主的东方蜜蜂微孢子虫体外培养体系。本研究中使用的东方蜜蜂微孢子虫纯化孢子是从患微孢子虫病的意大利蜜蜂外勤蜂中肠制备而来，其中可能既含有经增殖而致死宿主的孢子，也含有未侵入宿主细胞的孢子。鉴于目前没有区分此二种类型孢子的技术手段，我们参照前人研究中采用的接种西方蜜蜂所采用的孢子浓度和接种方法(Antúnezetal., 2009; Huang and Solter, 2013; Huangetal., 2016, 2018; Rodríguez-Garcíaetal., 2018)，使每头工蜂摄入等量的孢子。对于微孢子虫，相对于代谢和生命活动活跃的增殖阶段(如裂殖体)，休眠态孢子中仅存在较低水平的代谢和生命活动(Williamsetal., 2005; Corradietal., 2008)。本研究中，用于接种工蜂的孢子(其中已完成侵染的孢子仍然为休眠态孢子)来源于同一批制备的东方蜜蜂微孢子虫纯化孢子，通过对照组和处理组数据比较分析得到的基因差异表达信息应由病原侵染所致。本研究中，在NcCKvsNcT1, NcCKvsNcT2和NcT1vsNcT2中分别鉴定出1 397, 1 497和52个显著DEG，表明东方蜜蜂微孢子虫在其侵染意大利蜜蜂工蜂的过程中伴随着活跃的基因差异表达，暗示这些DEG中蕴含着病原侵染相关的重要信息。

Cornman等(2009)组装及注释了东方蜜蜂微孢子虫的基因组，分析发现东方蜜蜂微孢子虫的供能主要依赖碳水化合物代谢、糖酵解、氧化磷酸化和磷酸戊糖途径。本研究发现，有10个基因在NcCKvsNcT1, NcCKvsNcT2和NcT1vsNcT2 3个比较组均上调表达，其中4个上调基因涉及碳水化合物代谢，包括己糖激酶基因(GenBank登录号: 9423428)、甘油-3-磷酸脱氢酶基因(GenBank登录号: 9423035)、磷酸甘油酸激酶基因(GenBank登录号: 9423279)和丙酮酸脱氢酶e1-β亚基基因(GenBank登录号: 9423243)(表2)，表明东方蜜蜂微孢子虫在侵染过程中通过上调上述4个基因的表达量，满足自身增殖的能量需要。蓖麻毒素是一种可抑制蛋白质活性和改变宿主细胞激素编码基因的表达水平的细胞毒素，具有结合并粘附细胞的蓖麻毒素B凝集素的结构(Spooner and Lord, 2015)。有研究表明在东方蜜蜂微孢子虫侵染家蚕Bombyxmori卵巢BmN细胞的过程中，蓖麻毒素B凝集素能够增强家蚕微孢子虫对宿主细胞的粘附能力(Liuetal., 2016b)。Huang等(2016)研究发现东方蜜蜂微孢子虫感染西方蜜蜂的次日就能检测到蓖麻毒素编码基因的表达，进一步分析发现蓖麻毒素含有信号肽，同时宿主的核糖体蛋白受到抑制。本研究中，3个比较组共有的上调基因包含3个蓖麻毒素B凝集素编码基因(GenBank登录号分别为9422649, 9422645和9422719)(表2)，暗示东方蜜蜂微孢子虫在侵染意大利蜜蜂工蜂的过程中通过增强蓖麻毒素B凝集素的分泌以加强对中肠上皮细胞的粘附力，同时抑制宿主细胞蛋白质的合成。细胞膜蛋白糖基化是细胞识别和应答外界环境的重要方式之一(刘啸尘等, 2018)。本研究中，细胞表面糖基化蛋白编码基因(GenBank登录号: 9423191)在NcCKvsNcT1, NcCKvsNcT2和NcT1vsNcT2 3个比较组中分别上调1.74, 3.06和1.32倍(表2)，推测东方蜜蜂微孢子虫在侵染过程中通过上调细胞表面糖基化蛋白基因的表达对宿主细胞内环境产生应答。此外，共同上调基因中包含1个假定蛋白编码基因(GenBank登录号: 9422567)，共同下调的1个基因也为假定蛋白编码基因(GenBank登录号: 9422077)(表2)，说明东方蜜蜂微孢子虫的基因组功能注释信息尚不完全，有待通过分子克隆、基因和蛋白功能研究进一步完善。到目前为止，包括东方蜜蜂微孢子虫在内的绝大多数微孢子虫的基因功能未知，最大的制约因素是缺乏转基因操作技术平台。Guo等(2016)通过碱液刺激家蚕微孢子虫孢子体外发芽，继而感染家蚕BmN细胞系，3 d后利用脂质体包裹piggyBac转座子载体和pIZT/V5-His非转座子载体转染家蚕微孢子虫感染的细胞，通过连续的荧光观察和分子生物学检测证实将2种载体携带的gfp基因成功导入家蚕微孢子虫基因组；此外，作者还通过饲喂家蚕微孢子虫孢子感染家蚕个体，继而注射脂质体包裹的上述两种载体，进而也通过连续的荧光观察和分子生物学检测证实家蚕微孢子虫基因组插入了gfp基因。目前，笔者团队正在参考该方法尝试建立东方蜜蜂微孢子虫的替代细胞连续感染系和转基因操作技术体系。

NcCKvsNcT1和NcCKvsNcT2比较组包含的DEG富集在代谢进程、细胞进程、单组织进程等生物学进程相关功能条目，细胞、细胞组件和细胞器等细胞组分相关功能条目，结合、催化活性和转运活性等分子功能相关功能条目，表明东方蜜蜂微孢子虫通过改变新陈代谢、细胞生命活动相关基因的表达水平促进自身增殖。NcT1vsNcT2比较组包含的DEG主要富集在催化活性、代谢进程和细胞进程等功能条目且多数基因上调表达，说明东方蜜蜂微孢子虫在侵染意大利蜜蜂工蜂的过程中物质和能量代谢活跃。进一步对DEG进行KEGG通路分析，发现NcCKvsNcT1和NcCKvsNcT2比较组包含的DEG大量富集在新陈代谢相关通路，包括新陈代谢途径(分别为93和92个)、次级代谢产物合成(分别为37和39个)、抗生素合成(分别为29和31个)、嘌呤代谢(分别为21和21个)、嘧啶代谢(分别为20和19个)、碳代谢(分别为16和18个)、甘油磷脂代谢(分别为11和13个)、糖酵解/糖异生(分别为11和12个)等物质代谢通路，以及氧化磷酸化(分别为10和10个)和甲烷代谢(分别为4和6个)等能量代谢通路(表3)。上述结果再次表明东方蜜蜂微孢子虫侵染意大利蜜蜂工蜂的过程伴随着复杂的物质和能量代谢活动，病原本身需要提升物质和能量代谢满足增殖需要，但同时也面临宿主的抑制。真菌通过MAPK、cAMP和双组份信号通路应答外界环境，调控其生长发育和毒力 (侯彬彬等, 2015)。蜜蜂球囊菌是另一种常见的蜜蜂真菌病原，特异性侵染蜜蜂幼虫而导致白垩病，笔者团队在前期研究中发现，对于侵染意大利蜜蜂幼虫的蜜蜂球囊菌，有48个DEG富集在MAPK信号通路，它们的表达水平随蜜蜂球囊菌胁迫时间的延长而逐渐增强(陈大福等, 2017)。本研究中，对于NcCKvsNcT1和NcCKvsNcT2比较组，分别有7和7个DEG富集在MAPK信号通路；笔者前期研究发现，对于NcCK, NcT1和NcT2，分别有5, 7和7个涉及MAPK信号通路的基因高量表达(熊翠玲等, 2019)；深入分析发现本研究中MAPK信号通路上有2个DEG与前期研究的高表达基因重叠，其中1个DEG(GenBank登录号: 9424462)在NcCK, NcT1和NcT2均高量表达(FPKM值分别为966, 829和22)，另外1个DEG(GenBank登录号: 9422306)在NcT1和NcT2高量表达(FPKM值分别为209和205)；上述结果表明MAPK信号通路及其富集DEG在东方蜜蜂微孢子虫的侵染过程中发挥重要作用，下一步拟通过设计相应的siRNA对GenBank登录号分别为9424462和9422306的基因进行RNAi，以探究其功能。

东方蜜蜂微孢子虫侵染蜜蜂时通过中空的极管将具有侵染性的孢原质注入宿主的中肠上皮细胞，因此极管蛋白在东方蜜蜂微孢子虫的侵染过程扮演着关键角色(Rodríguez-Garcíaetal., 2018)。前人研究报道干扰极管蛋白3编码基因的表达可显著降低东方蜜蜂微孢子虫孢子的数量(Rodríguez-Garcíaetal., 2018)。本研究中，在NcCK, NcT1和NcT2中极管蛋白1编码基因(GenBank登录号: 9423198)的表达量分别为171, 8 047和13 854；极管蛋白2编码基因(GenBank登录号: 9423201)的表达量分别为353, 17 919和27 591；极管蛋白编码基因(GenBank登录号: 9422340)的表达量分别为190, 2 923和3 994(表4)，表明东方蜜蜂微孢子虫在侵染增殖过程中比孢子状态上调表达，并且在侵染过程持续上调，印证了极管蛋白对东方蜜蜂微孢子虫侵染的重要性(Rodríguez-Garcíaetal., 2018)。孢壁蛋白保护东方蜜蜂微孢子虫在孢子状态抵制外界不良环境，同时在应答外界环境的变化传递信号和附着等作用(Southernetal., 2007)。孢壁蛋白还可以与极管蛋白互相作用，例如家蚕微孢子虫的孢壁蛋白5与极管蛋白2和极管蛋白3相互作用，孢壁蛋白5对极管蛋白锚定在孢壁上起到非常重要的作用，孢壁蛋白9定位到极管上，并且可以与极管蛋白1和极管蛋白2互相作用(Yangetal., 2018)。本研究中，一共检测到5个孢壁蛋白编码基因(GenBank登录号分别为9422782, 9422311, 9422437, 9422653和9423533)、1个孢壁蛋白前体编码基因(GenBank登录号: 9422389)和1个孢壁和锚定盘复合蛋白编码基因(GenBank登录号: 9424403)，其中孢壁蛋白编码基因(GenBank登录号分别为9422653和9423533)在NcCKvsNcT1和NcCKvsNcT2上调表达(表4)，推测它们主要在孢子侵染和增殖过程中起到重要的作用，而孢壁和锚定盘复合蛋白编码基因(GenBank登录号: 9424403)在NcCKvsNcT1和NcCKvsNcT2上调表达，推测其在东方蜜蜂微孢子虫侵染增殖过程发挥重要的作用。此外，孢壁蛋白9编码基因在NcCK中的表达量很高(FPKM值为12 059)，而在NcT1和NcT2的表达量相对较低(FPKM值分别为2 807和2 885)，暗示其在东方蜜蜂微孢子虫孢子抵御外界不良环境中起特殊作用。本研究发现，另有4个蓖麻毒素B凝集素编码基因(GenBank登录号分别为9422642, 9423906, 9422644和9422647)分别在NcCKvsNcT1(分别为4.77, 2.32, 2.72和5.62倍), NcCKvsNcT2(分别为4.98, 3.23, 3.21和3.38倍)中上调表达(表4)，再次说明蓖麻毒素B凝集素在东方蜜蜂微孢子虫侵染过程中的重要性。

东方蜜蜂微孢子虫寄生蜜蜂中肠上皮细胞，由于缺少典型的线粒体，所以自身只能通过糖酵解途径产生少量ATP，绝大多数的ATP必需依赖宿主细胞提供(Cornmanetal., 2009)。Huang等(2016)通过对东方蜜蜂微孢子虫胁迫的西方蜜蜂工蜂1-6 d的肠道进行测序，发现东方蜜蜂微孢子虫催化糖酵解的关键酶的基因在感染后1 d就开始表达，己糖激酶在感染后2 d被检测到并且在2-6 d的表达量上调，糖酵解途径的其他核心酶在感染后3 d也被检测出来，同时发现ATP/ADP移位酶基因和ABC转运蛋白基因在东方蜜蜂微孢子虫侵染过程中也维持高量表达。本研究中，涉及东方蜜蜂微孢子虫糖酵解的几种关键酶基因如己糖激酶基因(GenBank登录号: 9423428)、丙酮酸激酶基因(GenBank登录号: 9424331)和6-磷酸果糖激酶基因(GenBank登录号: 9424162)在侵染过程较孢子状态表达量上调(表5)，与前人研究结果(Huangetal., 2016)一致，再次说明东方蜜蜂微孢子虫在其增殖过程中通过糖酵解途径提供部分能量。有研究表明家蚕微孢子虫和蝗虫微孢子虫Nosemalocustae的己糖激酶都包含信号肽，均为分泌蛋白(Timofeevetal., 2017; 陈红丽等, 2018)。本研究发现，东方蜜蜂微孢子虫的己糖激酶同样也包含信号肽，推测东方蜜蜂微孢子虫在侵染过程中将己糖激酶分泌到宿主细胞质以调控宿主的糖酵解。ATP/ADP移位酶是微孢子虫窃取宿主能量的关键酶之一，Paldi等(2010)通过RNAi对东方蜜蜂微孢子虫的ATP/ADP移位酶基因进行敲减，发现东方蜜蜂微孢子虫感染的西方蜜蜂体内的病原增殖被显著抑制。本研究中，有3个ATP/ADP移位酶的编码基因(GenBank登录号分别为9424363, 9422880和9424586)在NcCKvsNcT1 [log2(Fold change)值分别为0.72, 3.60和2.31]和NcCKvsNcT2[log2(Fold change)值分别为0.45, 3.44和2.34]中的表达水平上调(表5)，表明ATP/ADP移位酶对于东方蜜蜂微孢子虫窃取宿主细胞能量、促进自身增殖的重要性。ABC转运蛋白的主要功能是结合ATP，利用ATP水解产生能量将物质逆浓度梯度运输(冯振月等, 2018)。本研究中，GenBank登录号分别为9422518和9423041的ABC转运蛋白编码基因在NcCKvsNcT1 [log2(Fold change)值分别为0.93和1.83]和NcCKvsNcT2[log2(Fold change)值分别为0.64和1.85]中上调表达，体现出对东方蜜蜂微孢子虫侵染过程中物质运输的重要性；但另外4个GenBank登录号分别为9423050, 9423588, 9422298和9422297的ABC转运蛋白编码基因在NcCKvsNcT1 [log2(Fold change)值分别为-3.65, -1.72, -1.91和-0.19] 和NcCKvsNcT2[log2(Fold change)值分别为-3.34, -0.13, -1.56和-0.35]中的表达水平也为下调(表5)，鉴于东方蜜蜂微孢子虫与蜜蜂之间存在密切的相互作用，上述3个基因的下调表达可能是由于东方蜜蜂微孢子虫受到宿主的抑制，具体的机制有待于进一步深入研究。

本研究解析了东方蜜蜂微孢子虫在侵染意大利蜜蜂工蜂过程的差异基因表达谱，在mRNA组学层面揭示了东方蜜蜂微孢子虫的侵染机制，筛选出病原的毒力因子和侵染相关因子基因可为后续功能研究提供候选靶标，并为治疗蜜蜂微孢子虫病提供潜在的分子靶点。