东方蜜蜂微孢子虫孢子中微小RNA的鉴定与分析

张文德, 赵浩东, 孙明会, 余岢骏, 郭意龙, 朱乐冉,胡 颖, 赵 萧, 叶亚萍, 陈大福,2,*, 郭 睿,2,*

(1. 福建农林大学动物科学学院(蜂学学院),福州 350002; 2. 福建农林大学蜂疗研究所, 福州 350002)

东方蜜蜂微孢子虫Nosemaceranae是一种专性寄生蜜蜂中肠上皮细胞的单细胞真菌病原，最早在东方蜜蜂Apiscerana上被鉴定到，随后跨种感染西方蜜蜂Apismellifera，随着蜂产品贸易全球化和跨国引种日益频繁已传播至各养蜂国家(Papinietal., 2017)。该病原不仅对蜜蜂宿主的食物消化、营养吸收、免疫应答、能量代谢和细胞凋亡造成全方位的影响，还能导致蜜蜂的寿命缩短、生产力下降和采集日龄提前，严重危害蜜蜂健康和养蜂生产(Martín-Hernándezetal., 2011)。Cornman等(2009)利用第二代焦磷酸测序技术组装并公布了东方蜜蜂微孢子虫的参考基因组，为其组学及分子生物学研究奠定了基础。近期，Huang等(2021)利用第三代纳米孔(nanopore)长读段测序技术重新组装了东方蜜蜂微孢子虫基因组，基因组大小约为8.8 Mb，包含2 280个蛋白编码基因。

微小RNA(microRNA, miRNA)是一类长度约为19～24 nt的细胞内源性小非编码RNA(non-coding RNA, ncRNA)，在真核生物中广泛存在且高度保守。miRNA能通过切割靶mRNA或阻止蛋白翻译调控基因表达，进而参与调节胚胎发育﹑细胞分化和增殖、神经发生、能量代谢、细胞凋亡及昼夜节律等诸多生物学进程(Legeaietal., 2010)。目前，miRNA在小鼠Musmusculus、拟南芥Arabidopsisthaliana和黑腹果蝇Drosophilamelanogaster等模式生物中的研究较为深入(Babiarzetal., 2008; Enderetal., 2008; Zhan and Lukens, 2010)。但真菌miRNA的相关研究较为滞后，蜜蜂真菌病原的miRNA研究更为有限。前期研究中，笔者所在团队对蜜蜂的另一种主要真菌病原蜜蜂球囊菌Ascosphaeraapis的miRNA进行转录组测序和分析(郭睿等, 2018)，揭示了差异表达miRNA(differentially expressed miRNA, DEmiRNA)在菌丝生长、孢子萌发和有性生殖等方面的潜在作用(陈华枝等, 2020)，并全面解析了蜜蜂球囊菌在侵染意大利蜜蜂Apismelliferaligustica幼虫和中华蜜蜂Apisceranacerana幼虫发育过程的miRNA差异表达谱及DEmiRNA的调控作用(熊翠玲等, 2020; 杜宇, 2021)。此前，Huang和Evans(2016)利用二代测序技术对东方蜜蜂微孢子虫感染的西方蜜蜂工蜂中肠进行测序，并结合RT-qPCR鉴定到病原的6个miRNA-like RNA(milRNA)。近期，Shao等(2021)基于更新的东方蜜蜂微孢子虫基因组版本(Huangetal., 2021)又鉴定到3个新miRNA，分析发现3个新miRNA参与了病原增殖过程的自我调控。然而，东方蜜蜂微孢子虫的miRNA研究仍比较有限，对于病原孢子中是否存在miRNA以及miRNA发挥何种功能，相关研究仍然缺失。

