东方蜜蜂微孢子虫孢子中长链非编码RNA的竞争性内源RNA调...

奥鹏网院作业

东方蜜蜂微孢子虫孢子中长链非编码RNA的竞争性内源RNA调控网络及潜在功能
周丁丁，史小玉，王杰，范元婵，祝智威，蒋海宾，范小雪，熊翠玲，郑燕珍，付中民，徐国钧，陈大福，郭睿

（福建农林大学动物科学学院（蜂学学院），福州 350002）

摘要：【背景】长链非编码RNA（long non-coding RNA，lncRNA）是一类长度＞200 nt且不具有编码能力的转录本，能够在真核生物的转录水平和转录后水平通过顺式（cis）、反式（trans）或竞争性内源RNA（competing endogenous RNA，ceRNA）机制调控基因的表达，已被证明在生物体的生长和发育等生物学过程中发挥重要的调控作用。【目的】结合测序得到的东方蜜蜂微孢子虫（Nosema ceranae）纯化孢子的small RNA（sRNA）组学数据和前期获得的lncRNA组学数据对lncRNA的调控方式进行深入细致的分析和探讨，揭示lncRNA在东方蜜蜂微孢子虫孢子中的潜在功能。【方法】利用small RNA-seq（sRNA-seq）技术对东方蜜蜂微孢子虫纯化孢子进行测序，通过相关的生物信息学软件对测序数据进行质控。根据lncRNA与mRNA的位置关系预测lncRNA的上下游基因；使用Blast软件比对GO和KEGG数据库，对上下游基因进行功能和代谢通路注释；利用Target Finder软件预测lncRNA-miRNA及lncRNA-miRNA-mRNA调控网络，进而通过Cytoscape v3.7.2软件对调控网络进行可视化；利用RT-PCR对调控网络中的部分lncRNA、miRNA和mRNA的表达进行验证。【结果】东方蜜蜂微孢子虫纯化孢子样品的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.58%。共预测出lncRNA的上下游基因310个，它们可注释到代谢进程、细胞进程、催化活性、结合和细胞等相关的35个功能条目；此外还可注释到56条代谢通路，涉及嘌呤代谢、碳代谢和丙酮酸代谢等物质代谢通路，以及甲烷代谢、氧化磷酸化和糖酵解/糖异生等能量代谢通路。上述结果表明相应的lncRNA可能通过cis作用调节上下游基因的表达，从而参与调控东方蜜蜂微孢子虫孢子中的物质和能量代谢以及细胞生命活动。进一步构建lncRNA-miRNA调控网络，分析结果显示MSTRG.3636.1、MSTRG.4498.1和 MSTRG.4883.1可靶向结合 4个 miRNA（nce-miR-7729、nce-miR-7502、nce-miR-8639和nce-miR-8565），说明此3个lncRNA可作为ceRNA在东方蜜蜂微孢子虫孢子中发挥潜在作用。LncRNA-miRNA-mRNA调控网络分析结果显示，三者之间形成较为复杂的调控网络关系，miRNA处于调控网络的中心并联系 lncRNA和mRNA，nce-miR-7502和 nce-miR-8639结合的 mRNA最多，达到 28个；与 MSTRG.3636.1、MSTRG.4498.1和MSTRG.4883.1存在靶向关系的miRNA结合多个mRNA，表明上述3个lncRNA可能通过ceRNA机制发挥作用，从而影响东方蜜蜂微孢子虫孢子中的新陈代谢和生命活动。【结论】lncRNA可能通过cis作用调控上下游基因的表达，以及作为ceRNA通过吸附miRNA间接影响靶基因的表达，从而参与调节东方蜜蜂微孢子虫孢子中的物质和能量代谢等基本生命活动。

关键词：东方蜜蜂微孢子虫；非编码RNA；长链非编码RNA；上下游基因；竞争性内源RNA；调控网络

0 引言
【研究意义】东方蜜蜂微孢子虫（Nosema ceranae）是一种常见的蜜蜂真菌病原，特异性侵染成年蜜蜂的中肠上皮细胞，导致宿主产生诸多生理和病理变化[1-3]。目前，有关东方蜜蜂微孢子虫非编码 RNA（non-coding RNA，ncRNA）的信息尤为有限[4-5]。利用组学技术探究东方蜜蜂微孢子虫长链非编码 RNA（long non-coding RNA，lncRNA）的不同作用方式及潜在功能，可为探明东方蜜蜂微孢子虫孢子发芽以及在宿主细胞内的增殖相关机制提供依据，筛选出的关键lncRNA有望为蜜蜂微孢子虫病的诊断和治疗提供新型分子标记和靶点。【前人研究进展】真核细胞中存在大量的ncRNA，其中lncRNA是一类长度＞200 nt，没有长阅读框，由RNA聚合酶Ⅱ转录，同时 5′端加帽与 3′端加尾而形成的转录本，具有复杂而特定的二级结构，通过支架、向导和诱饵等方式参与调节细胞的增殖、分化、凋亡及代谢等生命活动[4-6]。有研究表明lncRNA可通过顺式（cis）作用调控同一染色体上邻近蛋白编码基因的表达，或通过反式（trans）作用调控距离较远基因的转录激活及表达[6-9]。此外，含有微小 RNA 反应元件（microRNA response element，MRE）的lncRNA还能够作为“分子海绵”吸附结合miRNA，通过发挥竞争性内源RNA（competing endogenous RNA，ceRNA）调控作用间接影响miRNA与mRNA的靶向结合[6,10-11]。基于链特异性建库策略的二代测序技术的发展与成熟为在全基因组水平鉴定lncRNA提供了强大工具[6,12]。目前，已在人类[13]、小鼠[14]、拟南芥[15]、果蝇[16]等模式生物中鉴定出大量lncRNA。但相比较而言，真菌的lncRNA研究起步较晚且发展滞后。FAUQUENOY等[9]运用生物信息学方法在裂殖酵母（Schizosaccharomyces pombe）中筛选出与mRNA邻近的68个lncRNA，进一步研究发现lncRNA对邻近蛋白编码基因具有调控作用；LIU等[17]利用RNA-seq技术对未发芽和发芽的家蚕微孢子虫（Nosema bombycis）进行转录组测序，基于测序数据分别组装出 2 756和 2 690个unigenes，进而通过比较分析筛选出66个显著差异表达基因。然而对于东方蜜蜂微孢子虫等蜜蜂病原，lncRNA的相关信息极其匮乏。GUO等[18]通过基于链特异性cDNA文库的lncRNA-seq技术对蜜蜂白垩病病原蜜蜂球囊菌（Ascosphaera apis，简称球囊菌）的菌丝和孢子混合样品进行测序，通过生物信息学方法鉴定出379个lncRNA，进一步分析发现这些lncRNA具有类似于其他物种已报道 lncRNA的结构特征；ATKINSON[19]利用RNA-seq技术在裂殖酵母中鉴定出5 775个lnRNA，进而通过结构优化发现其中的1 557个lnRNA已被注释，多数lncRNA表达与mRNA表达呈正相关，这些lncRNA的cis作用及核内的共转录功能介导了裂殖酵母的调控机制；GUO等[20]利用lncRNA的cDNA芯片技术对人的慢性骨髓白血病细胞的lncRNA进行了全面分析，发现lncRNA-BGL3作为一种 ceRNA靶向结合 miR-17家族的多个miRNA，进而调控肿瘤抑制因子 PTEN相关蛋白编码基因的表达。【本研究切入点】目前，东方蜜蜂微孢子虫的lncRNA研究非常滞后，相关信息匮乏。笔者团队前期已利用基于 cDNA链特异性文库的lncRNA-seq技术对东方蜜蜂微孢子虫纯化孢子进行深度测序，联用CPC、CNCI、CPAT和Pfam软件鉴定出83个lncRNA，进而对上述lncRNA的种类和结构特征进行了初步分析，并利用RT-PCR验证了其中13个lncRNA的真实表达[21]。然而，这些lncRNA在东方蜜蜂微孢子虫孢子中具有怎样的调控功能以及如何发挥作用仍未可知。【拟解决的关键问题】利用small RNA-seq（sRNA-seq）技术对东方蜜蜂微孢子虫纯化孢子进行测序，结合高质量的 sRNA数据和前期获得的 lncRNA数据进一步深入分析东方蜜蜂微孢子虫lncRNA的cis作用和ceRNA作用，为东方蜜蜂微孢子虫的lncRNA信息提供有益补充，并揭示 lncRNA在东方蜜蜂微孢子虫孢子中的潜在功能。

