微小RNA介導東方蜜蜂微孢子蟲侵染意大利蜜蜂工蜂的分子機制

耿四海，石彩云，范小雪，王杰，祝智威，蔣海賓，范元嬋，陳華枝，杜宇，王心蕊，熊翠玲，鄭燕珍，付中民，陳大福，郭睿

微小RNA介導東方蜜蜂微孢子蟲侵染意大利蜜蜂工蜂的分子機制

耿四海，石彩云，范小雪，王杰，祝智威，蔣海賓，范元嬋，陳華枝，杜宇，王心蕊，熊翠玲，鄭燕珍，付中民，陳大福，郭睿

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

【】通過對東方蜜蜂微孢子蟲（）純化孢子與侵染意大利蜜蜂（，簡稱意蜂）工蜂的東方蜜蜂微孢子蟲的差異表達miRNA（DEmiRNA）及其靶mRNA進行系統分析，篩選、分析和探討病原毒力因子和侵染因子相關的DEmiRNA及調控網絡，在miRNA組學層面揭示東方蜜蜂微孢子蟲對意蜂的侵染機制。利用small RNA-seq（sRNA-seq）技術對東方蜜蜂微孢子蟲感染7 d和10 d的意蜂工蜂中腸和東方蜜蜂微孢子蟲純化孢子（NcCK）進行深度測序，通過連續比對rRNA數據庫、西方蜜蜂（）基因組和東方蜜蜂微孢子蟲基因組篩濾出處于侵染過程的東方蜜蜂微孢子蟲（NcT1和NcT2）數據和東方蜜蜂微孢子蟲孢子的測序數據。根據≤0.05，｜log2fold change｜≥1的標準，通過比較分析篩選出各比較組中的差異表達miRNA（differentially expressed miRNA，DEmiRNA）。通過相關生物信息學軟件對DEmiRNA進行表達譜分析，靶mRNA預測及功能和代謝通路注釋，以及調控網絡的構建與分析。通過Stem-loop RT-qPCR驗證DEmiRNA的差異表達趨勢及測序數據的可靠性。NcCKNcT1、NcCKNcT2和NcT1NcT2比較組分別包含164、122和60個DEmiRNA。Venn分析結果顯示，3個比較組共有的上調和下調miRNA分別為5和6個。上述DEmiRNA分別預測出1 885、1 733和1 524個靶mRNA。這些靶mRNA分別注釋到27、25和26個功能條目，其中注釋數量最多的是新陳代謝進程、催化活性、細胞進程、結合和細胞。上述靶mRNA可分別注釋到84、84和84條代謝通路，其中注釋數量最多的是代謝途徑、核糖體和次級代謝產物生物合成。此外，對于NcCKNcT1、NcCKNcT2和NcT1NcT2中的DEmiRNA，分別有35、26和12個靶向結合MAPK信號通路相關靶mRNA，分別有49、40和17個DEmiRNA靶向結合糖酵解/糖異生通路相關靶mRNA。進一步分析發現，東方蜜蜂微孢子蟲的DEmiRNA參與調控蓖麻毒素B凝集素、細胞凋亡抑制因子、極管蛋白和孢壁蛋白等病原毒力因子的基因表達，以及己糖激酶、ATP/ADP移位酶、ABC轉運蛋白和轉錄因子ste12等侵染因子的基因表達。通過對東方蜜蜂微孢子蟲純化孢子與侵染意蜂工蜂的東方蜜蜂微孢子蟲進行深入細致的miRNA組學分析和探討，解析了病原侵染過程的miRNA差異表達譜，揭示了東方蜜蜂微孢子蟲可能通過調節相應miRNA的表達水平對蓖麻毒素B凝集素、細胞凋亡抑制因子、極管蛋白和孢壁蛋白等毒力因子及己糖激酶、ATP/ADP移位酶、ABC轉運蛋白和轉錄因子ste12等侵染因子基因表達進行調控，從而適應宿主細胞內的環境并促進自身的增殖與侵染。

東方蜜蜂微孢子蟲；意大利蜜蜂；微小RNA；侵染機制；調控網絡

0 引言

【研究意義】微孢子蟲是一類細胞內寄生的真菌，其寄主范圍非常廣泛，包括哺乳動物、魚類和昆蟲等[1]。自從1857[2]首次在家蠶（）中鑒定出家蠶微孢子蟲（）以來，目前已在脊椎動物和無脊椎動物中鑒定出超過1 400種微孢子蟲[3]。東方蜜蜂微孢子蟲（）專性寄生成年蜜蜂中腸上皮細胞[4]，可引起宿主腸道組織破壞、消化系統紊亂、免疫基因抑制、細胞凋亡抑制、壽命縮短以及工蜂進入采集工作的日齡提前等[5]。感染東方蜜蜂微孢子蟲的蜂群罹患蜜蜂微孢子蟲病，出現生產力下降和蜂群群勢衰弱等現象[6]，嚴重危害養蜂生產。1996年，Fries等[7]在北京附近地區的東方蜜蜂（）樣本中首次分離鑒定出東方蜜蜂微孢子蟲，隨后該微孢子蟲迅速蔓延至歐洲地區和中國臺灣等地的西方蜜蜂（）蜂群[8-9]，目前已廣泛感染世界主要養蜂國家的蜂群。然而對于蜜蜂微孢子蟲病，目前仍缺乏有效的治療手段，東方蜜蜂微孢子蟲侵染機制研究的滯后是最主要的限制因素。已有研究表明微小RNA（microRNA，miRNA）參與調控真菌的侵染、繁殖和毒力[10-11]。通過深度測序和組學分析解析東方蜜蜂微孢子蟲在侵染意大利蜜蜂（，簡稱意蜂）過程中的miRNA差異表達譜及調控網絡，可為闡明東方蜜蜂微孢子蟲的侵染機制提供數據基礎和理論依據，并為蜜蜂微孢子蟲病的治療提供潛在的分子靶點?！厩叭搜芯窟M展】東方蜜蜂微孢子蟲在體外以休眠態的孢子形式存在，孢子壁由幾丁質和孢壁蛋白組成，能夠抵御外界的不良環境，當孢子被蜜蜂攝取之后，在蜜蜂中腸的特殊理化環境的刺激下萌發，內部高度壓縮盤旋的極絲彈射而出，刺穿中腸上皮細胞并將有侵染性的孢原質注入，隨后進入繁殖期，利用宿主的物質與能量合成系統進行增殖，孢子數量逐漸增多，最終導致宿主細胞裂解，釋放出的成熟孢子繼續感染臨近的上皮細胞，或者隨糞便排出體外，感染新的宿主[12]。MiRNA是一類長度在18—25 nt的單鏈小RNA，由細胞內源產生的發卡結構轉錄本加工而來。1993年，Lee等[13]首次在線蟲中發現一種長度為22 nt的小RNA并命名為lin-4，lin-4可以調節線蟲的發育。let-7為第二個被發現的miRNA，let-7同樣可通過靶向結合靶mRNA調節線蟲的發育[14]。隨后在動植物中陸續鑒定出大量的miRNA[15-16]。越來越多的研究結果證實miRNA在生物體的新陳代謝、免疫應答、細胞分化、細胞周期等重要的生理生化過程中發揮著舉足輕重的作用[15-16]。然而，真菌miRNA的研究相對滯后。Lee等[17]于2010年第一次在粗糙脈孢菌（）中發現并報道類微小RNA（microRNA-like），此后在綠僵菌（）[18]、里氏木酶（）[11]和蜜蜂球囊菌（，簡稱球囊菌）[19]等中相繼鑒定出miRNA。有研究報道miRNA參與調控真菌的菌絲生長和毒力，但關于miRNA如何參與調控尚未明確。目前，東方蜜蜂微孢子蟲的miRNA研究極為有限，相關信息十分缺乏。近期，EVANS等[20]研究發現東方蜜蜂微孢子蟲的5個miRNA不僅能靶向其本身的1 545個基因，還能靶向宿主的947個基因。MiRNA的形成過程需要Dicer酶進行加工，多數微孢子蟲在進化過程中丟失了RNAi通路中、和依賴RNA為模板的RNA聚合酶（RNA-dependent RNA polymerase，RdRP）基因等關鍵基因，但東方蜜蜂微孢子蟲基因組中還保留有上述基因[21]。Huang等[22]利用RNAi敲減東方蜜蜂微孢子蟲的后，發現不但東方蜜蜂微孢子蟲的增殖水平顯著降低，而且侵染過程中病原基因的表達量也發生了變化，表明miRNA對東方蜜蜂微孢子蟲的侵染和增殖起到重要的調節作用。【本研究切入點】意蜂作為西方蜜蜂的亞種之一，是我國養蜂生產使用的主要蜂種，具有重要的經濟和生態價值。前期研究中，筆者團隊利用RNA-seq技術和生物信息學分析方法在mRNA組學水平初步探究了東方蜜蜂微孢子蟲對意蜂的侵染機制[23]。但東方蜜蜂微孢子蟲的侵染機制尚未闡明，miRNA在東方蜜蜂微孢子蟲侵染過程中的作用仍不清楚?！緮M解決的關鍵問題】利用small RNA-seq（sRNA-seq）技術對東方蜜蜂微孢子蟲純化孢子及東方蜜蜂微孢子蟲感染7 d和10 d的意蜂工蜂中腸進行深度測序，基于高質量的測序數據進行東方蜜蜂微孢子蟲miRNA的生物信息學預測，結構特征和表達譜分析，以及差異表達miRNA（differentially expressed miRNA，DEmiRNA）的靶mRNA預測、功能注釋、調控網絡構建及分析，進一步篩選和探討靶向病原毒力因子和侵染相關因子的DEmiRNA，以期在miRNA組學層面揭示DEmiRNA介導東方蜜蜂微孢子蟲的侵染機制。