本研究基于已获得的small RNA-seq(sRNA-seq)数据对东方蜜蜂微孢子虫孢子中的miRNA进行鉴定和分析，并利用分子生物学手段验证miRNA的表达和序列真实性，以期明确东方蜜蜂微孢子虫孢子中miRNA的数量、结构特征及潜在调控作用，丰富东方蜜蜂微孢子虫的miRNA信息，并为进一步探究miRNA在病原孢子及病原致病过程中的功能提供理论和实验依据。

1 材料与方法

1.1 东方蜜蜂微孢子虫孢子的sRNA-seq数据来源

笔者所在团队前期构建了东方蜜蜂微孢子虫纯净孢子3份平行样品Nc-1, Nc-2和Nc-3的cDNA文库，并委托北京百迈克生物科技有限公司进行单端测序，测序平台为Illumina HiSeq Xten。测序得到的原始数据已上传NCBI SRA数据库，BioProject号为PRJNA395137。 Nc-1, Nc-2和Nc-3的sRNA-seq分别产生16 597 883, 15 451 791和12 248 316条raw reads，经严格质控分别得到15 608 370, 14 249 255和11 440 684条clean reads，Q30分别达到98.61%, 98.65%和98.58%，分别有7 717 280, 7 801 135和6 482 422条clean reads能比对到东方蜜蜂微孢子虫参考基因组(周丁丁等, 2020)，说明测序数据质量较高，可满足本研究的相关分析。

1.2 sRNA分类注释及参考基因组比对

为获得包含miRNA的未注释读段(unannotated reads)，使用Bowtie软件(Langmeadetal., 2009)将clean reads分别比对至Silva(https:∥www.arb-silva.de/)、GtRNAdb(http:∥gtrnadb.ucsc.edu/)、Rfam(http:∥rfam.xfam.org/)和Repbase(https:∥www.girinst.org/repbase/)数据库，滤除核糖体RNA(rRNA)、转运RNA(tRNA)、核内小RNA(snRNA)、核仁小RNA(snoRNA)等ncRNA及重复序列。进一步采用Bowtie软件将上述未注释读段与东方蜜蜂微孢子虫参考基因组(https:∥www.ncbi.nlm.nih.gov/genome/?term=nosema%20ceranae)进行序列比对以获取未注释读段在参考基因组上的位置信息，从而得到比对上读段(mapped reads)。

1.3 miRNA鉴定

采用miRDeep2(Friedländeretal., 2012)软件包将1.2节的mapped reads与miRBase数据库(http:∥www.mirbase.org/)中的已知miRNA前体序列进行比对，以鉴定保守miRNA的表达。同时，通过reads比对到基因组上的位置信息得到可能的前体序列，基于reads在miRNA前体(mature, star和loop)上的分布信息及前体结构能量信息(RNAfold randfold)采用贝叶斯模型经打分最终实现新miRNA的鉴定。进一步对鉴定到的miRNA进行结构特征分析，包括长度分布和首位及每位碱基偏向性。采用TPM(transcripts per million)法(Fahlgrenetal., 2007)对鉴定到的miRNA进行表达量计算。

1.4 miRNA的茎环反转录PCR(stem-loop RT-PCR)检测和Sanger测序验证

根据1.3节鉴定到的10个miRNA的核酸序列，利用DNAMAN软件(https:∥www.lynnon.com/dnaman.html)设计特异性Stem-loop引物和上游引物(F)及通用下游引物(R)，委托上海生工生物工程股份有限公司进行合成(表1)。东方蜜蜂微孢子虫感染的意大利蜜蜂外勤蜂取自福州市闽侯县荆溪源安养蜂场。基于笔者所在团队前期已建立东方蜜蜂微孢子虫孢子纯化的技术流程(Chenetal., 2019; 耿四海, 2020)制备东方蜜蜂微孢子虫的纯净孢子，利用RNA抽提试剂盒(Promega，美国)提取东方蜜蜂微孢子虫孢子样品的总RNA，按照cDNA第1链合成试剂盒(诺唯赞，南京)说明书进行反转录，得到的cDNA作为模板进行常规PCR。PCR体系(20 μL): cDNA模板1 μL, 上下游引物(2.5 pmol/μL)各1 μL, PCR Mix 10 μL, 无菌水7 μL。PCR程序: 95℃ 5 min; 95℃ 50 s, 50℃ 30 s, 72℃ 50 s 34个循环; 72℃ 10 min。PCR产物经2%的琼脂糖凝胶电泳检测，目的片段经切胶回收和分子克隆后送至上海生工生物工程股份有限公司进行Sanger测序。