1 材料与方法
试验于2018年12月至2019年4月在福建农林大学完成。

1.1 生物材料
供试东方蜜蜂微孢子虫纯化孢子由福建农林大学动物科学学院（蜂学学院）蜜蜂保护实验室制备和保存[18,21]。

1.2 LncRNA组学数据来源
前期研究中，笔者团队通过Percoll不连续密度梯度离心获得了东方蜜蜂微孢子虫的纯化孢子，利用lncRNA-seq技术对纯化孢子（NcL-1、NcL-2、NcL-3）进行了深度测序，获得了高质量的 lncRNA组学数据[21]。链特异性cDNA文库构建及高通量测序委托北京百迈克生物科技有限公司完成，测序平台为Illumina HiSeq Xten。LncRNA-seq的原始数据已上传 NCBI SRA数据库，BioProject号：PRJNA395264。

1.3 sRNA-seq及数据质控
在东方蜜蜂微孢子虫的纯化孢子中加入适量液氮充分研磨至粉末，利用Trizol法分别抽提3份纯化孢子样品（NcS-1、NcS-2、NcS-3）的总 RNA，通过 Nanodrop超微量测定仪（Thermo公司，美国）检测上述RNA的浓度，通过Qubit 2.0（Life公司，美国）和Agilent 2100 bioanalyzer（Agilent公司，美国）分别检测上述RNA的纯度和完整性。对于检测合格的RNA，以1.5 μL的量作为RNA样本起始量，用无菌水补充体积至6 μL，使用small RNA Sample Pre Kit试剂盒（Illumina公司，美国）进行文库构建。由于 sRNA 的 5′端含磷酸基团，3′端含羟基，利用T4 RNA Ligase 1和T4 RNA Ligase 2分别在sRNA的3′端和5′端连接上接头序列，反转录合成cDNA，再进行PCR扩增，采用胶分离技术筛选目的片段，切胶回收得到的片段即为sRNA文库。文库构建完成后，使用Qubit 2.0对文库的浓度进行检测，将文库浓度稀释至1 ng·μL-1，使用Agilent 2100 bioanalyzer（Agilent公司，美国）对Insert Size进行检测，使用qPCR方法对文库的有效浓度进行准确定量，以保证文库质量。委托北京百迈克生物科技有限公司进行sRNA-seq，测序平台为Illumina MiSeq，sRNA的测序读长为single-end (SE) 50 nt。sRNA-seq的原始数据已上传 NCBI SRA数据库，BioProject号：PRJNA395137。

基于东方蜜蜂微孢子虫纯化孢子的 lncRNA-seq数据，笔者团队前期已联用 CPC、CNCI、CPAT和Pfam软件分别预测出162、295、342和206条lncRNA，四者的交集为83条lncRNA[21]。位于编码蛋白基因上下游的lncRNA可能在转录或转录后水平对其邻近蛋白编码基因的表达具有调控作用[8-9]。因此，将lncRNA上下游100 kb范围内的邻近基因作为其调控的潜在靶基因。利用 Blast软件对 lncRNA的上下游基因进行GO数据库（http://geneontology.org/）和 KEGG数据库（https://www.kegg.jp/）注释。

进一步通过 Bowtie软件将 clean reads分别与Silva数据库（https://www.arb-silva.de/）、GtRNAdb数据库（http://bigd.big.ac.cn/databasecommons/database/id/1306）、Rfam 数据库（http://rfam.xfam.org/）和Repbase数据库（http://www.repeatmasker.org/）进行序列比对，过滤掉其他 ncRNA及重复序列，得到含miRNA的未注释有效读段（unannotated clean reads）映射（mapping）到东方蜜蜂微孢子虫参考基因组（assembly ASM98816v1），获得clean reads在参考基因组上的位置信息。

1.4 LncRNA上下游基因的功能注释
不一会儿，门开了，罗丽进来了，她身后跟着一个高大的小伙子，他们手里提着摩托帽。看见家里坐着警察，罗丽的眼里没有恐惧，只有厌烦。她大模大样地拉着她的男朋友坐下，歪着脑袋等着提问。

1.5 LncRNA的靶miRNA预测与ceRNA调控网络构建及分析
参照郭睿等[22]的方法，利用Target Finder软件[23]预测东方蜜蜂微孢子虫孢子中 lncRNA 靶向结合的miRNA及miRNA靶向结合的mRNA，软件参数设置为：自由能≤-20 kcal·mol-1，空位开度=-0.9，空位扩展=-4.0。利用三者之间的靶向结合关系构建lncRNA-miRNA及lncRNA-miRNA-mRNA调控网络，并通过Cytoscape v3.7.2软件[24]对上述调控网络进行可视化。

1.6 LncRNA和靶mRNA的RT-PCR验证
根据 1.5中靶向预测的结果，随机选取 2个lncRNA（MSTRG.3636.1和 MSTRG.4883.1）和8个靶 mRNA（gene565、gene1719、gene1695、gene1563、gene1439、gene1366、gene1360和 gene1134）进行RT-PCR验证。根据上述lncRNA和mRNA的核酸序列利用DNAMAN软件设计相应的特异性引物，委托上海生工生物工程有限公司合成引物（表 1）。参照GUO等[21]的方法，利用RNA抽提试剂盒（TaKaRa，中国）提取东方蜜蜂微孢子虫孢子的总RNA，作为模板进行反转录，得到的cDNA作为模板进行PCR。PCR反应体系（20 μL）：cDNA 模板 1 μL，上游引物 1 μL，下游引物 1 μL，PCR mix 10 μL，无菌水 7 μL。将无菌水设为阴性对照，反应体系（20 μL）：无菌水模板1 μL，上述2个lncRNA和8个靶mRNA的混合上游引物 1 μL（每种 0.1 μL）、混合下游引物 1 μL（每种0.1 μL），PCR mix 10 μL，无菌水 7 μL。PCR 程序：95℃预变性 5 min；95℃变性 50 s，55℃退火 30 s，72℃延伸50 s，共35个循环；72℃再延伸10 min。PCR产物经 1.6%琼脂糖凝胶电泳和凝胶成像仪（上海培清，中国）检测。