1 材料與方法

試驗于2017年9月至2019年8月在福建農林大學動物科學學院（蜂學學院）蜜蜂保護實驗室完成。

1.1 生物材料

意蜂工蜂取自福建農林大學動物科學學院（蜂學學院）教學蜂場；東方蜜蜂微孢子蟲感染的意蜂外勤蜂取自福州市閩侯縣荊溪源安養蜂場。

1.2 東方蜜蜂微孢子蟲孢子的純化

參照Guo等[24-25]的方法進行東方蜜蜂微孢子蟲孢子的粗提和純化，流程簡述如下：（1）抓取東方蜜蜂微孢子蟲感染嚴重的意蜂外勤蜂若干只，拉取中腸放入研缽充分研磨，經濾網過濾后的濾液轉移到干凈的EP管，4℃，3 000 r/min離心15 min，棄上清；（2）加無菌水重懸底部白色沉淀，4℃，5 000 r/min離心15 min，棄上清，保留沉淀，如此重復兩次；（3）將25%、50%、75%、100%的percoll密度梯度溶液按密度從高到低的順序依次緩緩加入干凈的EP管，每種密度200 μL；用25%的percoll溶液重懸孢子，并緩緩加在上述密度梯度溶液的上方，4℃，14 000 r/min離心30 min；（4）利用無菌的注射器針頭小心吸取乳白色的孢子帶，并轉移到干凈的EP管，即得到東方蜜蜂微孢子蟲的純凈孢子。將純化得到的純凈孢子分為兩部分，一部分經液氮速凍后保存于-80℃冰箱，作為對照組（NcCK）；另一部分保存于4℃冰箱，用于飼喂接種意蜂工蜂。

1.3 東方蜜蜂微孢子蟲孢子的飼喂接種及意蜂工蜂中腸樣品的制備

按照本實驗室已建立的方法[26-28]進行意蜂工蜂的飼喂接種及中腸樣品的制備，流程簡述如下：（1）選取群勢較強且外觀健康的蜂群，參照文獻報道合成克什米爾蜜蜂病毒（KBV）[29]、黑蜂王臺病毒（BQCV）[30]、囊狀幼蟲病病毒（SBV）[31]、蜜蜂急性麻痹病毒（ABPV）[29]、慢性蜜蜂麻痹病病毒（CBPV）[32]、以色列急性麻痹病毒（IAPV）[31]的特異性引物，同時設計合成蜜蜂微孢子蟲（）[33]和東方蜜蜂微孢子蟲[33]的特異性引物進行PCR檢測，檢測結果顯示上述6種常見蜜蜂病毒、蜜蜂微孢子蟲和東方蜜蜂微孢子蟲均為陰性，將該蜂群作為本研究的實驗蜂群；（2）將實驗蜂群剛要出房的封蓋子脾迅速提至實驗室，放入（34±0.5）℃培養箱，將剛出房的工蜂（當天記為0 d）放入干凈的塑料盒（35只/盒），共計6盒，每盒上方插入一支裝有50%（w/v）無菌蔗糖溶液的飼喂器，（34±0.5）℃飼養過夜；（3）將上述工蜂饑餓處理2 h，然后其中3盒工蜂用于飼喂接種東方蜜蜂微孢子蟲，每只工蜂飼喂接種5 μL蔗糖溶液（含106個東方蜜蜂微孢子蟲孢子）；（4）每日清理死去的蜜蜂，于接種后7 d和10 d分別快速拉取工蜂中腸，放入干凈的1.5 mL RNA-Free離心管，每6只中腸置于1個離心管，經液氮速凍后保存于-80℃冰箱備用。將侵染意蜂后7 d和10 d工蜂中腸的東方蜜蜂微孢子蟲分別命名為NcT1（NcT1-1、NcT1-2、NcT1-3）和NcT2（NcT2-1、NcT2-2、NcT2-3）。

1.4 RNA提取、sRNA文庫構建及深度測序

參照筆者團隊已建立的流程[34]，利用Trizol法提取上述3個孢子樣品（NcCK1、NcCK2、NcCK3）和6個中腸樣品的總RNA；通過瓊脂糖凝膠電泳選擇18—30 nt的片段進行切膠回收，并連接3′接頭序列；連接產物通過15%的變性PAGE膠電泳進行分離，選擇36—44 nt目的片段進行切膠回收，再連接5′接頭序列，然后對已連接兩側接頭序列的小RNA進行反轉錄，產物通過3.5%的瓊脂糖凝膠電泳分離，進而切膠回收140—160 bp區域條帶。委托廣州基迪奧生物科技有限公司對建好的cDNA文庫進行上機測序，測序平臺為Illumina MiSeq。已將測序原始數據上傳NCBI SRA數據庫，BioProject號：PRJNA395264（NcCK）和PRJNA406998（NcT1、NcT2）。

1.5 測序數據質控及處理組東方蜜蜂微孢子蟲數據的篩濾

通過過濾原始數據中的低質量reads得到高質量reads，進而過濾掉不含3′接頭序列的reads，3′和5′接頭序列，短于18或長于30個核苷酸的reads及含ployA的reads，得到最終small RNA（sRNA）的tag序列。

按照以下流程篩濾處理組東方蜜蜂微孢子蟲（NcT1和NcT2）的數據：（1）利用Bowtie軟件將中腸測序得到的混合數據映射西方蜜蜂參考基因組（Amel4.5），去除比對上的數據；（2）繼而將未比對上的clean tags和NcCK組的tags分別比對到GeneBank和Pfam數據庫，過濾掉比對到核糖體RNA（rRNA）、細胞質小RNA（scRNA）、核仁小RNA（snoRNA）、核內小RNA（snRNA）和轉運RNA（tRNA）的tags；（3）進而將未比對上的clean tags映射到東方蜜蜂微孢子蟲參考基因組（Assembly ASM98816v1），比對上的數據即為NcT1和NcT2的數據。類似地，按照上述流程獲得比對上東方蜜蜂微孢子蟲參考基因組的NcCK數據。

1.6 MiRNA的預測與分析

首先將東方蜜蜂微孢子蟲的sRNA tags比對到miRBase中已知的miRNA，鑒定到保守miRNA。然后利用Mireap_v0.2軟件預測上述未比對上的sRNA tags的發夾結構，并結合其在參考基因組的位置鑒定新 miRNA。統計miRNA的長度和首位堿基偏好性。已知miRNA和新miRNA在各樣品中的表達量根據TPM（transcripts per million）算法計算。

1.7 DEmiRNA及其靶mRNA分析