1.5 miRNA靶基因预测、分析及调控网络构建

根据预测到的miRNA与东方蜜蜂微孢子虫的基因序列信息，利用TargetFinder(Allenetal., 2005)软件对1.3节鉴定到的miRNA进行靶向预测，得到靶基因集合。再利用BLAST软件将靶基因分别比对到Nr(ftp:∥ftp.ncbi.nih.gov/blast/db), Swiss-Prot(http:∥www.uniprot.org/uniprot/?query=*&fil=reviewed%3Ayes), eggNOG(http:∥eggnogdb.embl.de/#/app/emapper)和KOG(ftp:∥ftp.ncbi.nih.gov/pub/COG/KOG/kyva)数据库。利用基迪奥云平台(https:∥www.omicshare.com/tools)对靶基因进行GO(http:∥www.geneontology.org/)和KEGG(http:∥www.genome.jp/kegg/pathway.html)数据库注释。根据miRNA与靶基因的靶向结合关系构建调控网络，并通过Cytoscape(https:∥cytoscape.org/)软件进行可视化。

2 结果

2.1 东方蜜蜂微孢子虫孢子中miRNA的鉴定与特征

共鉴定到10个miRNA，长度分布介于21～25 nt，不同长度miRNA的首位碱基偏向性相似，多为U(图1: A)。对于特定长度的miRNA，其每一位碱基的偏向性差异明显(图1: B)。分布在24 nt和25 nt的miRNA数量较多，分别为3个和4个，分布在21, 22和23 nt的miRNA均为1个(表1)。 这10个>miRNA的TPM值介于844 910.04～274.06之间；其中相对表达量最高的为NW_003313949.1_26675，最低的为NW_003312818.1_15325(表1)。

2.2 东方蜜蜂微孢子虫孢子中miRNA的检测和测序验证

利用stem-loop RT-PCR对已鉴定的10个miRNA进行检测，电泳结果显示10个miRNA均能扩增出符合预期的目的片段(图2: A)；进一步对其中2个miRNA的扩增片段进行分子克隆与Sanger测序，结果显示基因组比对鉴定的miRNA序列与测序结果一致(图2: B)。结果表明本研究鉴定到的miRNA真实表达和存在。

2.3 东方蜜蜂微孢子虫孢子中miRNA的靶基因

图2 东方蜜蜂微孢子虫孢子中miRNA的Stem-loopRT-PCR检测(A)和Sanger测序验证(B)Fig. 2 Stem-loop RT-PCR detection (A) and Sangersequencing verification (B) of miRNAs inNosema ceranae sporesMarker: DNA分子量标准DNA molecular weight marker. miRNA信息见表1。For information of miRNAs, see Table 1. 下同The same below.

靶向预测结果显示，已鉴定的10个miRNA共靶向249个基因，其中分别有249, 118和3个靶基因可注释到Nr, Swiss-Prot和eggNOG数据库；有136个靶基因可注释到KOG数据库(图3)；有134个靶基因可注释到GO数据库中的30个功能条目，包括细胞进程(cellular process)(77)和代谢进程(metabolism process)(77)等15个生物学进程相关条目，细胞(cell)(51)和细胞组件(cell part)(51)等6个细胞组分相关条目，结合(binding)(84)和催化活性(catalytic activity)(78)等9个分子功能相关条目(图4: A)；有71个靶基因可注释到KEGG数据库中的54条通路，包括代谢进程(17)、次生代谢产物的生物合成(biosynthesis of secondary metabolites)(9)和氨酰基-tRNA的生物合成(biosynthesis of aminoacyl-tRNA)(7)、细胞周期-酵母(cell cycle-yeast)(6)和蛋白酶体(proteosome)(5)等(图4: B)。

图3 东方蜜蜂微孢子虫孢子中miRNA的靶基因的KOG数据库注释Fig. 3 KOG database annotation of target genes of miRNAs in Nosema ceranae spores

图4 东方蜜蜂微孢子虫孢子中miRNA的靶基因的GO(A)和KEGG(B)数据库注释Fig. 4 GO (A) and KEGG (B) database annotation of target genes of miRNAs in Nosema ceranae spores

2.4 东方蜜蜂微孢子虫孢子中miRNA与靶基因的调控网络

进一步根据已鉴定的miRNA和靶基因之间的靶向关系构建调控网络，分析结果显示东方蜜蜂微孢子虫的miRNA及其靶基因之间可形成较复杂的调控网络，其中NW_003314208.1_34537的靶基因数量最多，达到151个；NW_003311987.1_11248的靶基因最少，仅2个；NW_003311645.1_9769, NW_003312821.1_15338, NW_003312212.1_12220, NW_003313799.1_23928, NW_003313949.1_26675, NW_003312818.1_15325, NW_003311862.1_10660和NW_003311484.1_9172分别靶向28, 16, 15, 15, 13, 11, 9和5个基因(图5)。

图5 东方蜜蜂微孢子虫孢子中miRNA及其靶基因的调控网络Fig. 5 Regulatory network between miRNAs and their target genes in Nosema ceranae spores

3 讨论

本研究基于前期获得的东方蜜蜂微孢子虫纯净孢子的sRNA-seq测序数据，利用生物信息学方法共鉴定到10个miRNA，Stem-loop RT-PCR检测结果显示它们均真实表达，进一步的分子克隆与Sanger测序证实了其中2个miRNA的序列真实性(图2)。鉴于目前东方蜜蜂微孢子虫的研究较为缺乏，鉴定到的miRNA丰富了东方蜜蜂微孢子虫的miRNA信息，为进一步探究miRNA在病原孢子和病原侵染中的功能提供了候选分子。本研究为东方蜜蜂微孢子虫孢子中miRNA的首例报道，研究结果为进一步探究孢子中miRNA对基因表达的调控作用及对生命活动的调节作用提供了理论和实验依据。

结构特征分析结果显示东方蜜蜂微孢子虫的10个miRNA的首位碱基表现出U偏向性(图1: A)，与欧洲油菜Brassicanapus(Fuetal., 2019)、长牡蛎Crassostreagigas(王雪等, 2020)和家蚕微孢子虫Nosemabombycis(潘秋玲等, 2015)等物种的miRNA首位碱基偏向性一致。此外，特定长度miRNA的每一位碱基的偏向性表现出明显差异(图1: B)，这与其他动植物的研究报道(Fuetal., 2013; Saminathanetal., 2016)相似。以上结果表明不同物种的miRNA具有较强的结构保守性。此前，Huang和Evans(2016)以及Shao等(2021)共鉴定到东方蜜蜂微孢子虫的9个miRNA。通过比较发现，本研究鉴定到的10个miRNA与前人鉴定到的miRNA序列均不相同，这可能是由于本研究的测序数据来源于东方蜜蜂微孢子虫的孢子，而前人研究的测序数据主要来源于侵染西方蜜蜂工蜂状态下的东方蜜蜂微孢子虫。

本研究中，东方蜜蜂微孢子虫的10个miRNA共靶向249个基因，其中有118个靶基因可注释到Swiss-Prot数据库，注释数量最多的是染色体结构维持(structural maintenance of chromosomes, SMC)蛋白，涉及4个靶基因(gene2901, gene1851, gene1125和gene1649)及存在潜在靶向结合关系的3个miRNA(miRNA ID: NW_003311862.1_10660, NW_003314208.1_34537和NW_003312818.1_15325)(表1)。该蛋白与染色体结构的细胞周期性动态变化紧密相关，并参与有丝分裂染色体的集缩和分离，性染色体的剂量补偿效应，姐妹染色单体的内聚作用遗传重组以及DNA修复等过程(彭莉和张飞雄, 2001)。推测相应的miRNA通过调控相关基因表达影响SMC蛋白合成，进而参与维持东方蜜蜂微孢子虫孢子中的染色体结构。此外，有134个靶基因可注释到30个GO条目(图4: A)，包括应激反应(response to stimulus)、发育进程(developmental process)、生殖进程(reproductive process)、增殖(reproduction)和生长(growth)等。还发现有71个靶基因可注释到54条KEGG通路(图4: B)，包括代谢通路和次生代谢产物的生物合成等物质代谢通路，氨酰基-tRNA的生物合成和蛋白酶体等遗传信息处理通路，MAPK和磷脂酰肌醇(phosphatidylinositol)等信号通路。上述结果说明东方蜜蜂微孢子虫孢子中miRNA通过调控靶基因表达潜在参与生长发育、新陈代谢和信号转导等生物学进程。作为一种休眠态，真菌孢子中仅维持较低的代谢水平以维持基本的生命活动(Sephton-Clark and Voelz, 2018)。前期研究中，笔者所在团队在东方蜜蜂微孢子虫孢子中分别鉴定到83条长链非编码RNA(long non-coding RNA, lncRNA)和204条环状RNA(circular RNA, circRNA)，分析发现这些lncRNA和circRNA分别与4和3个miRNA存在靶向结合关系，并潜在通过竞争性内源RNA(competing endogenous RNA, ceRNA)调控网络影响诸多生物学过程(Guoetal., 2018a, 2018b)。以上结果表明东方蜜蜂微孢子虫孢子中存在活跃的基因转录活动，lncRNA和circRNA与miRNA可能通过ceRNA网络而密切互作，共同调节孢子中的生命活动。未来可联合蛋白组和翻译组技术进一步深入剖析东方蜜蜂微孢子虫孢子中的蛋白翻译及基因表达的表观调控机制。

微孢子虫是一类细胞内寄生的单细胞真菌病原，因缺乏线粒体，其增殖所需的物质和能量高度依赖宿主细胞提供(耿四海等, 2020)。 ABC转运蛋白和ATP/ADP移位酶在兔脑炎微孢子虫Encephalitozooncuniculi(Tsaousisetal., 2008)和家蚕微孢子虫(Panetal., 2017)等微孢子虫窃取宿主细胞中物质和能量过程中扮演关键角色。本研究发现，东方蜜蜂微孢子虫孢子中的2个miRNA(miRNA ID: NW_003313799.1_23928和NW_003314208.1_34537)可分别靶向ATP/ADP移位酶基因(gene1379)和ABC转运蛋白基因(gene2043)，说明这2个miRNA与东方蜜蜂微孢子虫的增殖与侵染具有潜在关联，值得进一步深入研究。鉴于目前来源于蜜蜂组织或器官的可传代培养的细胞系缺失，后续工作中，我们将根据已鉴定的miRNA的核酸序列设计合成相应的模拟物(mimic)和抑制物(inhibitor)，再通过饲喂东方蜜蜂微孢子虫感染的蜜蜂宿主对病原的miRNA进行过表达和敲减，以深入探究相关miRNA调控东方蜜蜂微孢子虫侵染的机制。

综上，本研究首次在东方蜜蜂微孢子虫孢子中鉴定到10个miRNA，并揭示它们的首位碱基具有U偏向性，且每一位碱基的偏向性差异明显，miRNA通过潜在靶向调控相关基因表达调节孢子中的细胞周期、代谢进程及次生代谢产物的生物合成等生物学过程。