東方蜜蜂微孢子蟲孢子中長鏈非編碼RNA的競爭性內源RNA調控網絡及潛在功能

周丁丁，史小玉，王杰，范元嬋，祝智威，蔣海賓，范小雪，熊翠玲，鄭燕珍，付中民，徐國鈞，陳大福，郭睿

（福建農林大學動物科學學院（蜂學學院），福州 350002）

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。

對于下機的原始數據，因含有接頭序列或低質量序列，首先以以下流程進行數據質控：（1）去除低質量值的讀段（reads）；（2）去除未知堿基（N）含量≥10%的reads；（3）去除沒有3′接頭序列的reads；（4）剪切掉3′接頭序列；（5）去除長度＜18 nt和＞30 nt的序列。經嚴格質控后得到的數據即為有效讀段（clean reads）。

進一步通過 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上下游基因的功能注釋

基于東方蜜蜂微孢子蟲純化孢子的 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/）注釋。

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

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）。括號內的數字代表注釋到該功能條目的上下游基因數。

圖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

表3 sRNA-seq數據比對東方蜜蜂微孢子蟲參考基因組的信息統計Table 3 Summary of mapping information of sRNA-seq data to the reference genome of N. ceranae

2.3 東方蜜蜂微孢子蟲孢子中lncRNA上下游基因的KEGG數據庫注釋

進一步對lncRNA的上下游基因進行KEGG數據庫注釋，結果顯示310個上下游基因可注釋到56條通路，注釋基因數最多的前10位分別是新陳代謝途徑（25）、抗生素的生物合成（12）、次生代謝產物的生物合成（11）、核糖體（10）、嘌呤代謝（7）、蛋白酶體（7）、RNA轉運（7）、真核生物中核糖體的生物發生（7）、碳代謝（6）和嘧啶代謝（6）（圖2）。括號內的數字代表注釋到該通路的上下游基因數。

圖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）。

圖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真實表達。

圖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

圖5 東方蜜蜂微孢子蟲孢子中lncRNA、靶miRNA和靶mRNA的RT-PCR驗證Fig. 5 RT-PCR verification of lncRNAs, target miRNAs and target mRNAs in N. ceranae spore

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作用影響靶基因的表達水平，從而對東方蜜蜂微孢子蟲孢子的新陳代謝和基本生命活動進行靈活調控。