根據各組miRNA的表達量統計NcCKNcT1、NcCKNcT2和NcT1NcT2比較組中的miRNA差異表達變化倍數（fold change，FC），按照 |log2fold change|≥1且≤0.05的標準篩選出顯著DEmiRNA。利用Graph Prism 7軟件統計DEmiRNA的個數并繪制柱狀圖。通過OmicShare在線平臺（www.omicshare. com）的相關工具進行Venn分析，采用默認參數。聯用RNAhybrid（v2.1.2）+ svm_light（v6.01）、Miranda（v3.3a）和TargetScan（Version: 7.0）軟件對東方蜜蜂微孢子蟲mRNA的3′-UTR進行靶向預測，取預測結果的交集作為可靠的預測結果。RNAhybrid軟件設置的參數為-f 2,8 -b 1 -v 3 -u 3 -e -10 -n 24，之后選取比對區域5′端 1—9 nt為seed region，用程序篩選過濾掉兩個連續GU matches的情況，再用svm_light和已經建立好的預測模型進行分類；Miranda（v3.3a）軟件設置的參數為-sc 140 -en -10 -strict -go -4.0 -ge -9.0；TargetScan軟件設置的參數：選取small RNA 5′端2—8 nt作為種子序列與轉錄本的3′-UTR區做預測。通過OmicShare在線平臺相關工具對DEmiRNA的靶mRNA進行GO（Gene Ontology）和KEGG（Kyoto Encyclopedia of Genes and Genomes）數據庫注釋，采用默認參數。利用Cytoscape軟件對DEmiRNA及其靶向結合的MAPK信號通路和糖酵解/糖異生通路相關mRNA的調控網絡進行可視化。

1.8 毒力因子和侵染相關因子相關DEmiRNA及靶mRNA調控網絡的構建及分析

東方蜜蜂微孢子蟲是專性寄生成年蜜蜂中腸上皮細胞的真菌病原，其增殖所需的物質和能量高度依賴宿主細胞提供。參考前人關于東方蜜蜂微孢子蟲的基因組學、轉錄組學和分子生物學的研究報道[22,35-38]，以及兔腦炎微孢子蟲（）[39]、斑馬魚微孢子蟲（）[40]和家蠶微孢子蟲[41]等其他微孢子蟲的相關研究報道，篩選出細胞凋亡抑制因子、孢壁蛋白、蓖麻毒素和極管蛋白等毒力因子，以及己糖激酶、ATP/ADP移位酶、ABC轉運蛋白和轉錄因子ste12等侵染因子。根據上述毒力因子和侵染因子相關mRNA及靶向結合關系預測出相應的DEmiRNA，利用Cytoscape軟件對DEmiRNA及靶mRNA的調控網絡進行可視化。

1.9 DEmiRNA的Stem-loop RT-qPCR驗證

隨機在NcCKNcT1、NcCKNcT2和NcT1NcT2中分別選取5、3和5個miRNA，根據成熟miRNA的序列，利用Premier 6軟件設計Stem-loop引物、特異性上游引物和通用下游引物。委托上海生工生物工程股份有限公司合成引物（表1）。利用RNA提取試劑盒（TaKaRa公司，日本）抽提各樣品的總RNA，利用Stem-loop引物和cDNA第一鏈合成試劑盒合成miRNA的cDNA，作為模板進行qPCR驗證。qPCR體系（20 μL）包括：SYBR Green Dye 10 μL，上下游引物（4 μmol·L-1）各1 μL，cDNA模板1 μL，Rox 0.44 μL，DEPC水補至20 μL。選擇東方蜜蜂微孢子蟲的5s rRNA作為內參。反應體系在ABI QuantStudio 3熒光定量PCR系統（ABI公司，美國）中運行。反應程序：95℃ 1 min；95℃ 15 s，49℃ 30 s，共40個循環。每個樣品進行3次重復。MiRNA的相對表達量根據2-ΔΔct算法計算得出。利用Graph Prism 7進行數據計算并繪圖。

表1 本研究使用的引物

2 結果

2.1 sRNA-seq的數據質控與過濾

東方蜜蜂微孢子蟲孢子（NcCK）測序得到 44 055 635條raw reads，經質控和過濾得到35 566 428條clean reads；經篩濾得到NcT1和NcT2的raw reads分別為40 126 400和42 145 184條，經質控和過濾得到的clean reads分別為21 525 444和30 316 756條。NcCK、NcT1和NcT2比對上參考基因組的平均clean reads數分別為7 188 407、2 103 266和4 926 712，平均比對率分別達到60.91%、29.26%和48.67%（表2）。

表2 sRNA-seq 數據總覽

2.2 東方蜜蜂微孢子蟲的miRNA預測與結構特征

基于質控后的高質量測序數據，共預測出545個miRNA，包括510個保守miRNA和35個新miRNA。結構特征分析結果顯示，上述miRNA的長度分布介于18—26 nt，分布在22 nt的miRNA數量最多（28.25%），其次為分布在18 nt和23 nt的miRNA（20.58%和14.77%）（圖1-A）；此外，miRNA的首位堿基偏向為U，其次為A和C（圖1-B）。

A：MiRNA長度分布Length distribution of miRNAs；B：MiRNA的首位堿基First nucleotide bias of miRNA

2.3 東方蜜蜂微孢子蟲的miRNA表達譜

表達量聚類分析結果顯示，miRNA的表達量在東方蜜蜂微孢子蟲孢子和侵染過程的東方蜜蜂微孢子蟲中差異明顯，一些miRNA（miR-44-x、miR-8797-x和miR-9277-y等）在孢子中的表達量較低，但隨著病原增殖表達量逐漸提高；另一些miRNA（miR-451-x、miR-486-x和miR-122-x等）在孢子中的表達量較高，但在病原的侵染過程呈下降趨勢；部分miRNA（miR-8-y、miR-12-z和miR-281-y等）在孢子中保持較低的表達水平，在侵染7 d的意蜂工蜂中腸中表達量較高，但在侵染10 d的宿主中腸中又降至較低的表達水平（圖2-A）。Venn分析結果顯示，有147個miRNA在3組樣品中共同表達，分別有86、139和41個僅在NcCK、NcT1和NcT2中特異性表達（圖2-B）。

A：各個樣品中miRNA表達量的熱圖分析Heat map analysis of miRNA expression in each sample；B：各個樣品中miRNA的Venn分析Venn analysis of miRNAs in each sample

上述共有miRNA的靶mRNA可注釋到27個功能條目，包括代謝過程（340）、催化活性（295）和細胞進程（292）；以及84條通路如代謝途徑（122）、核糖體（55）和次級代謝產物生物合成（46）。NcCK、NcT1和NcT2特有miRNA的靶mRNA可分別注釋到27、27和24個功能條目，注釋數目最多的條目均為代謝進程（333、329和242）、催化活性（294、291和212）和細胞進程（291、282和208）；這些靶mRNA還可注釋到84、84和83條通路，注釋數目最多的條目均為代謝途徑（121、119和87）、次級代謝產物生物合成（49、50和35）和核糖體（47、53和33）。括號內的數字代表注釋到該條目（或通路）的靶mRNA數。

2.4 東方蜜蜂微孢子蟲的miRNA差異表達譜

差異表達分析結果顯示，NcCKNcT1、NcCKNcT2和NcT1NcT2比較組中分別含有164、122和60個DEmiRNA，其中上調miRNA數分別為90、55和16個，下調miRNA數分別為74、67和44個。Venn分析結果顯示，上述3個比較組中共同上調表達的miRNA有5個，分別為miR-5106-y（log2FC= 15.76、19.69和3.93）、miR-301-x（log2FC=12.80、15.44和2.64）、novel-m0017-5p（log2FC=1.09、3.28和2.18）、miR-5108-y（log2FC =17.59、19.55和1.96）和miR-7465-x（log2FC=14.14、15.36和1.22）；共同下調表達的miRNA有6個，分別為miR-378-y（log2FC =-1.35、-2.61和-1.26）、miR-128-y（log2FC=-1.87、-3.87和-2.00）、miR-183-x（log2FC =-3.12、-16.80和-13.67）、miR-16-y（log2FC =-3.05、-16.85和-13.81）、miR-27-x（log2FC=-1.90、-16.05和-14.16）和miR-1307-y（log2FC =-1.46、-16.36和-14.90）。

2.5 東方蜜蜂微孢子蟲DEmiRNA的靶mRNA預測及功能和代謝通路注釋