1.7 MiRNA的Stem-loop RT-PCR验证
通过Stem-loop RT-PCR验证对ceRNA调控网络中 3个靶 miRNA（nce-miR-7502、nce-miR-7729和nce-miR-8565）的表达进行验证。利用RNA抽提试剂盒（TaKaRa公司，日本）提取东方蜜蜂微孢子虫孢子的总RNA。参照郭睿等[25]和杜宇等[26]的方法，利用DNAMAN软件（LynnonBiosoft公司，美国）设计miRNA的Stem-loop引物、上游和下游引物，委托上海生工生物工程有限公司合成引物（表 1）。利用Stem-loop引物经反转录得到 cDNA，作为模板进行PCR扩增。PCR体系（20 μL）：上下游引物（1.67 μmol·L-1）和 cDNA 模板各 1 μL，PCR mix 10 μL，无菌水7 μL。将DEPC水设为阴性对照，反应体系（20 μL）：DEPC水模板1 μL，上述3个靶miRNA的混合上游引物0.9 μL（每种0.3 μL）、混合下游引物0.9 μL（每种 0.3 μL），PCR mix 10 μL，无菌水 7.2 μL。PCR 程序：95℃ 5 min；95℃ 30 s，46℃ 30 s，72℃ 1 min，35个循环；72℃ 10 min。PCR产物经1.6%琼脂糖凝胶电泳检测。

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

pagenumber_ebook=202,pagenumber_book=2126
2 结果
2.1 LncRNA-seq及sRNA-seq的数据质控
前期研究中，东方蜜蜂微孢子虫纯化孢子样品（NcL-1、NcL-2、NcL-3）的lncRNA-seq得到的clean reads数分别为107 123 113、117 688 499和116 776 007条，占 raw reads的比例分别为 79.75%、78.84%和76.15%；Q30分别达到90.39%、90.49%和89.79%[21]。高质量的 lncRNA组学数据可用于本研究中 lncRNA的生物信息学分析。

东方蜜蜂微孢子虫纯化孢子样品的 sRNA-seq分别得到16 597 883、15 451 791和12 248 316条raw reads，经过严格过滤得到的clean reads数分别为15 608 370（94.04%）、14 249 255（92.22%）和11 440 684（93.41%），Q30分别为98.61%、98.65%和98.58%（表2），说明本研究得到的sRNA数据质量良好，可进行下一步分析。将各孢子样品的unannotated reads比对到东方蜜蜂微孢子虫的参考基因组，比对上的reads总数为22 000 837条，比对率≥52.87%（表3）。

2.2 东方蜜蜂微孢子虫孢子中lncRNA上下游基因的GO数据库注释
东方蜜蜂微孢子虫孢子中56个lncRNA共预测出 310个上下游基因。对这些上下游基因进行 GO数据库注释，结果显示它们可注释到生物学进程、分子功能和细胞组分相关的35个功能条目；在生物学进程中，注释基因数最多的条目分别为代谢进程（105）、细胞进程（100）、单组织进程（45）、生物学调控（19）和生物学进程调控（18）；在分子功能中，注释基因数最多的分别为催化活性（96）、结合（95）、结构分子活性（9）、核酸结合转录因子活性（3）和转运器活性（3）；在细胞组分中，注释基因数最多的分别为细胞（64）、细胞组件（64）、细胞器（46）、大分子复合物（34）和细胞器组件（13）（图1）。括号内的数字代表注释到该功能条目的上下游基因数。

pagenumber_ebook=203,pagenumber_book=2127
图1 东方蜜蜂微孢子虫孢子中lncRNA的上下游基因的GO数据库注释
Fig. 1 GO database annotation of upstream and downstream genes of lncRNAs in N. ceranae spore

表2 sRNA-seq测序数据信息统计
Table 2 Overview of sRNA-seq datasets

pagenumber_ebook=204,pagenumber_book=2128
表3 sRNA-seq数据比对东方蜜蜂微孢子虫参考基因组的信息统计
Table 3 Summary of mapping information of sRNA-seq data to the reference genome of N. ceranae

pagenumber_ebook=204,pagenumber_book=2128
2.3 东方蜜蜂微孢子虫孢子中lncRNA上下游基因的KEGG数据库注释
进一步对lncRNA的上下游基因进行KEGG数据库注释，结果显示310个上下游基因可注释到56条通路，注释基因数最多的前10位分别是新陈代谢途径（25）、抗生素的生物合成（12）、次生代谢产物的生物合成（11）、核糖体（10）、嘌呤代谢（7）、蛋白酶体（7）、RNA转运（7）、真核生物中核糖体的生物发生（7）、碳代谢（6）和嘧啶代谢（6）（图2）。括号内的数字代表注释到该通路的上下游基因数。

pagenumber_ebook=204,pagenumber_book=2128
图2 东方蜜蜂微孢子虫孢子中lncRNA的上下游基因的KEGG数据库注释
Fig. 2 KEGG database annotation of upstream and downstream genes of lncRNAs in N. ceranae spore

2.4 东方蜜蜂微孢子虫孢子中lncRNA-miRNA调控网络的构建及分析
对东方蜜蜂微孢子虫孢子中 lncRNA进行靶miRNA预测，发现3个lncRNA与4个miRNA存在靶向结合关系，其中 MSTRG.3636.1靶向 nce-miR-8565，MSTRG.4498.1靶向nce-miR-7502和nce-miR-8639，MSTRG.4883.1靶向nce-miR-7729（图3）。

2.5 东方蜜蜂微孢子虫孢子中 lncRNA-miRNA-mRNA调控网络的构建及分析
利用软件预测 lncRNA靶向结合 miRNA的靶mRNA，综合lncRNA、miRNA与mRNA的靶向结合关系构建调控网络，分析发现MSTRG.4883.1靶向的nce-miR-7729共结合16个mRNA，MSTRG.3636.1靶向的nce-miR-8565共结合15个mRNA，lncRNA、miRNA和mRNA之间形成密切的调控关系；3个lncRNA的4个靶miRNA可结合59个mRNA，其中MSTRG.4498.1靶向的nce-miR-7502和nce-miR-8639共结合28个mRNA，三者形成更为复杂的调控关系（图4）。

pagenumber_ebook=205,pagenumber_book=2129
图3 东方蜜蜂微孢子虫孢子中的lncRNA-miRNA调控网络
Fig. 3 LncRNA-miRNA regulatory network in N. ceranae spore

进一步分析发现，上述靶mRNA在Nr数据库中注释皆为假定蛋白；但在Swissprot数据库中均有功能注释信息，如gene2014、gene1662和gene1383分别注释为转录相关蛋白1、微管蛋白γ链和V型质子ATP酶催化亚基A（表4）。

