熊翠玲#, 杜宇#, 冯睿蓉, 蒋海宾, 史小玉, 王海朋, 范小雪, 王杰, 祝智威, 范元婵, 陈华枝, 周丁丁, 郑燕珍, 陈大福
, 郭睿 福建农林大学动物科学学院(蜂学学院), 福建 福州 350002
收稿日期：2019-09-10；修回日期：2019-11-06；网络出版日期：2020-01-20
基金项目：国家自然科学基金（31702190）；现代农业产业技术体系建设专项资金（CARS-44-KXJ7）；福建省教育厅中青年教师教育科研项目（JAT170158）；福建农林大学****科研人才计划（xjq201814）；福建农林大学科技创新专项基金（CXZX2017342，CXZX2017343）；福建省大学生创新创业训练计划（201910389011，201810389029，201810389082）
*通信作者：陈大福, E-mail:dfchen826@fafu.edu.cn;
郭睿, Tel/Fax:+86-591-83726835;E-mail:ruiguo@fafu.edu.cn.
#并列第一作者
摘要：[目的] 蜜蜂球囊菌（Ascosphaera apis，简称球囊菌）是专性侵染蜜蜂幼虫的致死性真菌病原。MicroRNA（miRNA）作为一类重要的基因表达调控因子，能够广泛参与真菌及其宿主的相互作用过程。本研究通过比较分析球囊菌孢子（AaCK）和侵染中华蜜蜂（Apis cerana cerana，简称中蜂）6日龄幼虫肠道内的球囊菌（AaT）的small RNA（sRNA）组学数据对球囊菌的差异表达miRNA（differentially expressed miRNA，DEmiRNA）、靶mRNA及二者间的调控网络进行全面解析，旨在揭示miRNA介导的球囊菌对中蜂幼虫的侵染机制。[方法] 对于球囊菌侵染的中蜂6日龄幼虫肠道的small RNA-seq（sRNA-seq）数据，利用BLAST工具连续比对东方蜜蜂（Apis cerana）和球囊菌的参考基因组筛滤得到AaT的sRNA组学数据。分别将AaCK和AaT的sRNA组学数据比对miRBase数据库，对球囊菌侵染宿主前后miRNA的数量和结构特征进行分析。联用RNAhybrid+svm_light、Miranda和TargetScan软件预测AaCK vs AaT比较组中DEmiRNA的靶mRNA，进而利用相关生物信息学软件对上述靶mRNA进行GO分类和KEGG代谢通路富集分析。通过Cytoscape软件对DEmiRNA-mRNA调控网络进行可视化。利用Stem-loop RT-PCR、RT-qPCR和分子克隆验证测序结果的可靠性。[结果] 在AaCK和AaT中分别鉴定到380和387个miRNA。结构特征分析结果显示，AaCK和AaT的miRNA皆集中分布在18-25 nt，且首位碱基主要偏向于U。AaCK vs AaT比较组共有270个DEmiRNA，包含155个上调miRNA和115个下调miRNA，分别靶向结合6091和6145个mRNA。GO分类结果显示，上述靶mRNA主要涉及代谢进程、细胞进程、应激反应等15个生物学进程；细胞、细胞组分、细胞器等12个细胞组分；催化活性、结合、转运子活性等11个分子功能。KEGG代谢通路富集分析结果显示，上述靶mRNA富集在123条代谢通路，参与对氨基酸代谢、碳水化合物代谢以及核苷酸代谢等物质代谢，氧化磷酸化、硫代谢、氮代谢等能量代谢，以及MAPK和Hippo等信号通路的调控。球囊菌DEmiRNA与靶mRNA之间存在复杂的调控关系，其中miR-29-x、miR-250-x、miR-4968-y、miR-11200-x、novel-m0023-5p、novel-m0130-5p和novel-m0135-5p等DEmiRNA可靶向结合与球囊菌的半胱氨酸蛋白酶、DNA甲基化转移酶以及几丁质酶相关的mRNA；此外，miR-7-x、miR-9-z、miR-319-y和miR-5951-y等同时参与调控MAPK信号通路；进一步分析发现，miR-250-x同时参与对DNA甲基化转移酶、MAPK信号通路及其他酶类合成与代谢途径的调控，并可能参与球囊菌与中蜂6日龄幼虫之间的跨界调控。通过Stem-loop RT-PCR和RT-qPCR验证了4个DEmiRNA的差异表达，并利用分子克隆和Sanger测序证实miR-7-x的序列与测序结果一致。[结论] 本研究解析了侵染中蜂6日龄幼虫的球囊菌的miRNA差异表达谱及DEmiRNA的调控网络，揭示了球囊菌DEmiRNA可能通过调控病原的物质和能量代谢、增殖、毒力、信号通路及相关mRNA参与对中蜂幼虫的侵染过程。miR-7-x、miR-250-x、novel-m0023-5p等关键DEmiRNA有望作为白垩病治疗的新型分子靶点。
关键词：蜜蜂球囊菌中华蜜蜂幼虫微小RNA调控网络侵染机制跨界调控
Differential expression pattern and regulation network of microRNAs in Ascosphaera apis invading Apis cerana cerana 6-day-old larvae
Cuiling Xiong#, Yu Du#, Ruirong Feng, Haibin Jiang, Xiaoyu Shi, Haipeng Wang, Xiaoxue Fan, Jie Wang, Zhiwei Zhu, Yuanchan Fan, Huazhi Chen, Dingding Zhou, Yanzhen Zheng, Dafu Chen
, Rui Guo College of Animal Sciences (College of Bee Science), Fujian Agriculture and Forestry University, Fuzhou 350002, Fujian Province, China
Received: 10 September 2019; Revised: 6 November 2019; Published online: 20 January 2020
*Corresponding author: Dafu Chen, dfchen826@fafu.edu.cn;
Tel/Fax: +86-591-83726835; E-mail: Rui Guo, ruiguo@fafu.edu.cn.
Foundation item: Supported by the National Natural Science Foundation of China (31702190), by the Earmarked Fund for China Agriculture Research System (CARS-44-KXJ7), by the Teaching and Scientific Research Fund of Education Department of Fujian Province (JAT170158), by the Outstanding Scientific Research Manpower Fund of Fujian Agriculture and Forestry University (xjq201814), by the Scientific and Technical Innovation Fund of Fujian Agriculture and Forestry University (CXZX2017342, CXZX2017343) and by the Undergraduate Innovation and Entrepreneurship Training Program of Fujian Province (201910389011, 201810389029, 201810389082)
#These authors contributed equally to this work.
Abstract: [Objective] This study aimed to reveal miRNA-mediated mechanism underlying Ascosphaera apis infection of Apis cerana cerana larvae. [Methods] Small RNA (sRNA) dataset of A. apis during infection (AaT) was screened out from sRNA-seq data from Ascosphaera apis-infected A. c. cerana 6-day-old larval guts. The filtered sRNA datasets from the purified spores (AaCK) and AaT were aligned against miRBase using Blast, followed by analyses of number and structural characteristics of pathogen miRNAs before and after Ascosphaera apis infection. Prediction, GO categorization and KEGG pathway enrichment analysis of targets of DEmiRNAs were conducted using related software. The regulation network between DEmiRNAs and corresponding targets was visualized using Cytoscape. Stem-loop RT-PCR, qPCR and molecular cloning were used to verify the reliability of our sequencing data. [Results] Totally, 380 and 387 miRNAs were identified in AaCK and AaT, respectively. The length of Ascosphaera apis miRNAs were mainly distributed between 18 nt and 25 nt; and the first base had a U bias. There were 155 up-regulated and 115 down-regulated miRNAs in AaCK vs AaT, targeting 6091 and 6145 mRNAs. Targets of DEmiRNAs were involved in 15 biological processes, 12 cell components and 11 molecular functions. Additionally, these targets were engaged in 123 pathways, regulating material metabolisms, energy metabolisms and signaling pathways. Moreover, complex regulation networks existed between DEmiRNA and corresponding targets, among them miR-29-x, miR-250-x, miR-4968-y, miR-11200-x, novel-m0023-5p, novel-m0130-5p and novel-m0135-5p can target mRNAs associated with cysteine proteinase, DNA methyltransferases and chitinase; miR-7-x, miR-9-z, miR-319-y and miR-5951-y can simultaneously regulate MAPK signaling pathway; miR-250-x may be involved in cross-kingdom regulation between A. apis and A. c. cerana larvae. [Conclusion] Our results revealed DEmiRNAs may participate in the infection process of A. apis via regulating targets associated with material and energy metabolisms, pathogen proliferation, virulence, and several signaling pathways; several key miRNAs including miR-7-x were potential targets for chalkbrood control.
Keywords: Ascosphaera apisApis cerana ceranalarvaemicroRNAregulation networkinfection mechanismcross-kingdom regulation
蜜蜂球囊菌(Ascosphaera apis，简称球囊菌)是导致蜜蜂幼虫罹患白垩病的致死性真菌病原。不同于其他昆虫病原菌主要依赖分泌几丁质酶等胞外酶从外界突破宿主表皮和肠腔，球囊菌与蜜蜂幼虫的长期协同进化过程中形成了高度适应蜜蜂幼虫肠道发育的特定侵染模式[1]。病原孢子通过内勤蜂的哺育行为经口饲喂进入幼虫中肠，在CO2的刺激下低水平萌发，并伴有菌丝的少量生长，6日龄末期随着中肠和后肠的隔膜消失，孢子随食物残渣涌入后肠，在O2的刺激下迅速萌发，此时菌丝大量生长，到了宿主免疫系统防御能力较低的预蛹期，菌丝先后穿透幼虫肠壁和体壁，导致幼虫死亡，最终从幼虫的尾部蔓延包裹整个幼虫，形成白垩病状虫尸[2]。在养蜂生产中，一般认为该病仅危害西方蜜蜂(Apis mellifera)。直至2015年，Maxfield-Taylor等[3]发现球囊菌能够侵染成年熊蜂蜂王。随后，笔者团队通过组织学和分子生物学手段证实球囊菌对中华蜜蜂(Apis cerana cerana，简称中蜂)的雄蜂幼虫具有侵染性，并通过交叉感染试验证实球囊菌对中蜂工蜂幼虫也具有侵染性[4]。前期研究中，为探究球囊菌的侵染机制，笔者团队开展了一系列的分子生物学和组学研究，基于高质量的测序数据de novo组装并注释了球囊菌的参考转录组，并大规模挖掘了球囊菌的SSR分子标记[5]；在mRNA组学层面解析了球囊菌在胁迫意大利蜜蜂(Apis mellifera ligustica，简称意蜂)幼虫和中华蜜蜂幼虫过程中的差异基因表达谱和转录组动态，初步揭示了病原差异表达基因(DEG)的功能[6-7]；此外，还在全基因组水平预测、分析和鉴定了球囊菌的微小RNA (microRNA，miRNA)[8]。相关研究为深入解析球囊菌对蜜蜂幼虫的侵染机制打下了良好基础。
MiRNA作为一类关键的基因表达调节因子，其种子序列在转录后水平能够特异性结合靶mRNA的3′ UTR，导致靶mRNA降解或抑制其翻译过程，从而调控生物体的各类生命活动[9]。2010年，Lee等[10]首次在粗糙脉胞菌(Neurospora crassa)中发现丝状真菌的miRNA-like RNA (milRNA)。鉴于真菌milRNA与动植物成熟miRNA的特征及对靶mRNA的调控方式高度相似[11]，以下统称真菌miRNA。此后，随着高通量测序技术的迅速发展与革新，人们在尖孢镰刀菌(Fusarium oxysporum)[12]、核盘菌(Sclerotinia sclerotiorum)[13]、金龟子绿僵菌(Metarhizium anisopliae)[14]、东方蜜蜂微孢子虫(Nosema ceranae)[15]等真菌中预测和鉴定出大量miRNA。李琼等[16]研究发现mro-miR-33可通过负向调节brlA促进罗伯茨绿僵菌(Metarhizium robertsii)的孢子形成。大量的研究表明，miRNA广泛参与真菌及其宿主的相互作用过程[17]。Shakeel等[18]研究发现小菜蛾(Plutella xylostella)接种M. anisopliae来源的僵菌素A后，miR-263、miR-279和miR-306等差异表达miRNA (differentially expressed miRNAs，DEmiRNA)参与调节Toll、IMD、Jak-STAT及细胞粘附分子等免疫相关信号通路。Evans等[19]研究推测N. ceranae来源的miRNA可分泌到蜜蜂的细胞质内，从而跨界调控宿主的新陈代谢与免疫过程。笔者团队曾利用二代测序技术和生物信息学方法在球囊菌菌丝和孢子的混合样品中预测和鉴定到118个miRNA，进而通过靶mRNA的功能分析及miRNA的调控网络分析揭示了球囊菌miRNA与致病性的潜在关联[8]。近期，笔者利用small RNA-seq (sRNA-seq)技术对正常及球囊菌胁迫的中蜂6日龄幼虫肠道进行了深度测序，获得了高质量的sRNA组学数据，并通过生物信息学和分子生物学方法深入分析了幼虫肠道DEmiRNA及其调控网络，研究结果表明球囊菌胁迫可导致宿主miR-30-x、miR-6052-x和miR-1277-x等54个miRNA差异表达，上述DEmiRNA可能参与调控宿主泛素介导的蛋白水解、细胞凋亡和Jak-STAT信号通路等免疫途径[20]。然而，对于病原miRNA在球囊菌-蜜蜂幼虫互作过程中的动态和功能，仍缺乏深入研究和相关报道。
因此，为进一步探究球囊菌对中蜂幼虫的侵染机制，本研究基于已获得的高质量sRNA组学数据，通过连续比对东方蜜蜂(Apis cerana)和球囊菌的参考基因组筛滤得到处理组球囊菌(AaT)的数据，结合已获得的对照组球囊菌孢子(AaCK)数据全面解析球囊菌侵染前后miRNA的数量和结构特征，DEmiRNA的靶mRNA功能注释，以及DEmiRNA-mRNA调控网络，并利用分子生物学技术对DEmiRNA进行验证。研究结果可为阐明球囊菌对中蜂幼虫的侵染机制及二者的互作机制提供新的见解和数据基础。
1 材料和方法 1.1 供试球囊菌及蜜蜂幼虫 球囊菌菌株由福建农林大学动物科学学院(蜂学学院)蜜蜂保护实验室分离和保藏；中蜂幼虫取自福建农林大学动物科学学院(蜂学学院)教学蜂场。
1.2 对照组和处理组球囊菌的sRNA-seq数据来源 球囊菌和中蜂幼虫肠道样本的制备参照笔者团队已建立的方法进行[20-21]。简述如下：将球囊菌菌株接种到PDA培养基，并置于33.0±0.5 ℃条件下培养10 d后，刮取黑色孢子囊，每100 mg置于1个RNA-free的EP管内，液氮速冻后置于–80 ℃冷藏。其中一份球囊菌孢子样品用于sRNA-seq，作为对照组。其余球囊菌利用血球计数板计数，并配置终浓度为1×107个孢子/mL的幼虫饲料，饲喂接种3日龄中蜂幼虫(50 μL/只)，此后每24 h更换正常饲料。实验幼虫置于35.0±0.5 ℃，90%相对湿度(RH)条件下饲养至6日龄，在超净工作台解剖肠道，每9只肠道放于1个RNA-Free的EP管中，液氮速冻后置于–80 ℃冷藏。上述接种感染实验进行3次生物学重复，制备的3份球囊菌感染幼虫肠道样品用于sRNA-seq。上述球囊菌孢子和球囊菌侵染的中蜂6日龄幼虫肠道样品委托广州基迪奥生物科技有限公司进行单端测序，测序平台为Illumina MiSeq。
对于球囊菌感染中蜂幼虫肠道测序得到的混合数据，首先将质控后的非注释标签序列映射到东方蜜蜂参考基因组(assembly ACSNU-2.0)，过滤映射上的序列标签(mapped tags) (即宿主数据)，并将unmapped tags比对到球囊菌参考基因组(assembly AAP 1.0)，由此得到的mapped tags即为病原本身的数据，从而得到球囊菌相关tags在基因组上的位置信息，获得球囊菌sRNA组学数据作为处理组数据(AaT-1、AaT-2、AaT-3)。球囊菌孢子测序得到的数据作为对照组数据(AaCK)。
根据笔者团队前期建立的数据质控方法对原始测序数据进行质控[22]。简述如下：首先过滤掉下机数据中质量值低于20的碱基数超过1个或含有N的低质量reads，剔除不带有3′接头的reads，并截取3′接头前面的序列，用作后续分析；进一步过滤掉含有5′接头的reads、不含插入片段的reads以及插入片段小于18 nt的reads、含有poly A的reads，最终得到小RNA的clean tags。
1.3 球囊菌miRNA的生物信息学预测及分析 目前，miRBase数据库内暂无球囊菌miRNA的相关信息，笔者参照郭睿等[8]的方法，利用miRDeep2软件将上述tags与miRBase数据库中其他物种已知miRNA的前体序列进行比对，并过滤掉比对位置为该miRNA左右各延伸2 bp的位置以外的tag，将最终比对得到的miRNA定义为known miRNA，其后缀“-x”、“-y”分别代表从miRNA前体的5′臂或3′臂加工而来，“-z”代表来源位置不确定。进一步过滤掉known miRNA的tag，并将剩余tag序列比对到球囊菌基因组位置，利用MIREAP_v0.2软件预测miRNA的二级结构，并将存在茎环结构、前体自由能大于–18 kcal/mol、且长度介于18–26 nt的成熟miRNA定义为novel miRNA。利用TPM (tags per million)算法对获得的所有miRNA的表达量进行归一化。
通过GraphPad Prism 7软件统计sRNA测序数据得到的包含rRNA、tRNA、miRNA、snRNA、snoRNA等大量的非编码RNA (non-coding RNA)片段数据，并分析AaCK和AaT组miRNA的长度分布和首位碱基偏向性。
1.4 DEmiRNA及其靶mRNA的预测和功能注释 AaCK vs AaT比较组的差异表达miRNA (DEmiRNA)的筛选标准为|log2 (Fold change)|≥1且P≤0.05。选取–10≤log2 (Fold change)≤10范围的DEmiRNA，利用在线heatmap工具(www.omicshare.com)进行表达量聚类分析。
分别利用RNAhybrid (v2.1.2)+svm_light (v6.01)、Miranda (v3.3a)、TargetScan (Version: 7.0)三种软件对球囊菌DEmiRNA的靶mRNA进行生物信息学预测，采用各软件的默认参数；将三者预测结果的交集作为可靠的靶mRNA集合。进一步筛选得到结合自由能小于–17 kcal/mol的靶向结合关系用以后续分析[23-24]。利用BLAST工具将预测得到的靶mRNA序列比对到GO和KEGG数据库，从而获得靶mRNA的功能注释信息。
1.5 DEmiRNA-mRNA调控网络的构建及分析 根据笔者团队的前期研究及上述靶mRNA的功能注释结果，筛选可能与球囊菌增殖、毒力和抗逆性等相关的靶mRNA，根据DEmiRNA与靶mRNA的靶向结合关系，利用Cytoscape软件构建和呈现DEmiRNA与靶mRNA的调控网络。
1.6 DEmiRNA的Stem-loop RT-PCR、Sanger测序及RT-qPCR验证 参照前期建立的验证体系[20]，随机挑选的4个DEmiRNA (miR-7-x、miR-9-z、miR-319-y和miR-5951-y)进行分子验证。根据相应的核酸序列，利用DNAMAN软件设计特异性Stem-loop引物和上游引物，以及通用下游引物(表 1)。利用RNA抽提试剂盒(TaKaRa，日本)分别抽取AaCK和AaT的总RNA，利用Stem-loop引物反转录得到cDNA模板进行PCR扩增，PCR产物经1.5%琼脂糖凝胶电泳检测。胶回收100 bp附近的片段，连接pMD-19T载体，转化DH5α大肠杆菌后挑斑摇菌，取菌液PCR呈阳性的菌液送生工生物工程(上海)股份有限公司单端测序。
表 1. 本研究使用的引物 Table 1. Primers used in this study
| Primer ID | Primer sequence (5′→3′) |
| miR-7-x-Loop | CTCAACTGGTGTCGTGGAGTCGGCAATTCAGTTGAGAACAACAA |
| miR-9-z-Loop | CTCAACTGGTGTCGTGGAGTCGGCAATTCAGTTGAGTCATACAG |
| miR-319-y-Loop | CTCAACTGGTGTCGTGGAGTCGGCAATTCAGTTGAGGGGAGCTC |
| miR-5951-y-Loop | CTCAACTGGTGTCGTGGAGTCGGCAATTCAGTTGAGCCGTCGTC |
| miR-7-x-F | ACACTCCAGCTGGGTGGAAGACTAGTGATT |
| miR-9-z-F | ACACTCCAGCTGGGTCTTTGGTTATCTAG |
| miR-319-y-F | ACACTCCAGCTGGGTTGGACTGAAGG |
| miR-5951-y-F | ACACTCCAGCTGGGTTTCGTCAGGAC |
| Universal R | CTCAACTGGTGTCGTGGA |
| 5.8S rRNA-F | GCAGCGAAATGCGATAAGTAA |
| 5.8S rRNA-R | CCCTCCTAAGACGGGACGAT |
表选项
利用ABI QuantStudio 3荧光定量PCR系统进行上述4个DEmiRNA的RT-qPCR检测。采用SYBR Green法，反应条件如下：95 ℃ 5 min；95 ℃ 10 s，59 ℃ 30 s，40个循环；熔解曲线程序默认系统设置。以球囊菌5.8S rRNA基因作为内参(GenBank：U68313.1)，利用2–△△Ct方法对上述4个DEmiRNA进行相对定量。利用GraphPad Prism 7软件处理3次生物学重复的数据并绘图，并利用t检验计算组间显著性。
2 结果和分析 2.1 数据质控与评估 AaCK的sRNA-seq产生12708832 raw reads，AaT测序得到的平均raw reads数为10773889；经质控后AaCK和AaT的clean reads数分别为12544551和10399091，占raw reads的比例均在96.09%以上(表 2)，说明测序数据质量良好。将AaCK和AaT的clean reads映射球囊菌基因组，mapping ratio分别为71.60%和6.43%。
表 2. sRNA-seq数据信息统计 Table 2. Summary of sRNA-seq datasets
| Samples | Raw reads | Clean reads | Mapping ratio/% |
| AaCK | 12708832 | 12544551 (98.71%) | 71.60 |
| AaT-1 | 10252283 | 9932883 (96.88%) | 8.02 |
| AaT-2 | 11367190 | 10922645 (96.09%) | 7.58 |
| AaT-3 | 10702194 | 10341746 (96.63%) | 3.69 |
表选项
2.2 球囊菌miRNA的预测与结构特征分析 通过深度测序得到较为复杂sRNA数据集，其中包含rRNA、tRNA、miRNA、snRNA、snoRNA等大量的非编码RNA (non-coding RNA)片段(图 1-A)。通过比对miRBase数据库，AaCK和AaT共鉴定到648个miRNA。其中，在AaCK中鉴定到380个miRNA，包含233个known miRNA和147个novel miRNA；在AaT中鉴定到387个miRNA，包含330个known miRNA和57个novel miRNA。结构特征分析结果显示，AaCK和AaT的miRNA皆集中分布在18–25 nt，miRNA分布最多的分别为18 nt和22 nt (图 1-B)；此外，AaCK的miRNA首位碱基主要偏向于U，其次偏向于G和C；AaT的miRNA首位碱基同样主要偏向于U，但其次偏向于A和C (图 1-C)。
 |