NcCKNcT1、NcCKNcT2和NcT1NcT2比較組中DEmiRNA分別預測出1 885、1 733和1 524個靶mRNA。GO數據庫注釋結果顯示，上述靶mRNA可分別注釋到27、25和26個功能條目，均涉及生物學調控、細胞組分和分子功能3類（圖3）。3個比較組中DEmiRNA的靶mRNA注釋數量最多的均為代謝進程、催化活性、細胞進程、結合和細胞。但每個比較組在這些功能條目上富集的靶mRNA數并不相同，NcCKNcT1組的功能條目上富集的靶mRNA最多，NcCKNcT2組次之，NcT1NcT2組的功能條目上富集的靶mRNA數最少。3個比較組中相差的功能條目為病毒、病毒組分和分子功能調節，但富集的靶mRNA數均僅為1個。NcCKNcT1組包含這3個功能條目，NcCKNcT2組僅有分子功能調節，NcT1NcT2組包含病毒和病毒組分2個功能條目。

A：NcCK vs NcT1靶向的mRNA Target mRNA in NcCK vs NcT1；B：NcCK vs NcT2靶向的mRNA Target mRNA in NcCK vs NcT2；C：NcT1 vs NcT2靶向的mRNA Target mRNA in NcT1 vs NcT2。1：代謝進程Metabolic process；2：細胞進程Cellular process；3：單一有機體進程Single-organism process；4：定位Localization；5：生物學調控Biological regulation；6：生物進程調控Regulation of biological process；7：細胞成分組織或生物合成Cellular component organization or biogenesis；8：應激反應Response to stimulus；9：信號Signaling；10：復制Reproduction；11：生長Growth；12：細胞Cell；13：細胞組分Cell part；14：細胞器Organelle；15：大分子復合物Macromolecular complex；16：細胞膜Membrane；17：細胞膜組分Membrane part；18：細胞器組分Organelle part；19：膜封閉腔Membrane-enclosed lumen；20：病毒Virion；21：病毒組分Virion part；22：催化活性Catalytic activity；23：結合Binding；24：運輸活性Transporter activity；25：分子結構活性Structural molecule activity；26：核苷酸結合轉錄因子活性Nucleic acid binding transcription factor activity；27：分子功能調節Molecular function regulator

KEGG數據庫注釋結果顯示，上述NcCKNcT1、NcCKNcT2和NcT1NcT2比較組中DEmiRNA的靶mRNA可分別注釋到84、84和84條通路。NcCKNcT1中DEmiRNA的靶mRNA注釋到最多的通路分別為代謝途徑、核糖體、次級代謝產物生物合成、真核生物核糖體生物合成以及抗生素的生物合成；NcCKNcT2中DEmiRNA的靶mRNA注釋到最多的通路分別為代謝途徑、核糖體、次級代謝產物生物合成、真核生物核糖體生物合成以及抗生素的生物合成；NcT1NcT2中DEmiRNA的靶mRNA注釋到最多的通路分別為代謝途徑、核糖體、次級代謝產物生物合成、抗生素的生物合成以及內質網蛋白質加工（表3）。

表3 各比較組中DEmiRNA的靶mRNA富集的前15條通路

2.6 東方蜜蜂微孢子蟲DEmiRNA的調控網絡構建及分析

因各比較組中DEmiRNA及其靶mRNA數量較多，形成的調控網絡十分復雜，故根據靶mRNA的KEGG數據庫注釋信息，選擇與東方蜜蜂微孢子蟲的增殖、環境應答和能量代謝密切相關的MAPK信號通路和糖酵解/糖異生通路富集靶mRNA及相應DEmiRNA構建調控網絡，分析結果顯示NcCKNcT1中35個DEmiRNA靶向結合10個MAPK信號通路相關靶mRNA，49個DEmiRNA靶向結合13個糖酵解/糖異生通路相關靶mRNA（圖4-A、4-B）；NcCKNcT2中26個DEmiRNA靶向結合8個MAPK信號通路相關靶mRNA，40個DEmiRNA靶向結合13個糖酵解/糖異生通路相關靶mRNA（圖4-C、4-D）；NcT1NcT2中12個DEmiRNA靶向結合9個MAPK信號通路相關靶mRNA，17個DEmiRNA靶向結合10個糖酵解/糖異生通路相關靶mRNA（圖4-E、4-F）。

進一步篩選出蓖麻毒素B凝集素、孢壁蛋白、孢壁蛋白前體、孢壁和錨定吸盤復合蛋白和極管蛋白等毒力因子以及ATP/ADP移位酶、ABC轉運蛋白、己糖激酶、丙酮酸激酶和6-磷酸果糖激酶等侵染相關因子相關mRNA存在靶向結合關系的DEmiRNA，并構建DEmiRNA-靶mRNA的調控網絡。3個比較組中，共有15個DEmiRNA靶向結合糖酵解途徑3個關鍵酶基因，其中有10個DEmiRNA靶向己糖激酶編碼基因（圖5-A）；2個DEmiRNA可靶向結合抑制細胞凋亡蛋白編碼基因（圖5-B）；78個DEmiRNA可靶向結合4個ATP/ADP移位酶編碼基因和7個ABC轉運蛋白編碼基因（圖5-C）；61個DEmiRNA可靶向結合3個極管蛋白編碼基因、5個孢壁蛋白編碼基因、孢壁蛋白前體編碼基因、孢壁和錨定吸盤復合蛋白編碼基因（圖5-D）；27個DEmiRNA可靶向結合7個蓖麻毒素B凝集素的蛋白編碼基因（圖5-E）。

圖4 MAPK信號通路和糖酵解/糖異生通路相關靶mRNA及相應DEmiRNA的調控網絡

圖5 東方蜜蜂微孢子蟲毒力因子和侵染相關因子相關mRNA及相應DEmiRNA的調控網絡

2.7 東方蜜蜂微孢子蟲DEmiRNA的RT-qPCR驗證

從NcCKNcT1、NcCKNcT2和NcT1NcT2比較組中分別隨機選取5、3和5個顯著性DEmiRNA進行RT-qPCR驗證，結果顯示其表達趨勢變化與測序數據中的變化趨勢一致，證實測序數據和miRNA差異表達趨勢的可靠性（圖6）。

圖6 DEmiRNA的RT-qPCR驗證

3 討論

3.1 東方蜜蜂微孢子蟲miRNA的數量和結構特征

HUANG等[35]利用二代測序技術對東方蜜蜂微孢子蟲感染1—6 d的西方蜜蜂中腸組織進行轉錄組測序，分析發現東方蜜蜂微孢子蟲的1 122個差異表達基因可歸入4種表達模式。隨后，該團隊研究報道東方蜜蜂微孢子蟲的5個miRNA不僅能靶向結合病原自身的mRNA，還可以靶向結合宿主的mRNA[20]，暗示miRNA可能介導西方蜜蜂對東方蜜蜂微孢子蟲的跨界調控。本研究通過對東方蜜蜂微孢子蟲感染7 d和10 d的意蜂工蜂中腸以及純化孢子進行測序，利用生物信息學軟件預測出545個miRNA，其長度主要分布在21—24 nt，首位堿基主要偏向U，上述結構特征類似于其他已報道的真菌miRNA[17]。本研究鑒定到的miRNA為東方蜜蜂微孢子蟲的miRNA信息提供了很好的補充。筆者團隊前期利用高通量測序技術和生物信息學手段從東方蜜蜂微孢子蟲孢子中分別鑒定出83個長鏈非編碼RNA（long non-coding RNA，lncRNA）和204個環狀RNA（circular RNA，circRNA）[24-25]。本研究鑒定出的miRNA和前期鑒定出的lncRNA、circRNA豐富了東方蜜蜂微孢子蟲的ncRNA信息，也為其他微孢子蟲的ncRNA研究提供了重要參考。

3.2 侵染意蜂工蜂的東方蜜蜂微孢子蟲的miRNA表達譜差異

本研究發現，共有147個miRNA在NcCK、NcT1和NcT2中共同表達，推測這些miRNA在東方蜜蜂微孢子蟲孢子和侵染狀態均具有一定的生物學功能。miRNA的表達具有時空特異性。東方蜜蜂微孢子蟲的生活史約為6 d[35]，感染意蜂后7 d和10 d的東方蜜蜂微孢子蟲理論上處于第二生活史時期，其中感染7 d的病原孢子可能處于彈射極絲侵染宿主中腸上皮細胞階段，而感染10 d的病原可能處于裂殖體階段。在本研究中分別有86、139和41個miRNA在NcCK、NcT1和NcT2中特異性表達，推測這些miRNA在東方蜜蜂微孢子蟲孢子感染工蜂的不同時期發揮特定的調控作用。從NcCKNcT1、NcCKNcT2和NcT1NcT2比較組中分別篩選出164、122和60個DEmiRNA，表明東方蜜蜂微孢子蟲的侵染和增殖過程伴隨著miRNA表達量的動態變化，上述DEmiRNA與病原的侵染和增殖具有密切的相關性。此外，miR-5106-y、miR-301-x、novel-m0017-5p、miR-5108-y和miR-7465-x 5個miRNA在各比較組中共同上調表達，而miR-378-y、miR-128-y、miR-183-x、miR-16-y、miR-27-x和miR-1307-y 6個miRNA共同下調表達，暗示這些DEmiRNA在東方蜜蜂微孢子蟲增殖和侵染過程具有特殊功能，值得進一步的功能研究。

