杜宇1#, 冯睿蓉1#, 王杰1, 祝智威1, 张文德1, 余岢骏1, 隆琦1, 蔡宗兵1, 解彦玲1, 熊翠玲1, 郑燕珍1,2, 陈大福1,2, 郭睿1,2
1. 福建农林大学动物科学学院(蜂学学院), 福建 福州 350002;
2. 福建农林大学蜂产品加工与应用教育部工程研究中心, 福建 福州 350002
收稿日期：2020-07-10；修回日期：2020-10-05；网络出版日期：2021-01-28
基金项目：国家自然科学基金(31702190);国家现代农业产业技术体系建设专项资金(CARS-44-KXJ7);福建省教育厅中青年教师教育科研项目(JAT170158);福建农林大学科技创新专项基金(CXZX2017342, CXZX2017343);福建农林大学硕士生导师团队项目(郭睿); 福建省病原真菌与真菌毒素重点实验室开放课题; 福建农林大学优秀硕士学位论文资助基金(杜宇)
*通信作者：郭睿, Tel/Fax: +86-591-83726835;E-mail: ruiguo@fafu.edu.cn.
#并列第一作者。
摘要：[目的] 本研究旨在探究长链非编码RNA（long non-coding RNA，lncRNA）在中华蜜蜂（Apis cerana cerana，简称中蜂）6日龄幼虫应答蜜蜂球囊菌（Ascosphaera apis，简称球囊菌）侵染过程中的差异表达谱及调控作用。[方法] 利用链特异性cDNA建库的RNA-seq技术对未被侵染及球囊菌侵染的中蜂6日龄幼虫肠道（AcCK和AcT）进行深度测序。通过相关生物信息学软件分析lncRNA的结构特征和表达谱。筛选并分析差异表达lncRNA（differentially expressed lncRNA，DElncRNA）的顺式（cis）作用及竞争性内源RNA（competing endogenous RNA，ceRNA）调控网络。采用RT-qPCR验证测序数据及DElncRNA差异变化趋势的可靠性。[结果] AcCK和AcT共预测出642个已知lncRNA和487个新lncRNA。与蛋白编码基因相比，上述中蜂lncRNA外显子数更少、长度更短且表达量更低。43个antisense lncRNA与40个正义链mRNA之间互补配对。AcCK和AcT比较组包含367个上调lncRNA和268个下调lncRNA。有194个DElncRNA潜在调控461个上下游基因，并涉及细胞进程、代谢进程和催化活性等38个功能条目以及氨基酸代谢、内吞作用和MAPK等191条通路。此外，有180个DElncRNA可靶向结合50个DEmiRNA，进而调控6365个mRNA；三者之间形成较为复杂的ceRNA调控网络。[结论] 中蜂6日龄幼虫肠道的部分lncRNA可作为antisense lncRNA参与应答球囊菌侵染；部分DElncRNA可通过cis作用调节物质代谢和免疫途径相关的上下游基因，从而介导宿主的侵染应答；TCONS_00010661和TCONS_00003104等DElncRNA可通过ceRNA网络调控Jak-STAT和氧化磷酸化等通路及富集基因，进而参与宿主的侵染应答。
关键词：中华蜜蜂幼虫蜜蜂球囊菌长链非编码RNA竞争性内源RNA免疫应答
Long non-coding RNA response of 6-day-old Apis cerana cerana larvae to Ascosphaera apis infection
Yu Du1#, Ruirong Feng1#, Jie Wang1, Zhiwei Zhu1, Wende Zhang1, Kejun Yu1, Qi Long1, Zongbing Cai1, Yanling Xie1, Cuiling Xiong1, Yanzhen Zheng1,2, Dafu Chen1,2, Rui Guo1,2
1. College of Animal Sciences(College of Bee Science), Fujian Agriculture and Forestry University, Fuzhou 350002, Fujian Province, China;
2. Engineering Research Center of Processing and Application of Bee Products of Ministry of Education, Fujian Agriculture and Forestry University, Fuzhou 350002, Fujian Province, China
Received: 10 July 2020; Revised: 5 October 2020; Published online: 28 January 2021
*Corresponding author: Rui Guo, Tel/Fax: +86-591-83726835; E-mail: 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 Scientific and Technical Innovation Fund of Fujian Agriculture and Forestry University (CXZX2017342, CXZX2017343), by the Master Supervisor Team Fund of Fujian Agriculture and Forestry University (Rui Guo), by the Open Project of Key Laboratory of Pathogenic Fungi and Mycotoxins of Fujian Province (Fujian Agriculture and Forestry University) and by the Fund for Excellent Master Dissertation of Fujian Agriculture and Forestry University (Yu Du)
#These authors contributed equally to this work.
Abstract: [Objective] This study aims to investigate the differential expression pattern of long non-coding RNAs (lncRNAs) and their regulatory function involved in Apis cerana cerana 6-day-old larval gut response to Ascosphaera apis infection. [Methods] Un-infected and Ascosphaera apis-infected 6-day-old larval guts of Apis cerana cerana (AcCK and AcT) were sequenced using strand-specific cDNA library-based RNA-seq technology. Structural characteristics and expression pattern of lncRNAs were analyzed using related bioinformatic softwares. DElncRNAs were screened followed by investigation of their cis-acting role and competitive endogenous RNA (ceRNA) network. RT-qPCR was conducted to verify the sequencing data and expression pattern of DElncRNAs. [Results] Here, 642 known lncRNAs and 487 novel lncRNAs were identified. Compared with mRNAs, these Apis cerana cerana lncRNAs were shorter in exon and intron length, fewer in exon number and lower in expression level. Additionally, 43 antisense lncRNAs were discovered to have complementary relationship with 40 sense-strand mRNAs. In AcCK vs. AcT comparison group, 367 up-regulated lncRNAs and 268 down-regulated ones were identified. In total, 194 DElncRNAs were found to potentially regulate 461 upstream and downstream genes, which were annotated to 38 functional terms such as cellular process, metabolic process and catalytic activity, as well as 191 pathways including amino acid metabolism, endocytosis and MAPK signaling pathways. Moreover, 180 DElncRNAs can target 50 DEmiRNAs, further regulating 6365 mRNAs; additionally, complex regulatory networks existed among them. [Conclusion] These results demonstrated that partial antisense lncRNAs may participate in host response to Ascosphaera apis infection; some DElncRNAs may regulate upstream and downstream genes relative to material metabolism and immune-associated pathways, thus mediating host Ascosphaera apis-response; a portion of DElncRNAs including TCONS_00010661 and TCONS_00003104 were likely to regulate Jak-STAT signaling pathway and oxidative phosphorylation and corresponding enriched genes via ceRNA networks, further participating in the response of host to Ascosphaera apis invasion.
Keywords: Apis cerana ceranalarvaAscosphaera apislong non-coding RNAcompetitive endogenous RNAimmune response
中华蜜蜂(Apis cerana cerana，简称中蜂)是东方蜜蜂(Apis cerana)的指名亚种，也是我国特有的蜜蜂种质资源，经过长期进化已高度适应我国的地理环境，具有特殊的生态价值和重要的经济价值[1]。近年来，由于西方蜜蜂(Apis mellifera)的大规模饲养环境污染、植被破坏、滥用杀虫剂和病原侵染等因素的影响，中蜂的抗逆性和抗病性降低。蜜蜂球囊菌(Ascosphaera apis，简称球囊菌)是一种主要导致蜜蜂幼虫罹患白垩病的致死性真菌病原[2]。近期研究结果表明，球囊菌对成年熊蜂[3]同样具有侵染性。在养蜂生产中该病常在意大利蜜蜂(Apis mellifera ligustica，简称意蜂)蜂群中暴发，仅在中蜂蜂群内偶尔发生，体现出二者在球囊菌抗性方面存在明显差异。笔者所在团队前期的转录组研究结果显示，意蜂幼虫和中蜂幼虫响应球囊菌侵染的细胞和体液免疫应答存在差异，可能是导致两个蜂种具有不同球囊菌抗性差异的重要原因[4-5]；进一步通过比较转录组分析揭示了免疫基因的数量差异可能是意蜂幼虫和中蜂幼虫产生球囊菌抗性差异的原因之一，在与球囊菌长期协同进化的过程中，中蜂幼虫可能通过提高能量利用限制病原增殖[6]。目前尚缺乏有效治疗白垩病的手段，而通过科学的饲养管理并配合长期的遗传品系筛选，虽能在一定程度上减轻该病的危害，但总体效果并不明显，也无法阻止该病的重复发生[7]。因此，探明中蜂幼虫响应球囊菌侵染的应答机制，对于深入理解二者的互作机制以及开发有效的白垩病治疗手段具有重要意义。
长链非编码RNA (long non-coding RNA，lncRNA)是一类长度大于200 nt，包含2个以上外显子，缺少有效开放性阅读框(ORF)且不直接编码蛋白的内源性非编码RNA (non-coding RNA，ncRNA)，可参与染色质修饰、核内运输及转录调控等各类生命进程[8]。LncRNA具有表达量相对较低，物种间保守性较差以及组织特异性较高等特点[9]。根据lncRNA基因在基因组上的位置和来源一般可分为正义lncRNA (sense lncRNA)、反义lncRNA (antisense lncRNA)、内含子lncRNA (intronic lncRNA)和基因间区lncRNA (intergenic lncRNA，lincRNA)[10]。此外，lncRNA还可作为微小RNA (microRNA，miRNA)的前体以及通过miRNA应答元件(miRNA response element，MRE)竞争性吸附结合miRNA，间接影响miRNA对靶mRNA的调控作用[11-12]。近年来，随着高通量测序技术的不断发展，人们已陆续在东方蜜蜂[13]、西方蜜蜂(Apis mellifera)[14]、黑腹果蝇(Drosophila melanogaster)[15]、家蚕(Bombyx mori)[16]、褐飞虱(Nilaparvata lugens)[17]和埃及伊蚊(Aedes aegypti)[18]等昆虫中鉴定出大量lncRNA，并证实lncRNA可参与调节昆虫的生殖发育、行为转变以及免疫应答等过程[13, 15, 18]。例如，Valanne等[15]发现藤黄微球菌(Micrococcus luteus)和布拉迪小环腹瘿蜂(Leptopilina boulardi)侵染导致果蝇体内lincRNA-IBIN的表达显著提高，继而影响碳水化合物代谢以及NF-κB、Toll信号通路和抗菌肽合成等先天性免疫途径的调控。Etebari等[18]研究发现埃及伊蚊被登革病毒感染后，lincRNA_1317、lincRNA_1613和lincRNA_2329等72个lincRNA在中肠内的表达显著提高，其中lincRNA_1317还与miR-1、miR-9a、miR-33、miR-252-3p和miR-278-5p等21个miRNA之间具有潜在结合关系，作者推测lincRNA_1317可作为竞争性内源RNA (competing endogenous RNA, ceRNA)参与埃及伊蚊的抗病毒免疫反应。此前，蜜蜂的lncRNA相关研究主要涉及蜜蜂的级型分化[19]、卵巢发育[14]及劳动分工[20]等方面，但lncRNA介导蜜蜂及其病原互作的研究十分有限[13, 21]。Jayakodi等[13]在基因组水平上分别从西方蜜蜂和东方蜜蜂体内鉴定到1514和2470个lincRNA，且二者之间的lincRNA具有较高的保守性(95%)，但与黑腹果蝇、秀丽隐杆线虫(Caenorhabditis elegans)以及其他哺乳动物之间的lincRNA保守性较低。此外，作者还发现11个lincRNA在东方蜜蜂感染囊状幼虫病毒(SBV)和残翅病毒(DWV)后显著差异表达，并推测上述lncRNA参与了东方蜜蜂的抗病毒免疫反应。近期，笔者课题组研究发现意蜂工蜂中肠的lncRNA可能通过顺式(cis)和反式(trans)作用调控上下游基因的表达，并可作为ceRNA参与影响宿主的物质和能量代谢以及细胞和体液免疫过程以响应东方蜜蜂微孢子虫(Nosema ceranae)的侵染过程，表明lncRNA在蜜蜂与病原互作过程中扮演重要角色[21]。
球囊菌孢子经蜜蜂幼虫摄食进入中肠后，在CO2刺激下处于低水平萌发状态，直至6日龄幼虫与预蛹过渡期，后肠与中肠隔膜消失，孢子随食物残渣进入后肠，接触O2后迅速萌发繁殖成大量菌丝并相继穿透肠壁和体壁，导致幼虫死亡，最终从幼虫尾部蔓延包裹至整个白垩病症虫尸[7]。因此，在球囊菌侵染后期(6日龄)，中蜂幼虫和球囊菌之间存在较强的相互作用。前期研究中，笔者课题组已从mRNA、miRNA和环状RNA (circular RNA，circRNA)的角度系统分析和探讨了中蜂6日龄幼虫响应球囊菌侵染的应答机制，揭示了差异表达基因(differentially expressed gene，DEG)、差异表达miRNA (DEmiRNA)和差异表达circRNA (DEcircRNA)在宿主胁迫应答中的潜在作用[22-24]。但对于lncRNA是否参与中蜂幼虫的球囊菌胁迫应答以及如何发挥作用，相关研究仍然缺失。鉴于宿主-病原互作是非常复杂的分子事件，不同RNA分子必然参与其中并扮演重要角色，本研究利用基于链特异性cDNA建库的RNA-seq技术对球囊菌侵染的中蜂6日龄幼虫肠道和未受侵染的6日龄幼虫肠道进行深度测序，全基因组鉴定lncRNA并分析其种类和结构特征，进而深入分析和探讨差异表达lncRNA (Differentially expressed lncRNA，DElncRNA)的作用方式和调控网络，以期解析中蜂6日龄幼虫响应球囊菌侵染的lncRNA应答，为阐明宿主应答的分子机制及开发白垩病的新型治疗手段打下基础。
1 材料和方法 1.1 供试蜜蜂幼虫及球囊菌 中蜂幼虫取自福建农林大学动物科学学院(蜂学学院)蜜蜂保护课题组的实验蜂群，球囊菌菌株[24-25]由福建农林大学动物科学学院(蜂学学院)蜜蜂保护实验室分离与活化。
1.2 中蜂幼虫肠道样品的制备及深度测序 根据笔者课题组前期已建立的技术流程制备球囊菌的纯化孢子及饲喂接种中蜂幼虫[23, 25-26]：将标记的空巢脾放至群势较强且无白垩病症状的中蜂蜂群内，供蜂王产卵。待幼虫发育至2日龄时，用移虫针将其移至已预置50 μL饲料的48孔细胞培养板，置于(35.0±0.5)℃、90%相对湿度(RH)条件下饲养。处理组的每只3日龄幼虫饲喂含球囊菌孢子的饲料(1×107个孢子/mL)，对照组的每只3日龄幼虫饲喂不含球囊菌孢子的饲料，待饲料食尽，每24 h更换1次不含孢子的饲料，并及时移除死亡幼虫。在超净工作台中分别剖取处理组中蜂6日龄幼虫肠道(AcT)和对照组中蜂6日龄幼虫肠道(AcCK)，每9只肠道置于1个RNA-Free EP管，经液氮速冻后立即转移到–80 ℃超低温冰箱保存备用。
参照笔者课题组前期研究的方法[21]提取AcCK和AcT的总RNA并构建cDNA链特异性文库。委托广州基迪奥生物科技有限公司对建好的文库进行双端测序，测序平台为Illumina HiSeqTM 4000。已将测序数据上传NCBI SRA数据库，BioProject号：PRJNA560730。
1.3 数据质控、lncRNA的预测及分析 参照已报道的方法[27-28]对下机的原始读段(raw reads)进行质控，利用Perl脚本剔除含有接头、未知碱基N > 10%及低质量的reads以获得有效读段(clean reads)。并进一步过滤可比对到rRNA数据库的clean reads。然后利用TopHat软件将过滤后的reads比对到东方蜜蜂参考基因组(assembly ACSNU-2.0，http://journals.im.ac.cn/html/actamicrocn/2021/5/www.ncbi.nlm.nih.gov/genome/12051?genome_assembly_id=258367)，并基于参考序列对AcCK和AcT组的转录本进行组装。
参照郭睿等和熊翠玲等的方法[27-28]对下机的原始读段(raw reads)，联用CPC和CNCI软件对转录本的编码能力进行预测，取二者的交集作为可靠的预测结果。其中，与其他物种已报道的lncRNA具有同源性的称之为已知lncRNA (known lncRNA)，在其他物种中不存在同源性的lncRNA且在东方蜜蜂中未报道的称之为新lncRNA (novel lncRNA)。采用FPKM (fragments per kilobase of transcript per million mapped reads)法计算lncRNA的表达量，并统计和比较lncRNA和mRNA的内含子数量、外显子数量、外显子长度和转录本长度。
反义lncRNA (antisense lncRNA)与正义链mRNA之间可通过碱基的互补配对作用在转录后水平参与调控基因沉默、转录及mRNA的稳定性等过程[29]。利用RNAplex软件预测antisense lncRNA与正义链mRNA之间的互补结合，并根据其热力学结构计算最小自由能来预测最佳碱基配对关系。
1.4 DELncRNA的上下游基因的功能及通路注释 根据edgeR方法(http://www.bioconductor.org/packages/release/bioc/html/edgeR.html)计算AcCK和AcT组间显著差异，并按照|log2(Fold change)|≥1且P≤0.05的标准筛选DElncRNA。LncRNA的功能与其坐标临近的蛋白编码基因相关，并通过cis作用对基因表达进行转录或转录后水平的调控[30]。根据lncRNA的基因位置筛选位于其上下游10 kb以内蛋白编码基因，并利用Omicshare平台(http://www.omicshare.com/tools/index.php/)的相关软件和超几何检验对DElncRNA的上下游基因进行GO和KEGG数据库注释。
1.5 DElncRNA的ceRNA调控网络构建和分析 前期研究中，笔者曾利用small RNA-seq技术对球囊菌侵染的中蜂6日龄幼虫和未受侵染的6日龄幼虫肠道进行测序，获得了高质量的miRNA组学数据，并解析了宿主响应病原侵染的miRNA差异表达谱及调控网络[23]。结合前期筛选出的DEmiRNA[23]，利用Target Finder软件预测DElncRNA靶向结合的DEmiRNA及DEmiRNA靶向结合的mRNA。根据靶mRNA的KEGG pathway注释信息进一步筛选注释到Jak-STAT信号通路和氧化磷酸化通路的靶mRNA及存在靶向结合关系的DEmiRNA及相应的DElncRNA，构建DElncRNA-DEmiRNA-mRNA调控网络并通过利用Cytoscape软件进行可视化，筛选靶向结合关系的标准为结合自由能小于-17 kcal/mol。
1.6 DElncRNA的RT-qPCR验证 为验证深度测序数据的可靠性，从DElncRNA中随机选取9个(TCONS_00013994、TCONS_00001628、TCONS_00013493、TCONS_00019559、TCONS_00024461、TCONS_00032169、XR_001765267、XR_001765889和XR_001766314)进行RT-qPCR验证。根据上述DElncRNA的碱基序列，利用DNAMAN软件设计相应的特异性上下游引物，并委托上海生工生物工程公司合成。利用RNA抽提试剂盒(TaKaRa公司，中国)分别提取AcCK和AcT的总RNA，反转录得到相应的cDNA。qPCR在ABI QuanStudio 3荧光定量PCR系统(Thermo Fisher公司，美国)中进行。反应程序设置为：95 ℃ 3 min；95 ℃ 10 s，60 ℃ 30 s，40个循环；溶解曲线程序默认系统设置。每个反应进行3次重复。以actin基因作为内参，采用2–△△Ct法[31]，对上述DElncRNA进行相对定量。利用Graphpad Prism 7软件进行数据处理及绘图，并利用t检验计算组间差异显著性。相关引物信息详见表 1。
表 1. 本研究使用的引物 Table 1. Primers used in this work
| Primer names | Primer sequences (5?→3?) |
| TCONS_00013994-F | CGGTGAAGGTAGAGGAGAGA |
| TCONS_00013994-R | CGATTGTTACGAAAGCACAG |
| TCONS_00001628-F | CGCTGGAATAAAGATGAAGG |
| TCONS_00001628-R | CAGATTGCCAAGCACCTCT |
| TCONS_00013493-F | CTGCGAGAAAGCACATTG |
| TCONS_00013493-R | GCGTAAGAACAGGATGGAA |
| TCONS_00019559-F | TGAAGATACAGGAAGTTGAGGC |
| TCONS_00019559-R | TATGGGTGACGCCAAATG |
| TCONS_00024461-F | AAACACGGTCACACTTCGC |
| TCONS_00024461-R | GAAATCCGATGATGGTTCG |
| TCONS_00032169-F | CGTGTGTATGGCGTATGAG |
| TCONS_00032169-R | GAGACATTGTGTTCTTCACCAC |
| XR_001765889-F | GTAGATACAACGCTTGGTGTGA |
| XR_001765889-R | TTTCCTTGATGACGCCAC |
| XR_001765267-F | CGTGAAGGTTAGGATTGACC |
| XR_001765267-R | TCAGTGAGCGAAGTGGTGA |
| XR_001766314-F | GCGTTGTAGTGGATAAGAATGG |
| XR_001766314-R | TGTCGTATGGTTGCTCTGTATG |
| actin-F | GCAGTTAGGCACTTTCTCTGA |
| actin-R | CGTTTCCCACTTCTTCGTA |
表选项
2 结果和分析 2.1 深度测序数据的质控 AcCK和AcT的测序分别得到213479408和249839658条raw reads，经严格质控分别得到213458846和249809780条clean reads，占raw reads的比例均达到99.99%，Q30分别达到96.65%和96.59%。分别有134564975和117418047条clean reads可比对到东方蜜蜂参考基因组。以上结果说明本研究的测序数据质量良好，可满足分析需要。
2.2 LncRNA的生物信息学预测及Venn分析 共鉴定到1129个lncRNA，包含642个已知(known) lncRNA和487个新(novel) lncRNA。在AcCK中鉴定到520个known lncRNA和464个novel lncRNA；在AcT中鉴定到529个known lncRNA和420个novel lncRNA。Venn分析结果显示，AcCK和AcT共有的known lncRNA和novel lncRNA分别为407和397个；特异性表达的lncRNA分别为180和145个(图 1)。
 |
| 图 1 球囊菌侵染及未侵染的中蜂6日龄幼虫肠道中lncRNA的Venn分析 Figure 1 Venn analysis of lncRNAs identified in A. apis-infected and uninfected guts of A. c. cerana 6-day-old larvae. |
| 图选项 |
2.3 LncRNA的结构特征分析及比较 结构特征的比较分析结果显示，中蜂lncRNA和mRNA的内含子长度介于0–10000 bp，其中多数lncRNA和mRNA的长度分布介于0–2500 bp；lncRNA的内含子长度较mRNA更短(图 2-A)。此外，lncRNA和mRNA的外显子长度介于0–15000 bp，多数lncRNA和mRNA的长度分布介于0–5000 bp；类似的，lncRNA的外显子长度短于mRNA (图 2-B)；lncRNA的外显子数量少于mRNA，含有3个外显子的lncRNA数量最多(35个)(图 2-C)。LncRNA的转录本长度总体较mRNA更短(图 2-D)。
 |
| 图 2 中蜂6日龄幼虫肠道中lncRNA与mRNA的结构特征比较 Figure 2 Comparison of structural property between lncRNAs and mRNAs identified in A. c. cerana 6-day-old larval gut. A: comparison of intron length; B: comparison of exon length; C: comparison of exon number; D: comparison of transcript length. |
| 图选项 |
2.4 Antisense lncRNA分析 从中蜂6日龄幼虫肠道中共预测出43个antisense lncRNA与40个正义链mRNA存在互补配对结合关系。例如，TCONS_00018864可靶向结合王浆主蛋白3编码基因(XM_017055601)的mRNA，TCONS_00018867可靶向结合王浆主蛋白8同种型X1编码基因(XM_017055589) 的mRNA，TCONS_00028357可靶向结合piopio-like蛋白编码基因(XM_017061795) 的mRNA，TCONS_00032421可靶向结合丝氨酸/苏氨酸蛋白激酶异构体编码基因(XM_017064272)的mRNA，TCONS_00033001可靶向结合微管相关蛋白FutsCH型异构体X9编码基因(XM_0170645)的mRNA。上述5个antisense lncRNA与相应的正义链mRNA的配对结合关系详见图 3。
 |
| 图 3 中蜂幼虫6日龄肠道中antisense lncRNA及相应正义链mRNA的配对预测 Figure 3 Prediction of pairing between antisense lncRNAs and corresponding sense-strand mRNAs in A. c. cerana 6-day-old larval gut. |
| 图选项 |
2.5 LncRNA的表达谱分析 AcCK中lncRNA的表达量介于0.001-142.60之间，表达量最高的lncRNA是TCONS_00032019 (FPKM=142.60)，其次是XR_001765880 (FPKM= 118.99)和XR_001765028 (FPKM=113.64)。AcT中lncRNA的表达量介于0.001与145.90之间，表达量最高的lncRNA为XR_001765880 (FPKM=145.90)，其次为XR_001765028 (FPKM=133.88)和TCONS_00009886 (FPKM=131.19)。
2.6 DElncRNA的筛选及cis作用分析 差异表达分析结果显示，AcCK vs. AcT比较组包含367个上调lncRNA和268个下调lncRNA。上调lncRNA的上调倍数介于1.00至12.47，上调倍数最大的lncRNA是TCONS_00017531[log2 (Fold change)=12.47]，其次是XR_001765250 [log2(Fold change)=10.94]和XR_001765500[log2 (Fold change)=10.91](表 2)；下调lncRNA的下调倍数介于–13.01至–1.07，下调幅度最大的是TCONS_00002633[log2(Fold change)=–13.01]，其次是TCONS_00019097[log2(Fold change)=–12.80]和TCONS_00011156[log2(Fold change)=–12.78](表 3)。
表 2. AcCK vs AcT比较组中上调倍数最高的前十位lncRNA Table 2. Top 10 lncRNAs with the most up-regulation times in AcCK vs AcT comparison group
| LncRNA ID | FPKM in AcCK | FPKM in AcT | Log2(Fold change) | P value | FDR |
| TCONS_00017531 | 0.001 | 5.69 | 12.47 | 2.59E-212 | 5.87E-210 |
| XR_001765250 | 0.001 | 1.96 | 10.94 | 6.15E-57 | 1.70E-55 |
| XR_001765500 | 0.001 | 1.92 | 10.91 | 2.01E-47 | 4.07E-46 |
| TCONS_00003662 | 0.001 | 1.52 | 10.57 | 7.59E-09 | 2.14E-08 |
| XR_001766516 | 0.001 | 1.48 | 10.53 | 1.51E-20 | 9.07E-20 |
| TCONS_00012588 | 0.001 | 1.47 | 10.52 | 9.64E-07 | 2.30E-06 |
| TCONS_00029241 | 0.001 | 1.34 | 10.39 | 2.88E-78 | 1.21E-76 |
| XR_001765462 | 0.001 | 0.76 | 9.57 | 2.99E-14 | 1.30E-13 |
| TCONS_00002261 | 0.001 | 0.56 | 9.13 | 9.50E-19 | 5.20E-18 |
| TCONS_00033491 | 0.001 | 0.50 | 8.97 | 0.0019575 | 0.0034385 |
表选项
表 3. AcCK vs AcT比较组中下调倍数最高的前十位lncRNA Table 3. Top 10 lncRNAs with the most down-regulation in AcCK vs. AcT comparison group
| LncRNA ID | FPKM in AcCK | FPKM in AcT | Log2(Fold change) | P value | FDR |
| TCONS_00002633 | 8.25 | 0.001 | –13.01 | 5.10E-48 | 1.09E-46 |
| TCONS_00019097 | 7.13 | 0.001 | –12.80 | 8.80E-39 | 1.33E-37 |
| TCONS_00011156 | 7.03 | 0.001 | –12.78 | 2.15E-36 | 2.79E-35 |
| TCONS_00023972 | 7.00 | 0.001 | –12.77 | 1.33E-34 | 1.58E-33 |
| TCONS_00005365 | 6.29 | 0.001 | –12.612 | 4.13E-33 | 4.63E-32 |
| TCONS_00029156 | 5.87 | 0.001 | –12.52 | 3.12E-46 | 5.89E-45 |
| TCONS_00017177 | 5.72 | 0.001 | –12.48 | 4.97E-28 | 4.40E-27 |
| TCONS_00034056 | 5.38 | 0.001 | –12.39 | 3.28E-78 | 1.33E-76 |
| TCONS_00015195 | 5.26 | 0.001 | –12.36 | 3.17E-29 | 2.99E-28 |
| TCONS_00026460 | 4.74 | 0.001 | –12.21 | 4.13E-33 | 4.63E-32 |
表选项
AcCK vs AcT比较组的194个DElncRNA共预测到461个上下游基因。GO分析结果表明这些上下游基因可注释到细胞进程(120)、单一有机体进程(109)和代谢进程(94)等17个生物学进程相关条目，细胞(60)、细胞组件(60)和细胞膜(55)等12个细胞组分相关条目，以及结合(125)、催化活性(89)和转运活性(26)等9个分子功能相关条目(图 4)。括号内的数字代表注释到该条目的上下游基因数量。
 |
| 图 4 DElncRNA的上下游基因的GO分类 Figure 4 GO categorization of upstream and downstream genes of DElncRNAs. |
| 图选项 |
KEGG数据库注释结果显示，上述上下游基因还可注释到191条通路，包括Hippo(6)、PI3K-Akt(3)、Notch(2)、MAPK(2)、FoxO(2)、mTOR(1)、Wnt(1)和Toll-like受体(1)等8条信号通路，碳水化合物代谢(10)、氨基酸代谢(8)和脂质代谢(2)等9条物质代谢通路。此外，部分上下游基因可注释到细胞凋亡(4)、内吞作用(2)、MAPK信号通路(2)、Toll-like受体信号通路(1)、吞噬体(1)及泛素介导的蛋白水解(1)等6条细胞和体液免疫途径。
2.7 DElncRNA-DEmiRNA-mRNA调控网络的构建及分析 靶向预测结果显示，50个DEmiRNA (包含29个上调和21个下调miRNA)与180个DElncRNA(包含112个上调和68个下调lncRNA)、6365个mRNA存在靶向结合关系，DElncRNA与DEmiRNA及mRNA之间形成复杂的ceRNA调控网络(图 5)。DElncRNA TCONS_00015199、TCONS_00030368和TCONS_00007304靶向结合的DEmiRNA数量最多，分别为13、8和7个。DEmiRNA novel-m0035-3p、miR-1277-x和miR-6039-x靶向的DElncRNA数量较多，分别为20、17和14个(图 5)。
 |
| 图 5 中蜂6日龄幼虫肠道响应球囊菌侵染的ceRNA调控网络 Figure 5 CeRNA regulation network engaged in A. apis-response of A. c. cerana 6-day-old larvae. |
| 图选项 |
进一步分析发现，在Jak-STAT信号通路相关的ceRNA调控网络中，23个上调lncRNA (TCONS_00019995、TCONS_00003104及TCONS_00032449等)和6个下调lncRNA (TCONS_00002228、TCONS_00002262及TCONS_00008538等)可靶向4个DEmiRNA (包含上调表达的miR-21-x、miR-6052-x和miR-101-y以及下调表达的miR-1277-x)，进而结合3个DEmRNA (包含上调表达的XM_017053583.1以及下调表达的XM_017053215.1和XM_017064778.1)。下调表达的TCONS_00010661靶向上调表达的miR-101-y，进而结合下调表达的XM_017064778.1；下调lncRNA (TCONS_00002228和TCONS_00026779)靶向上调miR-21-x，进而结合下调表达的XM_017053215.1(图 6-A)。
 |
| 图 6 中蜂6日龄幼虫肠道响应球囊菌侵染的Jak-STAT信号通路和氧化磷酸化通路的ceRNA调控网络 Figure 6 CeRNA regulation networks of Jak-STAT signaling pathway and oxidative phosphorylation engaged in A. apis-response of A. c. cerana 6-day-old larvae. A: Regulation network of Jak-STAT signaling pathway- associated target DEmRNAs and corresponding DEmiRNAs and DElncRNAs involved in ceRNA network; B: Regulation network of oxidative phosphorylation-associated target DEmRNAs and corresponding DEmiRNAs and DElncRNAs involved in ceRNA network. |
| 图选项 |
此外，在氧化磷酸化相关的ceRNA调控网络中，29个上调lncRNA (TCONS_00030368、TCONS_00019036及TCONS_00012533等)和16个下调lncRNA (TCONS_00010661、TCONS_00018603及TCONS_00018048等)靶向4个DEmiRNA (包含上调表达的miR-101-y和下调表达的miR-1277-x、novel-m0003-3p及novel-m0035- 3p)，进而结合4个DEmRNA (包含上调表达的XM_017057859.1和下调表达的XM_017061023.1、XM_017065759.1及XM_017067246.1)。下调TCONS_00010661靶向上调miR-101-y，进而结合下调表达的XM_017061023.1和XM_017065759.1；TCONS_00030368、TCONS_00019036和TCONS_00003104等14个上调lncRNA靶向下调miR-1277-x，进而结合上调表达的XM_017057859.1 (图 6-B)。
2.8 DElncRNA的RT-qPCR验证 RT-qPCR结果显示，随机挑取3个上调lncRNA和6个下调lncRNA的差异表达趋势与测序结果中的差异表达趋势一致(图 7)，说明本研究的测序数据和lncRNA的差异表达趋势真实可靠。
 |
| 图 7 DElncRNA的RT-qPCR验证 Figure 7 RT-qPCR verification of DElncRNAs. Three biological replicates were analyzed, and are represented as Mean±SD. ns: P > 0.05, *: P < 0.05, **: P < 0.01. |
| 图选项 |
2.9 LncRNA介导中蜂6日龄幼虫应答球囊菌胁迫的机制 中蜂6日龄幼虫肠道中的DElncRNA以多种方式参与宿主的胁迫应答，TCONS_00018864和TCONS_00032421等antisense lncRNA可能通过靶向结合编码MRJP3和丝氨酸/苏氨酸蛋白激酶等蛋白的正义链mRNA，调控mRNA的稳定性及相关基因的转录，进而参与宿主对球囊菌侵染的应答(图 8)；部分DElncRNA可能通过调节物质代谢通路和免疫通路相关上下游基因的转录介导中蜂6日龄幼虫应答球囊菌侵染(图 8)；TCONS_00010661和TCONS_00003104等DElncRNA可能作为“分子海绵”，与circRNA等通过竞争性结合相应的miRNA减轻miRNA对下游靶基因的抑制或降解，从而调节Jak-STAT信号通路和氧化磷酸化等通路及富集基因，并参与宿主的胁迫应答(图 8)。
 |
| 图 8 LncRNA介导中蜂6日龄幼虫应答球囊菌胁迫的机制模式图 Figure 8 A working model showing lncRNA-mediated mechanism underlying A. apis-response of A. c. cerana 6-day-old larvae. |
| 图选项 |
3 讨论 中蜂幼虫与球囊菌经过长期协同进化形成了特定的互作机制，因而是研究蜜蜂宿主及其病原互作的良好模型[7]。由于球囊菌基因组[32]和东方蜜蜂基因组[33]的序列和功能注释信息公布时间较晚，二者的组学和分子生物学研究较为滞后。近年来，笔者团队以中蜂幼虫-球囊菌的互作为研究模型进行了一系列探索，在组织学和分子水平首次证实球囊菌对中蜂雄蜂和工蜂幼虫具有侵染性[2]；在mRNA组学水平解析了中蜂幼虫对球囊菌的免疫应答[22]；全面分析和探讨了宿主miRNA和circRNA在响应球囊菌侵染过程中的潜在功能[23-24]。目前，越来越多的研究结果表明lncRNA广泛参与了果蝇[15]、家蚕[16]和埃及伊蚊[18]及其病原的相互作用过程，但迄今还没有lncRNA介导中蜂及其病原互作的研究报道。前期研究中，笔者团队在实验室条件下对保存的球囊菌菌株进行活化，利用纯化的病原孢子对中蜂幼虫进行饲喂接种，中蜂的死亡率在幼虫期很低，而在预蛹期爆发式上升，预蛹期球囊菌菌丝率先从宿主尾部穿出，进而蔓延包裹整个虫尸，此外通过PCR仅能在球囊菌接种的幼虫样品中检测到病原的特异性片段[5]，证实了球囊菌能够有效感染中蜂幼虫。本研究活化的球囊菌来源于同一保存菌株，中蜂幼虫取自相同实验蜂群，实验条件也与前期研究保持一致，理论上饲喂接种球囊菌孢子的中蜂幼虫可被有效感染。我们在人工饲养中蜂幼虫时也观察到了典型的白垩病症状。本研究通过对球囊菌侵染的中蜂6日龄幼虫肠道及未受侵染的6日龄幼虫肠道进行测序，共鉴定到642个known lncRNA和487个novel lncRNA；在AcCK和AcT中分别鉴定到984和949个lncRNA，二者的共有lncRNA为804个，特有lncRNA分别为180和145个，推测上述804个共有lncRNA在中蜂幼虫的生理和病理状态下均发挥作用，而145个特异性表达的lncRNA参与了中蜂幼虫的球囊菌侵染应答。结构特征分析结果显示，相较于mRNA，中蜂lncRNA转录本长度、内含子长度和外显子长度均更短，外显子数量更少，表达量更低，这与意蜂[21]和小菜蛾(Plutella xylostella)[34]等昆虫以及球囊菌[35]和东方蜜蜂微孢子虫[36]等真菌病原的lncRNA结构特征相似，说明上述结构特征是真核生物lncRNA的普遍结构特征。
Antisense lncRNA可作为潜在的调节器，并通过与正义链mRNA互补配对结合而调控mRNA的稳定性和基因的表达[37]。王浆主蛋白(major royal jelly protein，MRJP)作为蜂王浆蛋白中的主要成分，可为蜜蜂生长和繁殖提供营养成分。其中MRJP3富含可被幼虫吸收和利用的氮元素[38]，并具有显著的免疫调节效应[39]。本研究发现，lncRNA(TCONS_00018864)在AcT中显著上调表达[log2 (Fold Change)=1.53，P=0.03]，并可反向结合MRJP3编码基因(XM_017055601)，暗示在宿主的球囊菌侵染应答中，TCONS_00018864可作为antisense lncRNA影响mrjp3的表达水平，从而调节对球囊菌的免疫防御。蜜蜂等昆虫在与病原体的长期协同进化过程中，形成了以吞噬、包囊、黑化和蛋白的酶促水解为代表的细胞免疫，以及Jak-STAT、Toll-like受体和MAPK等信号通路为代表的体液免疫[25, 40]。Guo等[41]研究发现具有调节昆虫发育和免疫功能的蜕皮激素(Ec)和保幼激素(JH)在相互拮抗过程中能够通过共同作用激活体内MAPK信号途径，进而反式调控小菜蛾中肠内氨基肽酶N等多个受体基因的表达，并诱导宿主小菜蛾在不影响自身发育的前提下对苏云金芽胞杆菌(Bacillus thuringiensis)杀虫蛋白Cry1Ac产生高抗性。本研究中，lncRNA (TCONS_00032421)在AcT中下调表达[log2 (Fold Change)=–0.36，P=3.62E-06]且能反向结合注释在MAPK信号通路上的丝氨酸/苏氨酸蛋白激酶编码基因(XM_017064272)，推测TCONS_00032421可作为antisense lncRNA参与调控MAPK信号通路，进而介导宿主对球囊菌侵染的应答。
LncRNA已被证实可通过cis作用调控上下游基因的表达，进而影响昆虫对病原侵染的响应[16, 21]。Zhang等[16]研究发现家蚕核型多角体病毒在感染家蚕前后，宿主DElncRNA可能通过cis作用影响蛋白质转运、ncRNA代谢过程、细胞应激反应、泛素介导的蛋白水解、内吞作用和溶酶体等通路相关的上下游基因，从而抵御病毒侵染。Chen等[21]研究发现意蜂工蜂中肠的DElncRNA可能通过cis作用调控上下游基因的表达，间接影响宿主的物质和能量代谢及免疫途径以响应东方蜜蜂微孢子虫的侵染。本研究中，AcCK vs. AcT比较组的194个DElncRNA可通过cis作用调控461个上下游基因的表达，进而影响细胞、结合、催化活性、代谢进程、应激反应、发育进程及细胞杀伤等功能条目，碳水化合物代谢、氨基酸代谢和脂质代谢相关的诸多物质代谢通路，以及Hippo、Notch、FoxO、mTOR和Wnt等与昆虫生长发育密切相关的信号通路。以上结果表明球囊菌侵染通过改变lncRNA的表达水平对中蜂幼虫肠道的细胞进程、新陈代谢以及生长发育造成广泛影响。此外，分别有4、2、2、1和1个上下游基因涉及细胞凋亡、内吞作用、MAPK信号通路、Toll-like受体信号通路及泛素介导的蛋白水解等细胞和体液免疫途径，表明宿主的DElncRNA还可通过cis作用调控相关上下游基因的表达调节上述免疫途径，从而对球囊菌侵染产生免疫应答。但相关机制仍需进一步的深入研究。
近年来，随着实验技术手段的进步和研究的深入，ceRNA机制已在人类和一些模式生物中证实[42-43]。此外，已有一些研究报道lncRNA可作为ceRNA参与调节昆虫的抗药性、发育、生殖和免疫应答等过程[12, 14, 18]。Feng等[12]研究发现lincRNA_Tc13743.2通过竞争性结合miR-133-5p促进谷胱甘肽S-转移酶基因TcGSTm02的过表达，从而介导了朱砂叶螨(Tetranychus cinnabarinus)对丁氟螨酯的抗药性。Etebari等[18]研究发现埃及伊蚊被登革病毒感染后lincRNA_1317显著上调表达，并可作为ceRNA结合miR-1、miR-9a和miR-33等miRNA，参与调节宿主的抗病毒免疫反应。前期研究中，笔者课题组发现在意蜂工蜂中肠发育过程中，lncRNA可能作为ceRNA竞争性结合miRNA，进而影响物质和能量代谢及细胞和体液免疫[27, 44]。本研究中，宿主的180个DElncRNA与50个DEmiRNA之间存在靶向结合关系较为复杂，其中TCONS_00015199、TCONS_00030368和TCONS_00007304等DElncRNA靶向的DEmiRNA数量较多，分别达到13、8和7个。进一步分析发现DElncRNA可通过吸附结合DEmiRNA间接影响mRNA的表达，从而对Jak-STAT信号通路以及氧化磷酸化等与免疫和代谢相关通路进行调控。Jak-STAT信号通路作为昆虫体液免疫系统的重要组成部分，可参与西方蜜蜂[40]、果蝇[45]和家蚕[46]等对病原微生物的免疫应答。前期研究结果表明，miR-101-y和miR-21-x在中蜂6日龄幼虫肠道应答球囊菌侵染过程中显著上调表达[log2(Fold Change)=12.42，P=0.00015；log2(Fold Change)=11.96，P=4.70E-05]，并作为核心同时调控Jak-STAT信号通路和泛素介导的蛋白水解等免疫通路相关靶mRNA[23]。本研究中，上调表达的miR-101-y同时与被下调的lncRNA (TCONS_00010661)和mRNA (XM_017064778.1)靶向结合；上调表达的miR-21-x同时被2个下调的lncRNA (TCONS_00002228和TCONS_00026779)和1个下调的mRNA (XM_017053215.1)靶向结合。根据上述结果推测球囊菌通过与宿主的互作下调宿主肠道内TCONS_00010661、TCONS_00002228和TCONS_00026779的表达水平，减轻对miR-101-y和miR-21-x的结合，从而加强对Jak-STAT信号通路相关基因的抑制，进而促进自身的增殖。氧化磷酸化作为真核生物的主要能量代谢方式，在昆虫体内的能量代谢中同样居于核心地位[47]。前期研究中，笔者课题组通过比较转录组学分析揭示富集在氧化磷酸化通路的2个中蜂单拷贝基因受到强烈正向选择，中蜂幼虫可能通过提升能量利用限制病原增殖[6]。在中蜂6日龄幼虫肠道对球囊菌侵染的应答中，miR-1277-x显著下调表达[log2(Fold Change)= –7.82，P=1.35E-06]且参与调控Jak-STAT信号通路和细胞凋亡等免疫途径，同时与其他下调miRNA共同靶向结合氧化磷酸化通路相关的38个mRNA[23]。本研究发现在中蜂6日龄幼虫应答球囊菌侵染中，存在1对下调lncRNA与上调的miRNA (TCONS_00010661和miR-101-y)、28对上调lncRNA与下调的miRNA (如TCONS_00003104和miR-1277-x，TCONS_00032519和novel-m0003-3p以及TCONS_00032519和novel-m0035-3p等)的靶向结合关系，暗示宿主通过上调表达更多的lncRNA吸附结合更多的miRNA，进一步减弱对氧化磷酸化途径的抑制作用，以提高能量利用来限制球囊菌增殖，与前期研究结果一致；此外，球囊菌还可能通过下调TCONS_0001066的表达水平减弱对miR-101-y的吸附结合，一定程度上增强对氧化磷酸化途径的抑制，反映出中蜂幼虫与球囊菌相互作用的复杂性。
昆虫-病原的互作十分复杂，过程中伴随着ncRNA和mRNA之间复杂的相互作用。前期研究中，笔者团队通过解析中蜂6日龄幼虫响应球囊菌侵染的基因差异表达谱，发现部分上调基因广泛富集在氨基酸生物合成、碳代谢和酪氨酸代谢等多条物质代谢通路，细胞凋亡、溶酶体和内吞作用等细胞免疫途径，以及Toll-like和Jak-STAT信号通路等体液免疫途径，表明在球囊菌侵染后期，宿主的物质代谢和免疫相关通路被激活[22]。进一步研究发现，对于球囊菌侵染的中蜂6日龄幼虫，宿主的54个DEmiRNA及其靶mRNA存在靶向调控关系，miR-1277-x、miR-101-y、miR-21-x和miR-6052-x等差异表达且可能通过调控泛素介导的蛋白水解和Jak-STAT信号通路等相关基因参与对肠道免疫的调节[23]。近期，我们发现不同长度和环化类型的circRNA在中蜂6日龄幼虫响应球囊菌胁迫的过程中发生数量变化，并通过调控来源基因表达和作为ceRNA间接调控靶基因表达，参与调控宿主的一些物质和能量代谢通路以及细胞和体液免疫途径[24]。CeRNA机制认为凡具有miRNA反应元件的RNA如circRNA和lncRNA，都可以竞争性结合相同或类似的miRNA，从而间接调控基因表达[48]。本研究发现DElncRNA可能作为ceRNA，通过吸附结合miR-101-y、miR-1277-x和novel-m0003-3p等调控宿主侵染应答中的氧化磷酸化；通过吸附结合miR-101-y、miR-21-x和miR-1277-x等调控宿主的Jak-STAT信号通路。综上，系统的分析结果表明中蜂6日龄幼虫可通过差异表达lncRNA、circRNA、miRNA和mRNA对球囊菌侵染产生应答，lncRNA和circRNA可通过ceRNA机制经miRNA介导对下游的靶基因进行表达调控。目前，来源于蜜蜂组织或器官的可传代培养细胞系仍然缺失，成熟的转基因操作技术体系尚未建立。人工合成miRNA mimic和inhibitor，通过饲喂或注射对蜜蜂个体的miRNA进行过表达和敲减已被证明可行且有效[49-50]。下一步将针对免疫相关的关键DEmiRNA进行过表达和敲减，进而探究ceRNA调控网络介导中蜂幼虫应答球囊菌侵染的分子机制。
在我国的养蜂生产中，中蜂和意蜂为两个常用蜂种，其中前者很少受到白垩病的影响，而后者常暴发白垩病，体现出二者在群体水平具有不同的球囊菌抗性。笔者团队前期研究发现，中蜂幼虫和意蜂幼虫在个体水平表现出不同的细胞和体液免疫应答[4-5]；二者的细胞和体液免疫相关的特有基因的数量存在差异，中蜂幼虫的氧化磷酸化相关的2个单拷贝基因受到强烈正选择[6]。以上结果揭示中蜂幼虫和意蜂幼虫在个体水平同样表现出球囊菌抗性差异。LncRNA和其他ncRNA如miRNA和circRNA是否参与决定中蜂幼虫和意蜂幼虫的球囊菌抗性差异以及通过何种机制发挥作用？这是一个非常有趣且重要的科学问题。目前，笔者团队已制备球囊菌感染的意蜂6日龄幼虫肠道及未感染的6日龄幼虫肠道，下一步将利用链特异性cDNA建库的RNA-seq技术对上述肠道样品进行深度测序，进而对DElncRNA的调控方式和潜在作用进行深入分析，并结合本研究的结果进行比较和探讨。
综上所述，本研究首次在全基因组水平上对中蜂幼虫肠道的lncRNA进行了鉴定和结构特征分析，通过深入细致的组学解析揭示部分lncRNA在宿主应答球囊菌侵染的过程中可能作为antisense lncRNA反向结合并调控mrjp3，通过顺式调控上下游基因调节宿主的新陈代谢和免疫防御途径；一些DElncRNA还可作为ceRNA竞争性结合miR-1277-x等DEmiRNA，间接调控氧化磷酸化通路和Jak-STAT信号通路，从而介导宿主对球囊菌侵染的应答。
References
| [1] | Lin ZG, Page P, Li L, Qin Y, Zhang YY, Hu FL, Neumann P, Zheng HQ, Dietemann V. Go east for better honey bee health: Apis cerana is faster at hygienic behavior than A. mellifera. PLoS One, 2016, 11(9): e0162647. DOI:10.1371/journal.pone.0162647 |
| [2] | 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. DOI:10.1080/00218839.2018.1475943 |
| [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. DOI:10.1371/journal.pone.0124868 |
| [4] | 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. DOI:10.1016/j.gene.2017.04.022 |
| [5] | 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. DOI:10.1016/j.jip.2019.107210 |
| [6] | Xiong CL, Du Y, Wang HQ, Zheng YZ, Fu ZM, Wang HP, Zhang L, Chen DF, Guo R. Unraveling the mechanism regulating the Ascosphaera apis-resistance difference between Apis cerana cerana and Apis mellifera ligustica larvae based on comparative transcriptome analysis. Journal of China Agricultural University, 2019, 24(5): 106-114. (in Chinese) 熊翠玲, 杜宇, 王鸿权, 郑燕珍, 付中民, 王海朋, 张璐, 陈大福, 郭睿. 基于比较转录组学分析揭示中华蜜蜂及意大利蜜蜂幼虫的球囊菌抗性差异机制. 中国农业大学学报, 2019, 24(5): 106-114. |
| [7] | Aronstein KA, Murray KD. Chalkbrood disease in honey bees. Journal of Invertebrate Pathology, 2010, 103(Suppl 1): S20-S29. |
| [8] | Fatica A, Bozzoni I. Long non-coding RNAs: new players in cell differentiation and development. Nature Reviews Genetics, 2014, 15(1): 7-21. DOI:10.1038/nrg3606 |
| [9] | Wang TM, Qu LN, Li YH. Structures and functions of long non-coding RNAs and its roles in diseases. Chinese Journal of Biochemistry and Molecular Biology, 2015, 31(7): 659-666. (in Chinese) 王婷梅, 曲丽娜, 李莹辉. LncRNA的结构、功能及其与疾病的关系. 中国生物化学与分子生物学报, 2015, 31(7): 659-666. |
| [10] | Atianand MK, Fitzgerald KA. Long non-coding RNAs and control of gene expression in the immune system. Trends in Molecular Medicine, 2014, 20(11): 623-631. DOI:10.1016/j.molmed.2014.09.002 |
| [11] | Kopp F, Mendell JT. Functional classification and experimental dissection of long noncoding RNAs. Cell, 2018, 172(3): 393-407. DOI:10.1016/j.cell.2018.01.011 |
| [12] | Feng KY, Liu J, Wei P, Ou SY, Wen X, Shen GM, Xu ZF, Xu Q, He L. lincRNA_Tc13743.2-miR-133-5p-TcGSTm02 regulation pathway mediates cyflumetofen resistance in Tetranychus cinnabarinus. Insect Biochemistry and Molecular Biology, 2020, 123: 103413. DOI:10.1016/j.ibmb.2020.103413 |
| [13] | Jayakodi M, Jung JW, Park D, Ahn YJ, Lee SC, Shin SY, Shin C, Yang TJ, Kwon HW. Genome-wide characterization of long intergenic non-coding RNAs (lincRNAs) provides new insight into viral diseases in honey bees Apis cerana and Apis mellifera. BMC Genomics, 2015, 16(1): 680. DOI:10.1186/s12864-015-1868-7 |
| [14] | Chen X, Ma C, Chen C, Lu Q, Shi W, Liu ZG, Wang HH, Guo HK. Integration of lncRNA-miRNA-mRNA reveals novel insights into oviposition regulation in honey bees. PeerJ, 2017, 5: e3881. DOI:10.7717/peerj.3881 |
| [15] | Valanne S, Salminen TS, J?rvel?-St?lting M, Vesala L, R?met M. Immune-inducible non-coding RNA molecule lincRNA-IBIN connects immunity and metabolism in Drosophila melanogaster. PLoS Pathogens, 2019, 15(1): e1007504. DOI:10.1371/journal.ppat.1007504 |
| [16] | Zhang SL, Yin HT, Shen MM, Huang HL, Hou QR, Zhang ZD, Zhao WG, Guo XJ, Wu P. Analysis of lncRNA-mediated gene regulatory network of Bombyx mori in response to BmNPV infection. Journal of Invertebrate Pathology, 2020, 170: 107323. DOI:10.1016/j.jip.2020.107323 |
| [17] | Xiao HM, Yuan ZT, Guo DH, Hou BF, Yin CL, Zhang WQ, Li F. Genome-wide identification of long noncoding RNA genes and their potential association with fecundity and virulence in rice brown planthopper, Nilaparvata lugens. BMC Genomics, 2015, 16(1): 749. DOI:10.1186/s12864-015-1953-y |
| [18] | Etebari K, Asad S, Zhang GM, Asgari S. Identification of Aedes aegypti long intergenic non-coding RNAs and their association with Wolbachia and Dengue virus infection. PLoS Neglected Tropical Diseases, 2016, 10(10): e0005069. DOI:10.1371/journal.pntd.0005069 |
| [19] | Guo Y, Su SK, Chen SL, Zhang SW, Chen RS. The function of lncRNAs in the caste determination of the honeybee. Progress in Biochemistry and Biophysics, 2015, 42(8): 750-757. (in Chinese) 郭昱, 苏松坤, 陈盛禄, 张少吾, 陈润生. LncRNA在蜜蜂级型分化中的功能研究. 生物化学与生物物理进展, 2015, 42(8): 750-757. |
| [20] | Liu F, Shi TF, Qi L, Su X, Wang DQ, Dong J, Huang ZY. LncRNA profile of Apis mellifera and its possible role in behavioural transition from nurses to foragers. BMC Genomics, 2019, 20(1): 393. DOI:10.1186/s12864-019-5664-7 |
| [21] | Chen DF, Chen HZ, Du Y, Zhou DD, Geng SH, Wang HP, Wan JQ, Xiong CL, Zheng YZ, Guo R. Genome-wide identification of long non-coding RNAs and their regulatory networks involved in Apis mellifera ligustica response to Nosema ceranae infection. Insects, 2019, 10(8): 245. DOI:10.3390/insects10080245 |
| [22] | Guo R, Zhang L, Xu XJ, Shi XL, Xiong CL, Zheng YZ, Fu ZM, Huang ZJ, Wang HQ, Hou ZX, Chen DF. Analysis of the differentially expressed genes in the 6-day-old larval gut of Apis cerana cerana under the stress of Ascosphaera apis. Journal of Environmental Entomology, 2017, 39(3): 539-547. (in Chinese) 郭睿, 张璐, 徐细建, 史秀丽, 熊翠玲, 郑燕珍, 付中民, 黄枳腱, 王鸿权, 侯志贤, 陈大福. 中华蜜蜂6日龄幼虫肠道响应球囊菌胁迫的差异表达基因分析. 环境昆虫学报, 2017, 39(3): 539-547. |
| [23] | 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. |
| [24] | Chen HZ, Fu ZM, Wang J, Zhu ZW, Fan XX, Jiang HB, Fan YC, Zhou DD, Li WD, Xiong CL, Zheng YZ, Xu GJ, Chen DF, Guo R. Circular RNA response of Apis cerana cerana 6-day-old larvae to Ascosphaera apis stress. Acta Microbiologica Sinica, 2020, 60(10): 2292-2310. (in Chinese) 陈华枝, 付中民, 王杰, 祝智威, 范小雪, 蒋海宾, 范元婵, 周丁丁, 李汶东, 熊翠玲, 郑燕珍, 徐国钧, 陈大福, 郭睿. 中华蜜蜂6日龄幼虫响应蜜蜂球囊菌胁迫的环状RNA应答. 微生物学报, 2020, 60(10): 2292-2310. |
| [25] | Xiong CL, Du Y, Feng RR, Jiang HB, Shi XY, Wang HP, Fan XX, Wang J, Zhu ZW, Fan YC, Chen HZ, Zhou DD, Zheng YZ, Chen DF, Guo R. Differential expression pattern and regulation network of microRNAs in Ascosphaera apis invading Apis cerana cerana 6-day-old larvae. Acta Microbiologica Sinica, 2020, 60(5): 992-1009. (in Chinese) 熊翠玲, 杜宇, 冯睿蓉, 蒋海宾, 史小玉, 王海朋, 范小雪, 王杰, 祝智威, 范元婵, 陈华枝, 周丁丁, 郑燕珍, 陈大福, 郭睿. 侵染中华蜜蜂6日龄幼虫的蜜蜂球囊菌的微小RNA差异表达谱及调控网络. 微生物学报, 2020, 60(5): 992-1009. |
| [26] | 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. |
| [27] | Guo R, Geng SH, Xiong CL, Zheng YZ, Fu ZM, Wang HP, Du Y, Tong XY, Zhao HX, Chen DF. Differential expression analysis of long non-coding RNAs during the developmental process of Apis mellifera ligustica worker's midgut. Scientia Agricultura Sinica, 2018, 51(18): 3600-3613. (in Chinese) 郭睿, 耿四海, 熊翠玲, 郑燕珍, 付中民, 王海朋, 杜宇, 童新宇, 赵红霞, 陈大福. 意大利蜜蜂工蜂中肠发育过程中长链非编码RNA的差异表达分析. 中国农业科学, 2018, 51(18): 3600-3613. DOI:10.3864/j.issn.0578-1752.2018.18.016 |
| [28] | Xiong CL, Geng SH, Wang XR, Liu SY, Chen DF, Zheng YZ, Fu ZM, Du Y, Wang HP, Chen HZ, Zhou DD, Guo R. Prediction, analysis and identification of long non-coding RNA in the midguts of Apis mellifera ligustica workers. Chinese Journal of Applied Entomology, 2018, 55(6): 1034-1044. (in Chinese) 熊翠玲, 耿四海, 王心蕊, 刘思亚, 陈大福, 郑燕珍, 付中民, 杜宇, 王海朋, 陈华枝, 周丁丁, 郭睿. 意大利蜜蜂工蜂中肠的长链非编码RNA的预测、分析及鉴定. 应用昆虫学报, 2018, 55(6): 1034-1044. |
| [29] | Carrieri C, Cimatti L, Biagioli M, Beugnet A, Zucchelli S, Fedele S, Pesce E, Ferrer I, Collavin L, Santoro C, Forrest A R R, Carninci P, Biffo S, Stupka E, Gustincich S. Long non-coding antisense RNA controls Uchl1 translation through an embedded SINEB2 repeat. Nature, 2012, 491(7424): 454-457. DOI:10.1038/nature11508 |
| [30] | Rinn JL, Chang HY. Genome regulation by long non-coding RNAs. Annual Review of Biochemistry, 2012, 81: 145-166. DOI:10.1146/annurev-biochem-051410-092902 |
| [31] | Livak KJ, Schmittgen TD. Analysis of relative gene expression data using real-time quantitative PCR and the method. Methods, 2001, 25(4): 402-408. DOI:10.1006/meth.2001.1262 |
| [32] | Shang YF, Xiao GH, Zheng P, Cen K, Zhan S, Wang CS. Divergent and convergent evolution of fungal pathogenicity. Genome Biology and Evolution, 2016, 8(5): 1374-1387. DOI:10.1093/gbe/evw082 |
| [33] | Park D, Jung JW, Choi BS, Jayakodi M, Lee J, Lim J, Yu Y, Choi YS, Lee ML, Park Y, Choi IY, Yang TJ, Edwards O R, Nah G, Kwon HW. Uncovering the novel characteristics of Asian honey bee, Apis cerana, by whole genome sequencing. BMC Genomics, 2015, 16(1): 1. DOI:10.1186/1471-2164-16-1 |
| [34] | Zhu B, Xu MY, Shi HY, Gao XW, Liang P. Genome-wide identification of lncRNAs associated with chlorantraniliprole resistance in diamondback moth Plutella xylostella (L.). BMC Genomics, 2017, 18(1): 380. DOI:10.1186/s12864-017-3748-9 |
| [35] | 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. DOI:10.1016/j.jip.2018.06.001 |
| [36] | Guo R, Chen DF, Xiong CL, Hou CS, Zheng YZ, Fu ZM, Liang Q, Diao QY, Zhang L, Wang HQ, Hou ZX, Kumar D. First identification of long non-coding RNAs in fungal parasite Nosema ceranae. Apidologie, 2018, 49(5): 660-670. DOI:10.1007/s13592-018-0593-z |
| [37] | Amit-Avraham I, Pozner G, Eshar S, Fastman Y, Kolevzon N, Yavin E, Dzikowski R. Antisense long non-coding RNAs regulate var gene activation in the malaria parasite Plasmodium falciparum. Proceedings of the National Academy of Sciences of the United States of America, 2015, 112(9): E982-E991. DOI:10.1073/pnas.1420855112 |
| [38] | Albertová V, Su SK, Brockmann A, Gadau J, Albert S. Organization and potential function of the mrjp3 locus in four honeybee species. Journal of Agricultural and Food Chemistry, 2005, 53(20): 8075-8081. DOI:10.1021/jf051417x |
| [39] | Peixoto LG, Calábria LK, Garcia L, Capparelli FE, Goulart LR, de Sousa MV, Espindola FS. Identification of major royal jelly proteins in the brain of the honeybee Apis mellifera. Journal of Insect Physiology, 2009, 55(8): 671-677. DOI:10.1016/j.jinsphys.2009.05.005 |
| [40] | 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. DOI:10.1111/j.1365-2583.2006.00682.x |
| [41] | Guo ZJ, Kang S, Sun D, Gong LJ, Zhou JL, Qin JY, Guo L, Zhu LH, Bai Y, Ye F, Wu QJ, Wang SL, Crickmore N, Zhou XG, Zhang YJ. MAPK-dependent hormonal signaling plasticity contributes to overcoming Bacillus thuringiensis toxin action in an insect host. Nature Communications, 2020, 11(1): 3003. DOI:10.1038/s41467-020-16608-8 |
| [42] | Cheng EC, Lin HF. Repressing the repressor: a lincRNA as a microRNA sponge in embryonic stem cell self-renewal. Developmental Cell, 2013, 25(1): 1-2. DOI:10.1016/j.devcel.2013.03.020 |
| [43] | Xie CR, Wang F, Zhang S, Wang FQ, Zheng S, Li Z, Lv J, Qi HQ, Fang QL, Wang XM, Yin ZY. Long non-coding RNA HCAL facilitates the growth and metastasis of hepatocellular carcinoma by acting as a ceRNA of LAPTM4B. Molecular Therapy: Nucleic Acids, 2017, 9: 440-451. DOI:10.1016/j.omtn.2017.10.018 |
| [44] | Guo R, Chen HZ, Xiong CL, Zheng YZ, Fu ZM, Xu GJ, Du Y, Wang HP, Geng SH, Zhou DD, Liu SY, Chen DF. Analysis of differentially expressed circular RNAs and their regulation networks during the developmental process of Apis mellifera ligustica worker's midgut. Scientia Agricultura Sinica, 2018, 51(23): 4575-4590. (in Chinese) 郭睿, 陈华枝, 熊翠玲, 郑燕珍, 付中民, 徐国钧, 杜宇, 王海朋, 耿四海, 周丁丁, 刘思亚, 陈大福. 意大利蜜蜂工蜂中肠发育过程中的差异表达环状RNA及其调控网络分析. 中国农业科学, 2018, 51(23): 4575-4590. DOI:10.3864/j.issn.0578-1752.2018.23.015 |
| [45] | Myllym?ki H, R?met M. JAK/STAT pathway in Drosophila immunity. Scandinavian Journal of Immunology, 2014, 79(6): 377-385. DOI:10.1111/sji.12170 |
| [46] | Geng T, Lv DD, Huang YX, Hou CX, Qin GX, Guo XJ. JAK/STAT signaling pathway-mediated immune response in silkworm (Bombyx mori) challenged by Beauveria bassiana. Gene, 2016, 595(1): 69-76. DOI:10.1016/j.gene.2016.09.043 |
| [47] | Li YY, Zhang R, Liu SL, Donath A, Peters RS, Ware J, Misof B, Niehuis O, Pfrender ME, Zhou X. The molecular evolutionary dynamics of oxidative phosphorylation (OXPHOS) genes in Hymenoptera. BMC Evolutionary Biology, 2017, 17(1): 269. DOI:10.1186/s12862-017-1111-z |
| [48] | Thomson DW, Dinger ME. Endogenous microRNA sponges: evidence and controversy. Nature Reviews Genetics, 2016, 17(5): 272-283. DOI:10.1038/nrg.2016.20 |
| [49] | Liu F, Shi TF, Yin W, Su X, Qi L, Huang ZY, Zhang SW, Yu LS. The microRNA ame-miR-279a regulates sucrose responsiveness of forager honey bees (Apis mellifera). Insect Biochemistry and Molecular Biology, 2017, 90: 34-42. DOI:10.1016/j.ibmb.2017.09.008 |
| [50] | 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. DOI:10.1371/journal.pgen.1006946 |