| 图 1 球囊菌miRNA的预测及特征分析 Figure 1 Prediction and characterization of A. apis miRNA. A: Composition of various sRNA populations generated from AaCK and AaT; B: Length distribution of A. apis miRNAs; C: First nucleotide bias of A. apis miRNAs |
| 图选项 |
2.3 球囊菌miRNA的差异表达分析 AaCK vs AaT比较组经筛选共得到270个DEmiRNA，上调和下调表达的miRNA数分别为155和115个(图 2-A)。其中miR-31-x、miR-9-z、miR-1-y以及miR-375-y等显著上调，log2(Fold change)分别达到6.28、4.80、4.72和3.96；novel-m0072-3p、novel-m0003-3p、novel-m0046-5p和miR-4968-y等显著下调，log2(Fold change)分别达到–9.49、–9.60、–9.88和–9.92 (图 2-B)。
 |
| 图 2 球囊菌miRNA的差异表达(A)及表达量聚类(B) Figure 2 Differential expression (A) and expression clustering (B) of A. apis miRNAs |
| 图选项 |
2.4 球囊菌DEmiRNA的靶mRNA预测及分析 靶mRNA预测结果显示，上述270个DEmiRNA共结合6348个靶mRNA，其中上调miRNA和下调miRNA分别靶向结合6091和6145个mRNA。GO分类结果显示，上述靶mRNA主要涉及代谢进程(1328)、细胞进程(1309)、单组织进程(1028)、定位(355)、生物调节(310)、细胞成分组织或生物发生(294)、应激反应(169)等15个生物学进程；细胞(1134)、细胞组分(1134)、细胞器(770)、大分子复合物(436)、细胞器组件(356)、细胞膜(352)、细胞膜组件(238)、膜关闭内腔(76)等12个细胞组分；催化活性(1264)、结合(920)、转运子活性(129)等11个分子功能(表 3)。
表 3. 球囊菌DEmiRNA的靶mRNA富集数前20位的GO条目 Table 3. Top 20 GO terms enriched by target mRNAs of A. apis DEmiRNAs
| GO category | GO term | NR description | Number of target mRNAs |
| BP | GO:0009987 | Cellular process | 1328 |
| BP | GO:0008152 | Metabolic process | 1309 |
| MF | GO:0003824 | Catalytic activity | 1264 |
| CC | GO:0005623 | Cell | 1134 |
| CC | GO:0044464 | Cell part | 1134 |
| BP | GO:0044699 | Single-organism process | 1028 |
| MF | GO:0005488 | Binding | 920 |
| CC | GO:0043226 | Organelle | 770 |
| CC | GO:0032991 | Macromolecular complex | 436 |
| CC | GO:0044422 | Organelle part | 356 |
| BP | GO:0051179 | Localization | 355 |
| CC | GO:0016020 | Membrane | 352 |
| BP | GO:0065007 | Biological regulation | 310 |
| BP | GO:0071840 | Cellular component organization or biogenesis | 294 |
| CC | GO:0044425 | Membrane part | 238 |
| BP | GO:0050896 | Response to stimulus | 169 |
| MF | GO:0005215 | Transporter activity | 129 |
| CC | GO:0031974 | Membrane-enclosed lumen | 76 |
| BP | GO:0023052 | Signaling | 52 |
| MF | GO:0005198 | Structural molecule activity | 44 |
| BP: Biological process; MF: Molecular function; CC: Cellular component. | |||
表选项
进一步对上述靶mRNA进行KEGG代谢通路富集分析，发现它们富集到123条KEGG代谢通路，涉及新陈代谢(2872)、遗传信息加工(854)、细胞进程(385)、环境信息加工(57)以及有机系统(19)等5大类，富集mRNA数最多的前10位分别是新陈代谢总览(708)，翻译(295)，折叠、分类与降解(237)，碳水化合物代谢(235)，氨基酸代谢(207)，转运和分解代谢(193)，转录(138)，脂质代谢(130)，能量代谢(115)，核苷酸代谢(110) (图 3)。
 |
