蜜蜂球囊菌中长链非编码RNA的调控作用

周丁丁，范元婵，王杰，蒋海宾，祝智威，范小雪，陈华枝，杜宇，周紫彧，熊翠玲，郑燕珍，付中民，陈大福，郭睿

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

0 引言

【研究意义】蜜蜂球囊菌（Ascosphaera apis，简称球囊菌）是一种专性侵染蜜蜂幼虫的真菌性病原[1]，对西方蜜蜂（Apis mellifera）幼虫、中华蜜蜂（Apis cerana cerana，简称中蜂）雄蜂和工蜂幼虫以及成年熊蜂均具有侵染性[2-3]。球囊菌可导致蜜蜂幼虫罹患白垩病，能够引起成年蜜蜂数量和蜂群生产力大幅下降[2-5]。长链非编码 RNA（long non-coding RNA，lncRNA）是一类主要由RNA聚合酶II转录，不具备蛋白编码能力，缺少完整开放阅读框（ORF）且长度＞200 nt的转录本[6]，已被证明在真核生物的基因组印记和细胞周期等诸多生物学过程中扮演关键角色[7-8]。然而，球囊菌的lncRNA研究十分滞后。本研究基于前期获得的球囊菌菌丝和孢子混合样品的高质量lncRNA组学数据开展相关生物信息学分析和分子生物学验证，对球囊菌lncRNA的顺式（cis）作用、反义lncRNA（antisense lncRNA）作用和竞争性内源RNA（competing endogenous RNA，ceRNA）作用进行全面分析和深入探讨。研究结果可丰富球囊菌的lncRNA信息，为其他真菌的lncRNA提供参考，并揭示lncRNA在球囊菌中的潜在调控功能。【前人研究进展】LncRNA可通过顺式作用调控同一染色体上邻近蛋白编码基因的表达[9]；通过反式作用调控距离较远基因的转录激活及表达[9-10]；含有微小 RNA反应元件（microRNA response element，MRE）的lncRNA还能够通过发挥 ceRNA作用吸附结合miRNA，从而间接调节下游靶mRNA的表达[11]；还能与特定蛋白质结合形成核酸蛋白复合物，调控蛋白活性或改变蛋白在细胞中的定位，以及作为某些小RNA的前体分子[12-13]。VAN WERVEN等[14]研究发现在酿酒酵母（Saccharomyces cerevisiae）的孢子形成过程中，lncRNAIRT1通过顺式作用抑制转录因子Rme1与启动子IME1的结合，从而抑制IME1的表达；孢子形成过程中另一个lncRNAIME4-AS的表达能够抑制IME4的表达。杜庆国[15]通过对正常磷水培和低磷水培2 d和8 d的地下和地上的玉米组织进行链特异性建库的RNA-seq测序，基于测序数据预测出65 099个lncRNA，进一步利用psRobot软件预测6 787个差异表达lncRNA（differentially expressed lncRNA，DElncRNA）可靶向结合miR399，推测这些DElncRNA通过影响miR399的表达介导玉米的低磷胁迫响应。SHAO等[16]通过对生理盐水和吗啡耐受两组小鼠的比较分析，鉴定到 136个 DElncRNA，进而通过构建和分析 DElncRNA-miRNA-mRNA调控网络发现部分DElncRNA可作为潜在的ceRNA竞争性结合miRNA，从而间接影响泛素化途径、GRCP和TLR信号通路。真菌lncRNA的研究起步较晚且进展缓慢，有限的 lncRNA研究主要集中于芽殖酵母（Schizosaccharomyces cerevisiae）、裂殖酵母（Schizosaccharomyces pombe）和粗糙脉孢菌（Neurospora crassa）等少数真菌物种[14,17-19]，昆虫病原真菌的lncRNA研究未见报道。球囊菌的基因组序列信息早在2006年就已公布[20]；但其完整的基因序列和功能注释信息直到 2016年才正式公布[21]。因此，球囊菌的分子生物学及组学研究进展缓慢。前期研究中，笔者团队组装和注释了球囊菌的参考转录组[22]；分析和探讨了侵染意大利蜜蜂（Apis mellifera ligustica，简称意蜂）幼虫和中蜂幼虫的球囊菌的转录组动态和致病机理[23-24]；近期又利用PacBio单分子实时测序技术（SMRT）对实验室条件下纯培养的球囊菌进行了测序，重新组装和注释了更高质量的全长转录组，包含394 142条高可信度的全长转录本以及11 623条参考基因组未注释的全长转录本[25]，可为参考基因组的完善、新基因的全长克隆以及可变剪接和可变腺苷酸化的研究提供数据基础。【本研究切入点】笔者团队前期利用链特异性建库的 lncRNA-seq技术分别对纯化的球囊菌菌丝和孢子的混合样品进行了深度测序，对球囊菌lncRNA的数量、种类和保守性进行了分析，并全面比较了lncRNA和mRNA的结构特点[26]。但上述lncRNA在球囊菌中的作用方式及潜在调控功能尚不清楚。【拟解决的关键问题】基于已获得的高质量 lncRNA组学数据对球囊菌 lncRNA的顺式作用、反义lncRNA作用和ceRNA作用进行深入细致的分析和探讨，以期丰富球囊菌lncRNA的相关信息，并在组学水平揭示球囊菌lncRNA的潜在调控功能。