3.3 東方蜜蜂微孢子蟲DEmiRNA通過調控MAPK信號通路調節病原的環境應答和增殖

MAPK、cAMP和雙組分信號傳導系統是真菌應答外界環境的3個關鍵信號通路，均涉及調控真菌生長發育和毒力因子[42]。真菌的MAPK級聯是由外界刺激觸發，調節細胞周期、繁殖、形態發生、應激反應和毒性等一系列過程[43]。MAPK信號通路參與調節釀酒酵母（）的細胞壁完整性、配合、菌絲生長、孢子形成以及細胞應激反應[43]。筆者團隊前期研究發現，球囊菌可能通過激活MAPK信號通路促進對意蜂幼蟲的侵染[44]；但對于侵染中蜂幼蟲的球囊菌，該信號通路受到明顯抑制[45]；此外還發現在侵染意蜂工蜂的過程中，東方蜜蜂微孢子蟲的MAPK信號通路富集基因多數上調表達[23]。本研究中，NcCKNcT1、NcCKNcT2和NcT1NcT2比較組中DEmiRNA的靶mRNA分別有10、8和9個注釋到MAPK信號通路，表明相應的DEmiRNA可能調控MAPK信號通路相關基因的表達水平，從而調節東方蜜蜂微孢子蟲對外界環境的應答以及病原的增殖；進一步分析發現，上述DEmiRNA包含的下調miRNA（20、16和7個）多于上調DEmiRNA（15、10和5個），體現出意蜂工蜂與東方蜜蜂微孢子蟲互作的復雜性，宿主可能通過跨界調控等機制對病原的MAPK信號通路進行調節。轉錄因子ste12基因（）是MAPK信號通路下游的關鍵基因，有研究報道灰霉菌（）、白色念球菌（）、新型隱球菌（）、粗糙脈胞菌和構巢曲霉（）等真菌如缺少將導致生殖周期產生缺陷[46]。前期研究發現，東方蜜蜂微孢子蟲的在NcCKNcT1和NcCKNcT2比較組中均上調表達[23]。本研究中，靶向mRNA的miR-223-x（log2FC=-17.58和-17.58）、miR-148-y（log2FC=-4.54和-5.29）、miR-224-x（log2FC=-1.98和-3.53）和Novel-m0028-5p（log2FC=-14.19和-1.95）在NcCKNcT1、NcCKNcT2中表達量均為下調，表明東方蜜蜂微孢子蟲在侵染過程中通過下調上述miRNA的表達水平減弱對mRNA的抑制，從而上調的表達量，以促進自身的增殖和侵染。

3.4 東方蜜蜂微孢子蟲DEmiRNA通過調控糖酵解/糖異生、ATP/ADP移位酶和ABC轉運蛋白滿足增殖過程的能量需求

東方蜜蜂微孢子蟲感染會對西方蜜蜂工蜂造成能量脅迫，引起后者較未感染工蜂攝取更多的糖水[6]。有研究表明受東方蜜蜂微孢子蟲感染的西方蜜蜂的糖酵解相關基因表達量并無明顯變化，但東方蜜蜂微孢子蟲的催化糖酵解的核心酶基因在感染后24 h即開始表達，己糖激酶基因在感染后48 h也開始表達[35]。DOLGIKH等[47]研究發現東方蜜蜂微孢子蟲的己糖激酶基因含有信號肽，合成并釋放的己糖激酶在宿主的細胞核中積累。筆者團隊前期發現，NcCKNcT1和NcCKNcT2比較組中分別有11和12個差異表達基因富集在糖酵解/糖異生通路，其中包括己糖激酶基因在內的3個關鍵酶基因在東方蜜蜂微孢子蟲侵染意蜂工蜂的過程中均上調表達[23]。本研究中，NcCKNcT1和NcCKNcT2比較組中DEmiRNA的靶mRNA分別有13和13個注釋到糖酵解/糖異生通路；共有10個DEmiRNA可靶向結合己糖激酶基因，其中miR-139-x（log2FC= -15.66和-15.66）、miR-1-z（log2FC=-15.16和-15.16）、miR-223-x（log2FC=-17.58和-17.58）、miR-23-y（log2FC=-2.67和-2.67）在NcCKNcT1和NcCKNcT2中均為下調表達，推測東方蜜蜂微孢子蟲在侵染過程中通過下調上述DEmiRNA的表達量減弱對己糖激酶基因的抑制，從而增強己糖激酶的合成，以調節宿主的能量代謝并滿足自身增殖的能量需求。

微孢子蟲生長與繁殖所需的能量高度依賴宿主細胞提供。ATP/ADP移位酶和ABC轉運蛋白參與微孢子蟲對宿主物質和能量的竊取[37]。Paldi等研究發現通過RNAi敲減東方蜜蜂微孢子蟲的ATP/ADP移位酶基因可導致病原的增殖水平下降[36]。前期研究發現3個ATP/ADP移位酶基因（XM_002996497.1、XM_002995682.1和XM_002996655.1）在NcCKNcT1和NcCK vs NcT2中表達量上調[23]；本研究中，NcCKNcT1、NcCKNcT2和NcT1NcT2中分別有4、5和2個DEmiRNA與上述3個基因的mRNA存在靶向結合關系且差異變化趨勢相反，表明東方蜜蜂微孢子蟲侵染過程中相應的DEmiRNA參與調節ATP/ADP移位酶基因的表達，通過加強ATP/ADP移位酶的合成滿足自身能量需求并促進增殖。此外，另有1個ATP/ADP移位酶基因（XM_002996538.1）在NcCKNcT1和NcCKNcT2中下調表達[23]；但與其mRNA靶向結合的miR-1332-y（log2FC=15.51和16.19）、miR-216-x（log2FC=16.33和15.48）、miR-283-x（log2FC=5.50和4.62）、miR-301-y（log2FC=15.58和14.83）、miR-8-y（log2FC=7.03和5.53）5個DEmiRNA的差異變化趨勢與該基因的表達趨勢恰好相反，體現出病原增殖過程能量代謝的活躍性和相關DEmiRNA參與調控能量代謝的廣泛性。前期還發現東方蜜蜂微孢子蟲的7個ABC轉運蛋白編碼基因在NcCKNcT1和NcCKNcT2中差異表達，其中2個ABC轉運蛋白編碼基因（XM_002995835.1和XM_002996245.1）上調表達，5個ABC轉運蛋白基因（XM_002996253.1、XM_002995069.1、XM_ 002996674.1、TCONS_00003043和XM_002996675.1）下調表達[23]，說明東方蜜蜂微孢子蟲與意蜂工蜂之間存在復雜的互作，后者可能通過調節病原的部分ABC轉運蛋白基因的表達水平限制其增殖。本研究中，NcCKNcT1和NcCKNcT2比較組中分別有6和12個DEmiRNA（包含3和7個下調miRNA）靶向結合XM_002995835.1和XM_ 002996245.1。NcCKNcT1和NcCKNcT2比較組中分別有11、17、12和12個DEmiRNA（包含3、5、6和2個上調miRNA）靶向結合XM_002996253.1、XM_002995069.1、XM_002996674.1、TCONS_00003043和XM_002996675.1。上述結果表明東方蜜蜂微孢子蟲通過差異表達相應的miRNA調控ABC轉運蛋白編碼基因的表達。

3.5 東方蜜蜂微孢子蟲DEmiRNA通過調控極管蛋白和孢壁蛋白基因的表達調節病原的侵染過程