| 图 3 球囊菌DEmiRNA的靶mRNA的KEGG代谢通路富集分析 Figure 3 KEGG pathway enrichment analysis of A. apis DEmiRNA target mRNAs |
| 图选项 |
进一步分析发现富集通路还包括氧化磷酸化(64)、硫代谢(17)、氮代谢(9)等能量代谢相关通路；糖酵解/糖异生(46)、氨基糖和核苷酸糖代谢(40)、丙酮酸代谢(31)等碳水化合物代谢相关通路；甘氨酸、丝氨酸和苏氨酸代谢(37)、缬氨酸、亮氨酸和异亮氨酸降解(30)、丙氨酸、天冬氨酸和谷氨酸代谢(28)等氨基酸代谢通路，以及MAPK信号通路(29)、磷脂酰肌醇信号系统(9)、Hippo信号通路(3)等信号传导相关通路。
2.5 球囊菌DEmiRNA的调控网络构建及分析 鉴于球囊菌DEmiRNA数量达到270个，构建出的调控网络极为复杂，因此根据Nr数据库注释信息，筛选与半胱氨酸蛋白酶、DNA甲基化转移酶、几丁质酶以及MAPK信号通路相关的靶mRNA，并构建它们与靶向结合的球囊菌DEmiRNA之间的调控网络(图 4)，分析发现分别有11个上调miRNA和18个下调miRNA结合2个与半胱氨酸蛋白酶相关的靶mRNA；分别有6个上调miRNA和7个下调miRNA结合1个与DNA甲基化转移酶相关的靶mRNA；分别有3个上调miRNA和6个下调miRNA结合1个与几丁质酶相关的靶mRNA；分别有114个上调miRNA和88个下调miRNA结合29个MAPK信号通路相关的靶mRNA。
 |
