参与调控意大利蜜蜂工蜂中肠基因表达的东方蜜蜂微孢子虫miRNA的组学解析及其调控网络

范小雪, 杜 宇, 张文德, 王 杰, 蒋海宾, 范元婵, 冯睿蓉,万洁琦, 周紫彧, 熊翠玲,2, 郑燕珍,2, 陈大福,2,3, 郭 睿,2,3,*

(1. 福建农林大学动物科学学院(蜂学学院), 福州 350002; 2. 福建农林大学蜂疗研究所, 福州 350002;3. 福建农林大学蜂产品加工与应用教育部工程研究中心, 福州 350002)

东方蜜蜂微孢子虫Nosemaceranae是一种专性侵染成年蜜蜂中肠上皮细胞的单细胞真菌病原，可对蜜蜂宿主造成肠道破坏、能量胁迫、免疫抑制、消化紊乱及寿命缩短等不良影响(Higesetal., 2007; Mayack and Naug, 2009; Alauxetal., 2010; Botíasetal., 2013; Chenetal., 2013)，但到目前为止二者互作的分子机制尚未阐明。

微小RNA(microRNA, miRNA)是一类长度约为18～25 nt且高度保守的内源性单链非编码RNA(non-coding RNA, ncRNA)(Bartel, 2004)，可通过靶向结合mRNA使其降解或抑制其翻译，从而在转录后水平调控基因表达，进而发挥广泛而重要的生物学功能(Pillai, 2005)。此外，miRNA已被证实能够作为媒介介导真菌病原及其宿主的互作(Linetal., 2016; Crostonetal., 2018)。2012年，张辰宇团队研究发现植物来源的miR-168a可在动物体内稳定积累表达，并通过靶向抑制小鼠低密度脂蛋白LDLRAP1的表达水平，降低小鼠血浆LDL的去除率，首次证实miRNA的跨界调控作用(Zhangetal., 2012)。较多的研究证实结果表明小RNA(small RNA, sRNA)可在真菌病原及其宿主之间双向传播，真菌病原可将sRNA作为效应因子转移到宿主细胞内，参与宿主的基因表达调控，从而影响自身毒力、宿主免疫力以及疾病进程等。例如，灰霉菌Botrytiscinerea和球孢白僵菌Beauveriabassiana来源的sRNA在细胞囊泡的包裹和保护作用下转移到宿主体内，通过与宿主Argonaute 1蛋白结合并劫持宿主RNAi机制，进而选择性沉默宿主免疫基因，以达到有效感染的目的(Weibergetal., 2013; Cuietal., 2019)。对于miRNA介导的西方蜜蜂Apismellifera及东方蜜蜂微孢子虫互作，前人开展了少量研究(Evans and Huang, 2018; Huangetal., 2019)，但总体进展缓慢。前期研究中，笔者所在团队联用链特异性建库的RNA-seq和small RNA-seq(sRNA-seq)技术对东方蜜蜂微孢子虫感染7 d和10 d的意大利蜜蜂Apismelliferaligustica工蜂中肠及未感染中肠进行测序，基于高质量的全转录组数据系统解析了lncRNA, miRNA及ceRNA调控网络介导的宿主免疫应答(Chenetal., 2019; 付中民等, 2019)；此外，还在lncRNA, miRNA和mRNA组学水平深入分析和探讨了东方蜜蜂微孢子虫侵染意大利蜜蜂工蜂的分子机制(耿四海等, 2020a, 2020b)。近期，笔者所在团队基于已获得的东方蜜蜂微孢子虫感染及未感染的意大利蜜蜂工蜂中肠转录组数据(付中民等, 2019)以及纯净的病原孢子全转录组数据(耿四海等, 2020a)筛选出宿主的差异表达miRNA(differentially expressed miRNA, DEmiRNA)和病原的差异表达mRNA(differentially expressed mRNA, DEmRNA)，通过生物信息学手段预测宿主DEmiRNA和病原DEmRNA的靶向关系，并对二者间的调控网络进行了构建和分析，进一步结合本团队的前期研究结果和前人的相关研究结果探讨了意大利蜜蜂工蜂的DEmiRNA对东方蜜蜂微孢子虫基因和通路的调控作用(未发表数据)。在长期的协同进化中，东方蜜蜂微孢子虫的基因组大大简化，多数物质和能量代谢的途径丢失，其增殖高度依赖宿主细胞提供物质和能量(Burrietal., 2006)。同时东方蜜蜂微孢子虫也进化出一些特殊的能力，能够对宿主的生理生化功能进行操纵，如东方蜜蜂微孢子虫侵染可对西方蜜蜂工蜂造成较强的能量胁迫，致使宿主的能量代谢速率加快且糖水的摄入量提高(Antúnezetal., 2009; Chenetal., 2013)。然而，对于东方蜜蜂微孢子虫通过何种方式和机制对西方蜜蜂的基因表达和生理活动进行调控，相关研究仍然匮乏，迄今仅有一例报道(Evans and Huang, 2018)。此前，Evans和Huang(2018)通过高通量测序和生物信息学分析方法在东方蜜蜂微孢子虫感染的西方蜜蜂工蜂中肠中预测出病原的6条miRNA，进一步分析发现其中5条miRNA不仅能靶向病原本身的1 545条mRNA，还可能被分泌到宿主细胞的胞质内，通过靶向918条宿主mRNA调控宿主的新陈代谢、细胞凋亡和免疫防御等过程。

本研究拟通过比较分析筛选出侵染意大利蜜蜂工蜂的东方蜜蜂微孢子虫的DEmiRNA以及意大利蜜蜂工蜂响应东方蜜蜂微孢子虫侵染的意大利蜜蜂DEmRNA，通过生物信息学方法预测病原DEmiRNA靶向的宿主mRNA和DEmRNA，并通过注释数据库获得相应的功能和通路信息，进一步结合东方蜜蜂微孢子虫侵染西方蜜蜂的生物学背景，笔者所在课题组的前期研究结果以及前人的相关研究结果深入探讨东方蜜蜂微孢子虫的DEmiRNA对意大利蜜蜂工蜂中肠基因表达的调控作用。研究结果可为探明东方蜜蜂微孢子虫的侵染机制、东方蜜蜂微孢子虫-意大利蜜蜂工蜂互作机制提供重要的参考信息和基础。