東方蜜蜂微孢子蟲的孢壁蛋白能夠保護孢子抵御外界不良環境，同時在病原應答外界環境變化和附著宿主細胞等方面發揮關鍵作用[12]。東方蜜蜂微孢子蟲孢子被蜜蜂攝入后，在中腸的堿性環境和離子濃度刺激下，內部極管快速彈射并刺入上皮細胞，然后經中空管道將有感染性的孢原質注入細胞[7-9]。在細胞識別和吸附以及極絲彈射的過程中，孢壁蛋白和極管蛋白發揮著非常重要的作用[48]。有研究表明家蠶微孢子蟲的孢壁蛋白5與極管蛋白2、極管蛋白3存在相互作用，孢壁蛋白5對于極管蛋白錨定在孢壁上具有重要作用[48]；孢壁蛋白9定位于極管上且能與極管蛋白1和極管蛋白2互作[49]。前期研究發現，東方蜜蜂微孢子蟲的極管蛋白1、極管蛋白2和極管蛋白基因在NcCK、NcT1和NcT2中的表達量持續上調，印證了極管蛋白對東方蜜蜂微孢子蟲侵染的重要性[23]。RODRíGUEZ-GARCíA等[38]研究報道干擾極管蛋白3基因的表達可顯著降低東方蜜蜂微孢子蟲孢子的數量。本研究發現，靶向結合極管蛋白1基因的miR-101-y（log2FC=-0.45和-2.75）、miR-451-x（log2FC=-7.03和-8.14）和miR-8117-y（log2FC=-4.63和-3.07）在NcCKNcT1和NcCKNcT2中下調表達；此外靶向結合極管蛋白2基因的novel-m0004-5p（log2FC=-3.81和-1.41）和novel-m0021-5p（log2FC=-4.41和-1.52）在NcCKNcT1和NcCKNcT2中也下調表達；還發現多達20個DEmiRNA靶向結合極管蛋白，其中有13個miRNA下調表達。以上結果表明東方蜜蜂微孢子蟲在侵染過程中通過下調相關miRNA的表達量提高極管蛋白基因的表達水平，從而增強極管蛋白的合成，以促進病原對宿主細胞的侵染。前期研究還發現，孢壁蛋白8基因在NcCKNcT1中表達量下調，在NcCKNcT2中表達量未發生變化[23]；本研究中，靶向結合孢壁蛋白8蛋白編碼基因的miR-315-x在NcCKNcT1中上調表達（log2FC=12.80），在NcCKNcT2中沒有發生差異表達（log2FC=0），說明miR-315-x通過調控孢壁蛋白8蛋白基因的表達參與東方蜜蜂微孢子蟲對意蜂工蜂的侵染過程。此外，靶向結合孢壁蛋白9基因mRNA的miR-1788-y（log2FC=13.75和13.94）和miR-424-x（log2FC=13.19和13.86）在NcCKNcT1和NcCKNcT2比較組中上調表達，與該基因的差異變化趨勢相反，推測這是病原和宿主之間的復雜互作的結果。

3.6 東方蜜蜂微孢子蟲DEmiRNA通過調控病原蓖麻毒素B凝集素和宿主細胞凋亡抑制基因促進侵染

蓖麻毒素B凝集素作為蓖麻毒素的B鏈，主要通過與細胞表面受體結合起到凝集素的作用，從而介導蓖麻毒素A鏈進入細胞發揮功能[50]。LIU等[41]對家蠶微孢子蟲休眠狀態和發芽狀態的孢子進行轉錄組測序和分析，發現蓖麻毒素B凝集素在發芽孢子中的表達量下調；Liu等[51]發現在家蠶微孢子蟲侵染家蠶卵巢BmN細胞的研究過程中加入蓖麻毒素B凝集素，能夠增強病原對宿主細胞的黏附能力。前期研究發現，東方蜜蜂微孢子蟲的7個蓖麻毒素凝集素基因在NcCKNcT1和NcCKNcT2比較組中表達量上調[23]；本研究中，NcCKNcT1和NcCKNcT2比較組中的DEmiRNA可以靶向結合其中的6個蓖麻毒素凝集素基因的mRNA，其中有10個DEmiRNA在NcCKNcT1和NcCKNcT2比較組中的差異變化趨勢與該基因相反，說明東方蜜蜂微孢子蟲在侵染意蜂工蜂的過程中通過調節部分miRNA的表達水平對蓖麻毒素基因進行表達調控，從而增強對宿主細胞的黏附性和侵染，但宿主同樣能夠通過互作對病原的蓖麻毒素凝集素基因進行一定程度的抑制，再次體現出二者相互作用的復雜性；該基因或許能作為一個蜜蜂微孢子蟲病治療的有效靶點。東方蜜蜂微孢子蟲能抑制西方蜜蜂的細胞凋亡[5]。細胞凋亡抑制基因通過結合活化半胱天冬酶抑制細胞的凋亡[35]。HUANG等[35]研究發現東方蜜蜂微孢子蟲的細胞凋亡抑制基因在病原侵染西方蜜蜂后24 h即開始表達。前期研究發現，細胞凋亡抑制基因（XM_002995997.1）在NcCKNcT1和NcCKNcT2中表達量上調；本研究中，NcCKNcT1和NcCKNcT2比較組中有2個DEmiRNA可以靶向結合細胞凋亡抑制基因，其中novel-m0017-3p（log2FC =-14.05和-14.05）在NcCKNcT1和NcCKNcT2比較組中的差異變化趨勢與該基因相反；而miR-186-x（log2FC=1.02和0.23）在NcCKNcT1和NcCKNcT2比較組中的差異變化趨勢與該基因變化趨勢一致，但差異變化倍數較小。以上結果表明東方蜜蜂微孢子蟲主要通過下調novel-m0017-3p增強對宿主細胞凋亡抑制基因的抑制作用，從而延長宿主細胞的存活時間，更好地促進自身的繁殖。

4 結論

在miRNA組學層面對侵染意蜂工蜂的東方蜜蜂微孢子蟲的DEmiRNA、調控網絡及潛在作用進行了深入解析，揭示了DEmiRNA在東方蜜蜂微孢子蟲侵染過程中廣泛參與調控病原的物質和能量代謝、細胞生命活動，部分DEmiRNA通過調節ste12、己糖激酶、孢壁蛋白、極管蛋白、ATP/ADP移位酶、ABC轉運蛋白、蓖麻毒素等毒力及侵染相關因子基因的表達促進東方蜜蜂微孢子蟲在宿主細胞內的增殖和侵染。

[1] CUOMO C A, DESJARDINS C A, BAKOWSKI M A, GOLDBERG J, MA A T, BECNEL J J, DIDIER E S, FAN L, HEIMAN D I, LEVIN J Z, YOUNG S, ZENG Q, TROEMEL E R. Microsporidian genome analysis reveals evolutionary strategies for obligate intracellular growth., 2012, 22(12): 2478-2488.

[2] NAGELI K W. Uber die neue krankheit der seidenraupe und verwandte organismen., 1857, 15: 760-761.

[3] 吳杰, 項勛, 趙屹欽, 鄧紫艷, 王世玉, 常華. 微孢子蟲研究進展. 動物醫學進展, 2017, 38(6): 78-81.

WU J, XIANG X, ZHAO Y Q, DENG Z Y, WANG S Y, CHANG H. Progress on microsporidia., 2017, 38(6): 78-81. (in Chinese)

[4] FORSGREN E, FRIES I. Comparative virulence ofandin individual European honey bees., 2010, 170(3/4): 212-217.

[5] MARTíN-HERNáNDEZ R, BARTOLOMé C, CHEJANOVSKY N, LE CONTE Y, DALMON A, DUSSAUBAT C, GARCíA-PALENCIA P, MEANA A, PINTO M A, SOROKER V, HIGES M.in: a 12 years post-detection perspective., 2018, 20(4): 1302-1329.

[6] BOTíAS C, MARTíN-HERNáNDEZ R, BARRIOS L, MEANA A, HIGES M.spp. infection and its negative effects on honey bees () at the colony level., 2013, 44(1): 25.

[7] FRIES I, FENG F, DA SILVA A, SLEMENDA S B, PIENIAZEK N J.n. sp. (Microspora, Nosematidae), morphological and molecular characterization of a microsporidian parasite of the Asian honey bee(Hymenoptera, Apidae)., 1996, 32(3): 356-365.

[8] HIGES M, MARTíN R, MEANA A., a new microsporidian parasite in honeybees in Europe., 2006, 92(2): 93-95.

[9] HUANG W F, JIANG J H, CHEN Y W, WANG C H. Aisolate from the honeybee., 2007, 38: 30-37.

[10] WANG L L, XU X Y, YANG J, CHEN L H, LIU B, LIU T, JIN Q. Integrated microRNA and mRNA analysis in the pathogenic filamentous fungus., 2018, 19(1): Article number 933.