| 图 4 球囊菌的DEmiRNA-mRNA调控网络 Figure 4 DEmiRNA-mRNA regulation networks in A. apis. A: Regulation networks of cysteine proteinase-associated mRNAs and corresponding DEmiRNAs; B: Regulation networks of chitinase-associated mRNAs and corresponding DEmiRNAs; C: Regulation networks of DNA methyltransferases-associated mRNAs and corresponding DEmiRNAs; D: Regulation networks of MAPK signaling pathway-associated mRNAs and corresponding DEmiRNAs. Orange circles indicate up-regulated miRNAs; green circles indicate down-regulated miRNAs. Purple diamond indicate target mRNAs |
| 图选项 |
进一步分析结果显示，显著下调的miR-4968-y同时参与对DNA甲基化转移酶、半胱氨酸蛋白酶和MAPK信号通路的调控；显著下调的miR-11200-x和novel-m0130-5p均同时参与对DNA甲基化转移酶、几丁质酶和MAPK信号通路的调控；显著下调的novel-m0023-5p和novel-m0135-5p均同时参与对几丁质酶、半胱氨酸蛋白酶和MAPK信号通路的调控；显著上调的miR-6037-y、miR-250-x、miR-8503-x、miR-3786-x、miR-1344-x、miR-3718-y等6个miRNA和显著下调novel-m0123-3p、novel-m0028-5p、miR-29-x、miR-5658-x等4个显著miRNA均同时参与对DNA甲基化转移酶和MAPK信号通路的调控。显著上调的miR-2-x、miR-10-y、miR-100-y、miR-190-x、miR-971-y、miR-2788-y、miR-3715-x、miR-6005-x、miR-6497-x、miR-29b3-y等10个miRNA和显著下调的miR-29-y、miR-1546-x、miR-5951-y、miR-7977-x、miR-8577-x、novel-m0019-3p、novel-m0060-3p、novel-m0062-3p、novel-m0077-3p、novel-m0108-5p、novel-m0112-3p等11个miRNA均同时参与对半胱氨酸蛋白酶和MAPK信号通路的调控。
2.6 球囊菌DEmiRNA的Stem-loop RT-PCR和RT-qPCR验证 利用Stem-loop RT-PCR和RT-qPCR对随机选择的2个上调miRNA和2个下调miRNA在AaCK和AaT中的真实表达与相对表达变化趋势进行验证，结果显示上述4个DEmiRNA在AaCK vs AaT比较组内均真实表达，且与sRNA-seq中相应的表达量变化趋势一致(图 5 A–D)。对miR-7-x的Stem-loop RT-PCR扩增片段进行TA克隆和Sanger测序，结果与其测序序列一致(图 5-E)，进一步证明了测序结果的可靠性。
 |