2.6 东方蜜蜂微孢子虫孢子中lncRNA、靶miRNA和靶mRNA的分子验证
通过RT-PCR对上述ceRNA网络中lncRNA和靶mRNA的进行扩增，扩增产物的电泳结果显示8个靶mRNA均能扩增出预期大小的目的条带，阴性对照未扩增出片段；2个 lncRNA也能扩增出预期片段，阴性对照未扩增出片段；利用Stem-loop RT-PCR对上述ceRNA网络中靶miRNA进行扩增，扩增产物的电泳结果显示3个靶miRNA均扩增出大小约100 bp的目的片段，阴性对照未扩增出片段（图 5）。上述结果表明东方蜜蜂微孢子虫孢子中lncRNA、靶miRNA和靶mRNA真实表达。

pagenumber_ebook=205,pagenumber_book=2129
图4 东方蜜蜂微孢子虫孢子中3个lncRNA的竞争性内源RNA调控网络
Fig. 4 Competing endogenous network of three lncRNAs in N. ceranae spore

表4 东方蜜蜂微孢子虫孢子中lncRNA靶向miRNA的靶mRNA的Swissprot和Nr数据库注释
Table 4 Swissprot and Nr database annotations of target mRNAs of miRNAs targeted by lncRNAs in N. ceranae spore

pagenumber_ebook=206,pagenumber_book=2130
pagenumber_ebook=207,pagenumber_book=2131
图5 东方蜜蜂微孢子虫孢子中lncRNA、靶miRNA和靶mRNA的RT-PCR验证
Fig. 5 RT-PCR verification of lncRNAs, target miRNAs and target mRNAs in N. ceranae spore

M：DNA marker；1：假定蛋白NCER_101243编码基因Hypothetical protein NCER_101243 coding gene；2：假定蛋白NCER_102113编码基因Hypothetical protein NCER_102113 coding gene；3：假定蛋白NCER_100958编码基因Hypothetical protein NCER_100958 coding gene；4：假定蛋白NCER_101315编码基因Hypothetical protein NCER_101315 coding gene；5：假定蛋白NCER_101322编码基因Hypothetical protein NCER_101322 coding gene；6：假定蛋白NCER_101544编码基因Hypothetical protein NCER_101544 coding gene；7：假定蛋白NCER_100991编码基因Hypothetical protein NCER_100991 coding gene；8：假定蛋白NCER_1011114编码基因Hypothetical protein NCER_101114 coding gene；9：MSTRG.3636.1；10：MSTRG.4883.1；11：nce-miR-8565；12：nce-miR-7502；13：nce-miR-7729；N1、N2：阴性对照（无菌水）Negative control (sterile water)；N3：阴性对照（DEPC水）Negative control (DEPC water)

3 讨论
东方蜜蜂微孢子虫孢子处于休眠态，人们对于孢子中是否存在转录、翻译和新陈代谢等生命活动的认识十分有限。前期研究中，笔者团队证实东方蜜蜂微孢子虫孢子同样存在一定水平的转录，并基于高质量的lncRNA组学数据鉴定出83条lncRNA，同时对这些lncRNA的种类和结构特征进行了全面分析[21]，为东方蜜蜂微孢子虫 lncRNA的深入研究提供了必要的数据基础。本研究结合前期获得的东方蜜蜂微孢子虫lncRNA组学数据和测序得到的sRNA组学数据，进一步深入探究东方蜜蜂微孢子虫孢子中 lncRNA的ceRNA调控网络及潜在功能。目前，动植物和微生物等物种的sRNA-seq数据在相应参考基因组上的比对率一般介于50%—70%[27-29]，原因在于 sRNA的长度很短，比对参考基因组时要求每个碱基均要匹配上，因此导致比对率偏低。本研究中，NcS-1、NcS-2和NcS-3的clean reads在东方蜜蜂微孢子虫参考基因组上的比对率介于 52.87%—60.36%，与其他物种的研究报道类似，属于正常的比对率范围。目前，NONCODE[30]、LncRNAdb[31]和 LncRNome[32]等lncRNA数据库中有限的功能注释信息局限于人类和其他少数模式物种（如小鼠、果蝇和拟南芥等）。GUO等[33]构建了家蚕微孢子虫感染家蚕（Bombyx mori）卵巢BmN细胞的连续感染系，进而利用脂质体介导将含 gfp的转座子载体 piggyBac和非转座子载体pIZT/V5-His转染感染细胞，通过分子生物学手段证实外源gfp被成功导入家蚕微孢子虫基因组；此外，GUO等[33]还通过饲喂家蚕微孢子虫孢子感染家蚕个体，然后通过无菌操作将蚕血加入正常 BmN细胞，再进行脂质体转染的方法也成功将gfp导入了家蚕微孢子虫基因组；这是对微孢子虫直接进行转基因操作的首例报道，为东方蜜蜂微孢子虫等其他微孢子虫的转基因研究提供了重要的思路和方法借鉴。目前，东方蜜蜂微孢子虫的基因功能研究仍处于初步阶段。PALDI等[34]利用饲喂的方法将人工合成的靶向 ADP/ATP转运器蛋白编码基因的双链RNA（dsRNA）导入被东方蜜蜂微孢子虫感染的西方蜜蜂体内，研究发现病原的ADP/ATP转运器蛋白转录本受到干扰而沉默，该沉默影响了病原增殖水平和宿主生理状态；HUANG等[35]通过饲喂小干扰RNA（siRNA）对东方蜜蜂微孢子虫感染的西方蜜蜂的Dicer进行敲减，发现病原孢子数显著下降，此外超过10%的病原蛋白编码基因呈现显著性差异表达。对于包括东方蜜蜂微孢子虫在内的绝大多数物种，通过分析和验证 lncRNA的 cis、trans和ceRNA作用，从而间接推测lncRNA的潜在功能仍是常用方法[6-11]。较多的研究表明lncRNA可以与一些较远基因上的增强子或启动子结合，通过 trans作用行使其对编码蛋白基因表达的调控[7]。然而，本研究并没有预测到表达量存在正相关或负相关关系的lncRNA和mRNA，暗示lncRNA在东方蜜蜂微孢子虫的休眠态孢子中不行使trans作用。鉴于lncRNA的表达具有时空特异性[36]，东方蜜蜂微孢子虫在孢子发芽和侵染宿主细胞的过程中是否存在 trans作用或通过哪些方式发挥主要作用，值得进一步深入研究。利用链特异性建库的 lncRNA-seq技术对东方蜜蜂微孢子虫感染的蜜蜂中肠组织进行测序，通过连续比对宿主蜜蜂参考基因组和东方蜜蜂微孢子虫参考基因组过滤得到病原的 lncRNA组学数据，进而分析病原lncRNA的trans作用是具有可行性的方法。为探究东方蜜蜂微孢子虫孢子中lncRNA的潜在功能，本研究利用生物信息学方法对lncRNA的上下游基因和靶向结合的miRNA进行预测，进而对lncRNA的cis作用方式进行分析，发现东方蜜蜂微孢子虫孢子中 56条lncRNA可能调控310个上下游基因，它们涉及代谢进程、细胞进程、催化活性、结合和细胞组件等 35条功能条目，说明lncRNA可能广泛参与孢子休眠状态的基础生命活动，但其具体功能仍需要进一步实验验证。