[11] KANG K, ZHONG J, JIANG L, LIU G, GOU C Y, WU Q, WANG Y, LUO J, GOU D. Identification of microRNA-like RNAs in thefilamentous fungusby solexa sequencing.

, 2013, 8(10): e76288.

[12] FRIES I, GRANADOS R, MORSE R. Intracellular germination of spores ofZ., 1992, 23: 61-71.

[13] LEE R C, FEINBAUM R L, AMBROS V. Theheterochronic gene lin-4 encodes small RNAs with antisense complementarity to lin-14., 1993, 75(5): 843-854.

[14] REINHART B J, SLACK F J, BASSON M, PASQUINELLI A E, BETTINGER J C, ROUGVIE A E, HORVITZ H R, RUVKUN G. The 21-nucleotide let-7 RNA regulates developmental timing in., 2000, 403(6772): 901-906.

[15] BARTEL D P. MicroRNAs: Target recognition and regulatory functions., 2009, 136(2): 215-233.

[16] GRIFFITHS-JONES S, SAINI H K, VAN DONGEN S, ENRIGHT A J. MiRBase: tools for microRNA genomics., 2008, 36(Database issue): D154-D158.

[17] LEE H C, LI L, GU W, XUE Z, CROSTHWAITE S K, PERTSEMLIDIS A, LEWIS Z A, FREITAG M, SELKER E U, MELLO C C, LIU Y. Diverse pathways generate microRNA-like RNAs and dicer-independent small interfering RNAs in fungi., 2010, 38(6): 803-814.

[18] ZHOU Q, WANG Z, ZHANG J, MENG H, HUANG B. Genome- wide identification and profiling of microRNA-like RNAs fromduring development., 2012, 116(11): 1156-1162.

[19] 郭睿, 王海朋, 陳華枝, 熊翠玲, 鄭燕珍, 付中民, 趙紅霞, 陳大福. 蜜蜂球囊菌的microRNA鑒定及其調控網絡分析. 微生物學報, 2018, 58(6): 1077-1089.

GUO R, WANG H P, CHEN H Z, XIONG C L, ZHENG Y Z, FU Z M, ZHAO H X, CHEN D F. Identification ofmicroRNAs and investigation of their regulation networks., 2018, 58(6): 1077-1089. (in Chinese)

[20] EVANS J D, HUANG Q. Interactions among host-parasite microRNAs duringproliferation in., 2018, 9: 698.

[21] HUANG Q. Evolution of Dicer and Argonaute orthologs in microsporidian parasites., 2018, 65: 329-332.

[22] HUANG Q, LI W, CHEN Y, RETSCHNIG-TANNER G, YANEZ O, NEUMANN P, EVANS J D. Dicer regulatesproliferation in honeybees., 2019, 28(1): 74-85.

[23] 耿四海, 周丁丁, 范小雪, 蔣海賓, 祝智威, 王杰, 范元嬋, 王心蕊, 熊翠玲, 鄭燕珍, 付中民, 陳大福, 郭睿. 轉錄組學分析揭示東方蜜蜂微孢子蟲侵染意大利蜜蜂的分子機制. 昆蟲學報, 2020, 63(3): 294-308.

GENG S H, ZHOU D D, FAN X J, JIANG H B, ZHU Z W, WANG J, FAN Y C, WANG X R, XIONG C L, ZHENG Y Z, FU Z M, CHEN D F, GUO R. transcriptomic analysis reveals the molecular mechanism underlyinginfection of., 2020, 63(3): 294-308. (in Chinese)

[24] GUO R, CHEN D F, CHEN H Z, XIONG C L, ZHENG Y Z, HOU C S, DU Y, GENG S H, WANG H P, ZHOU D D, GUO Y L. Genome-wide identification of circular RNAs in fungal parasite., 2018, 75(12): 1655-1660.

[25] 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 X, DHIRAJ K. First identification of long non-coding RNAs in fungal parasite, 2018, 49: 660-670.

[26] GUO R, CHEN H Z, DU Y, ZHOU D D, GENG S H, WANG H P, ZHU Z W, SHI C Y, WAN J Q, XIONG C L, ZHENG Y Z, CHEN D F. Genome-wide identification of long non-coding RNAs and their regulatory networks involved inresponse toinfection., 2019, 10(8): 245.

[27] 付中民, 陳華枝, 劉思亞, 祝智威, 范小雪, 范元嬋, 萬潔琦, 張璐, 熊翠玲, 徐國鈞, 陳大福, 郭睿. 意大利蜜蜂響應東方蜜蜂微孢子蟲脅迫的免疫應答. 中國農業科學, 2019, 52(17): 3069-3082.

FU Z M, CHEN H Z, LIU S Y, ZHU Z W, FAN X X, FAN Y C, WAN J Q, ZHANG L, XIONG C L, XU G J, CHEN D F, GUO R. Immune responses oftostress., 2019, 52(17): 3069-3082. (in Chinese)

[28] 周倪紅, 王海朋, 周丁丁, 付中民, 祝智威, 范元嬋, 張曦, 熊翠玲, 鄭燕珍, 陳大福, 郭睿. 意大利蜜蜂工蜂中腸響應東方蜜蜂微孢子蟲脅迫的可變剪接基因分析. 福建農林大學學報 (自然科學版), 2020, 49(3): 372-379.

ZHOU N H, WANG H P, ZHOU D D, FU Z M, ZHU Z W, FAN Y C, ZHANG X, XIONG C L, ZHENG Y Z, CHEN D F, GUO R. Analysis on the response of alternatively splicing genes inworkers’ midguts tostress., 2020, 49(3): 372-379. (in Chinese)

[29] STOLTZ D, SHEN X R, BOGGIS C, SISSON G. Molecular diagnosis of Kashmir bee virus infection., 1995, 34(3): 153-160.

[30] BENJEDDOU M, LEAT N, ALLSOPP M, DAVISON S. Detection of Acute bee paralysis virus and Black queen cell virus from honeybees by reverse transcriptase PCR., 2001, 67(5): 2384-2387.

[31] SINGH R, LEVITT A L, RAJOTTE E G, HOLMES E C, OSTIGUY N, VAN ENGSLEDORP D, LIPKIN W L, DE PAMPHILIS C W, TOTH A L, COX-FOSTER D L. RNA viruses in hymenopteran pollinators: evidence of inter-taxa virus transmission via pollen and potential impact on non-hymenopteran species., 2010, 5(12): e14357.

[32] RIBIèRE M, TRIBOULOT C, MATHIEU L, AURIèRES C, FAUCON J P, PéPIN M. Molecular diagnosis of chronic bee paralysis virus infection., 2002, 33(3): 339-351.

[33] CHEN Y, EVANS J D, SMITH I B, PETTIS J S.is a long-present and wide-spread microsporidian infection of the European honey bee () in the United States., 2008, 97(2): 186-188.

[34] 郭睿, 杜宇, 童新宇, 熊翠玲, 鄭燕珍, 徐國鈞, 王海朋, 耿四海, 周丁丁, 郭意龍, 吳素珍, 陳大福. 意大利蜜蜂幼蟲腸道在球囊菌侵染前期的差異表達microRNA及其調控網絡. 中國農業科學, 2019, 52(1): 166-180.

GUO R, DU Y, TONG X Y, XIONG C L, ZHENG Y Z, XU G J, WANG H P, GENG S H, ZHOU D D, GUO Y L, WU S Z, CHEN D F. Differentially expressed microRNAs and their regulation networks inlarval gut during the early stage ofinfection., 2019, 52(1): 166-180. (in Chinese)

[35] HUANG Q, CHEN Y P, WANG R W, CHENG S, EVANS J D. Host-parasite interactions and purifying selection in a microsporidian parasite of honey bees., 2016, 11(2): e0147549.

[36] 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 () colony declines., 2010, 76(17): 5960-5964.

[37] PELIN A, SELMAN M, ARIS-BROSOU S, FARINELLI L, CORRADI N. Genome analyses suggest the presence of polyploidy and recent human-driven expansions in eight global populations of the honeybee pathogen., 2015, 17(11): 4443-4458.

[38] RODRíGUEZ-GARCíA C, EVANS J D, LI W, BRANCHICCELA B, LI J H, HEERMAN M C, BANMEKE O, ZHAO Y, HAMILTON M, HIGES M, MARTíN-HERNáNDEZ R, CHEN Y P. Nosemosis control in European honey bees,, by silencing the gene encodingpolar tube protein 3., 2018, 221(19): jeb184606.