1 材料与方法

试验于2019年在福建农林大学动物科学学院（蜂学学院）蜜蜂保护实验室完成。

1.1 生物材料

球囊菌菌株由福建农林大学动物科学学院（蜂学学院）蜜蜂保护实验室分离、培养和保存[22-24]。

1.2 LncRNA组学数据来源

笔者团队前期已利用链特异性建库的 lncRNA-seq技术对球囊菌菌丝和孢子混合样品进行测序，共获得了70 612 299条原始读段（raw reads），质控后共获得48 268 696条有效读段（clean reads），平均Q30达到90.21%，测序数据质量良好[26]；利用软件CNCI、CPC、Pfam和CAPT分别预测出578、1 848、1 460和 1 662条 lncRNA，四者的交集为 379条lncRNA，即高可信度的球囊菌lncRNA[26]。高质量的lncRNA组学数据可用于本研究中lncRNA顺式作用、反义lncRNA作用和ceRNA作用的生物信息学分析。原始数据已上传NCBI SRA数据库，BioProject号：PRJNA395108。

1.3 sRNA组学数据来源

笔者团队前期已利用small RNA-seq（sRNA-seq）技术对球囊菌菌丝和孢子的混合样品进行测序，获得了70 612 299条raw reads，质控后得到48 268 696条clean reads，平均Q30达到98.75%[27]。高质量的sRNA组学数据为本研究中lncRNA与miRNA以及miRNA和mRNA靶向结合关系的预测提供数据支撑。原始数据已上传NCBI SRA 数据库，BioProject号：PRJNA560456。

1.4 LncRNA上下游基因的功能及通路注释

LncRNA位于编码蛋白基因上下游，可能与启动子、增强子和可诱导元件等产生部分重叠，在转录或转录后水平对其邻近蛋白编码基因的表达具有调控作用[8-9]。参照周丁丁等[28]的方法，将lncRNA上下游100 kb范围内的邻近基因作为其调控的靶基因，利用Blast软件（采用默认参数）将lncRNA的上下游基因比对到GO数据库（http://geneontology.org/）和KEGG数据库（https://www.kegg.jp/），获得上下游基因的功能和通路注释。

1.5 LncRNA序列互补的mRNA预测及分析

LncRNA序列可以与对应的靶mRNA序列碱基互补配对，在配对的过程中，可以携带某些作用因子到靶mRNA，进而调控相关编码蛋白基因的表达[11]。参照周丁丁等[28]的方法，利用 LncTar软件[29]对lncRNA的靶 mRNA进行预测，基于 lncRNA和mRNA之间的互补序列，计算配对位点自由能和标准化自由能，标准化自由能阈值以下的则认为是lncRNA的靶mRNA。利用Blast软件（采用默认参数）将lncRNA的互补序列基因比对到eggNOG数据库（http://eggnogdb.embl.de/）和KEGG数据库，获得相应的蛋白功能注释信息。

1.6 LncRNA的ceRNA调控网络构建及分析

参照郭睿等[30]的方法，利用Target Finder软件[31]预测球囊菌lncRNA靶向结合的miRNA以及miRNA靶向结合的mRNA，利用三者之间的靶向结合关系构建lncRNA-miRNA及lncRNA-miRNA-mRNA调控网络，并通过Cytoscape v 3.7.1软件[28,31]对上述调控网络进行可视化。