目前，lncRNA的相关研究主要集中在少数模式生物，如人类[13]、小鼠[14]、果蝇[16]、酵母[7,9]等；其中多数研究集中在lncRNA的生物信息学预测和分析，仅少数研究涉及lncRNA的功能验证，有限的技术手段是最大的限制因素。ZHAN等[37]研究发现山羊骨骼肌细胞中的1 153个lncRNA潜在调控1 455个上游及下游基因，部分上下游基因可能参与了骨骼肌细胞的发育，说明相应的lncRNA在山羊骨骼肌细胞发育过程中具有潜在的cis调控作用；FAUQUENOY等[9]通过对野生型和缺失一个拷贝rse1等位基因的杂合二倍体突变型裂殖酵母菌株进行研究，发现lncRNA通过抑制复合物的支架结构特异性抑制邻近的rse1等位基因表达，表明lncRNA通过cis作用抑制裂殖酵母的细胞分化。本研究中，对于东方蜜蜂微孢子虫的lncRNA，分别有105、100和45个上下游基因注释在代谢进程、细胞进程和单组织进程，说明相应的lncRNA通过cis作用维持东方蜜蜂微孢子虫孢子处于休眠态期间的基础新陈代谢。此外，分别有64、64和46个上下游基因涉及细胞组件、细胞和细胞器，表明lncRNA对东方蜜蜂微孢子虫孢子的细胞生命活动中具有潜在的调节功能。还发现有96个上下游基因涉及催化活性，进一步表明lncRNA可能参与东方蜜蜂微孢子虫孢子内新陈代谢过程的生化反应。LncRNA通过与 DNA、RNA或蛋白质的相互作用而发挥调控作用[4,36]。东方蜜蜂微孢子虫在进行发芽和极丝弹射之前，首先需要识别和接触宿主细胞，进而与宿主细胞相互作用[38-41]，上述过程必然也伴随着不同生物大分子之间的结合。本研究中，有 95个上下游基因涉及分子功能中的结合，暗示lncRNA在东方蜜蜂微孢子虫孢子内可能通过与其他生物大分子的相互结合发挥特定的调控功能，也可能在病原识别和接触宿主细胞中扮演特定角色。微孢子虫具有一层坚硬的孢子壁，能够保护内部结构免受外界不良环境破坏[41]。东方蜜蜂微孢子虫在宿主细胞内增殖的过程需要不断适应外部环境并作出应激反应。本研究中，有14个上下游基因涉及应激反应，说明相应的 lncRNA（MSTRG.4856.1、MSTRG.5126.1、MSTRG.5617.1、MSTRG.6243.1、MSTRG.6290.1、MSTRG.6328.1、MSTRG.6473.1、MSTRG.6633.1、MSTRG.6654.2、MSTRG.6688.1、MSTRG.6683.1、MSTRG.6692.1和 MSTRG.6703.1）在孢子的环境应激方面发挥一定作用。

在糖酵解途径中，2-磷酸甘油酸可在烯醇化酶作用下转变为磷酸烯醇式丙酮酸，进而在丙酮酸激酶作用下转变为烯醇式丙酮酸，过程中产生两分子的H2O，同时两分子 ADP转化为两分子 ATP[42]。DOLGIKH等[43]检测了黄斑黑蟋蟀微孢子虫（Nosema grylli）孢子中6-磷酸葡萄糖脱氢酶、葡萄糖磷酸变位酶、磷酸葡萄糖异构酶、6-磷酸果糖激酶及3-磷酸甘油酸激酶等14种酶的活性，并分析了相关能量代谢途径，发现黄斑黑蟋蟀微孢子虫可能利用自身储备的碳水化合物而不是外源物质，推测糖酵解是黄斑黑蟋蟀微孢子虫碳水化合物分解代谢的主要途径。本研究发现，gene1402、gene2290、gene2665和gene926这4个上下游基因涉及糖酵解/糖异生途径，其中 gene926和gene2290分别注释到 Swissprot数据库中的胃瘤真菌的烯醇化酶基因和兔脑炎微孢子虫（Encephalitozoon cuniculi）的丙酮酸激酶基因，表明相应的 lncRNA（MSTRG.6328.1和MSTRG.4700.1）具有调节糖酵解/糖异生途径的潜在功能，可能在东方蜜蜂微孢子虫孢子的能量代谢方面发挥作用。此外，共有85个上下游基因注释在10条碳水化合物代谢相关通路，包括代谢途径（25）、抗生素的生物合成（12）和次生代谢产物的生物合成（11）等，进一步说明相应的 lncRNA参与调节东方蜜蜂微孢子虫孢子的能量代谢。核苷酸是组成DNA和RNA的重要部分，并作为电子信号和能量载体在生物体内发挥重要作用[44]。HEINZ等[45]研究了来自 HIV/AIDS患者的人气管普孢子虫（Trachipleistophora hominis）中的4种核苷酸转运蛋白，验证了这些蛋白可以转运有放射性标记的2种嘌呤核苷酸，并未发现嘧啶核苷酸；作者进一步将人工合成酶与细胞内人气管普孢子虫、兔脑炎微孢子虫和东方蜜蜂微孢子虫酶比较，发现存在嘌呤和嘧啶核苷酸在磷酸化和氧化作用之间相互转换的一类核心酶，表明微孢子虫中的此类酶能够保持嘌呤和嘧啶处于不同的激活状态满足生命活动过程中的相关代谢途径。本研究中，分别有7和6个上下游基因注释到嘌呤代谢和嘧啶代谢，其中的共有基因为 6个，分别为gene2676、 gene2395、gene1404、gene2450、gene1439和gene1058，说明相应的lncRNA（MSTRG.6703.1、MSTRG.6459.2、MSTRG.6459.1、MSTRG.5334.1、MSTRG.6512.1、MSTRG.5375.1、MSTRG.4883.1 和MSTRG.4883.1）具有通过调控东方蜜蜂微孢子虫孢子中的嘌呤代谢和嘧啶代谢影响遗传信息传递、电子传递、信号通路及能量代谢的潜在作用。还发现分别有6、2、1、1和1个上下游基因分别注释到碳代谢、丙酮酸代谢、磷酸肌醇代谢、脂肪酸代谢和果糖甘露糖代谢等物质代谢通路；有 2、1和1个基因分别注释到甲烷代谢、氧化磷酸化和过氧化物酶体等能量代谢通路；上述结果进一步表明相应的 lncRNA（MSTRG.4700.1、MSTRG.5292.1、MSTRG.5334.1、MSTRG.6328.1、MSTRG.6383.1、MSTRG.6688.1、MSTRG.6692.1和MSTRG.6703.1）通过cis作用参与东方蜜蜂微孢子虫孢子中的基本物质和能量代谢的调控。