1 材料与方法

1.1 东方蜜蜂微孢子虫的miRNA组学数据来源

前期研究中，笔者所在课题组分别提取了3个东方蜜蜂微孢子虫的纯净孢子样品NcCK(NcCK1, NcCK2, NcCK3)的总RNA，并按照已建立的方法(耿四海等, 2020a)构建了相应的cDNA文库，单端测序由广州基迪奥生物科技有限公司完成，测序平台为Illumina MiSeq。通过连续比对西方蜜蜂参考基因组(Assembly Amel_4.5)、GenBank数据库、Rfam数据库及东方蜜蜂微孢子虫参考基因组(Assembly ASM98816v1)筛滤得到比对上东方蜜蜂微孢子虫参考基因组的东方蜜蜂微孢子虫纯净孢子miRNA组学数据(NcCK)，侵染意大利蜜蜂工蜂7 d(NcT1)和10 d(NcT2)的东方蜜蜂微孢子虫miRNA组学数据(耿四海等, 2020a)。测序数据已上传至NCBI SRA数据库，NcCK的Bioproject登记号为PRJNA395264。通过对下机的原始数据进行质控得到了最终的有效标签序列(clean tags)(耿四海等, 2020a)。

1.2 东方蜜蜂微孢子虫的DEmiRNA筛选

按照耿四海等(2020a)的方法进行东方蜜蜂微孢子虫miRNA的鉴定和表达量计算以及NcCKvsNcT1和NcCKvsNcT2比较组显著性DEmiRNA的筛选，筛选出的显著性DEmiRNA(|log2fold change|≥1且P≤0.05)用于本研究中靶向意大利蜜蜂工蜂中肠mRNA和DEmRNA的预测和分析。

1.3 意大利蜜蜂工蜂中肠的mRNA组学数据来源

笔者所在团队前期已对接种50%(w/v)蔗糖溶液7 d(AmCK1)和10 d(AmCK2)的意大利蜜蜂工蜂中肠样品以及接种含东方蜜蜂微孢子虫孢子的50%(w/v)蔗糖溶液7 d (AmT1)和10 d(AmT2)的意大利蜜蜂工蜂中肠样品进行RNA抽提、cDNA文库构建、转录组测序以及测序数据质控(付中民等, 2019)。测序数据已上传至NCBI SRA数据库，Bioproject号为PRJNA395264。质控后的意大利蜜蜂工蜂中肠mRNA组学数据可用于本研究中东方蜜蜂微孢子虫的显著性DEmiRNA的靶标预测及分析。

1.4 东方蜜蜂微孢子虫显著性DEmiRNA靶向意大利蜜蜂工蜂中肠mRNA的预测

为探究1.2节筛选出的东方蜜蜂微孢子虫DEmiRNA调控意大利蜜蜂工蜂中肠基因表达的广泛性，利用TargetFinder软件(Allenetal., 2005)预测NcCKvsNcT1和NcCKvsNcT2比较组中显著性DEmiRNA靶向的1.3节的意大利蜜蜂工蜂中肠的全部mRNA，采用默认参数。

1.5 意大利蜜蜂工蜂中肠响应东方蜜蜂微孢子虫侵染的DEmRNA的筛选及分析

miRNA一般对基因进行负调控(Pillai, 2005)，故差异变化趋势相反的miRNA-mRNA关系有重要意义。因此，本研究对东方蜜蜂微孢子虫上调miRNA-意大利蜜蜂工蜂下调mRNA、东方蜜蜂微孢子虫下调miRNA-意大利蜜蜂工蜂上调mRNA关系进行重点分析和探讨。采用FPKM(fragments per kilobase of transcript per million fragments mapped)法计算和归一化基因表达量。利用edgeR软件(Robinsonetal., 2010)筛选AmCK1vsAmT1和AmCK2vsAmT2比较组的显著性DEmRNA，筛选标准为：|log2fold change|≥1且P≤0.05。

利用OmicShare云平台(www.omicshare.com)的GO分析软件和KEGG通路分析软件对1.4节的靶mRNA和上述靶DEmRNA进行GO数据库(George, 2008)和KEGG(Kanehisa and Goto, 2006)通路注释，采用默认参数。

1.6 东方蜜蜂微孢子虫显著性DEmiRNA与意大利蜜蜂工蜂中肠DEmRNA的调控网络构建及分析

根据1.4和1.5节预测出的东方蜜蜂微孢子虫显著性DEmiRNA与意大利蜜蜂工蜂中肠显著性DEmRNA的靶向关系，构建NcCKvsNcT1比较组中DEmiRNA与AmCK1vsAmT1比较组中DEmRNA以及NcCKvsNcT2比较组中DEmiRNA与AmCK2vsAmT2比较组中DEmRNA的调控网络，并利用Cytoscape软件(Smootetal., 2011)进行可视化，采用默认参数。

2 结果

2.1 东方蜜蜂微孢子虫DEmiRNA靶向意大利蜜蜂工蜂中肠mRNA和DEmRNA的预测及分析

NcCKvsNcT1比较组中东方蜜蜂微孢子虫的165条显著性DEmiRNA整体上可靶向意大利蜜蜂工蜂中肠的11 711条mRNA，其中东方蜜蜂微孢子虫的91条显著上调和55条显著下调miRNA分别靶向意大利蜜蜂工蜂中肠AmCK1vsAmT1的9 920和8 391条mRNA，其余19条DEmiRNA与意大利蜜蜂工蜂中肠mRNA之间无靶向关系。进一步分析发现东方蜜蜂微孢子虫的77条显著上调miRNA和52条显著下调miRNA可分别靶向意大利蜜蜂工蜂中肠的118条显著下调mRNA和135条显著上调mRNA(图1)。