1.7 LncRNA和mRNA的RT-PCR验证

根据1.6中靶向预测的结果，随机选取9个调控网络中的lncRNA（MSTRG.6443.1、MSTRG.2135.1、MSTRG.2134.1、MSTRG.3422.1、MSTRG.2302.1、MSTRG.5023.1、MSTRG.5584.1、MSTRG.1870.1 和MSTRG.1614.1）和8个靶mRNA（gene2674、gene4126、gene4125、gene1602、gene5384、gene1986、gene989、和gene1970）进行RT-PCR验证。根据上述lncRNA和mRNA的核酸序列，利用Primer Primer 5.0软件设计相应的特异性引物，委托上海生工生物工程有限公司合成引物（表 1）。选取actin（gene6001）作为 lncRNA和mRNA的内参基因（阳性对照）。将9对lncRNA（或8对mRNA）特异性上游引物（1 μL/种）和下游引物（1 μL/种）混合，以无菌水为模板进行RT-PCR，作为阴性对照。利用RNA抽提试剂盒（TaKaRa公司，日本）分别提取球囊菌菌丝的总 RNA和孢子的总RNA，将1 μg菌丝RNA与1 μg孢子RNA等量混合，作为模板进行反转录，得到的 cDNA作为模板进行PCR。PCR反应体系（20 μL）：cDNA模板1 μL，上游引物 1 μL，下游引物 1 μL，Mixture 10 μL，无菌水7 μL。PCR程序：95℃预变性 5 min；95℃变性 30 s，55℃退火30 s，72℃延伸30 s，共34个循环；72℃再延伸5 min。PCR产物经1.5%琼脂糖凝胶电泳和凝胶成像仪（上海培清，中国）检测。

1.8 MiRNA的Stem-loop RT-PCR验证

根据1.6中靶向预测的结果，随机选取7个调控网络中的靶miRNA（miR-281-y、miR-13-y、miR-11980-y、miR-6057-y、miR-1-z、miR-9-z 和miR-13-x）参照郭睿等[27]和杜宇等[32]的方法，利用DNAMAN软件（Lynnon Biosoft 公司，美国）设计miRNA的Stem-loop引物、特异性上游引物和通用下游引物，委托上海生工生物工程有限公司合成引物（表1）。选取actin（gene6001）作为miRNA的内参基因（阳性对照）。将7对miRNA特异性上游引物（1 μL/种）和通用下游引物（1 μL/种）混合，以无菌水为模板进行RT-PCR，作为阴性对照。利用RNA抽提试剂盒（TaKaRa公司，日本）分别提取球囊菌菌丝的总RNA和孢子的总RNA，将1 μg菌丝RNA与1 μg孢子RNA等量混合作为反转录的模板。利用Stem-loop引物反转录得到 cDNA，作为模板进行 PCR扩增。PCR 体系（20 μL）：上下游引物（1.67 mol·L-1）各 1 μL，cDNA 模板 1 μL，PCR mix 10 μL，无菌水 7 μL。PCR程序：95℃预变性5 min；95℃变性30 s，梯度退火30 s，72℃延伸30 s，共34个循环；72℃再延伸5 min。PCR产物经1.5%琼脂糖凝胶电泳和凝胶成像仪（上海培清，中国）检测。

2 结果

2.1 球囊菌lncRNA的顺式作用分析

球囊菌的371个lncRNA共预测出5 852个上下游基因。GO数据库注释结果显示，这些上下游基因可注释到48个功能条目，涉及生物学进程、分子功能和细胞组分 3大类；其中生物学进程相关的前 5个条目分别是细胞进程（1 707）、代谢进程（1 705）、单组织进程（1 231）、定位（473）和生物学调控（325）；细胞组分相关的前5个条目分别是细胞（1 163）、细胞组件（1 163）、细胞器（795）、细胞膜（532）和大分子复合物（501）；分子功能相关的前5个条目分别是催化活性（1 796）、结合（1 388）、转运活性（180）、分子结构活性（106）和核酸结合转化因素（40）（图1）。括号内的数字表示注释到该条目的上下游基因数。