SALMENA等[46]于 2011年提出了“ceRNA 假说”，该假说认为含有 MRE的 mRNA、假基因和lncRNA等RNA，能够竞争性结合miRNA，从而间接影响靶基因的表达量[46-47]。随后，越来越多的研究结果证实了这一主流假说[22,46-48]。WANG 等[48]利用荧光素酶RNA免疫沉淀技术探究在低氧或炎症因子处理下 lncRNA的ceRNA机制，结果显示 lncRNA H19和HULC可作为ceRNA靶向结合let-7a/let-7b和miR-372/miR-373，从而参与炎症反应促进胆管癌细胞的迁移和侵袭。前期研究中，笔者团队曾对意大利蜜蜂（Apis mellifera ligustica）工蜂7 d和10 d中肠发育过程中lncRNA和circRNA差异表达谱及调控网络进行研究，分析发现此2种ncRNA均可作为潜在的ceRNA吸附miRNA，从而减少对靶mRNA的抑制或降解[22,49]。本研究发现，MSTRG.4498.1、MSTRG.4883.1和MSTRG.3636.1与4个miRNA（nce-miR-8565、nce-miR-7502、nce-miR-8639和nce-miR-7729）之间存在靶向结合关系；此外，上述4个miRNA与59个mRNA也存在靶向结合关系，表明lncRNA可能通过竞争性结合miRNA调控mRNA的表达量，进而影响东方蜜蜂微孢子虫孢子休眠态下的诸多生命活动。对于 lncRNA-miRNA-mRNA调控网络，通过 RT-PCR证实了2个 lncRNA和8个靶mRNA的真实表达，通过Stem-loop RT-PCR证实了3个靶miRNA的真实表达（图5）。但本研究中的靶向结合关系是基于生物信息学软件预测出来的，上述靶向结合关系是否真实存在仍需通过双荧光素酶报告系统加以验证。对调控网络中的靶mRNA进行功能和代谢通路注释，结果显示它们涉及19条功能条目，包括代谢进程、细胞进程和催化活性等，进一步分析发现对于调控网络中间接调控靶mRNA的lncRNA，其中部分lncRNA也参与了对上下游基因的调控；此外，有24个靶mRNA富集在18条代谢通路，涉及嘌呤代谢、嘧啶代谢、脂肪酸代谢、过氧化物酶体及氧化磷酸化等。上述结果表明部分 lncRNA能同时通过 cis和ceRNA作用影响靶基因的表达水平，从而对东方蜜蜂微孢子虫孢子的新陈代谢和基本生命活动进行灵活调控。

4 结论
部分lncRNA可能通过cis作用调控上下游基因的表达，进而调节东方蜜蜂微孢子虫孢子中的物质代谢、能量代谢和环境应激等生物学过程；部分lncRNA可能作为 ceRNA参与调节东方蜜蜂微孢子虫孢子休眠态下的诸多生命活动；此外，少数lncRNA能同时通过cis和ceRNA作用影响靶基因的表达水平，从而对东方蜜蜂微孢子虫孢子的新陈代谢和基本生命活动进行灵活调控。

References

[1] KEELING P J, FAST N M. Microsporidia: Biology and evolution of highly reduced intracellular parasites. Annual Review of Microbiology,2002, 56: 93-116.

[2] HOLT H L, VILLAR G, CHENG W, SONG J, GROZINGER C M.Molecular, physiological and behavioral responses of honey bee (Apis mellifera) drones to infection with microsporidian parasites. Journal of Invertebrate Pathology, 2018, 155: 14-24.

[3] KURZE C, CONTE Y L, KRYGER P, LEWKOWSKI O, MÜLLER T,MORITZ R F A. Infection dynamics of Nosema ceranae in honey bee midgut and host cell apoptosis. Journal of Invertebrate Pathology,2018, 154: 1-4.

[4] ARNDT K M, REINES D. Termination of transcription of short noncoding RNAs by RNA polymerase II. Annual Review of Biochemistry, 2015, 84: 381-404.

[5] WU Z M, LIU X X, LIU L, DENG H L, ZHANG J J, XU Q, CEN B H, JI A M. Regulation of lncRNA expression. Cellular and Molecular Biology Letters, 2014, 19(4): 561-575.

[6] KOPP F, MENDELL J T. Functional classification and experimental dissection of long noncoding RNAs. Cell, 2018, 172(3): 393-407.

[7] DE BOER P, VOS H R, FABER A W, VOS J C, RAUÉ H A. Rrp5p, a trans-acting factor in yeast ribosome biogenesis, is an RNA-binding protein with a pronounced preference for U-rich sequences. RNA,2006, 12(2): 263-271.

[8] YAN P, LUO S, LU J Y, SHEN X. Cis- and trans-acting lncRNAs in pluripotency and reprogramming. Current Opinion in Genetics and Development, 2017, 46: 170-178.

[9] FAUQUENOY S, MIGEOT V, FINET O, YAGUE-SANZ C,KHOROSJUTINA O, EKWALL K, HERMAND D. Repression of cell differentiation by a cis-acting lincRNA in fission yeast. Current Biology, 2018, 28(3): 383-391.

[10] SARVER A L, SUBRAMANIAN S. Competing endogenous RNA database. Bioinformation, 2012, 8(15): 731-733.

[11] CESANA M, CACCHIARELLI D, LEGNINI I, SANTINI T,STHANDIER O, CHINAPPI M, TRAMONTANO A, BOZZONI I. A long noncoding RNA controls muscle differentiation by functioning as a competing endogenous RNA. Cell, 2011, 147(2): 358-369.

[12] CHAN W L, HUANG H D, CHANG J G. LncRNAMap: A map of putative regulatory functions in the long non-coding transcriptome.Computational Biology and Chemistry, 2014, 50: 41-49.

[13] KRISHNAN J, MISHRA R K. Emerging trends of long non-coding RNAs in gene activation. The FEBS Journal, 2014, 281(1): 34-45.

[14] GUTTMAN M, AMIT I, GARBER M, FRENCH C, LIN M F,FELDSER D, HUARTE M, ZUK O, CAREY B W, CASSADY J P,CABILI M N, JAENISCH R, MIKKELSEN T S, JACKS T,HACOHEN N, BERNSTEIN B E, KELLIS M, REGEV A, RINN J L,LANDER E S. Chromatin signature reveals over a thousand highly conserved large non-coding RNAs in mammals. Nature, 2009,458(7235): 223-227.

[15] LIU J, JUNG C, XU J, WANG H, DENG S, BERNAD L,ARENAS-HUERTERO C, CHUA N H. Genome-wide analysis uncovers regulation of long intergenic noncoding RNAs in Arabidopsis. The Plant Cell, 2012, 24(11): 4333-4345.

[16] CHEN B, ZHANG Y, ZHANG X, JIA S L, CHEN S, KANG L.Genome-wide identification and developmental expression profiling of long noncoding RNAs during Drosophila metamorphosis. Scientific Reports, 2016, 6: 23330.

[17] LIU H, LI M Q, HE X Y, CAI S F, HE X K, LU X M. Transcriptome sequencing and characterization of ungerminated and germinated spores of Nosema bombycis. Acta Biochimica et Biophysica Sinica,2016, 48(3): 246-256.

[18] GUO R, CHEN D F, XIONG C L, HOU C S, ZHENG Y Z, FU Z M,DIAO Q Y, ZHANG L, WANG H Q, HOU Z X, LI W D, 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.

[19] ATKINSON S R. The fission yeast non-coding transcriptome[D].London: University College London, 2014.

[20] GUO G J, KANG Q Z, CHEN Q H, JING O Y, CHEN J L.Abstract 3543: A long noncoding RNA plays a critical role in Bcr-Abl-mediated cellular transformation. Cancer Research, 2014,74(19): 3543.

[21] GUO R, CHEN D F, XIONG C L, HOU C S, ZHENG Y Z, FU Z M,LIANG Q, DIAO Q Y, ZHANG L, WANG H Q, HOU Z , KUMAR D.First identification of long non-coding RNAs in fungal parasite Nosema ceranae. Apidologie, 2018, 49(5): 660-670.