NcCKvsNcT2比较组中东方蜜蜂微孢子虫的124条显著性DEmiRNA整体上可靶向意大利蜜蜂工蜂中肠的11 295条mRNA，其中东方蜜蜂微孢子虫的56条显著上调和49条显著下调miRNA可分别靶向AmCK2vsAmT2比较组中意大利蜜蜂工蜂中肠的8 138和9 055条mRNA，其余19条东方蜜蜂微孢子虫DEmiRNA与意大利蜜蜂工蜂中肠mRNA之间无靶向关系。进一步分析结果显示，东方蜜蜂微孢子虫的52条显著上调miRNA可靶向AmCK2vsAmT2中意大利蜜蜂工蜂中肠的97条显著下调mRNA，东方蜜蜂微孢子虫的49条显著下调miRNA可靶向意大利蜜蜂工蜂中肠的210条显著上调mRNA(图2)。

图1 NcCK vs NcT1比较组中东方蜜蜂微孢子虫的显著性DEmiRNA与AmCK1 vs AmT1比较组中意大利蜜蜂工蜂中肠的显著性DEmRNA之间的调控网络

图2 NcCK vs NcT2比较组中东方蜜蜂微孢子虫的显著性DEmiRNA与AmCK2 vs AmT2比较组中意大利蜜蜂工蜂中肠的显著性DEmRNA之间的调控网络

进一步分析发现，NcCKvsNcT1和NcCKvsNcT2中的11条共有的显著上调miRNA(miR-3968-y, miR-5119-y和miR-317-y等)可靶向AmCK1vsAmT1和AmCK2vsAmT2中的6条皆显著下调的mRNA (GenBank登录号: XM_006564187.2, XM_006567501.2和XM_006558063.2等)；19条共有的显著下调miRNA(novel-m0016-3p, let-7-x和miR-122-x等)靶向AmCK1vsAmT1和AmCK2vsAmT2中14条皆显著上调的mRNA (GenBank登录号: XM_016912123.1, XM_006558673.2和XM_016916036.1等)。

NcCKvsNcT1中共63条特有的DEmiRNA，其中36条显著上调miRNA和18条显著下调miRNA可分别靶向AmCK1vsAmT1中的97条显著下调mRNA和91条显著上调mRNA，其余9条东方蜜蜂微孢子虫DEmiRNA与意大利蜜蜂工蜂中肠mRNA之间不存在靶向关系。NcCKvsNcT2中共22条特有的DEmiRNA，其中8条显著上调miRNA和14条显著下调miRNA可分别靶向AmCK2vsAmT2中51条显著下调mRNA和139条显著上调mRNA。

2.2 东方蜜蜂微孢子虫显著性DEmiRNA靶向意大利蜜蜂工蜂中肠显著性DEmRNA的GO数据库注释

NcCKvsNcT1中的显著上调miRNA靶向AmCK1vsAmT1中意大利蜜蜂工蜂中肠的显著下调mRNA涉及细胞(25条下调mRNA)、细胞组分(25条下调mRNA)和细胞膜(24条下调mRNA)等10个细胞组分相关条目，细胞进程(75条下调mRNA)、代谢进程(60条下调mRNA)和单一组织进程(49条下调mRNA)等15个生物学进程相关条目，以及结合(59条下调mRNA)、催化活性(55条下调mRNA)和转运活性(11条下调mRNA)等6个分子功能相关条目(图3: A)。东方蜜蜂微孢子虫的显著下调miRNA靶向的意大利蜜蜂工蜂中肠显著上调mRNA涉及细胞膜(15条上调mRNA)、细胞膜组分(15条上调mRNA)和细胞(12条上调mRNA)等7个细胞组分相关条目，细胞进程(34条上调mRNA)、代谢进程(30条上调mRNA)和单一组织进程(21条上调mRNA)等10个生物学进程相关条目，以及结合(34条上调mRNA)、催化活性(24条上调mRNA)和核酸结合转录因子活性(5条上调mRNA)等8个分子功能相关条目(图3: B)。

图3 NcCK vs NcT1比较组中东方蜜蜂微孢子虫显著性DEmiRNA靶向AmCK1 vs AmT1比较组中意大利蜜蜂工蜂中肠显著性DEmRNA的GO数据库注释

NcCKvsNcT2中东方蜜蜂微孢子虫的显著上调miRNA靶向AmCK2vsAmT2中意大利蜜蜂工蜂中肠的显著下调mRNA涉及细胞膜组分(20条下调mRNA)、细胞膜(20条下调mRNA)和细胞(14条下调mRNA)等11个细胞组分相关条目，细胞进程(52条下调mRNA)、代谢进程(42条下调mRNA)和单一组织进程(35条下调mRNA)等10个生物学进程相关条目，以及结合(49条下调mRNA)、催化活性(41条下调mRNA)和转运活性(9条下调mRNA)等6个分子功能相关条目(图4: A)。东方蜜蜂微孢子虫的显著下调miRNA靶向的意大利蜜蜂工蜂中肠显著上调mRNA涉及细胞(24条上调mRNA)、细胞组分(24条上调mRNA)和细胞膜组分(22条上调mRNA)等9个细胞组分相关条目，细胞进程(59条上调mRNA)、代谢进程(47条上调mRNA)和单一组织进程(42条上调mRNA)等12个生物学进程相关条目，以及结合(57条上调mRNA)、催化活性(35条上调mRNA)和转运活性(14条上调mRNA)等9个分子功能相关条目(图4: B)。

此外，NcCKvsNcT1和NcCKvsNcT2中共同的显著上调miRNA靶向AmCK1vsAmT1和AmCK2vsAmT2中共同的显著下调mRNA可注释到单一组织进程(2条下调mRNA)、细胞进程(2条下调mRNA)和结合(2条下调mRNA)等7个条目；NcCKvsNcT1和NcCKvsNcT2中共同的显著下调miRNA靶向AmCK1vsAmT1和AmCK2vsAmT2中共同的显著上调mRNA可注释到结合(3条上调mRNA)、代谢进程(2条上调mRNA)和细胞(2条上调mRNA)等10个条目。