KEGG数据库注释结果显示，这些上下游基因可注释到121条通路，注释基因数最多的前10位分别为新陈代谢途径（564）、次生代谢产物的生物合成（233）、抗生素的生物合成（176）、氨基酸的生物合成（89）、碳代谢（81）、核糖体（78）、嘌呤代谢（75）、剪接体（66）、RNA转运（65）和细胞周期-酵母（62）。括号内的数字表示注释到该通路的上下游基因数。

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

2.2 球囊菌中反义lncRNA分析

共预测出球囊菌的7个lncRNA与7个mRNA序列互补，它们的结合关系如图2所示。其中，有1个靶mRNA（gene3444）在KEGG数据库注释为核孔复合体蛋白An-Nup120和假定蛋白（K14303），另有1个靶mRNA（gene1554）注释为假定蛋白（K07117）；共有 5个靶基因（gene1316、gene1194、gene4650、gene5233和gene1554）在eggNOG数据库中注释为假定蛋白（表2）。

2.3 球囊菌lncRNA的ceRNA作用分析

共预测出227个lncRNA靶向结合73个miRNA，二者之间形成较为复杂的调控网络，其中MSTRG.2597.1、MSTRG.4497.1和MSTRG.1789.9靶向结合的miRNA数量较多，分别为7、7和6个；但多数 lncRNA（79.02%）仅能结合 1—2个 miRNA；部分 miRNA可靶向多个 lncRNA，其中 miR-34-x、miR-10-x和miR-375-y结合的lncRNA数量最多，分别为31、27和24个（图3）；但多数miRNA（80.82%）可靶向的lncRNA数少于10个。

表2 球囊菌中与lncRNA序列互补靶标的功能注释Table 2 Functional annotation of complementary targets with lncRNA sequence in A. apis

2.4 氧化磷酸化和MAPK信号通路相关的ceRNA调控网络

根据lncRNA与miRNA、miRNA与mRNA的靶向结合关系构建 lncRNA-miRNA-mRNA调控网络，发现三者之间具有更为复杂的调控关系，227个lncRNA靶向73个miRNA，这些miRNA可结合50个靶mRNA。对上述靶mRNA进行KEGG数据库注释，结果显示它们共注释到28条通路，注释最多的通路是新陈代谢途径（9），其次为次生代谢产物的生物合成（4）、自噬-其他真核生物（3）、线粒体自噬-酵母（3）和自噬-酵母（3）（表3）。进一步分析发现有10个靶mRNA可注释到9条物质能量代谢通路，分别为谷胱甘肽代谢（2）、脂肪酸延长（1）、不饱和脂肪酸的生物合成（1）、鞘脂类代谢（1）、脂肪酸代谢（1）、碳代谢（1）、氧化磷酸化（1）、丙酸代谢（1）、糖酵解和糖质生（1）。括号内的数字代表注释在该通路的靶mRNA数。

根据靶mRNA对应的Nr数据库蛋白注释信息，筛选与氧化磷酸化通路和 MAPK信号通路相关的靶mRNA及与上述靶 mRNA存在靶向结合关系的miRNA及相应lncRNA，分析发现氧化磷酸化通路相关调控网络包含222个lncRNA、78个靶miRNA和50个靶mRNA，其中MSTRG.2597.1、MSTRG.4497.1及MSTRG.1789.9结合的miRNA数最多，分别为7、7和6个（图4-A）；MAPK信号通路相关调控网络包含222个lncRNA、76个靶miRNA和46个靶mRNA，其中MSTRG.4497.1、MSTRG.2597.1及MSTRG.1789.9结合的miRNA数最多，分别为7、7和6个（图4-B）；共有217个lncRNA和66个miRNA同时涉及上述两条通路的调控网络。

2.5 球囊菌ceRNA调控网络中lncRNA、miRNA和mRNA的表达验证

为验证ceRNA网络中 lncRNA、靶miRNA和靶mRNA的表达，通过RT-PCR对上述ceRNA网络中lncRNA和靶mRNA的进行扩增，扩增产物的电泳结果显示9个lncRNA均能扩增出预期大小的目的条带，阳性对照能扩增出目的片段而阴性对照未扩增出片段；8个靶mRNA能扩增出预期片段，阳性对照能扩增出目的片段而阴性对照未扩增出片段；上述结果证实了lncRNA和靶mRNA的真实表达。利用Stem-loop RT-PCR对7个miRNA进行检测，电泳结果显示它们均能扩增出预期大小的目的片段，阳性对照能扩增出目的片段而阴性对照未扩增出片段（图5），证实了miRNA的真实表达。上述结果表明球囊菌中lncRNA、靶miRNA和靶mRNA真实表达。