[22] 郭睿, 耿四海, 熊翠玲, 郑燕珍, 付中民, 王海朋, 杜宇, 童新宇,赵红霞, 陈大福. 意大利蜜蜂工蜂中肠发育过程中长链非编码RNA的差异表达分析. 中国农业科学, 2018, 51(18): 3600-3613.GUO R, GENG S H, XIONG C L, ZHENG Y Z, FU Z M, WANG H P,DU Y, TONG X Y, ZHAO H X, CHEN D F. 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)

[23] LI J, MA W, ZENG P, WANG J, GENG B, YANG J, CUI Q. LncTar:A tool for predicting the RNA targets of long noncoding RNAs.Briefings in Bioinformatics, 2015, 16(5): 806-812.

[24] SMOOT M E, ONO K, RUSCHEINSKI J, WANG P L, IDEKER T.Cytoscape 2.8: New features for data integration and network visualization. Bioinformatics, 2011, 27(3): 431-432.

[25] 郭睿, 杜宇, 熊翠玲, 郑燕珍, 付中民, 徐国钧, 王海朋, 陈华枝,耿四海, 周丁丁, 石彩云, 赵红霞, 陈大福. 意大利蜜蜂幼虫肠道发育过程中的差异表达 microRNA及其调控网络. 中国农业科学,2018, 51(21): 4197-4209.GUO R, DU Y, XIONG C L, ZHENG Y Z, FU Z M, XU G J, WANG H P, CHEN H Z, GENG S H, ZHOU D D, SHI C Y, ZHAO H X,CHEN D F. Differentially expressed microRNA and their regulation networks during the developmental process of Apis mellifera ligustica larval gut. Scientia Agricultura Sinica, 2018, 51(21): 4197-4209. (in Chinese)

[26] 杜宇, 童新宇, 周丁丁, 陈大福, 熊翠玲, 郑燕珍, 徐国钧, 王海朋, 陈华枝, 郭意龙, 隆琦, 郭睿. 中华蜜蜂幼虫肠道响应球囊菌胁迫的 microRNA应答分析. 微生物学报, 2019, 59(9):1747-1764.DU Y, TONG X Y, ZHOU D D, CHEN D F, XIONG C L, ZHENG Y Z, XU G J, WANG H P, CHEN H Z, GUO Y L, 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)

[27] BI F C, MENG X C, MA C, YI G J. Identification of miRNAs involved in fruit ripening in Cavendish bananas by deep sequencing.BMC Genomics, 2015, 16: 776.

[28] CASSELLI T, QURESHI H, PETERSON E, PERLEY D, BLAKE E,JOKINEN B, ABBAS A, NECHAEV S, WATT J A, DHASARATHY A, BRISSETTE C A. MicroRNA and mRNA transcriptome profiling in primary human astrocytes infected with Borrelia burgdorferi. PLoS ONE, 2017, 12(1): e0170961.

[29] JIANG X F, QIAO F, LONG Y L ,CONG H Q, SUN H P. MicroRNA-like RNAs in plant pathogenic fungus Fusarium oxysporum f. sp.niveum are involved in toxin gene expression fine tuning. 3 Biotech,2017, 7(5): Article 354.

[30] FANG S S, ZHANG L L, GUO J C, NIU Y W, WU Y, LI H, ZHAO L H, LI X Y, TENG X Y, SUN X H, SUN L, ZHANG M Q, CHEN R S,ZHAO Y. NONCODEV5: A comprehensive annotation database for long non-coding RNAs. Nucleic Acids Research, 2018, 46(Database issue): D308-D314.

[31] QUEK X C, THOMSON D W, MAAG J L V, BARTONICEK N,SIGNAL B, CLARK M B, GLOSS B S, DINGER M E. LncRNAdb v2.0: Expanding the reference database for functional long noncoding RNAs. Nucleic Acids Research, 2015, 43(Database issue):D168-D173.

[32] HARROW J, FRANKISH A, GONZALEZ J M, TAPANARI E,DIEKHANS M, KOKOCINSKI F, AKEN B L, BARRELL D,ZADISSA A, SEARLE S, BARNES I, BIGNELL A, BOYCHENKO V, HUNT T, KAY M, MUKHERJEE G, RAJAN J, DESPACIOREYES G, SAUNDERS G, STEWARD C, HARTE R, LIN M,HOWALD C, TANZER A, DERRIEN T, CHRAST J, WALTERS N,BALASUBRAMANIAN S, PEI B, TRESS M, RODRIGUEZ J M,EZKURDIA I, VAN BAREN J, BRENT M, HAUSSLER D, KELLIS M, VALENCIA A, REYMOND A, GERSTEIN M, GUIGÓ R,HUBBARD T J. GENCODE: The reference human genome annotation for the ENCODE Project. Genome Research, 2012, 22(9):1760-1774.

[33] GUO R, CAO G L, LU Y H, XUE R Y, KUMAR D, HU X L, GONG C L. Exogenous gene can be integrated into Nosema bombycis genome by mediating with a non-transposon vector. Parasitology Research, 2016, 115(8): 3093-3098.

[34] PALDI N, GLICK E, OLIVA M, ZILBERBERG Y, AUBIN L,PETTIS J, CHEN Y P, EVANS J D. Effective gene silencing in a microsporidian parasite associated with honeybee (Apis mellifera)colony declines. Applied and Environmental Microbiology, 2010,76(17): 5960-5964.

[35] HUANG Q, LI W, CHEN Y P, RETSCHNIG-TANNER G, YANEZ O,NEUMANN P, EVANS J D. Dicer regulates Nosema ceranae proliferation in honeybees. Insect Molecular Biology, 2019, 28(1):74-85.

[36] WANG K C, CHANG H Y. Molecular mechanisms of long noncoding RNAs. Molecular Cell, 2011, 43(6): 904-914.

[37] ZHAN S Y, DONG Y, ZHAO W, GUO J Z, ZHONG T, WANG L J, LI L, ZHANG H P. Genome-wide identification and characterization of long non-coding RNAs in developmental skeletal muscle of fetal goat.BMC Genomics, 2016, 17: 666.

[38] FRANZEN C, MULLER A, HARTMANN P, SALZBERGER B. Cell invasion and intracellular fate of Encephalitozoon cuniculi(Microsporidia). Parasitology, 2005, 130(3): 285-292.

[39] CHEN D F, CHEN H Z, DU Y, ZHOU D D, GENG S H, WANG H P,WAN J Q, XIONG C L, ZHENG Y Z, 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.

[40] YANG D L, PAN L X, CHEN Z Z, DU H H, LUO B, LUO J, PAN G Q. The roles of microsporidia spore wall proteins in the spore wall formation and polar tube anchorage to spore wall during development and infection processes. Experimental Parasitology, 2018, 187:93-100.

[41] GISDER S, MÖCKEL N, LINDE A, GENERSCH E. A cell culture model for Nosema ceranae and Nosema apis allows new insights into the life cycle of these important honey bee-pathogenic microsporidia.Environmental Microbiology, 2011, 13(2): 404-413.

[42] NAKANISHI N, OZAWA K, YAMADA S. Enzymes of the glycolytic pathway-phosphofructokinase, pyruvate kinase and lactate dehydrogenase//Dynamic Aspects of Dental Pulp, 1990: 203-220.