NcCKvsNcT1中特有的显著上调miRNA靶向AmCK1vsAmT1中显著下调mRNA涉及26个条目，包括细胞膜组分(9条下调mRNA)、细胞膜(9条下调mRNA)和细胞(8条下调mRNA)等7个细胞组分相关条目，细胞进程(27条下调mRNA)、单一组织进程(21条下调mRNA)和代谢进程(21条下调mRNA)等13个生物学进程相关条目，结合(24条下调mRNA)、催化活性(19条下调mRNA)、分子转导活性(5条下调mRNA)等6个分子功能相关的条目。NcCKvsNcT1特有的东方蜜蜂微孢子虫显著下调miRNA靶向AmCK1vsAmT1中显著上调mRNA涉及22个功能条目，包括9个生物学进程相关条目，7个细胞组分相关条目，以及6个分子功能相关条目。NcCKvsNcT2中特有的显著上调miRNA靶向AmCK2vsAmT2中显著下调mRNA可注释到22个功能条目，包括细胞(10条下调mRNA)、细胞组分(10条下调mRNA)和细胞器(9条下调mRNA)等7个细胞组分相关条目，细胞进程(19条下调mRNA)、代谢进程(16条下调mRNA)和单一组织进程(10条下调mRNA)等8个生物学进程相关条目，结合(24条下调mRNA)、催化活性(13条下调mRNA)和核酸结合转录因子活性(4条下调mRNA)等7个分子功能相关条目；NcCKvsNcT2中特有的东方蜜蜂微孢子虫显著下调miRNA靶向AmCK2vsAmT2中意大利蜜蜂工蜂显著上调mRNA涉及28个功能条目，包括9个细胞组分相关条目，11个生物学进程相关条目，以及8个分子功能相关条目。

2.3 东方蜜蜂微孢子虫显著性DEmiRNA靶向意大利蜜蜂工蜂中肠显著性DEmRNA的KEGG通路注释

NcCKvsNcT1中东方蜜蜂微孢子虫的显著上调miRNA靶向AmCK1vsAmT1中意大利蜜蜂工蜂中肠的显著下调mRNA可注释到剪接体(4条下调mRNA)、光传导-果蝇(3条下调mRNA)和磷酸肌醇代谢(3条下调mRNA)等113条通路，其中包括胞吞作用(2条下调mRNA)等细胞免疫通路，以及泛素介导的蛋白水解(2条下调mRNA)、黑色素生成(3条下调mRNA)等体液免疫通路(图5: A)。NcCKvsNcT1中东方蜜蜂微孢子虫的显著下调miRNA靶向意大利蜜蜂工蜂中肠的显著上调mRNA可注释到Hedgehog信号通路-果蝇(3条上调mRNA)、粘着斑(3条上调mRNA)和生理节律-果蝇(2条上调mRNA)等107条通路，其中包括胞吞作用(4条上调mRNA)、泛素介导的蛋白水解(1条上调mRNA)等细胞免疫通路，黑色素生成(1条上调mRNA)等体液免疫通路，以及色氨酸代谢(1条上调mRNA)等3条氨基酸代谢通路，柠檬酸循环(TCA循环)(2条上调mRNA)等4条碳水化合物代谢通路和氧化磷酸化(1条上调mRNA)等2条能量代谢通路(图5: B)。注释数量前20位的通路信息详见图5。

图5 NcCK vs NcT1比较组中东方蜜蜂微孢子虫显著性DEmiRNA靶向AmCK1 vs AmT1比较组中意大利蜜蜂工蜂中肠显著性DEmRNA注释的前20位KEGG通路

NcCKvsNcT2中东方蜜蜂微孢子虫的显著上调miRNA靶向AmCK2vsAmT2中意大利蜜蜂工蜂中肠的显著下调mRNA涉及代谢通路(11条下调mRNA)、鞘脂信号通路(5条下调mRNA)及硫代谢(1条下调mRNA)等97条通路，其中包括溶酶体(1条下调mRNA)、黑色素生成(1条下调mRNA)、胞吞作用(1条下调mRNA)等细胞免疫通路。NcCKvsNcT2中东方蜜蜂微孢子虫的显著下调miRNA靶向意大利蜜蜂工蜂中肠的显著上调mRNA涉及代谢通路(9条上调mRNA)、Hippo信号通路(5条上调mRNA)和其他聚糖降解(4条上调mRNA)等127条通路，其中包括溶酶体(4条上调mRNA)、胞吞作用(2条上调mRNA)等细胞免疫通路，以及Toll/Imd信号通路(1条上调mRNA)、Toll样受体信号通路(1条上调mRNA)等体液免疫通路，嘧啶代谢(1条上调mRNA)等2条核苷酸代谢通路，角质、亚黄素和蜡的生物合成(1条上调mRNA)等3条脂质代谢通路，谷胱甘肽代谢(1条上调mRNA)等2条其他氨基酸的代谢通路，1条萜类和聚酮类化合物的代谢通路。注释数量前20位的通路信息详见图6。

图6 NcCK vs NcT2比较组中东方蜜蜂微孢子虫显著性DEmiRNA靶向AmCK2 vs AmT2比较组中意大利蜜蜂工蜂中肠显著性DEmRNA注释的前20位KEGG通路

此外，NcCKvsNcT1和NcCKvsNcT2中共同的显著上调miRNA靶向AmCK1vsAmT1和AmCK2vsAmT2中共同的显著下调mRNA未注释到任何通路；共同的东方蜜蜂微孢子虫显著下调miRNA靶向共同的意大利蜜蜂工蜂中肠显著上调mRNA可注释到9条通路，包括生理节律(1条上调mRNA)、生理节律-果蝇(1条上调mRNA)、Hippo信号通路-多物种(1条上调mRNA)、Hedgehog信号通路(1条上调mRNA)、Hedgehog信号通路-果蝇(1条上调mRNA)、FoxO信号通路(1条上调mRNA)、Hippo信号通路(1条上调mRNA)、Hippo信号通路-果蝇(1条上调mRNA)和Wnt信号通路(1条上调mRNA)等。