表3 球囊菌ceRNA网络中靶mRNA注释数前12位的KEGG通路Table 3 Top 12 KEGG pathways annotated by target mRNAs involved in A. apis ceRNA networks

3 讨论

球囊菌侵染蜜蜂幼虫导致白垩病，该病危害蜜蜂健康并长期困扰养蜂生产。相对于酵母等模式真菌，球囊菌的分子生物学及组学研究较为滞后，miRNA和lncRNA等ncRNA的相关研究更为匮乏。目前，较多的研究表明 lncRNA在人类、黑腹果蝇（Drosophila melanogaster）、秀丽隐杆线虫（Caenorhabditis elegans）、拟南芥（Arabidopsis thaliana）和酵母（Saccharomyces）等模式生物中发挥特殊而重要的调控功能[8,10-14,33-35]。笔者团队前期已对球囊菌的miRNA进行了全基因组鉴定、分析和验证[27]；并对侵染中华蜜蜂6日龄幼虫的球囊菌的miRNA差异表达谱及调控网络进行了解析[36]。近期，为最大限度鉴定球囊菌的 lncRNA，笔者团队利用链特异性建库的 lncRNA-seq技术对球囊菌菌丝和孢子混合样品进行了深度测序，通过生物信息学方法对lncRNA的数量、种类、结构特征进行了初步分析，鉴定到379个lncRNA，包括242个反义lncRNA、123个基因间区lncRNA（lincRNA）、13个正义链 lncRNA（sense lncRNA）和 1个内含子 lncRNA（intronic transcript lncRNA）[26]；并对部分lncRNA的真实表达进行了分子验证[26]；此为球囊菌lncRNA的首次研究报道。基于已获得的高质量lncRNA组学数据，本研究进一步探究球囊菌lncRNA的顺式作用、反义lncRNA作用和ceRNA作用，发现lncRNA在球囊菌的物质和能量代谢、遗传信息传递、环境刺激应答、信号通路、细胞生长和细胞周期等生物学过程中具有潜在的调控功能。

3.1 LncRNA对球囊菌生长发育以及物质能量代谢具有潜在的顺式调控作用

LncRNA发挥顺式作用主要是通过影响同一染色体上的邻近蛋白编码基因的表达发挥调控功能[9]。CAI等[37]利用RNAi技术和细胞实验对7周龄的兴化雌鸡各组织和器官进行研究，发现位于Six1基因上游432 bp处的lncRNA-Six1可编码一个约7.26 kD的小肽，进而通过顺式作用影响Six1的表达以促进鸡肉细胞增殖及诱导细胞分裂。FENG等[38]利用lncRNA-seq技术对长牡蛎（Crassostrea gigas）的白色、黑色、金色和正常外壳进行深度测序，预测出1 157个lncRNA及潜在调控的24 057个上下游基因，进一步分析结果显示427个DElncRNA可靶向2 088个上下游基因，这些上下游基因涉及RNA甲基转移酶和tRNA甲基转移酶活性，以及ECM-受体相互作用、Jak-STAT信号通路和Notch信号通路，作者推测相关DElncRNA通过顺式作用调节牡蛎壳的色素沉淀。本研究中，371个lncRNA共预测出5 852个上下游基因，其中分别有1 705、64、31和61个上下游基因注释到代谢进程、发育进程、生长和繁殖。还发现分别有1 163、1 163和795个上下游基因注释到细胞、细胞组件和细胞器，表明lncRNA可能通过顺式作用调节球囊菌的细胞生命活动。