[43] DOLGIKH V V, SOKOLOVA J J, ISSI I V. Activities of enzymes of carbohydrate and energy metabolism of the spores of the microsporidian,Nosema grylli. The Journal of Eukaryotic Microbiology, 1997, 44(3):246-249.

[44] ZIMMERMANN H. Extracellular ATP and other nucleotides—ubiquitous triggers of intercellular messenger release. Purinergic Signalling, 2015, 12(1): 25-57.

[45] HEINZ E, HACKER C, DEAN P, MIFSUD J, GOLDBERG A V,WILLIAMS T A, NAKJANG S, GREGORY A, HIRT R P, LUCOCQ J M, KUNJI E R S, EMBLEY T M. Plasma membrane-located purine nucleotide transport proteins are key components for host exploitation by microsporidian intracellular parasites. PLoS Pathogens, 2014,10(12): e1004547.

[46] SALMENA L, POLISENO L, TAY Y, KATS L, PANDOLFI P P. A ceRNA hypothesis: The Rosetta Stone of a hidden RNA language?Cell, 2011, 146(3): 353-358.

[47] POLISENO L, SALMENA L, ZHANG J, CARVER B, HAVEMAN W J, PANDOLFI P P. A coding-independent function of gene and pseudogene mRNAs regulates tumour biology. Nature, 2010,465(7301): 1033-1038.

[48] WANG W T, YE H, WEI P P, HAN B W, HE B, CHEN Z H, CHEN Y Q. LncRNAs H19 and HULC, activated by oxidative stress, promote cell migration and invasion in cholangiocarcinoma through a ceRNA manner. Journal of Hematology and Oncology, 2016, 9: 117.

[49] 郭睿, 陈华枝, 熊翠玲, 郑燕珍, 付中民, 徐国钧, 杜宇, 王海朋,耿四海, 周丁丁, 刘思亚, 陈大福. 意大利蜜蜂工蜂中肠发育过程中的差异表达环状RNA及其调控网络分析. 中国农业科学, 2018,51(23): 4575-4590.GUO R, CHEN H Z, XIONG C L, ZHENG Y Z, FU Z M, XU G J,DU Y, WANG H P, GENG S H, ZHOU D D, LIU S Y, CHEN D F.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)

Investigation of Competing Endogenous RNA Regulatory Network and Putative Function of Long Non-Coding RNAs in Nosema ceranae Spore

ZHOU DingDing, SHI XiaoYu, WANG Jie, FAN YuanChan, ZHU ZhiWei, JIANG HaiBin, FAN XiaoXue,XIONG CuiLing, ZHENG YanZhen, FU ZhongMin, XU GuoJun, CHEN DaFu, GUO Rui
(College of Animal Sciences (College of Bee Science), Fujian Agriculture and Forestry University, Fuzhou 350002)

Abstract: 【Background】Long non-coding RNAs (lncRNAs) are a kind of non-coding transcripts with a length of more than 200 nt. LncRNAs can regulate gene expression at transcriptional level and post-transcriptional level via cis effect, trans effect or competing endogenous RNA (ceRNA) mechanism. LncRNAs have been proved to play pivotal roles in many biological processes such as growth and development.【Objective】The objective of this study is to perform detailed investigation and discussion of the regulatory manner of lncRNAs by combining small RNA (sRNA) omics dataset obtained in this work and previously gained lncRNA omics dataset, and to reveal the putative function of lncRNAs in Nosema ceranae spores.【Method】Here, purified spores of N.ceranae were sequenced using small RNA-seq (sRNA-seq), related bioinformatics software was used to conduct quality control of sequencing data. The upstream and downstream genes of N. ceranae lncRNAs were predicted based on the positional relationship between lncRNAs and mRNAs. Annotations of these upstream and downstream genes in GO and KEGG databases were carried out using Blast software. The lncRNA-miRNA and lncRNA-miRNA-mRNA regulatory networks were predicted using Target Finder software and then visualized with Cytoscape v3.7.2 software. The expression of partial lncRNAs, miRNAs and mRNAs within regulatory networks was validated with RT-PCR.【Result】In total, 16 597 883, 15 451 791 and 12 248 316 raw reads were obtained from sRNA-seq of N. ceranae spore samples, respectively, and after quality control, 15 608 370, 14 249 255 and 11 440 684 clean reads with a mean Q30 above 98.58% were gained, respectively. A total of 310 upstream and downstream genes of lncRNAs were predicted. These genes could be annotated to 35 functional terms associated with metabolic process, cell process, catalytic activity,binding and cells, etc. Additionally, these genes could be annotated to 56 metabolic pathways, including material metabolic pathways such as purine metabolism, carbon metabolism, pyruvate metabolism; and energy metabolic pathways including methane metabolism,oxidative phosphorylation, glycolysis/gluconeogenesis and so forth. The results indicated that corresponding lncRNAs could regulate the expression level of upstream and downstream genes through cis effect, thus participating in regulation of material and energy metabolisms as well as cell activities in N. ceranae spore. Furthermore, lncRNA-miRNA regulatory networks were constructed and analyzed, the result showed that MSTRG.3636.1, MSTRG.4498.1 and MSTRG.4883.1 could bind to four miRNAs including nce-miR-7729, nce-miR-7502, nce-miR-8639 and nce-miR-8565, suggesting that these three lncRNAs as ceRNAs could exert potential function in N. ceranae spore. LncRNA-miRNA-mRNA regulatory network analysis demonstrated that complex networks existed among them, miRNAs were located in the center, connecting lncRNAs and mRNAs, nce-miR-7502 and nce-miR-8639 could bind to the most mRNAs (28). MiRNAs targeting MSTRG.3636.1, MSTRG.4498.1 and MSTRG.4883.1 could target several mRNAs,indicating that these three lncRNAs might play a role via ceRNA mechanism, thus affecting the mechanism and vital activity in N.ceranae spore.【Conclusion】LncRNAs are likely to regulate the expression of upstream and downstream genes via cis effect, and indirectly affect the expression of target genes by absorbing miRNAs as ceRNAs, thus controlling basic activities in N. ceranae spore such as material metabolism and energy metabolism.

Key words: Nosema ceranae; non-coding RNA (ncRNA); long non-coding RNA (lncRNA); upstream and downstream genes;competing endogenous RNA (ceRNA); regulatory network

开放科学（资源服务）标识码（OSID）：pagenumber_ebook=198,pagenumber_book=2122

收稿日期：2019-10-23；接受日期：2019-11-26

基金项目：国家现代农业产业技术体系建设专项资金（CARS-44-KXJ7）、福建省自然科学基金（2018J05042）、福建省教育厅中青年教师教育科研项目（JAT170158）、福建农林大学杰出青年科研人才计划（xjq201814）、福建农林大学科技创新专项基金（CXZX2017342，CXZX2017343）、福建省大学生创新创业训练计划（3165602032，3155006018）

联系方式：周丁丁，E-mail：zdd03569981@163.com。史小玉，E-mail：2395965008@qq.com。周丁丁和史小玉为同等贡献作者。

通信作者郭睿，E-mail：ruiguo@fafu.edu.cn

（责任编辑 岳梅）