NcCKvsNcT1中特有的显著上调miRNA靶向AmCK1vsAmT1中特有的显著下调mRNA涉及代谢途径(5条下调mRNA)、剪接体(4条下调mRNA)和光传导-果蝇(3条下调mRNA)等106条通路。NcCKvsNcT1东方蜜蜂微孢子虫特有的显著下调miRNA靶向意大利蜜蜂工蜂中肠特有的显著上调mRNA涉及代谢途径(4条上调mRNA)、Hippo信号通路(3条上调mRNA)和生理节律-果蝇(2条上调mRNA)等92条通路。

NcCKvsNcT2中特有的显著上调miRNA靶向AmCK2vsAmT2中显著下调mRNA涉及代谢通路(5条下调mRNA)、糖胺聚糖降解(1条下调mRNA)及Notch信号通路(1条下调mRNA)等36条通路。NcCKvsNcT2中东方蜜蜂微孢子虫特有的显著下调miRNA靶向意大利蜜蜂工蜂中肠特有的显著上调mRNA涉及各种类型的N-聚糖的生物合成(4条上调mRNA)、溶酶体(4条上调mRNA)及cAMP信号通路(4条上调mRNA)等112条通路。

2.4 东方蜜蜂微孢子虫显著性DEmiRNA及其靶向意大利蜜蜂工蜂中肠的免疫防御相关DEmRNA的调控网络

根据2.3节通路注释的结果，东方蜜蜂微孢子虫显著性DEmiRNA靶向意大利蜜蜂工蜂中肠DEmRNA涉及的免疫防御相关通路包括：胞吞作用、黑色素生成、溶酶体、自噬、Toll样受体信号通路、细胞凋亡、Ras信号通路、泛素介导的蛋白水解和MAPK信号通路。

图7 东方蜜蜂微孢子虫的显著性DEmiRNA靶向意大利蜜蜂工蜂中肠免疫防御和能量代谢相关通路注释的显著性DEmRNA的调控网络

NcCKvsNcT1中显著上调的miR-14-y [log2(fold change)=16.41,P=3.76E-05]靶向意大利蜜蜂工蜂中肠的mRNA (GenBank登录号: XM_395448.6) [log2(fold change)=-10.02,P=0.0076]； miR-283-x [log2(fold change)=5.50,P=0.049]和miR-216-x [log2(fold change)=16.33,P=1.43E-07]共同靶向意大利蜜蜂工蜂中肠的E3泛素蛋白连接酶Nedd-4亚型X8编码基因的mRNA (GenBank登录号: XM_016917135.1) [log2(fold change)=-10.57,P=0.0020]。NcCKvsNcT1中显著下调的miR-181-x[log2(fold change)=-2.22,P=0.0051]和miR-128-y[log2(fold change)=-1.87,P=5.55E-04]共同靶向的意大利蜜蜂工蜂中肠mRNA (GenBank登录号: XM_016915020.1) [log2(fold change)=1.04,P=0.019], miR-340-x[log2(fold change)=-1.12,P=0.024]靶向的mRNA(GenBank登录号: XM_006570195.2)[log2(fold change)=9.52,P=0.043], miR-484-z[log2(fold change)=-1.64,P=0.029]靶向的mRNA (GenBank登录号: XM_006571076.2)[log2(fold change)=9.40,P=0.033]和miR-320-y[log2(fold change)=-4.73,P=1.18E-04]靶向的mRNA (GenBank登录号: XM_394836.6) [log2(fold change)=4.01,P=0.0086]均可注释到胞吞作用(图7: A)。显著上调的miR-3202-x[log2(fold change)=12.47,P=1.75E-07], miR-317-y[log2(fold change)=19.06,P=1.72E-06]和miR-252-x[log2(fold change)=14.99,P=9.90E-05]等6条东方蜜蜂微孢子虫DEmiRNA共同靶向的意大利蜜蜂工蜂中肠mRNA (GenBank登录号: XM_016911563.1) [log2(fold change)=-8.28,P=0.0093]， miR-283-x, miR-8-y [log2(fold change)=7.03,P=0.032]和miR-216-x等7条DEmiRNA共同靶向的1-磷脂酰肌醇-4,5-二磷酸磷酸二酯酶I类和II类亚型X1编码基因mRNA (GenBank登录号: XM_003250942.3) [log2(fold change)=-3.41,P=0.0084]，miR-339-y [log2(fold change)=12.47,P=1.75E-07], miR-9226-y [log2(fold change)=21.32,P=1.98E-10]和miR-317-y 4条DEmiRNA共同靶向mRNA (GenBank登录号: XM_006559228.1) [log2(fold change)=-2.45,P=0.045]以及显著下调的novel-m0016-3p [log2(fold change)=-7.07,P=1.82E-05]和novel-m0028-5p [log2(fold change)=-14.19,P=3.86E-07]共同靶向mRNA (GenBank登录号: XM_006567354.2) [log2(fold change)=2.08,P=9.21E-04]均可注释到黑色素生成(图7: A)。显著上调的miR-14-y靶向的意大利蜜蜂工蜂中肠mRNA (GenBank登录号: XM_395448.6) [log2(fold change)=-10.02,P=0.0076], miR-283-x和miR-216-x共同靶向的mRNA (GenBank登录号: XM_016917135.1)以及显著下调的miR-106-x [log2(fold change)=-2.30,P=0.026], miR-17-x [log2(fold change)=-1.26,P=0.036], let-7-y [log2(fold change)=-3.76,P=0.0070], miR-20-x [log2(fold change)=-3.70,P=9.29E-04)和miR-146-x [log2(fold change)=-2.86,P=0.037)共同靶向的mRNA (GenBank登录号: XM_006564220.2) [log2(fold change)=10.41,P=0.010]均可注释到泛素介导的蛋白水解(图7: A)。显著上调的miR-283-x, miR-216-x和miR-4796-y [log2(fold change)=15.51,P=1.21E-06]共同靶向的mRNA (GenBank登录号: XM_016914371.1)，显著下调的miR-126-y [log2(fold change)=-3.77,P=0.0098]和miR-128-y共同靶向的mRNA (GenBank登录号: XM_394432.6) [log2(fold change)=7.97,P=0.035]均可注释到MAPK信号通路(图7: A)。