碳代谢在真菌生长和发育中扮演关键角色，其中一碳代谢在机体代谢活动发挥重要作用；而中心碳代谢是生物体所需能量的主要来源，并为其他代谢提供前体反应物，糖酵解途径和磷酸戊糖途径都属于中心碳代谢[39]。WANG等[40]利用转录组测序技术和分子生物学手段对野生型和突变型的谷氨酸棒杆菌（Corynebacterium glutamicum）进行深入分析，通过13C代谢通量分析，发现两个转录因子gntR1-E70K和ramA-A52V可促进突变型谷氨酸棒杆菌的糖酵解和磷酸戊糖途径相关基因上调表达，说明这两个转录因子加快了谷氨酸棒杆菌的中心碳代谢。本研究中，分别有81、39和21个上下游基因注释到碳代谢、糖酵解/糖异生途径和磷酸戊糖途径，表明相应的 lncRNA可能通过顺式作用调控球囊菌的一碳代谢和中心碳代谢。此外，分别有57、21和11个上下游基因注释到氧化磷酸化、甲烷代谢和硫代谢等能量代谢途径，说明相应的lncRNA可能通过顺式作用调控上述能量代谢通路。

3.2 LncRNA对球囊菌次级代谢产物的合成与降解以及过氧化物酶体具有潜在的顺式调控作用

生物体能够以糖类、氨基酸和核苷酸等初级代谢产物为前体，合成一些对于生长发育不可或缺的次级代谢产物，例如苯丙素类、萜类和黄酮类化合物[41-42]。本研究中，分别有30、27、20、16、13、12和7个上下游基因涉及氨基糖和核苷糖代谢、缬氨酸和亮氨酸及异亮氨酸的降解、果糖和甘露糖代谢、烟酸和烟酰胺代谢、赖氨酸降解、鞘脂类代谢和维生素B6代谢，表明相应的lncRNA通过调控上述初级代谢产物的加工，合成必要的次级代谢产物供球囊菌的生长和发育；此外，分别有233、15、7、4、1和1个lncRNA上下游基因涉及次级代谢物生物合成、萜类化合物生物合成、泛醌和其他烯萜类生物合成、酮体的合成与降解、碳青霉烯类生物合成和倍半萜三萜类生物合成，再次说明相关lncRNA通过调控次级代谢产物的合成与降解参与对球囊菌生长和发育的调节。HUARTE-BONNET等[43]利用含烷烃培养基培养球孢白僵菌（Beauveria bassiana），发现过氧化物酶体生物合成途径相关的Bbhyd1和Bbhyd2两个基因被激活，在过氧化物酶活性高的菌丝团形成菌丝球过程中，球孢白僵菌表现出过氧化物酶体增多和表面明显增厚的现象。本研究中，有32个上下游基因注释到过氧化物酶体，说明相关lncRNA（MSTRG.4299.3、MSTRG.6298.3和MSTRG.6220.4等）对该通路具有潜在的顺式调控作用。

3.3 LncRNA对球囊菌的环境适应及细胞生化反应具有潜在的顺式调控作用

应激反应是生物体对外界刺激的一种适应性机制，真菌在面对营养物质缺乏、温度胁迫、物理化学因子刺激和机械损伤时可在一定环境压力范围内通过自身应激调节来适应和生存[44]。金属伴侣蛋白是一类胞内可溶性金属结合蛋白，可以将金属离子准确地交付给目标蛋白，促使机体激活金属酶，同时也保证细胞免受金属危害[45]。本研究发现，分别有1 796、247、106、19、3和3个上下游基因涉及催化活性、应激反应、结构性分子活性、信号传感器活性、分子传感器活性和金属伴侣蛋白活性，说明相应的lncRNA可能通过顺式作用参与调控球囊菌的环境适应及细胞内的生化反应，以及保护细胞抵御外界不良因素的影响。

3.4 球囊菌反义lncRNA对核孔复合体蛋白的合成等生物学过程具有潜在调控作用