| 图 5 球囊菌DEmiRNAs的Stem-loop RT-PCR及RT-qPCR验证 Figure 5 Stem-loop RT-PCR and RT-qPCR confirmation of A. apis DEmiRNAs. A, B, C, D: Agarose gel electrophoresis of Stem-loop RT-PCR products and RT-qPCR examination of miR-319-y, miR-5951-y, miR-9-z and miR-7-x; E: Sanger sequencing of Stem-loop RT-PCR product from miR-7-x. Three biological replicates were analyzed, and are represented as Mean±SD. *: P < 0.05 compared with AaCK (t-test) |
| 图选项 |
3 讨论 球囊菌专性侵染蜜蜂幼虫，在与宿主的长期协同进化过程中形成特定的侵染模式，但目前对球囊菌的致病机理的认识还停留在组织学和生物化学水平[1]，缺乏分子生物学和组学水平的深入阐释。前期研究中，笔者团队利用二代测序技术和生物信息学方法在mRNA组学层面初步解析了球囊菌与意蜂及中蜂幼虫的复杂互作机制[6-7, 25-26]，发现在侵染中蜂幼虫肠道过程中，部分与物质和能量代谢相关的球囊菌基因表达受到不同程度的抑制[4]；中蜂幼虫肠道的免疫相关基因被球囊菌大量激活[25]；此外，在临近白垩病暴发的中蜂6日龄幼虫肠道鉴定到537个DEmiRNA，它们通过靶向调控各类代谢活动和免疫途径对球囊菌胁迫产生应答[20]。本研究在前期研究基础上进一步探究侵染中蜂6日龄幼虫肠道的球囊菌的miRNA差异表达谱及其潜在调控作用，在AaCK和AaT中分别鉴定出380和387个miRNA，它们的长度主要介于18–25 nt，首位碱基偏向性主要偏向于U。前人研究发现，番茄在接种厚垣孢普可尼亚菌(Pochonia chlamydosporia)后，其体内用于调控基因沉默的24 nt的sRNA丰度显著降低[27]。本研究发现AaCK中miRNA在18 nt分布的数量最多，而AaT中miRNA在22 nt分布的数量最多，这种差异体现出球囊菌侵染中蜂幼虫前后，其miRNA的结构特征发生改变，暗示球囊菌可能通过改变miRNA结构参与对中蜂幼虫的致病过程。本研究鉴定到的miRNA与前期研究[8]鉴定到的miRNA互为补充，丰富了球囊菌miRNA的信息，为其功能研究提供了数据基础。进一步的差异分析发现AaCK vs AaT比较组的上调和下调miRNA数分别为155和115个，分别靶向结合6091和6145个mRNA，这些靶mRNA参与调控代谢进程、细胞进程、催化活性、应激反应等功能条目，核苷酸代谢、氨基酸代谢和碳水化合物代谢等物质代谢途径，氧化磷酸化、硫代谢和氮代谢等能量代谢途径，以及MAPK信号通路、Hippo信号通路和磷脂酰肌醇信号系统等信号转导通路，上述结果表明球囊菌在侵染中蜂幼虫时通过改变自身miRNA的表达量，以影响靶mRNA的表达水平，进而适应幼虫肠道内的环境并调节自身的增殖。
DNA甲基化是真核生物的重要表观遗传标记，通过影响基因的表达调控遗传组分与环境胁迫的应答反应[28]。DNA甲基化转移酶(DNMTase)是甲基化的关键酶，与真菌病原的生长、抗逆性及毒力密切相关[29]。Jeon等[30]研究发现缺失DNMTase编码基因DIM-2的稻瘟病菌(Magnaporthe oryzae)，其菌落形态、径向生长和孢子形成过程发生异常，导致其致病力明显下降。Wang等[31]发现MrDIM-2和MrRID可影响M. robertsii的甲基化过程和基因组甲基化特异性，缺失MrDIM-2或RID/DIM-2的菌株，其富集在细胞进程调控和胁迫应答等功能条目上基因及菌株毒力相关的半胱氨酸蛋白酶和几丁质酶等基因均下调表达，导致M. robertsii在经过辐照和热激等不良环境胁迫后的孢子活性降低；此外，作者还发现缺失MrDIM-2或RID/DIM-2的M. obertsii感染的大蜡螟(Galleria mellonella)死亡率和体表分生孢子数也明显降低。Lewsey等[32]研究发现拟南芥(Arabidopsis thaliana)的siRNA在植株体内传播过程中可介导全基因组DNA甲基化(RdDM)。目前，相关miRNA主要通过靶向结合DNMTase mRNA的3′UTR并抑制其表达，进而影响整个机体的甲基化模式和生物学功能[33-34]。前期研究中，我们发现侵染中蜂6日龄肠道的球囊菌的DNMTase编码基因(ID：Unigene0012984)显著上调[log2(Fold Change)=2.52，P=0.02]。本研究中，分别有6个上调miRNA (miR-250-x、miR-1344-x、miR-3718-y、miR-3786-x、miR-6037-y、miR-8503-x)和7个下调miRNA (miR-29-x、miR-4968-y、miR-5658-x、miR-11200-x、novel-m0028-5p、novel-m0123-3p、novel-m0130-5p)共同靶向1个与DNA甲基化转移酶相关的mRNA (DMAP1，ID：KZZ94679.1)。推测侵染中蜂6日龄幼虫的球囊菌通过下调上述7个miRNA的表达水平，部分解除对DMAP1的抑制作用，从而维持正常的甲基化水平，保证与其抗逆性和致病性相关基因的表达。中蜂幼虫与球囊菌存在复杂互作，宿主的免疫系统必然产生相应的免疫防御，以限制病原的增殖造成部分相关miRNA上调表达。目前，笔者已对球囊菌菌丝和孢子进行全基因组DNA甲基化测序，未来可结合本研究中的miRNA组学数据，进一步探究球囊菌侵染中蜂幼虫过程中病原miRNA与DNA甲基化之间的调控关系。
在长期协同进化过程中，以蜜蜂为代表的昆虫进化出以吞噬、包囊和黑化作用为代表的细胞免疫和以Jak-STAT、Imd、JNK和Toll信号通路为代表的体液免疫[35]，与复杂而稳定的肠道菌群[36]和不饱和脂肪酸[2]等生物化学环境协同抵御外源病原菌的入侵。对于真菌病原，其致病性主要依赖于利用宿主营养满足自身繁殖，另一方面抑制和削弱宿主先天性免疫的能力，例如分泌以半胱氨酸蛋白酶和几丁质酶等蛋白酶和其他次生代谢产物抑制宿主免疫反应并降解宿主防御分子[37]。前期研究中，笔者团队发现球囊菌的12个几丁质酶相关基因在胁迫意蜂6日龄幼虫肠道过程中差异表达，说明这些差异表达基因参与了球囊菌的致病过程[38]。本研究发现，有3个上调miRNA (bantam-y、miR-81-y、miR-282-y)和6个下调miRNA (miR-11200-x、novel-m0023-5p、novel-m0076-5p、novel-m0096-5p、novel-m0130-5p和novel-m0135-5p)靶向结合1个几丁质酶相关mRNA (ID：KZZ95064.1)；此外，miR-2-x、miR-10-y、miR-100-y、miR-200-x和miR-971-y等11个上调miRNA以及miR-1546-x、miR-4968-y、miR-8577-x、novel-m0023-5p和novel-m0135-5p等18个下调miRNA靶向结合2个半胱氨酸蛋白酶相关mRNA (ID：KZZ90947.1；KZZ98031.1)；进一步分析发现novel-m0023-5p和novel-m0135-5p均显著下调且同时参与对几丁质酶和半胱氨酸蛋白酶的调控。上述结果表明球囊菌在中蜂幼虫肠道的特定环境下，通过调节miRNA的表达水平，动态调控几丁质酶和半胱氨酸蛋白酶基因的表达，影响自身毒力并协助穿透幼虫肠道围食膜和肠壁，可能在球囊菌侵染后期引起白垩病暴发的过程中发挥关键作用，但需要进一步的实验验证。
病原真菌受到外界环境刺激后，可通过复杂而保守的信号级联反应引起下游靶标基因的表达变化[39]。其中，丝裂原活化蛋白激酶(mitogen- activated protein kinase，MAPK)级联反应在真菌信号通路中极为关键，参与到真菌交配繁殖、渗透压调节、孢子形成以及毒力水平的调控[40]，其依赖的MAP激酶可协助病原真菌侵染宿主[41]。笔者团队前期在mRNA组学层面探究了球囊菌在胁迫中蜂幼虫和意蜂幼虫过程中的基因表达谱，发现对于侵染中蜂6日龄幼虫的球囊菌，富集在MAPK信号通路上的高表达基因数量显著增加[42]，另有64个富集基因在球囊菌侵染意蜂6日龄幼虫肠道后显著上调[38]，暗示球囊菌MAPK信号通路上的大量基因被显著激活。本研究中，共有29个MAPK信号通路相关的mRNA (例如MAP激酶编码基因，ID：KZZ97981.1；真菌信息素交配因子编码基因，ID：KZZ87749.1，丝氨酸/苏氨酸蛋白激酶编码基因，ID：KZZ87358.1等)与114个上调miRNA (bantam-y、miR-2-x、miR-7-x、miR-9-z、miR-81-y、miR-100-y和miR-250-x等)和88个下调miRNA (miR-29-x、miR-319-y、miR-4968-y、miR-5951-y、miR-11200-x、novel-m0023-5p和novel-m0135-5p等)存在靶向结合关系。RT-qPCR结果显示其中的miR-7-x、miR-9-z、miR-319-y和miR-5951-y的相对表达变化水平与测序数据一致。以上结果表明上述DEmiRNA参与了对侵染中蜂幼虫的球囊菌MAPK信号通路的调控，病原可通过激活或抑制该通路富集靶mRNA的表达，影响球囊菌对宿主肠道环境的适应和增殖。
球囊菌DEmiRNA与靶mRNA之间存在复杂的调控关系。一个miRNA可调控多个靶mRNA，反之亦然[8]。本研究发现，结合262个靶mRNA的miR-250-x显著上调[log2(Fold Change)=11.52，P=0.01]，同时参与DNA甲基化转移酶、MAPK信号通路以及其他酶类合成与代谢途径的调控，推测球囊菌在侵染中蜂6日龄幼虫的过程中，宿主通过与球囊菌的互作上调后者miR-250-x的表达水平，增强对相关靶mRNA的抑制，从而影响病原的孢子萌发和菌丝生长。前人研究发现miRNA在不同物种间具有一定的同源性和保守性[43]，并参与宿主与病原相互之间的跨界调控[19, 44-45]。前期研究发现，被球囊菌侵染的中蜂6日龄幼虫肠道的miR-250显著上调表达，并参与调控泛素介导的蛋白水解和Jak-STAT信号通路等细胞和体液免疫。结合本研究的结果进一步分析发现，中蜂幼虫肠道和球囊菌的miR-250的5′端种子区序列一致，具有较高的同源性，推测宿主来源的miR-250对球囊菌的毒力和增殖产生跨界调控，而病原来源的miR-250可能分泌到宿主细胞质内，对包括免疫防御在内的生物学进程产生影响，从而促进自身的增殖。下一步将结合病原与宿主的miRNA数据，预测中蜂幼虫肠道miRNA与球囊菌mRNA、球囊菌miRNA与中蜂幼虫肠道mRNA之间的靶向结合关系，并通过生物信息学和分子生物学方法解析二者的跨界调控作用，深入探究二者的互作机制。近年来，非编码RNA (non-coding RNA，ncRNA)日益受到国内外****的关注，长链非编码RNA (long non-coding RNA，lncRNA)和环状RNA (circular RNA，circRNA)已成为全世界的研究热点。笔者团队前期已利用链特异性建库的二代测序技术对球囊菌进行测序，在全基因组水平分析和鉴定了lncRNA[46]和circRNA[47]。竞争性内源RNA (competitive endogenous RNA，ceRNA)假说认为lncRNA和circRNA等ncRNA可充当miRNA的分子海绵，通过miRNA应答元件竞争性结合miRNA，进而影响其对靶mRNA以及各类生命进程的调控[48]。球囊菌对蜜蜂幼虫的侵染机制必然伴随着复杂分子事件的发生，其中ncRNA扮演着什么角色，本研究鉴定出DEmiRNA是否与球囊菌的lncRNA和circRNA存在cross talk，是否存在DEmiRNA介导的更为复杂的ceRNA调控网络，是非常值得关注的研究方向。
综上所述，本研究对胁迫中蜂6日龄幼虫的球囊菌的miRNA差异表达谱及调控网络进行全面分析，研究结果为深入探究球囊菌对中蜂幼虫的侵染机制以及二者之间的互作机制提供了新的思路和线索，也为球囊菌关键致病因子的筛选和功能研究提供了基础。
References
| [1] | Li JH, Zheng ZY, Chen DF, Liang Q. Factors influencing Ascosphaera apis infection on honeybee larvae and observation on the infection process. Acta Entomologica Sinica, 2012, 55(7): 790-797. (in Chinese) 李江红, 郑志阳, 陈大福, 梁勤. 影响蜜蜂球囊菌侵染蜜蜂幼虫的因素及侵染过程观察. 昆虫学报, 2012, 55(7): 790-797. |
| [2] | Aronstein KA, Murray KD. Chalkbrood disease in honey bees. Journal of Invertebrate Pathology, 2010, 103(S1): S20-S29. |
| [3] | Maxfield-Taylor SA, Mujic AB, Rao S. First detection of the larval chalkbrood disease pathogen Ascosphaera apis (Ascomycota:Eurotiomycetes:Ascosphaerales) in adult Bumble bees. PLoS One, 2015, 10(4): e0124868. |
| [4] | Chen DF, Guo R, Xiong CL, Zheng YZ, Hou CS, Fu ZM. Morphological and molecular identification of chalkbrood disease pathogen Ascosphaera apis in Apis cerana cerana. Journal of Apicultural Research, 2018, 57(4): 516-521. |
| [5] | Zhang ZN, Xiong CL, Xu XJ, Huang ZJ, Zheng YZ, Luo Q, Liu M, Li WD, Tong XY, Zhang Q, Liang Q, Guo R, Chen DF. De novo assembly of a reference transcriptome and development of SSR markers for Ascosphaera apis. Acta Entomologica Sinica, 2017, 60(1): 34-44. (in Chinese) 张曌楠, 熊翠玲, 徐细建, 黄枳腱, 郑燕珍, 骆群, 刘敏, 李汶东, 童新宇, 张琦, 梁勤, 郭睿, 陈大福. 蜜蜂球囊菌的参考转录组de novo组装及SSR分子标记开发. 昆虫学报, 2017, 60(1): 34-44. |
| [6] | Guo R, Chen DF, Huang ZJ, Liang Q, Xiong CL, Xu XJ, Zheng YZ, Zhang ZN, Xie YL, Tong XY, Hou ZX, Jiang LL, Dao C. Transcriptome analysis of Ascosphaera apis stressing larval gut of Apis cerana cerana. Acta Microbiologica Sinica, 2017, 57(12): 1865-1878. (in Chinese) 郭睿, 陈大福, 黄枳腱, 梁勤, 熊翠玲, 徐细建, 郑燕珍, 张曌楠, 解彦玲, 童新宇, 侯志贤, 江亮亮, 刀晨. 球囊菌胁迫中华蜜蜂幼虫肠道过程中病原的转录组学研究. 微生物学报, 2017, 57(12): 1865-1878. |
| [7] | Chen DF, Guo R, Xiong CL, Liang Q, Zheng YZ, Xu XJ, Huang ZJ, Zhang ZN, Zhang L, Li WD, Tong XY, Xi WJ. Transcriptomic analysis of Ascosphaera apis stressing larval gut of Apis mellifera ligustica (Hyemenoptera:Apidae). Acta Entomologica Sinica, 2017, 60(4): 401-411. (in Chinese) 陈大福, 郭睿, 熊翠玲, 梁勤, 郑燕珍, 徐细建, 黄枳腱, 张曌楠, 张璐, 李汶东, 童新宇, 席伟军. 胁迫意大利蜜蜂幼虫肠道的球囊菌的转录组分析. 昆虫学报, 2017, 60(4): 401-411. |
| [8] | Guo R, Wang HP, Chen HZ, Xiong CL, Zheng YZ, Fu ZM, Zhao HX, Chen DF. Identification of Ascosphaera apis microRNAs and investigation of their regulation networks. Acta Microbiologica Sinica, 2018, 58(6): 1077-1089. (in Chinese) 郭睿, 王海朋, 陈华枝, 熊翠玲, 郑燕珍, 付中民, 赵红霞, 陈大福. 蜜蜂球囊菌的microRNA鉴定及其调控网络分析. 微生物学报, 2018, 58(6): 1077-1089. |
| [9] | Hussain M, Asgari S. MicroRNAs as mediators of insect host-pathogen interactions and immunity. Journal of Insect Physiology, 2014, 70: 151-158. |
| [10] | Lee HC, Li LD, Gu WF, Xue ZH, Crosthwaite SK, Pertsemlidis A, Lewis ZA, Freitag M, Selker EU, Mello CC, Liu Y. Diverse pathways generate microRNA-like RNAs and Dicer-independent small interfering RNAs in fungi. Molecular Cell, 2010, 38(6): 803-814. |
| [11] | Nunes CC, Sailsbery JK, Dean RA. Characterization and application of small RNAs and RNA silencing mechanisms in fungi. Fungal Biology Reviews, 2011, 25(4): 172-180. |
| [12] | Chen R, Jiang N, Jiang QY, Sun XJ, Wang Y, Zhang H, Hu Z. Exploring microRNA-like small RNAs in the filamentous fungus Fusarium oxysporum. PLoS One, 2014, 9(8): e104956. |
| [13] | Zhou JH, Fu YP, Xie JT, Li B, Jiang DH, Li GQ, Cheng JS. Identification of microRNA-like RNAs in a plant pathogenic fungus Sclerotinia sclerotiorum by high-throughput sequencing. Molecular Genetics and Genomics, 2012, 287(4): 275-282. |
| [14] | Zhou Q, Wang ZX, Zhang J, Meng HM, Huang B. Genome-wide identification and profiling of microRNA-like RNAs from Metarhizium anisopliae during development. Fungal Biology, 2012, 116(11): 1156-1162. |
| [15] | Huang Q, Evans JD. Identification of microRNA-like small RNAs from fungal parasite Nosema ceranae. Journal of Invertebrate Pathology, 2016, 133: 107-109. |
| [16] | Li Q, Cui CL, Song HS, Wang SB. The effects of mro-miR-33 on the conidial production in Metarhizium robertsii. Mycosystema, 2017, 36(6): 671-678. (in Chinese) 李琼, 崔春来, 宋红生, 王四宝. mro-miR-33在绿僵菌产孢中的作用. 菌物学报, 2017, 36(6): 671-678. |
| [17] | Croston TL, Lemons AR, Beezhold DH, Beezhold DH, Green BJ. MicroRNA regulation of host immune responses following fungal exposure. Frontiers in Immunology, 2018, 9: 170. |
| [18] | Shakeel M, Xu XX, Xu J, Li SZ, Yu JL, Zhou XQ, Xu XJ, Hu QB, Yu XQ, Jin FL. Genome-wide identification of destruxin A-responsive immunity-related microRNAs in diamondback moth, Plutella xylostella. Frontiers in Immunology, 2018, 9: 185. |
| [19] | Evans JD, Huang Q. Interactions among host-parasite microRNAs during Nosema ceranae proliferation in Apis mellifera. Frontiers in Microbiology, 2018, 9: 698. |
| [20] | Du Y, Tong XY, Zhou DD, Chen DF, Xiong CL, Zheng YZ, Xu GJ, Wang HP, Chen HZ, Guo YL, Long Q, Guo R. MicroRNA responses in the larval gut of Apis cerana cerana to Ascosphaera apis stress. Acta Microbiologica Sinica, 2019, 59(9): 1747-1764. (in Chinese) 杜宇, 童新宇, 周丁丁, 陈大福, 熊翠玲, 郑燕珍, 徐国钧, 王海朋, 陈华枝, 郭意龙, 隆琦, 郭睿. 中华蜜蜂幼虫肠道响应球囊菌胁迫的microRNA应答分析. 微生物学报, 2019, 59(9): 1747-1764. |
| [21] | Guo R, Du Y, Zhou NH, Liu SY, Xiong CL, Zheng YZ, Fu ZM, Xu GJ, Wang HP, Geng SH, Zhou DD, Chen DF. Comprehensive analysis of differentially expressed microRNAs and their target genes in the larval gut of Apis mellifera ligustica during the late stage of Ascosphaera apis stress. Acta Entomologica Sinica, 2019, 62(1): 49-60. (in Chinese) 郭睿, 杜宇, 周倪红, 刘思亚, 熊翠玲, 郑燕珍, 付中民, 徐国钧, 王海朋, 耿四海, 周丁丁, 陈大福. 意大利蜜蜂幼虫肠道在球囊菌胁迫后期的差异表达微小RNA及其靶基因分析. 昆虫学报, 2019, 62(1): 49-60. |
| [22] | Xiong CL, Du Y, Chen DF, Zheng YZ, Fu ZM, Wang HP, Geng SH, Chen HZ, Zhou DD, Wu SZ, Shi CY, Guo R. Bioinformatic prediction and analysis of miRNAs in the Apis mellifera ligustica larval gut. Chinese Journal of Applied Entomology, 2018, 55(6): 1023-1033. (in Chinese) 熊翠玲, 杜宇, 陈大福, 郑燕珍, 付中民, 王海朋, 耿四海, 陈华枝, 周丁丁, 吴素珍, 石彩云, 郭睿. 意大利蜜蜂幼虫肠道的miRNAs的生物信息学预测及分析. 应用昆虫学报, 2018, 55(6): 1023-1033. |
| [23] | Rajewsky N, Socci ND. Computational identification of microRNA targets. Developmental Biology, 2004, 267(2): 529-535. |
| [24] | Zhu KG, Liu MH, Fu Z, Zhou Z, Kong Y, Liang HW, Lin ZG, Luo J, Zheng HQ, Wan P, Zhang JF, Zen K, Chen J, Hu FL, Zhang CY, Ren J, Chen X. Plant microRNAs in larval food regulate honeybee caste development. PLoS Genetics, 2017, 13(8): e1006946. |
| [25] | Guo R, Chen DF, Diao QY, Xiong CL, Zheng YZ, Hou CS. Transcriptomic investigation of immune responses of the Apis cerana cerana larval gut infected by Ascosphaera apis. Journal of Invertebrate Pathology, 2019, 166: 107210. |
| [26] | Chen DF, Guo R, Xu XJ, Xiong CL, Liang Q, Zheng YZ, Luo Q, Zhang ZN, Huang ZJ, Kumar D, Xi WJ, Zou X, Liu M. Uncovering the immune responses of Apis mellifera ligustica larval gut to Ascosphaera apis infection utilizing transcriptome sequencing. Gene, 2017, 621: 40-50. |
| [27] | Pentimone I, Lebrón R, Hackenberg M, Rosso LC, Colagiero M, Nigro F, Ciancio A. Identification of tomato miRNAs responsive to root colonization by endophytic Pochonia chlamydosporia. Applied Microbiology and Biotechnology, 2018, 102(2): 907-919. |
| [28] | Suzuki MM, Bird A. DNA methylation landscapes:provocative insights from epigenomics. Nature Reviews Genetics, 2008, 9(6): 465-476. |
| [29] | Kouzminova E, Selker EU. Dim-2 encodes a DNA methyltransferase responsible for all known cytosine methylation in Neurospora. The EMBO Journal, 2001, 20(15): 4309-4323. |
| [30] | Jeon J, Choi J, Lee GW, Park SY, Huh A, Dean RA, Lee YH. Genome-wide profiling of DNA methylation provides insights into epigenetic regulation of fungal development in a plant pathogenic fungus, Magnaporthe oryzae. Scientific Reports, 2015, 5(1): 8567. |
| [31] | Wang YL, Wang TT, Qiao LT, Zhu JY, Fan JR, Zhang TT, Wang ZX, Li WZ, Chen AH, Huang B. DNA methyltransferases contribute to the fungal development, stress tolerance and virulence of the entomopathogenic fungus Metarhizium robertsii. Applied Microbiology and Biotechnology, 2017, 101(10): 4215-4226. |
| [32] | Lewsey MG, Hardcastle TJ, Melnyk CW, Molnar A, Valli A, Urich MA, Nery JR, Baulcombe DC, Ecker JR. Mobile small RNAs regulate genome-wide DNA methylation. Proceedings of the National Academy of Science of the United States of America, 2016, 113(6): E801-E810. |
| [33] | Xu ZB, Song LX. Research progress in miRNAs and methylation regulation. Chemistry of Life, 2016, 36(5): 677-682. (in Chinese) 徐振彪, 宋林霞. miRNA与甲基化调控研究进展. 生命的化学, 2016, 36(5): 677-682. |
| [34] | Fabbri M, Garzon R, Cimmino A, Liu ZF, Zanesi N, Callegari E, Liu SJ, Alder H, Costinean S, Fernandez-Cymering C, Volinia S, Guler G, Morrison CD, Chan KK, Marcucci G, Calin GA, Huebner K, Croce CM. MicroRNA-29 family reverts aberrant methylation in lung cancer by targeting DNA methyltransferases 3A and 3B. Proceedings of the National Academy of Science of the United States of America, 2007, 104(40): 15805-15810. |
| [35] | Evans JD, Aronstein K, Chen YP, Hetru C, Imler JL, Jiang H, Kanost M, Thompson GJ, Zou Z, Hultmark D. Immune pathways and defence mechanisms in honey bees Apis mellifera. Insect Molecular Biology, 2006, 15(5): 645-656. |
| [36] | Huang SK, Ye KT, Huang WF, Ying BH, Su X, Lin LH, Li JH, Chen YP, Li JL, Bao XL, Hu JZ. Influence of feeding type and Nosema ceranae infection on the gut microbiota of Apis cerana workers. mSystems, 2018, 3(6): e00177-18. |
| [37] | Vilcinskas A. Coevolution between pathogen-derived proteinases and proteinase inhibitors of host insects. Virulence, 2010, 1(3): 206-214. |
| [38] | Du Y, Xiong CL, Shi XL, Zheng YZ, Fu ZM, Xu XJ, Chen DF, Guo R. Transcriptome analysis of differentially expressed genes in Ascosphaera apis stressing the 6-day-old larval gut of Apis mellifera ligustica. Acta Agriculturae Zhejiangensis, 2017, 29(7): 1119-1128. (in Chinese) 杜宇, 熊翠玲, 史秀丽, 郑燕珍, 付中民, 徐细建, 陈大福, 郭睿. 意大利蜜蜂6日龄幼虫肠道内球囊菌的差异表达基因分析. 浙江农业学报, 2017, 29(7): 1119-1128. |
| [39] | Xu XG, Zhang ZP, Cheng B. Progression on the research of two-component signal transduction system in fungus and its inhibitors. Chinese Journal of Mycology, 2011, 6(6): 370-372, 380. (in Chinese) 徐西光, 张子平, 程波. 真菌双组分信号转导系统及其抑制剂研究进展. 中国真菌学杂志, 2011, 6(6): 370-372, 380. |
| [40] | Jiang C, Zhang X, Liu HQ, Xu JR. Mitogen-activated protein kinase signaling in plant pathogenic fungi. PLoS Pathogens, 2018, 14(3): e1006875. |
| [41] | Turrà D, Segorbe D, Di Pietro A. Protein kinases in plant-pathogenic fungi:conserved regulators of infection. Annual Review of Phytopathology, 2014, 52: 267-288. |
| [42] | Chen DF, Wang HQ, Li WD, Xiong CL, Zheng YZ, Fu ZM, Xu XJ, Huang ZJ, Guo R. Analysis of highly expressed genes of Ascosphaera apis infecting the gut of Apis cerana cerana larvae and its in vitro culture. Journal of Fujian Agriculture and Forestry University (Natural Science Edition), 2017, 46(5): 562-568. (in Chinese) 陈大福, 王鸿权, 李汶东, 熊翠玲, 郑燕珍, 付中民, 徐细建, 黄枳腱, 郭睿. 胁迫中华蜜蜂幼虫肠道的球囊菌及其体外培养的高表达基因分析. 福建农林大学学报(自然科学版), 2017, 46(5): 562-568. |
| [43] | Lau NC, Lim LP, Weinstein EG, Bartel DP. An abundant class of tiny RNAs with probable regulatory roles in Caenorhabditis elegans. Science, 2001, 294(5543): 858-862. |
| [44] | McClure LV, Sullivan CS. Kaposi's sarcoma herpes virus taps into a host microRNA regulatory network. Cell Host & Microbe, 2008, 3(1): 1-3. |
| [45] | LaMonte G, Philip N, Reardon J, Lacsina JR, Majoros W, Chapman L, Thornburg CD, Telen MJ, Ohler U, Nicchitta CV, Haystead T, Chi JT. Translocation of sickle cell erythrocyte microRNAs into Plasmodium falciparum inhibits parasite translation and contributes to malaria resistance. Cell Host & Microbe, 2012, 12(2): 187-199. |
| [46] | Guo R, Chen DF, Xiong CL, Hou CS, Zheng YZ, Fu ZM, Diao QY, Zhang L, Wang HQ, Hou ZX, Li WD, Kumar D, Liang Q. Identification of long non-coding RNAs in the chalkbrood disease pathogen Ascospheara apis. Journal of Invertebrate Pathology, 2018, 156: 1-5. |
| [47] | Guo R, Chen DF, Chen HZ, Fu ZM, Xiong CL, Hou CS, Zheng YZ, Guo YL, Wang HP, Du Y, Diao QY. Systematic investigation of circular RNAs in Ascosphaera apis, a fungal pathogen of honeybee larvae. Gene, 2018, 678: 17-22. |
| [48] | Karreth FA, Tay Y, Perna D, Ala U, Tan SM, Rust AG, DeNicola G, Webster KA, Weiss D, Perez-Mancera PA, Krauthammer M, Halaban R, Provero P, Adams DJ, Tuveson DA, Pandolfi PP. In vivo identification of tumor-suppressive PTEN ceRNAs in an oncogenic BRAF-induced mouse model of melanoma. Cell, 2011, 147(2): 382-395. |