在NcCKvsNcT2中，东方蜜蜂微孢子虫显著下调的novel-m0016-3p [log2(fold change)=-4.37,P=3.24E-04]靶向意大利蜜蜂工蜂中肠的mRNA (GenBank登录号: XM_016917619.1) [log2(fold change)=12.12,P=4.14E-05]和miR-320-y [log2(fold change)=-4.87,P=0.0016]靶向的mRNA (GenBank登录号: XM_006572368.2) [log2(fold change)=8.59,P=0.0041]以及显著上调的miR-5106-y [log2(fold change)=19.69,P=9.27E-09]靶向的mRNA (GenBank登录号: XM_006561133.2)[log2(fold change)=-1.13,P=0.045]皆可注释到胞吞作用(图7: B)。显著下调的novel-m0021-5p [log2(fold change)=-1.52,P=0.036]靶向意大利蜜蜂工蜂中肠的4条mRNAs (GenBank登录号: XM_006567428.2) [log2(fold change)=1.52,P=0.0042]，(GenBank登录号: XM_006567431.2) [log2(fold change)=1.34,P=0.039]，(GenBank登录号: XM_016914010.1) [log2(fold change)=1.22,P=1.47E-04]和(GenBank登录号: XM_016914011.1) [log2(fold change)=3.09,P=0.042]以及显著上调的novel-m0022-5p [log2(fold change)=14.13,P=7.97E-06]和miR-216-x [log2(fold change)=16.33,P=1.43E-07]共同靶向的类α-乙酰氨基葡萄糖苷酶编码基因mRNA (GenBank登录号: XM_016913157.1) [log2(fold change)=-9.33,P=0.010]皆可注释到溶酶体(图7: B)。显著下调的novel-m0002-3p [log2(fold change)=-5.85,P=6.42E-05], novel-m0021-5p和miR-374-x[log2(fold change)=-15.48,P=1.23E-07]共同靶向的mRNA (GenBank登录号: XM_016911751.1) [log2(fold change)=1.65,P=0.0040]以及miR-26-x [log2(fold change)=-1.81,P=0.012]和miR-486-x [log2(fold change)=-8.86,P=0.028]共同靶向的mRNA (GenBank登录号: XM_006557238.2) [log2(fold change)=3.73,P=0.0047]皆可注释到自噬-动物(图7: B)。显著上调的bantam-y [log2(fold change)=16.57,P=1.84E-06], miR-8-y [log2(fold change)=5.53,P=0.044], miR-216-x和miR-5119-y [log2(fold change)=9.19,P=0.022]共同靶向的1-磷脂酰肌醇-4,5-二磷酸磷酸二酯酶I类和II类亚型X3编码基因mRNA (GenBank登录号: XM_006567225.2) [log2(fold change)=-8.75,P=0.016]可注释到黑色素生成(图7: B)。miR-424-x [log2(fold change)=13.86,P=9.14E-06], miR-7465-x [log2(fold change)=15.36,P=4.75E-07]和miR-547-x [log2(fold change)=13.98,P=5.96E-08]共同靶向的mRNA (GenBank登录号: XM_006564132.2) [log2(fold change)=-8.92,P=0.038]可注释到MAPK信号通路-果蝇。显著下调的miR-26-x和miR-486-x共同靶向的mRNA (GenBank登录号: XM_006557238.2)可同时注释到细胞凋亡、MAPK信号通路和Toll样受体信号通路(图7: B)。显著上调的miR-5106-y靶向的mRNA (GenBank登录号: XM_006561133.2),显著下调的miR-26-x和miR-486-x共同靶向的mRNA (GenBank登录号: XM_006557238.2), novel-m0021-5p, miR-221-z [log2(fold change)=-2.89,P=0.0039], miR-374-x, miR-423-y [log2(fold change)=-2.04,P=0.012], miR-17-x [log2(fold change)=-1.92,P=0.039]和miR-16-y [log2(fold change)=-16.85,P=5.38E-08]共同靶向的mRNA (GenBank登录号: XM_006569471.2) [log2(fold change)=5.64,P=0.027]以及显著上调的miR-8-y, bantam-y和novel-m0022-5p [log2(fold change)=14.13,P=7.97E-06]等10条miRNA共同靶向的1-磷脂酰肌醇-4,5-二磷酸磷酸二酯酶ε-1编码基因mRNA (GenBank登录号: XM_016911813.1) [log2(fold change)=-7.89，P=0.023]皆可注释到Ras信号通路(图7: B)。

此外，根据2.3节通路注释的结果，东方蜜蜂微孢子虫显著性DEmiRNA靶向的意大利蜜蜂工蜂中肠DEmRNA可注释到氧化磷酸化和硫代谢等能量代谢相关通路。NcCKvsNcT1中显著下调的miR-584-x [log2(fold change)=-1.32,P=0.042]和miR-151-y [log2(fold change)=-2.21,P=6.92E-04]等2条miRNA共同靶向的细胞色素c氧化酶亚基4亚型X1，线粒体亚型X1编码基因mRNA (GenBank登录号: XM_006572254.2) [log2(fold change)=-2.21,P=6.92E-04],显著上调的miR-459-x[log2(fold change)=16.24,P=1.76E-07], miR-315-x[log2(fold change)=12.80,P=4.70E-06]和miR-5119-y分别靶向的mRNAs (GenBank登录号: XM_001122063.4)[log2(fold change)=-1.07,P=0.022]、(GenBank登录号: XM_006567193.2)[log2(fold change)=-1.02,P=8.31E-04]和(GenBank登录号: XM_623229.5) [log2(fold change)=-2.86,P=0.0095]皆可注释到氧化磷酸化(图7: C)。NcCKvsNcT2中显著上调的miR-301-y [log2(fold change)=15.58,P=6.07E-09]靶向的mRNA (GenBank登录号: XM_396499.6) [log2(fold change)=-2.09,P=0.0048]可注释到硫代谢(图7: D)。