有些 lncRNA能够与一些距离较远的基因序列互补配对，也可以携带某些调控因子到互补的靶mRNA上，通过影响基因的表达发挥调控功能[9,28]。CHACKO等[46]基于Northern blot、单分子荧光原位杂交和RNA-seq等技术对新型隐球酵母（Cryptococcus neoformans）进行深入研究，发现lncRNARZE1影响新型隐球酵母的关键基因ZNF2的转录，RZE1可以调控新型隐球酵母菌丝的形成。本研究预测到 7个lncRNA（MSTRG.3576.1、MSTRG.2795.1、MSTRG.1238.1、MSTRG.5393.1、MSTRG.1387.1、MSTRG.1672.1和 MSTRG.3829.1）与7个靶 mRNA（gene1316、gene548、gene1194、gene3444、gene4650、gene5233和gene1554）之间存在序列互补；分析发现其中5个靶mRNA对应的蛋白在eggNOG数据库中注释为假定蛋白（表 2），说明球囊菌的基因注释信息尚不完善，基因组质量偏低和分子生物学研究滞后是主要原因。仅有gene3444对应的蛋白注释为核孔复合体蛋白（表2），表明相应的lncRNA（MSTRG.5393.1）通过与相应 mRNA序列互补调控球囊菌核孔复合体蛋白的合成，可能会影响核苷酸、酶和转录因子进出细胞核，以及遗传信息的传递。然而，本研究并未发现可以作为miRNA前体的lncRNA。鉴于本研究的组学数据来源于实验室条件下的球囊菌纯培养，而lncRNA的表达具有细胞、组织、发育阶段和胁迫时期特异性[47]，因此对于球囊菌 lncRNA是否能通过作为miRNA前体发挥更为广泛而灵活的调控作用，仍需要进一步研究。下一步将利用链特异性建库的lncRNA-seq技术对球囊菌侵染的不同日龄蜜蜂幼虫肠道进行测序，通过连续比对核糖体数据库、宿主基因组和球囊菌基因组筛滤得到纯净的病原lncRNA数据，进而开展处于侵染过程的球囊菌的lncRNA调控作用研究。

3.5 球囊菌lncRNA对物质能量代谢以及生长发育具有潜在的ceRNA调控作用

ceRNA假说认为，具有 MRE的不同种类 RNA如lncRNA、circRNA和假基因等，可通过竞争性结合miRNA间接调节基因表达[48]。目前，该假说已被越来越多的研究证实。周丁丁等[28]研究发现东方蜜蜂微孢子虫（Nosema ceranae）的lncRNA通过ceRNA作用调控孢子中的物质和能量代谢等基本生命活动。笔者团队前期系统分析和探讨了意蜂工蜂中肠发育过程的lncRNA差异表达谱及调控网络[30]、意蜂工蜂中肠响应东方蜜蜂微孢子虫胁迫的lncRNA差异表达谱及调控网络[49]，发现相关DElncRNA可能通过ceRNA作用参与意蜂工蜂中肠、意蜂幼虫肠道的发育及真菌病原胁迫应答过程。氧化磷酸化是提供真菌生命活动所需能量的基础代谢途径[50]。MAPK信号通路在酵母等真菌的信号转导、生长发育、胁迫应答和致病性等方面起到关键作用[51]。本研究中，球囊菌的 227个lncRNA与92个miRNA及96个mRNA之间形成的调控网络十分复杂，因此针对氧化磷酸化通路和MAPK信号通路相关的 lncRNA、靶 miRNA和靶mRNA构建调控网络，分析发现 MSTRG.2597.1、MSTRG.4497.1和MSTRG.1789.9等球囊菌lncRNA靶向结合 miRNA数量较多；共有 217个 lncRNA（MSTRG.1135.4、MSTRG.2904.2和 MSTRG.5757.1等）同时参与上述两个通路的调控；表明相关的球囊菌lncRNA可能作为ceRNA间接调节靶mRNA的表达。进一步分析发现，调控网络中的靶mRNA可注释到碳代谢和脂肪酸代谢等18条物质能量代谢通路，以及细胞周期和减数分裂。上述结果说明相应的球囊菌lncRNA可能通过ceRNA作用对物质和能量代谢以及生殖进行调控。本研究中，MSTRG.3576.1、MSTRG.1238.1和MSTRG.3829.1同时参与了顺式作用、反义lncRNA作用和ceRNA作用，说明这3个lncRNA作为多功能的调控因子，可通过多种方式对球囊菌的生物学过程进行灵活调控。

4 结论

球囊菌的371个lncRNA通过顺式作用潜在调控5 852个上下游基因的表达，227个lncRNA可通过92个 miRNA间接调控 96个靶 mRNA的表达，部分lncRNA可能通过顺式作用和ceRNA作用调节球囊菌的生长发育及物质能量代谢；球囊菌的7个 lncRNA可作为反义RNA调控7个互补配对的mRNA的表达，MSTRG.5393.1潜在参与调控球囊菌的核孔复合体的蛋白合成。