[39] KATINKA M D, DUPRAT S, CORNILLOT E, MéTéNIER G, THOMARAT F, PRENSIER G, BARBE V, PEYRETAILLADE E, BROTTIER P, WINCKER P, DELBAC F, ALAOUI H E, PEYRET P, SAURIN W, GOUY M, WEISSENBACH J, VIVARèS C P. Genome sequence and gene compaction of the eukaryote parasite., 2001, 414(6862): 450-453.

[40] NDIKUMANA S, PELIN A, WILLIOT A, SANDERS J L, KENT M, CORRADI N. Genome analysis of: a microsporidian parasite of Zebrafish ()., 2017, 64(1): 18-30.

[41] LIU H, LI M, HE X, CAI S, HE X, LU X. Transcriptome sequencing and characterization of ungerminated and germinated spores of., 2016, 48(3): 246-256.

[42] 侯彬彬, 劉霞, 張振穎, 劉曉明. 參與真菌形態及毒力形成的信號轉導通路的研究進展. 中國真菌學雜志, 2015, 10(4): 241-244.

HOU B B, LIU X, ZHANG Z Y, LIU X M. Advances in signal transduction pathways involved in fungal morphology and virulence., 2015, 10(4): 241-244. (in Chinese)

[43] LENGELER K B, DAVIDSON R C, D’SOUZA C, HARASHIMA T, SHEN W C, WANG P, PAN X, WAUGH M, HEITMAN J. Signal transduction cascades regulating fungal development and virulence., 2000, 64(4): 746-785.

[44] 陳大福, 郭睿, 熊翠玲, 梁勤, 鄭燕珍, 徐細建, 黃枳腱, 張瞾楠, 張璐, 李汶東, 童新宇, 席偉軍. 脅迫意大利蜜蜂幼蟲腸道的球囊菌的轉錄組分析.昆蟲學報, 2017, 60(4): 401-411.

CHEN D F, GUO R, XIONG C L, LIANG Q, ZHENG Y Z, XU X J, HUANG Z J, ZHANG Z N, ZHANG L, LI W D, TONG X Y, XI W J. Transcriptomic analysis ofstressing larval gut of(Hyemenoptera: Apidae)., 2017, 60(4): 401-411. (in Chinese)

[45] 郭睿, 陳大福, 黃枳腱, 梁勤, 熊翠玲, 徐細建, 鄭艷珍, 張瞾楠, 解彥玲, 童新宇, 侯志賢, 江亮亮, 刀晨. 球囊菌脅迫中華蜜蜂幼蟲腸道過程中病原的轉錄組學研究. 微生物學報, 2017, 57(12): 1865-1878.

GUO R, CHEN D F, HUANG Z J, LIANG Q, XIONG C L, XU X J, ZHENG Y Z, ZHANG Z N, XIE Y L, TONG X Y, HOU Z X, JIANG L L, DAO C. Transcriptome analysis ofstressing larval gut of., 2017, 57(12): 1865-1878. (in Chinese)

[46] RISPAIL N, DI PIETRO A. The homeodomain transcription factor Ste12: Connecting fungal MAPK signalling to plant pathogenicity., 2010, 3(4): 327-332.

[47] DOLGIKH V V, TSAREV A A, TIMOFEEV S A, ZHURAVLYOV V S. Heterologous overexpression of active hexokinases from microsporidiaandconfirms their ability to phosphorylate host glucose., 2019, 118(5): 1511-1518.

[48] SOUTHERN T R, JOLLY C E, LESTER M E, HAYMAN J R. Enp1, a microsporidian spore wall protein that enables spores to adhere to and infect host cells., 2007, 6(8): 1354-1362.

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

[50] ENDO Y, TSURUGI K. The RNA N-glycosidase activity of ricin A-chain. The characteristics of the enzymatic activity of ricin A-chain with ribosomes and with rRNA., 1988, 263(18): 8735-8739.

[51] LIU H, LI M, CAI S, HE X, SHAO Y, LU X. Ricin-B-lectin enhances microsporidiainfection in BmN cells from silkworm., 2016, 48(11): 1050-1057.

The mechanism underlying microRNAs-mediatedinfection toworker

Geng SiHai, Shi CaiYun, Fan XiaoXue, Wang Jie, Zhu ZhiWei, Jiang HaiBin, FAN YuanChan, Chen HuaZhi, Du Yu, Wang XinRui, XIONG CuiLing, ZHENG YanZhen, FU ZhongMin, CHEN DaFu, GUO Rui

(College of Animal Sciences (College of Bee Science), Fujian Agriculture and Forestry University, Fuzhou 350002)

【】Differentially expressed miRNAs (DEmiRNAs) and their target mRNAs in purified spores ofandinfectingworker were systematically analyzed, followed by screening, investigation and exploration of DEmiRNAs and corresponding regulatory networks associated with pathogenic virulence factors and infection factors to reveal the mechanism ofinfection toworker.【】workers’ midguts at 7 d and 10 d post infection (dpi) withand purified sporesof(NcCK) were deeply sequenced using small RNA-seq (sRNA-seq), followed by screening out datasets ofduring the infection process (NcT1 and NcT2) and purified spores through sequential mapping to rRNA database,genome andgenome. Based on the criteria of≤0.05, |log2fold change|≥1, DEmiRNAs in each comparison group were filtered out by comparative analysis. Expression profiling of DEmiRNAs, prediction of target mRNAs, function and pathway annotations of target mRNAs, and construction and survey of regulatory networks were performed by relevant bioinformatic software.The differential expression trend of DEmiRNAs and the reliability of sequencing data were verified by RT-qPCR.【】There were 164, 122 and 60 DEmiRNAs in NcCKNcT1, NcCKNcT2 and NcT1NcT2 comparison groups. Venn analysis showed that 5 up-regulated miRNAs and 6 down-regulated ones were shared by these three comparison groups. In total, 1 885, 1 733 and 1 524 target mRNAs of DEmiRNAs were respectively predicted. These targets were annotated to 27, 25, and 26 functional terms, respectively, with the most abundant annotations being metabolic processes, catalytic activity, cellular processes, binding and cell. These targets were as well respectively annotated to 84, 84, and 84 pathways, among them the largest groups were metabolic pathway, ribosome, and secondary metabolite biosynthesis. In addition, 35, 26 and 12 target mRNAs associated with MAPK signaling pathway were targeted by DEmiRNAs in NcCKNcT1, NcCKNcT2 and NcT1NcT2, while 49, 40 and 17 target mRNAs related to glycolysis/gluconeogenesis pathway were targeted by DEmiRNAs in the aforementioned three comparison groups. Further analysis revealed thatDEmiRNAs ofwere involved in regulating the expression of genes encoding virulence factors such as ricin B lectin, apoptosis inhibitor, polar tube protein and spore wall protein; and genes encoding infection factors such as hexokinase and ATP/ADP transferase, ABC transporter and transcription factor ste12.【】The differential expression profile of miRNAs ininfectingworker was parsed via in-depth and detailed investigation and exploration.may regulate the expression level of genes encoding virulence factors including ricin B lectin, apoptosis inhibitor, polar tube protein and spore wall protein, as well as genes encoding infection factors such as hexokinase, ATP/ADP transferase, ABC transporter and transcription factor ste12 corresponding miRNAs, by controlling the expression of corresponding miRNAs, thus adapting to the environment within the host cell and promoting its proliferation and infection.

;; microRNA (miRNA); infection mechanism; regulatory network

10.3864/j.issn.0578-1752.2020.15.018

2019-12-02；

2019-12-28

國家現代農業產業技術體系建設專項（CARS-44-KXJ7）、福建省自然科學基金（2018J05042）、福建省教育廳中青年教師教育科研項目（JAT170158）、福建農林大學杰出青年科研人才計劃（xjq201814）、福建農林大學碩士生導師團隊項目、福建農林大學科技創新專項基金（CXZX2017342，CXZX2017343）、福建省大學生創新創業訓練計劃（3165602032，3155006018）

耿四海，E-mail：15737313592@163.com。石彩云，E-mail：1901638174@qq.com。耿四海和石彩云為同等貢獻作者。通信作者郭睿，Tel： 0591-87640197；E-mail：ruiguo@fafu.edu.cn

（責任編輯 岳梅）