NcCKvsNcT1中特有的miR-186-x [log2(fold change)=1.02,P=0.047767213]靶向意大利蜜蜂工蜂中肠的mRNAs (GenBank登录号: XM_006559228.1)和(GenBank登录号: NM_001011611.2) [log2(fold change)=-4.44,P=0.0048]可分别注释到黑色素生成通路和免疫系统进程条目，NcCKvsNcT2中特有的novel-m0022-5p靶向意大利蜜蜂工蜂中肠的mRNA (GenBank登录号: XM_006559395.1) [log2(fold change)=-1.23,P=0.0088], mRNA(GenBank登录号: XM_006560586.2) [log2(fold change)=-1.13,P=1.03E-07]和mRNA(GenBank登录号: XM_016911813.1)可注释到代谢途径。

NcCKvsNcT1中特有的miR-22-x [log2(fold change)=-13.99,P=6.33E-06]靶向意大利蜜蜂工蜂中肠的mRNAs (GenBank登录号: XM_006563625.2) [log2(fold change)=8.42,P=0.027]和(GenBank登录号: XM_006559419.2) [log2(fold change)=1.30,P=0.0078]可分别注释到生物附着和应激反应条目。NcCKvsNcT2中特有的miR-26-x靶向意大利蜜蜂工蜂中肠的mRNA (GenBank登录号: XM_006557238.2)可同时注释到Toll样受体信号通路、Hippo信号通路、MAPK信号通路及生物进程调节功能条目。

3 讨论

3.1 东方蜜蜂微孢子虫DEmiRNA对意大利蜜蜂工蜂中肠的基因表达具有潜在的广泛影响

本研究中，NcCKvsNcT1中东方蜜蜂微孢子虫的显著下调miRNA靶向的意大利蜜蜂工蜂中肠显著上调mRNA涉及7个细胞组分相关条目，10个生物学进程相关条目，8个分子功能相关条目(图3: B)；以及3条氨基酸代谢通路，4条碳水化合物代谢通路，2条能量代谢通路(图5: B)。NcCKvsNcT2中东方蜜蜂微孢子虫的显著下调miRNA靶向意大利蜜蜂工蜂中肠的显著上调mRNA涉及9个细胞组分相关条目，12个生物学进程相关条目，9个分子功能相关条目(图4: B)；以及2条核苷酸代谢通路，3条脂质代谢通路，2条其他氨基酸的代谢通路，1条萜类和聚酮类化合物的代谢通路(图6: B)。这些结果暗示东方蜜蜂微孢子虫在侵染意大利蜜蜂工蜂中肠的不同时间点通过差异表达部分miRNA对相应的意大利蜜蜂工蜂基因表达进行调控，进而影响上述与意大利蜜蜂工蜂物质和能量代谢相关的多个方面。

此外，NcCKvsNcT1和NcCKvsNcT2比较组共同的显著下调miRNA靶向AmCK1vsAmT1和AmCK2vsAmT2比较组共同的显著上调mRNA可注释到10个GO条目，以及9条KEGG通路。上述结果表明东方蜜蜂微孢子虫在侵染意大利蜜蜂工蜂中肠的过程中可能通过持续差异表达一些miRNA对宿主的细胞活动、生化反应、生理活动和信号转导产生广泛的影响。

3.2 东方蜜蜂微孢子虫可能通过下调表达部分miRNA对意大利蜜蜂工蜂中肠的能量代谢进行调控

在长期的协同进化中，东方蜜蜂微孢子虫的基因组大大简化，多数物质和能量代谢的途径丢失，其增殖高度依赖宿主细胞提供物质和能量(Burrietal., 2006)。Chen等(2013)研究发现东方蜜蜂微孢子虫侵染可对西方蜜蜂工蜂造成较强的能量胁迫，致使宿主的能量代谢速率加快且糖水的摄入量提高。ATP/ADP转移酶和ABC转运蛋白参与微孢子虫对宿主物质和能量的窃取(Paldietal., 2010)。前期研究中，笔者所在课题组发现东方蜜蜂微孢子虫的7个ABC转运蛋白编码基因在NcCKvsNcT1和NcCKvsNcT2中差异表达(耿四海等, 2020b)。本研究发现，NcCKvsNcT1中miR-151-y和miR-584-x显著下调表达，且能靶向意大利蜜蜂工蜂中肠的mRNA (GenBank登录号: XM_006572254.2)，该基因可注释到能量代谢相关的氧化磷酸化信号通路。上述结果暗示东方蜜蜂微孢子虫的DEmiRNA参与意大利蜜蜂工蜂中肠的关键能量代谢过程，控制宿主能量代谢增强以产生更多的ATP，从而掠夺更多的ATP利于自身增殖。

3.3 东方蜜蜂微孢子虫可能通过上调表达部分miRNA调控意大利蜜蜂工蜂中肠的免疫防御

昆虫的天然免疫系统包括细胞免疫和体液免疫，前者主要是指东方蜜蜂微孢子虫体被血淋巴细胞等识别、包被和吞噬，而后者主要包括各种识别受体、胞内信号因子、抗微生物肽、蛋白酶以及酶级联反应导致血淋巴凝结、黑化等(Leulieretal., 2000)。较多的研究证明，东方蜜蜂微孢子虫在侵染宿主的过程中可通过多种方式跨界调控宿主的基因表达，并抑制宿主免疫防御(Wangetal., 2017; Shahidetal., 2018; Cuietal., 2019; Huangetal., 2019)。菟丝子将Bc-siR3.2装载到AGO1蛋白，靶向拟南芥Arabidopsisthaliana的丝裂原活化蛋白激酶2基因(MPK2)和MPK1基因的转录本，从而抑制宿主的免疫反应(Shahidetal., 2018)。球孢白僵菌将bba-milR-1载入囊泡并转运到斯氏按蚊Anophelesstephensi的细胞内，进而跨界调控宿主基因CLIPB9和Spz4的表达，抑制Toll样受体信号通路和逃避黑化反应(Cuietal., 2019)。本研究发现，NcCKvsNcT1和NcCKvsNcT2比较组中miR-216-x, miR-8-y, miR-4796-y, bantam-y和miR-5119-y等多个共同显著上调的miRNA可靶向注释到黑色素生成、MAPK信号通路、Ras信号通路、泛素介导的蛋白水解和溶酶体等通路的多条显著下调mRNA。其中，NcCKvsNcT1中的miR-216-x同时靶向AmCK1vsAmT1中可注释到宿主免疫相关的黑色素生成、MAPK信号通路及泛素介导的蛋白水解通路的mRNAs (GenBank登录号分别为XM_003250942.3, XM_016914371.1和XM_016917135.1)(图7: A)；NcCKvsNcT2中的miR-216-x还可同时靶向AmCK2vsAmT2中可注释到意大利蜜蜂工蜂免疫相关的黑色素生成、溶酶体及Ras信号通路通路的mRNAs (GenBank登录号分别为XM_006567225.2, XM_016913157.1和XM_016911813.1)(图7: B)。NcCKvsNcT1中miR-8-y靶向AmCK1vsAmT1中的mRNA可注释到黑色素生成通路(图7: A)；NcCKvsNcT2中miR-8-y靶向AmCK2vsAmT2中的2个mRNA可分别注释到黑色素生成和Ras信号通路(图7: B)。NcCKvsNcT1中miR-4796-y靶向AmCK1vsAmT1中的mRNA可注释到MAPK信号通路(图7: A)；NcCKvsNcT2中miR-4796-y靶向AmCK2vsAmT2中的mRNA可注释到Ras信号通路(图7: B)。NcCKvsNcT1中bantam-y靶向AmCK1vsAmT1中的mRNA可注释到黑色素生成(图7: A)；NcCKvsNcT2中bantam-y同时靶向AmCK2vsAmT2中2条mRNA可分别注释到Ras信号通路和黑色素生成(图7: B)。NcCKvsNcT1中miR-5119-y靶向AmCK1vsAmT1中的mRNA可注释到黑色素生成(图7: A)；NcCKvsNcT2中miR-5119-y还同时靶向AmCK2vsAmT2中的2条mRNA可分别注释到黑色素生成和MAPK信号通路(图7: B)。

值得注意的是，笔者所在课题组前期研究发现，NcCKvsNcT2中bantam-y可靶向东方蜜蜂微孢子虫的己糖激酶编码基因mRNA (GenBank登录号: XM_002995838.1)，该基因可注释到东方蜜蜂微孢子虫能量代谢相关的糖酵解/糖异生途径；并且bantam-y和miR-5119-y可分别靶向东方蜜蜂微孢子虫自身的ABC转运蛋白编码基因mRNA (GenBank登录号分别为XM_002995835.1和XM_002996675.1)(耿四海等, 2020a)。据此猜测东方蜜蜂微孢子虫在侵染意大利蜜蜂工蜂中肠的过程中通过不同程度地差异表达miR-216-x和miR-8-y等miRNA对意大利蜜蜂工蜂中肠的黑色素生成、溶酶体和Ras信号通路进行调控，此外东方蜜蜂微孢子虫还能通过差异表达bantam-y和miR-5119-y等miRNA一方面调控意大利蜜蜂工蜂中肠的黑色素生成、Ras信号通路和MAPK信号通路，另一方面调节自身的物质和能量运输，从而促进增殖，体现了miRNA作为具有多面性的关键调控因子在东方蜜蜂微孢子虫和意大利蜜蜂工蜂中肠的互作中扮演着特殊角色。

3.4 东方蜜蜂微孢子虫在不同侵染时间点通过特异性差异表达部分miRNA对意大利蜜蜂工蜂中肠基因进行表达调控

NcCKvsNcT1中特有的miR-186-x靶向意大利蜜蜂工蜂中肠的mRNAs (GenBank登录号分别为XM_006559228.1和NM_001011611.2)可分别注释到黑色素生成通路和免疫系统进程条目，NcCKvsNcT2中特有的novel-m0022-5p靶向意大利蜜蜂工蜂可注释到代谢途径的mRNAs (GenBank登录号分别为XM_006559395.1, XM_006560586.2和XM_016911813.1)。我们推测东方蜜蜂微孢子虫在感染意大利蜜蜂工蜂7 d通过上调miR-186-x抑制意大利蜜蜂工蜂中肠的免疫防御相关基因的表达，在感染10 d时通过上调novel-m0022-5p对意大利蜜蜂工蜂中肠的能量代谢相关基因进行抑制，从而加强对宿主细胞内能量的掠夺，供其增殖所需。此外，NcCKvsNcT1中特有的miR-22-x靶向意大利蜜蜂工蜂可分别注释到生物附着和应激反应条目的mRNAs (GenBank登录号分别为XM_006563625.2和XM_006559419.2)，NcCKvsNcT2中特有的miR-26-x靶向意大利蜜蜂工蜂的mRNA (GenBank登录号: XM_006557238.2)可同时注释到Toll样受体信号通路、MAPK信号通路及生物进程调节功能条目。这些结果表明东方蜜蜂微孢子虫在感染意大利蜜蜂工蜂7 d和10 d可能通过分别下调miR-22-x和miR-26-x对上述重要的功能条目和通路相关基因表达进行抑制，从而影响宿主的正常生命活动。

综上所述，本研究结合前期获得的高质量miRNA和mRNA组学数据对东方蜜蜂微孢子虫DEmiRNA靶向的意大利蜜蜂工蜂中肠mRNA和DEmRNA进行预测、分析和探讨，揭示病原的DEmiRNA对宿主的基因表达调控具有广泛性，病原可能通过上调部分miRNA调控宿主的免疫防御以促进增殖，通过下调部分miRNA跨界调控宿主的能量代谢以加强能量窃取并促进增殖。研究结果为阐明东方蜜蜂微孢子虫跨界调控意大利蜜蜂工蜂中肠的分子机制提供了理论依据和数据基础，也为深入理解二者的互作机制提供了有价值的参考信息。