陳華枝，祝智威，蔣海賓，王杰，范元嬋，范小雪，萬潔琦，盧家軒，熊翠玲，鄭燕珍，付中民，陳大福，郭睿

蜜蜂球囊菌菌絲和孢子中微小RNA及其靶mRNA的比較分析

陳華枝，祝智威，蔣海賓，王杰，范元嬋，范小雪，萬潔琦，盧家軒，熊翠玲，鄭燕珍，付中民，陳大福，郭睿

（福建農(nóng)林大學動物科學學院（蜂學學院），福州 350002）

【】蜜蜂球囊菌(，球囊菌)專性侵染蜜蜂幼蟲而導致白堊病。本研究旨在通過small RNA-seq（sRNA-seq）技術和生物信息學方法對球囊菌純化菌絲（AaM）和純化孢子（AaS）進行深度測序和比較分析，明確球囊菌菌絲miRNA和孢子miRNA的數(shù)量、結構和表達譜差異，并揭示菌絲和孢子共有miRNA、特有miRNA和差異表達miRNA（differentially expressed miRNA，DEmiRNA）及其靶mRNA與球囊菌菌絲和孢子生長、發(fā)育和病原致病性的潛在關系。實驗室條件下獲得純培養(yǎng)的球囊菌，利用sRNA-seq技術對AaM和AaS分別進行測序，通過對原始讀段（raw reads）進行過濾和質控獲得有效標簽序列（clean tags）。通過Venn分析篩選菌絲和孢子共有miRNA和特有miRNA。根據(jù)≤0.05且|log2fold change|≥1的標準篩選AaM vs AaS的DEmiRNA。對上述共有miRNA、特有miRNA和DEmiRNA的靶mRNA進行預測，并對靶mRNA進行GO及KEGG數(shù)據(jù)庫注釋。根據(jù)靶向結合關系構建DEmiRNA和靶mRNA的調控網(wǎng)絡。利用RT-qPCR驗證測序數(shù)據(jù)的可靠性。AaM和AaS中分別得到12 982 320和12 708 832條raw reads，經(jīng)過濾和質控分別得到10 800 101和9 888 848條clean tags。AaM中miRNA的長度介于18—26 nt，AaS中miRNA的長度介于18—24 nt，分布miRNA數(shù)量最多的長度均為18 nt，AaM和AaS中首位堿基為U的miRNA數(shù)量最多。AaM和AaS中表達量最高的miRNA均為miR6478-x、miR10516-x和miR482-x。菌絲和孢子共有miRNA靶向結合5 946個mRNA，二者特有miRNA分別靶向結合6 141和6 346個mRNA。共有miRNA的靶mRNA主要參與代謝進程、細胞進程和催化活性等42個功能條目，以及翻譯、碳水化合物代謝和能量代謝等120條通路。AaM vs AaS比較組包含93個DEmiRNA，可靶向結合6 090個mRNA，這些靶mRNA可注釋到38個功能條目和120條通路。DEmiRNA與靶mRNA之間形成較為復雜的調控網(wǎng)絡，miR-4968-y位于調控網(wǎng)絡的中心且能夠靶向結合多達118個mRNA。RT-qPCR結果顯示5個DEmiRNA的表達趨勢與測序數(shù)據(jù)一致，證實了本研究中測序數(shù)據(jù)的可靠性。球囊菌菌絲和孢子中的miRNA具有類似的結構特征，但表達譜表現(xiàn)出明顯差異；菌絲和孢子可能通過特異性表達和差異表達部分miRNA對其生長、發(fā)育和生殖進行調控。

蜜蜂球囊菌；菌絲；孢子；微小RNA；信使RNA；靶向結合

0 引言

【研究意義】蜜蜂球囊菌（，球囊菌）是一種特異性侵染蜜蜂幼蟲的致死性真菌病原，能夠使蜜蜂幼蟲罹患白堊病，可導致成年蜜蜂數(shù)量和群勢的急劇下降[1]。近期的研究表明，球囊菌對中華蜜蜂（，中蜂）幼蟲和成蟲、成年熊蜂也具有侵染性[2]。食物連同球囊菌孢子被蜜蜂幼蟲攝入后進入中腸，幼蟲期中腸與后腸之間存在一層隔膜，至預蛹期隔膜消失，孢子隨食物殘渣涌入后腸，在O2的刺激下孢子劇烈萌發(fā)、菌絲大量生長，先穿透腸壁再從尾部穿透體壁，進而包裹整個蟲尸[3]。前人在球囊菌的分類鑒定[4-5]、增殖方式[6]和病理學[7-8]等方面開展了較多研究，但其分子生物學及組學研究較少。球囊菌微小RNA（microRNA，miRNA）等非編碼RNA的相關研究尤為滯后。利用small RNA-seq（sRNA-seq）技術和生物信息學方法對球囊菌菌絲和孢子分別進行深度測序和組學分析，明確二者miRNA的數(shù)量和結構特征差異，揭示miRNA與菌絲生長、孢子萌發(fā)、雜交產(chǎn)孢、毒力因子合成與分泌等生物學過程的關系，可為明晰相關分子機理提供參考信息和數(shù)據(jù)基礎。【前人研究進展】MiRNA是一類長度約為23 nt的細胞內(nèi)源性小RNA分子，在真核生物中廣泛存在且高度保守[9]，可通過切割靶mRNA或阻止蛋白質的翻譯過程來調控基因表達[10]，從而參與調控生長發(fā)育及信號轉導等相關進程[11]。自從在秀麗隱桿線蟲（）中發(fā)現(xiàn)第一個miRNA以來[12]，人們已在動物[13]、植物[14]和微生物中[15-16]陸續(xù)發(fā)現(xiàn)了大量miRNA，但真菌miRNA的相關研究總體仍較為滯后。前期研究中，筆者團隊已在mRNA組學水平對球囊菌侵染意大利蜜蜂（，意蜂）幼蟲和中蜂幼蟲過程的宿主與病原相互作用進行了系統(tǒng)探究[17-21]，組裝并注釋了球囊菌的參考轉錄組[21]；基于球囊菌菌絲和孢子混合樣品的高質量組學數(shù)據(jù)鑒定出379個長鏈非編碼RNA和551個環(huán)狀RNA，并通過分子生物學手段驗證了它們的真實表達[22-23]?！颈狙芯壳腥朦c】筆者團隊前期基于球囊菌菌絲和孢子混合樣品的sRNA組學數(shù)據(jù)鑒定出118個miRNA，并預測和分析了miRNA的靶mRNA及調控網(wǎng)絡[24]。然而，球囊菌菌絲和孢子miRNA的數(shù)量、結構特征、表達譜差異仍未明確，miRNA與菌絲和孢子生長、發(fā)育、生殖和致病性的關系還不清楚?！緮M解決的關鍵問題】明確球囊菌菌絲和孢子miRNA的數(shù)量、結構特征和表達譜差異，揭示共有miRNA、特有miRNA和差異表達miRNA（differentially expressed miRNA，DEmiRNA）與菌絲和孢子生長、發(fā)育、生殖及致病性的潛在關系。

1 材料與方法

試驗于2019年在福建農(nóng)林大學動物科學學院（蜂學學院）蜜蜂保護實驗室完成。

1.1 供試真菌

球囊菌菌株由福建農(nóng)林大學動物科學學院（蜂學學院）蜜蜂保護實驗室分離和保存[22-24]。

1.2 球囊菌培養(yǎng)及測序樣品制備

參照筆者團隊前期已建立的方法[19-20,22-23]，將實驗室保存的球囊菌孢子接種于10塊馬鈴薯葡萄糖瓊脂（potato dextrose agar，PDA）固體培養(yǎng)基上，置于33℃生化培養(yǎng)箱中培養(yǎng)，培養(yǎng)7 d的PDA固體培養(yǎng)基上層為蓬松的白色菌絲，覆蓋著下層的黑色孢子囊，為避免菌絲和孢子囊之間的交叉污染，首先在超凈臺中用干凈的接種環(huán)刮取最上層的白色菌絲，避免接觸與孢子囊接觸的菌絲，將刮取的菌絲（AaM）集中于一個RNA Free的EP管，經(jīng)液氮速凍后迅速移至-80℃超低溫冰箱保存?zhèn)溆茫辉诔瑑襞_中通過無菌操作將覆蓋在孢子囊上的薄層菌絲小心刮去，再用另一個干凈的接種環(huán)將孢子囊刮至一個RNA Free的EP管；繼而按照Jensen等[25]的方法進行差速離心，球囊菌孢子密度較大，沉淀在EP管底部，棄去上清；重復上述差速離心3次，棄去上清，保留孢子沉淀；取少量孢子制備懸液，取少量懸液進行顯微制片及觀察，視野中可見較多的游離孢子，未見菌絲，因而得到的孢子沉淀即為球囊菌的純凈孢子（AaS）；將上述EP管投入液氮速凍，然后迅速移至-80℃超低溫冰箱保存?zhèn)溆谩?/p>

1.3 cDNA文庫構建及深度測序

（1）用Trizol法分別提取AaM和AaS的總RNA，瓊脂糖凝膠電泳切膠選擇18—30 nt的片段，分別連接3′接頭和5′接頭；（2）對連接了兩側接頭的small RNA進行反轉錄和PCR擴增；（3）瓊脂糖凝膠電泳回收并純化約140 bp的條帶，完成文庫構建。委托廣州基迪奧生物科技有限公司對上述樣品進行單端測序，測序平臺為Illumina MiSeq。測序數(shù)據(jù)已上傳到NCBI SRA數(shù)據(jù)庫，BioProject號：PRJNA560456。

1.4 測序數(shù)據(jù)質控

按照前期已建立的方法[24,26-27]對下機的原始讀段（raw reads）進行質量控制：（1）過濾掉質量值低于20的堿基數(shù)超過1個的reads；（2）過濾除掉含有未知堿基（N）的reads；（3）過濾3′或5′接頭的reads，并去除長度＜18 bp的reads；（4）過濾包含poly A的reads。過濾后得到的clean reads用于后續(xù)分析。

1.5 菌絲和孢子miRNA的預測與Venn分析

將嚴格質控后得到的small RNA有效標簽序列（clean tags）分別比對到GenBank數(shù)據(jù)庫（http://www. ncbi.nlm.nih.gov/Web/Genbank/）、Rfam數(shù)據(jù)庫（http://rfam.xfam.org/）、核糖體數(shù)據(jù)庫（http://oberon. fvms.ugent.be:8080/rRNA/lsu/index.html）和球囊菌參考基因組（assembly AAP 1.0，www.ncbi.nlm.nih.gov/ genome/656?genome_assembly_id=274809），以去除可能的rRNA、scRNA、snoRNA、snRNA和tRNA，以及比對上的基因組外顯子區(qū)域、內(nèi)含子區(qū)域和重復序列區(qū)域的clean tags，將未比對上的clean tags與miRBase數(shù)據(jù)庫（www.mirbase.org/）中的miRNA前體序列進行比對，從而鑒定已知miRNA序列。利用Mireap_V 0.2軟件將剩余的未比對上的clean tags比對球囊菌參考基因組（assembly AAP 1.0），得到可能的前體序列，根據(jù)clean tags在前體序列的分布信息和前體結構能量信息，采用貝葉斯模型打分預測新miRNA。按照TPM=T×106/N（T表示miRNA的tags，N表示總miRNA的tags）算法對AaM和AaS中每個miRNA的表達量進行歸一化處理。利用基迪奧在線工具集合（www.omicshare.com）對AaM和AaS的miRNA進行Venn分析，采用默認參數(shù)，分別篩選出AaM和AaS的共有miRNA及特有miRNA。

1.6 共有與特有miRNA的靶mRNA的GO和KEGG數(shù)據(jù)庫注釋分析

采用RNA hybrid（v2.1.2）+svm_light（v6.01）、Miranda（v3.3a）和TargetScan（Version 7.0）軟件[28-31]分別預測共有miRNA和特有miRNA的靶mRNA，將預測結果的交集作為可信度高的靶標集合。利用BLAST工具將靶mRNA比對到GO數(shù)據(jù)庫（http:// geneontology.org/）和KEGG數(shù)據(jù)庫（https://www. kegg.jp/），獲得相應的功能和通路注釋信息，并統(tǒng)計比對上各功能條目或通路（pathway）的靶mRNA數(shù)量。

1.7 DEmiRNA的篩選和靶mRNA的預測及分析

以|log2fold change (FC)|≥1且≤0.05為標準，利用edgeR軟件[32]篩選AaM vs AaS的DEmiRNA。前期已利用鏈特異性建庫的RNA-seq技術對球囊菌菌絲和孢子分別進行測序，其中高質量的mRNA組學數(shù)據(jù)可作為本研究中靶mRNA的數(shù)據(jù)來源（未發(fā)表數(shù)據(jù)）。分別采用RNA hybrid+svm_light、Miranda和TargetScan軟件[28-31]預測DEmiRNA的靶mRNA，將預測結果的交集作為可信度高的靶標集合。利用BLAST工具將DEmiRNA的靶mRNA比對到GO和KEGG數(shù)據(jù)庫，獲得相應的功能和通路注釋信息，并統(tǒng)計比對上的功能條目和通路的靶mRNA數(shù)量。根據(jù)上述DEmiRNA與mRNA的靶向結合關系，按照自由能＜-40 kcal·mol-1和＜0.05的閾值篩選靶向結合關系并構建調控網(wǎng)絡，通過Cytoscape軟件進行調控網(wǎng)絡的可視化。

1.8 DEmiRNA的Stem-loop RT-qPCR驗證

隨機選取5個DEmiRNA（miR5658-x、miR-10285- y、miR-3245-y、miR4404-x、miR-9-z）進行Stem-loop RT-qPCR驗證。參照Chen等[33]和郭睿等[34]的方法，利用DNAMAN軟件（Lynnon Biosoft公司，美國）設計上述DEmiRNA的Stem-loop引物、特異性上游引物和通用下游引物。委托上海生工生物工程股份有限公司合成引物（表1）。利用總RNA抽提試劑盒（Promega公司，中國）分別抽提菌絲和孢子的總RNA，然后用RNase-free DNase I去除各自基因組DNA殘留。吸取1 μL RNA，用NanoDrop One（Thermo Scientific公司，美國）分別對菌絲、孢子的RNA濃度進行測定。利用Stem-loop引物，按照cDNA第1鏈合成試劑盒（TaKaRa公司，日本）說明書進行總RNA的反轉錄，得到的cDNA作為模板進行qPCR。選用球囊菌的（5417）作為內(nèi)參基因。反應體系為20 μL：SYBR Green Dye 10 μL，特異性上游引物1 μL，通用下游引物1 μL，cDNA模板1 μL，RNA-Free H2O 7 μL。每個反應進行3次技術重復。采用2-ΔΔCt法計算DEmiRNA的相對表達量，然后利用GraphPad Prism 7軟件進行Student’s-test檢驗及繪圖。

表1 本研究使用的引物

2 結果

2.1 測序數(shù)據(jù)的質控與評估

AaM和AaS分別測得12 982 320和12 708 832條raw reads，經(jīng)嚴格過濾和質控后分別得到10 800 101和9 888 848條clean tags，占raw reads的比例分別為84.34%和78.83%；此外，比對上參考基因組的clean tags比例分別為80.59%和71.60%（表2）。上述結果表明本研究的sRNA-seq數(shù)據(jù)質量良好，可滿足后續(xù)分析。

2.2 球囊菌菌絲和孢子miRNA的預測及比較分析

基于AaM和AaS的高質量數(shù)據(jù)分別預測出193和275個miRNA。其中，AaM中miRNA的長度介于18—26 nt（圖1-A），AaS中miRNA的長度介于18—24 nt（圖1-C）；AaM和AaS中首位堿基偏向于U的miRNA數(shù)量最多，占比分別達到30.77%和45.73%（圖1-B、1-D）。

AaM中表達量最高的是miR6478-x、miR10516-x和miR482-x，TPM值分別達到294 585.5787、200 923.7875和99 563.767；AaS中表達量最高的同樣為miR6478-x、miR482-x和miR10516-x，TPM值分別為209 814.8692、71 995.2983和71 407.5815。AaM和AaS中表達量最高的前10位miRNA的詳細信息如表3和表4所示。

2.3 球囊菌菌絲和孢子共有和特有miRNA的篩選及靶mRNA分析

Venn分析結果顯示，AaM和AaS的共有miRNA為76個（19.4%），特有miRNA分別為117個（29.8%）和199個（50.8%）。AaM和AaS的共有miRNA可靶向5 946個mRNA，二者的特有miRNA可分別靶向6 141和6 346個mRNA。GO數(shù)據(jù)庫注釋結果顯示，上述共有miRNA的靶mRNA可注釋到42個功能條目，其中細胞組分大類注釋靶mRNA數(shù)最多的條目是細胞組件（1 027）、細胞（1 027）和細胞器（692），分子功能大類注釋靶mRNA數(shù)最多的條目是催化活性（1 199）、結合（871）和轉運活性（129），生物學進程大類注釋靶mRNA數(shù)最多的條目是代謝進程（1 303）、細胞進程（1 301）和單一組織進程（975）；AaM的特有miRNA的靶mRNA可注釋到42個功能條目，其中細胞組分大類注釋靶mRNA數(shù)最多的條目是細胞組件（1 095）、細胞（1 095）和細胞器（736），分子功能大類注釋靶mRNA數(shù)最多的條目是催化活性（1 241）、結合（900）和轉運活性（128），生物學進程大類注釋靶mRNA數(shù)最多的條目是代謝進程（1 362）、細胞進程（1 350）和單一組織進程（1 021）；AaS的特有miRNA的靶mRNA可注釋到42個功能條目，其中細胞組分大類注釋靶mRNA數(shù)最多的條目是細胞組件（1 121）、細胞（1 121）和細胞器（762），分子功能大類注釋靶mRNA數(shù)最多的條目是催化活性（1 256）、結合（915）和轉運活性（126），生物學進程大類注釋靶mRNA數(shù)最多的條目是代謝進程（1 383）、細胞進程（1 379）和單一組織進程（1 039）。括號內(nèi)的數(shù)字表示注釋到該條目的靶mRNA數(shù)。

表2 sRNA-seq數(shù)據(jù)概覽

A：AaM中miRNA的長度分布Length distribution of miRNAs in AaM；B：AaM中miRNA的首位堿基偏向性 First base bias of miRNAs in AaM；C：AaS中miRNA的長度分布Length distribution of miRNAs in AaS；D：AaS中miRNA的首位堿基偏向性 First base bias of miRNAs in AaS

表3 AaM中表達量最高的前10位miRNA

表4 AaS中表達量最高的前10位miRNA

KEGG數(shù)據(jù)庫注釋結果顯示，上述共有miRNA和二者特有miRNA的靶mRNA均可注釋到120條通路，涉及代謝、遺傳信息處理、環(huán)境信息處理和細胞進程4個大類。其中，共有miRNA的靶mRNA注釋數(shù)量最多的通路是新陳代謝總覽（690）、翻譯（264）、碳水化合物代謝（232）、折疊、分類與降解（223）、運輸與代謝（186）、傳染性疾?。?56）、信號轉導（138）、能量代謝（133）、細胞生長與死亡（133）和脂質代謝（124）；AaM的特有miRNA的靶mRNA注釋數(shù)量最多的是新陳代謝總覽（712）、翻譯（283）、碳水化合物代謝（232）、折疊與分類降解（231）、氨基酸代謝（202）、運輸與代謝（189）、傳染性疾病（164）、信號轉導（144）、能量代謝（137）和細胞生長與死亡（134）；AaS的特有miRNA的靶mRNA注釋數(shù)量最多的是新陳代謝總覽（725）、翻譯（292）、碳水化合物代謝（236）、折疊、分類與降解（236），氨基酸代謝（209）、運輸與代謝（194）、傳染性疾?。?70）、信號轉導（148）、細胞生長與死亡（138）和能量代謝（135）。進一步分析發(fā)現(xiàn)，AaM的116個特有miRNA靶向的287個mRNA可注釋到次級代謝物的生物合成通路；AaS的197個特有miRNA靶向的288個mRNA可注釋到次級代謝通路的生物合成通路；44個共有miRNA靶向的14個mRNA可注釋到自噬通路；56個共有miRNA靶向的6個mRNA，以及AaM的95個特有miRNA靶向的28個mRNA可注釋到MAPK信號通路。括號內(nèi)的數(shù)字表示注釋到該通路的靶mRNA數(shù)。

2.4 球囊菌菌絲和孢子中miRNA的差異表達分析

從AaM vs AaS中篩選出93個DEmiRNA，包括65個上調miRNA和28個下調miRNA；其中上調幅度最大的是novel-m0040-3p（log2fc=20.65019，=5.98E-17），其次是novel-m0016-3p（log2fc= 20.45754，=3.77E-15）和miR319-y（log2fc=20.23515，=4.75E-13）；下調幅度最大的為miR-4028-y（log2fc= -19.6472765，=4.77E-10），其次為miR-4171-x（log2fc=-18.2322391，=0.00049）和miR7787-y（log2fc=-18.1067082，=0.000979）（表5）。上述DEmiRNA能靶向結合6 090個mRNA。

DEmiRNA的靶mRNA可注釋到38個功能條目，包括細胞進程（1 268）、代謝進程（1 254）、單一組織進程（981）、定位（332）和生物學調控（290）等15條生物學進程相關條目；催化活性（1 230）、結合（879）、轉運活性（125）、結構分子活性（35）及分子功能調節(jié)器（18）等11條分子功能相關條目；細胞（1058）、細胞組件（1 058）、細胞器（713）、大分子復合物（390）和細胞膜（332）等12條細胞組分相關條目（圖2）。此外，這些靶mRNA可注釋到120條通路，注釋數(shù)量最多的是新陳代謝通路（652）、次級代謝物的生物合成（285）、抗生素的生物合成（212）、微生物在不同環(huán)境中的代謝（172）、氨基酸的生物合成（115）、碳代謝（101）、嘌呤代謝（81）、核糖體（78）、剪接體（77）及RNA轉運（76）（表6）。括號內(nèi)的數(shù)字代表注釋到該條目（通路）的靶mRNA數(shù)。

表5 AaM vs AaS比較組中前10位上調和下調miRNA

2.5 DEmiRNA的調控網(wǎng)絡構建及分析

根據(jù)靶向結合關系構建調控網(wǎng)絡，分析結果顯示13個DEmiRNA可靶向結合131個mRNA；同一個DEmiRNA可靶向結合多個mRNA，如miR-4968-y可靶向結合多達118個mRNA；同時，部分mRNA可靶向結合1—2個DEmiRNA，如KZZ97168.1分別靶向結合miR-4968-y和miR-6769-y，KZZ91616靶向結合miR5782-y（圖3）。

2.6 DEmiRNA的RT-qPCR驗證

RT-qPCR結果顯示，miR5658-x、miR-10285-y、miR-3245-y、miR4404-x和miR-9-z的表達變化趨勢與測序數(shù)據(jù)相符（圖4），證實了本研究中sRNA-seq數(shù)據(jù)的可靠性。

圖2 AaM vs AaS比較組中DEmiRNA的靶mRNA的GO數(shù)據(jù)庫注釋

表6 AaM vs AaS比較組中DEmiRNA的靶mRNA注釋數(shù)前10位通路

3 討論

MiRNA已被證實能夠參與調控真菌的菌絲生長和孢子形成過程[35-38]。例如，Shao等通過比較分析篩選出冬蟲夏草（）無性生殖階段和有性生殖階段的19個DEmiRNA，進而通過過表達和敲除實驗證實milR4和milR16參與了菌絲生長過程的調控[36]。前期研究中，為最大限度鑒定miRNA，筆者利用sRNA-seq技術對球囊菌菌絲和孢子的混合樣品進行測序，利用miRDeep2軟件鑒定到118個novel miRNA，這是關于球囊菌miRNA的首例報道[24]。然而，由于測序得到的混合數(shù)據(jù)無法區(qū)分來源于菌絲的數(shù)據(jù)和來源于孢子的數(shù)據(jù)，導致難以進一步對菌絲和孢子中的miRNA進行數(shù)量、結構特征、表達譜、靶mRNA及調控網(wǎng)絡的比較分析和深入挖掘。因此，本研究首先在實驗室條件下獲得純培養(yǎng)的球囊菌，利用sRNA-seq技術對純凈的球囊菌菌絲樣品、孢子樣品分別進行測序，基于二者的高質量sRNA組學數(shù)據(jù)分別鑒定到193和275個miRNA，它們的長度分別介于18—26和18—24 nt。分布在18 nt長度的miRNA數(shù)量最多，且首位堿基主要偏向U，其結構特征與灰蓋鬼傘菌（）[37]、新月彎孢（）[39]和馬爾尼菲青霉（）[40]等真菌以及蜜蜂和棉花等動植物[41-42]的miRNA結構高度相似。本研究中，有76個miRNA同時在球囊菌菌絲和孢子中表達，占二者全部miRNA的比例分別為39.38%和27.64%，推測這些共有miRNA在球囊菌的不同生長發(fā)育時期均具有一定的調控功能；此外，分別有117和199個miRNA在球囊菌菌絲和孢子中特異性表達，鑒于孢子是球囊菌的休眠態(tài)，新陳代謝等生命活動較之菌絲更低，該結果一定程度說明較多的miRNA通過在孢子中特異性表達發(fā)揮更強的基因表達抑制（或降解）作用。

圖3 DEmiRNA-mRNA的調控網(wǎng)絡

RT-qPCR 組中， *表示P＜0.05， **表示P＜0.01 In RT-qPCR group, * indicates P＜0.05, ** indicates P＜0.01

共有miRNA的靶mRNA可注釋到代謝進程、生殖、生殖進程、生長等42個功能條目，以及代謝途徑、次級代謝物的生物合成、不同環(huán)境下微生物的代謝、氨基酸的生物合成、碳代謝及嘌呤代謝和氧化磷酸化等120條通路，表明共有miRNA在球囊菌菌絲和孢子的生長發(fā)育、物質和能量代謝、環(huán)境適應等方面具有廣泛的調控功能。菌絲的特有miRNA可靶向結合6 141個mRNA，這些靶標可注釋到42個功能條目和120條通路；孢子的特有miRNA可靶向結合6 346個mRNA，同樣可注釋到42個功能條目和120條通路。對于真菌孢子中是否存在轉錄和翻譯等生命活動，相關研究報道很少。筆者團隊前期通過分子生物學和組學手段證實另一種廣泛存在的蜜蜂真菌病原東方蜜蜂微孢子蟲（）的孢子中存在基因轉錄[43]。本研究中，在球囊菌孢子中鑒定到199個特有miRNA，上述結果一定程度表明球囊菌孢子中同樣存在基因轉錄以及miRNA介導的基因表達調控現(xiàn)象。此外，共有93個miRNA在菌絲和孢子中差異表達，其中上調miRNA（65）的數(shù)量明顯多于下調miRNA（28），進一步說明對于球囊菌孢子，除了具有較多的特異性表達miRNA外，其部分miRNA還能通過上調表達量增強對靶基因的抑制（或降解）作用，從而維持較低的新陳代謝水平。推測這對于球囊菌孢子抵抗外界不良環(huán)境及長期存活具有重要意義。

真菌在侵染宿主時會分泌一些次級代謝物促進自身增殖并使宿主致死[44-45]。例如，球孢白僵菌（）和卵孢白僵菌（）合成及分泌的草酸、類草酸晶體和檸檬酸等次級代謝物可協(xié)同致死宿主[44]。有研究表明真菌在菌絲狀態(tài)產(chǎn)生的次級代謝物主要與真菌毒素有關[45]。本研究發(fā)現(xiàn)，球囊菌菌絲特有的116個miRNA（miR-11971-y、miR-14-x、miR-12227-y等）的287個靶mRNA（KZZ86592.1、KZZ86645.1和KZZ86652.1等）可注釋到次級代謝物的生物合成通路，表明相應的菌絲特有miRNA參與調控次級代謝物的生物合成過程，進而影響球囊菌毒素的合成。真菌病原對昆蟲寄主的侵染能力取決于蛋白酶、幾丁質酶及脂酶等毒力因子的協(xié)同作用，以及菌絲對圍食膜、腸壁和表皮的機械破壞力[46]。Cornman等[46]通過同源性比較發(fā)現(xiàn)球囊菌包含61個脂酶基因、51個蛋白酶基因和4個幾丁質酶基因。本研究分別有16個（miR-14-x、miR11173-y、miR-1002-x等）和9個（miR-12227-y、miR-4171-x、miR-1002-x等）菌絲特有miRNA靶向幾丁質酶合成相關的mRNA（KZZ93915.1和KZZ93066.1），暗示這些miRNA參與調控幾丁質酶合成，在球囊菌突破宿主幾丁質體表過程中發(fā)揮重要調控功能。次級代謝物還能影響真菌的孢子萌發(fā)[45,47]，例如玉米赤霉烯酮和禾谷鐮孢（）的次級代謝物可誘導禾谷鐮孢分生孢子和菌落的產(chǎn)生[47]。本研究發(fā)現(xiàn)，197個孢子特有miRNA調控的288個靶mRNA（KZZ86592.1、KZZ86645.1和KZZ86652.1等）注釋到了次級代謝物的生物合成，筆者推測孢子特有miRNA可能通過調控次級代謝物合成的相關mRNA影響孢子的萌發(fā)。李瓊等[48]研究發(fā)現(xiàn)，mro-miR-33負調控羅伯茨綠僵菌（）的產(chǎn)孢關鍵基因的表達，敲除mro-miR-33后的表達量明顯上調，同時病原的產(chǎn)孢量增加。本研究中，miR-33-x與mro-miR-33高度同源且在孢子中特異性表達，可能參與調控球囊菌的雜交產(chǎn)孢。細胞自噬與絲狀真菌的產(chǎn)孢、程序性細胞死亡及致病力緊密相關[49]。在米曲霉（）中，、及等自噬基因的功能缺失將影響其產(chǎn)孢過程[50]。此外，有研究表明煙曲霉（）也依賴自噬來調控孢子的形成[51]。本研究中，菌絲和孢子的44個共有miRNA的14個靶mRNA（KZZ87046.1、KZZ87260.1和KZZ87338.1等）可注釋到自噬通路，推測這些共有miRNA可通過負調控相關靶mRNA調節(jié)球囊菌的雜交產(chǎn)孢。

絲裂原活化蛋白激酶（mitogen-activated protein kinase，MAPK）級聯(lián)反應參與真菌的交配繁殖、致病性、孢子形成以及毒力水平的調控[52-53]，其依賴的MAPK激酶可協(xié)助病原真菌侵染宿主[54]。前期研究發(fā)現(xiàn)，對于侵染意蜂幼蟲的球囊菌，有48個差異表達基因富集在MAPK信號通路且表達水平隨著球囊菌脅迫時間的延長而顯著增強[19]；但對于侵染中蜂幼蟲的球囊菌，富集在MAPK信號通路的11個差異表達基因均表現(xiàn)為下調趨勢[20]。本研究中，菌絲和孢子的56個共有miRNA（let-7-x、miR-10285-y和miR-11980-x等）的靶mRNA，以及菌絲的95個特有miRNA（miR-11971-y、miR-12227-y和miR-14-x等）的28個靶mRNA均可注釋到MAPK信號通路，表明上述miRNA與MAPK信號通路具有潛在的調控關系。目前，筆者團隊已利用sRNA-seq技術對正常及球囊菌侵染的意蜂幼蟲腸道、中蜂幼蟲腸道進行測序，下一步將過濾得到侵染兩種蜜蜂幼蟲的球囊菌的sRNA組學數(shù)據(jù)，并結合本研究中的球囊菌孢子的sRNA組學數(shù)據(jù)進行比較分析，更深入地探討球囊菌對不同抗性蜜蜂幼蟲的侵染機制。

本研究中，DEmiRNA與mRNA之間存在較為復雜的調控關系，miR-4968-y可同時被118個mRNA同時靶向結合，表明miR-4968-y處于調控網(wǎng)絡的核心位置，可能在球囊菌菌絲和孢子的生長和發(fā)育過程發(fā)揮關鍵調控功能，值得進一步深入研究。根癌農(nóng)桿菌介導（AtMT）的真菌遺傳轉化體系已成功用于球孢白僵菌、金龜子綠僵菌（）和蠟蚧輪枝菌（）等昆蟲病原真菌的基因功能研究[55-57]。目前，球囊菌的基因功能研究未見報道。利用AtMT技術對本研究篩選出的菌絲和孢子的關鍵miRNA進行轉基因操作，進而探究其在菌絲和孢子生長發(fā)育以及病原致病性方面的功能是下一步的工作重點。

4 結論

分別在球囊菌菌絲和孢子中鑒定出193和275個miRNA，二者中特異性表達的miRNA分別為117和199個；菌絲和孢子中的miRNA具有類似的結構特征，但表達譜表現(xiàn)出明顯差異；菌絲和孢子可能通過特異性表達和差異表達部分miRNA對其生長、發(fā)育和生殖進行調控。

[1] ARONSTEIN K A, MURRAY K D. Chalkbrood disease in honey bees.,2010, 103(Suppl. 1): S20-S29.

[2] MAXFIELD-TAYLOR S A, MUJIC A B, RAO S. First detection of the larval chalkbrood disease pathogen(Ascomycota: Eurotiomycetes: Ascosphaerales) in adult bumble bees., 2015, 10(4): e0124868.

[3] EVISON S E. Chalkbrood: epidemiological perspectives from the host-parasite relationship., 2015, 10: 65-70.

[4] SPILTOIR C F. Life cycle of()., 1955, 42(6): 501-518.

[5] WIJAYAWARDENE N N, HYDE K D, LUMBSCH H T, LIU J K, MAHARACHCHIKUMBURA S S N, EKANAYAKA A H, TIAN Q, PHOOKAMSAK R. Outline of: 2017., 2018, 88: 167-263.

[6] P?GGELER S. Mating-type genes for classical strain improvements of ascomycetes., 2001, 56(5/6): 589-601.

[7] FLORES J M, SPIVAK M, GUTIERREZ I. Spores ofcontained in wax foundation can infect honeybee brood., 2005, 108(1/2): 141-144.

[8] FLORES J M, GUTIéRREZ I, ESPEJO R. The role of pollen in chalkbrood disease in: transmission and predisposing conditions., 2005, 97(6): 1171-1176.

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

[10] BURKLEW C E, XIE F U, AshlockJ, ZhangB h. Expression of microRNAs and their targets regulates floral development in tobacco ()., 2014, 14(2): 299-306.

[11] KIDNER C A, MARTIENSSEN R A. Spatially restricted microRNA directs leaf polarity through ARGONAUTE1., 2004, 428(6978): 81-84.

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

[13] GRIMSON A, SRIVASTAVA M, FAHEY B, WOODCROFT B J, CHIANG H R, KING N, DEGNAN B M, ROKHSAR D S, BARTEL D P. Early origins and evolution of microRNAs and Piwi-interacting RNAs in animals., 2008, 455(7217): 1193-1197.

[14] LLAVE C, KASSCHAU K D, RECTOR M A, CARRINGTON J C. Endogenous and silencing-associated small RNAs in plants., 2002, 14(7): 1605-1619.

[15] MOLNáR A, SCHWACH F, STUDHOLME D J, THUENEMANN E C, BAULCOMBE D C. MiRNAs control gene expression in the single-cell alga., 2007, 447(7148): 1126-1129.

[16] ZHAO T, LI G, MI S, LI S, HANNON G J, WANG X J, QI Y. A complex system of small RNAs in the unicellular green alga., 2007, 21(10): 1190-1203.

[17] CHEN D F, GUO R, XU X J, XIONG C L, LIANG Q, ZHENG Y Z, LUO Q, ZHANG Z N, HUANG Z J, KUMAR D, XI W J, ZOU X, LIU M. Uncovering the immune responses oflarval gut toinfection utilizing transcriptome sequencing., 2017, 621: 40-50.

[18] GUO R, CHEN D F, DIAO Q Y, XIONG C L, ZHENG Y Z, HOU C S. Transcriptomic investigation of immune responses of thelarval gut infected by., 2019, 166: 107210.

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

[20] 郭睿, 陳大福, 黃枳腱, 梁勤, 熊翠玲, 徐細建, 鄭燕珍, 張曌楠, 解彥玲, 童新宇, 侯志賢, 江亮亮, 刀晨. 球囊菌脅迫中華蜜蜂幼蟲腸道過程中病原的轉錄組學研究. 微生物學報, 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)

[21] 張曌楠, 熊翠玲, 徐細建, 黃枳腱, 鄭燕珍, 駱群, 劉敏, 李汶東, 童新宇, 張琦, 梁勤, 郭睿, 陳大福. 蜜蜂球囊菌的參考轉錄組組裝及SSR分子標記開發(fā). 昆蟲學報, 2017, 60(1): 34-44.

ZHANG Z N, XIONG C L, XU X J, HUANG Z J, ZHENG Y Z, LUO Q, LIU M, LI W D, TONG X Y, ZHANG Q, LIANG Q, GUO R, CHEN D F.assembly of a reference transcriptome and development of SSR markers for, 2017, 60(1): 34-44. (in Chinese)

[22] 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., 2018, 156: 1-5.

[23] GUO R, CHEN D F, CHEN H Z, FU Z M, XIONG C L, HOU C S, ZHENG Y Z, GUO Y L, WANG H P, DU Y, DIAO Q Y. Systematic investigation of circular RNAs in, a fungal pathogen of honeybee larvae., 2018, 678: 17-22.

[24] 郭睿, 王海朋, 陳華枝, 熊翠玲, 鄭燕珍, 付中民, 趙紅霞, 陳大福. 蜜蜂球囊菌的microRNA鑒定及其調控網(wǎng)絡分析. 微生物學報, 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)

[25] JENSEN A B, ARONSTEIN K, FLORES J M, VOJVODIC S, PALACIO M A, SPIVAK M. Standard methods for fungal brood disease research., 2013, 52(1): DOI: 10.3896/IBRA.1.52.1.13.

[26] 郭睿, 杜宇, 熊翠玲, 鄭燕珍, 付中民, 徐國鈞, 王海朋, 陳華枝, 耿四海, 周丁丁, 石彩云, 趙紅霞, 陳大福. 意大利蜜蜂幼蟲腸道發(fā)育過程中的差異表達microRNA及其調控網(wǎng)絡. 中國農(nóng)業(yè)科學, 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 oflarval gut., 2018, 51(21): 4197-4209. (in Chinese)

[27] 杜宇, 范小雪, 蔣海賓, 王杰, 范元嬋, 祝智威, 周丁丁, 萬潔琦, 盧家軒, 熊翠玲, 鄭燕珍, 陳大福, 郭睿. 微小RNA及其介導的競爭性內(nèi)源RNA調控網(wǎng)絡在意大利蜜蜂工蜂中腸發(fā)育過程中的潛在作用. 中國農(nóng)業(yè)科學, 2020, 53(12): 2512-2526.

DU Y, FAN X X, JIANG H B, WANG J, FAN Y C, ZHU Z Z, ZHOU D D, WAN J Q, LU J X, XIONG C L, ZHENG Y Z, CHEN D F, GUO R. The potential role of microRNAs and microRNA-mediated competitive endogenous networks during the developmental process ofa worker’s midgut., 2020, 53(12): 2512-2526. (in Chinese)

[28] RITCHIE W. microRNA target prediction., 2017, 1513: 193-200.

[29] REHMSMEIER M, STEFFEN P, HOCHSMANN M, GIEGERICH R. Fast and effective prediction of microRNA/target duplexes., 2004, 10(10): 1507-1517.

[30] KRüGER J, REHMSMEIER M. RNAhybrid: microRNA target prediction easy, fast and flexible., 2006, 34: W451-W454.

[31] ALLEN E, XIE Z x, GUSTAFSON A M, CARRINGTON J C. microRNA-directed phasing during-acting siRNA biogenesis in plants., 2005, 121(2): 207-221.

[32] ROBINSON M D, MCCARTHY D J, SMYTH G K. EdgeR: a Bioconductor package for differential expression analysis of digital gene expression data., 2010, 26(1): 139-140.

[33] CHEN C F, RIDZON D A, BROOMER A J, ZHOU Z H, LEE D H, NGUYEN J T, BARBISIN M, XU N L, MAHUVAKAR V R, ANDERSEN M R, LAO K Q, LIVAK K J, GUEGLER K J. Real-time quantification of microRNAs by stem-loop RT-PCR., 2005, 33(20): e179.

[34] 郭睿, 杜宇, 周倪紅, 劉思亞, 熊翠玲, 鄭燕珍, 付中民, 徐國鈞, 王海朋, 耿四海, 周丁丁, 陳大福. 意大利蜜蜂幼蟲腸道在球囊菌脅迫后期的差異表達微小RNA及其靶基因分析. 昆蟲學報, 2019, 62(1): 49-60.

GUO R, DU Y, ZHOU N H, LIU S Y, XIONG C L, ZHENG Y Z, FU Z M, XU G J, WANG H P, GENG S H, ZHOU D D, CHEN D F. Comprehensive analysis of differentially expressed microRNAs and their target genes in the larval gut ofduring the late stage ofstress., 2019, 62(1): 49-60. (in Chinese)

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

[36] SHAO Y, TANG J, CHEN S L, WU Y H, WANG K, MA B, ZHOU Q M, CHEN A H, WANG Y L. milR4 and milR16 mediated fruiting body development in the medicinal fungus., 2019, 10: 83.

[37] LAU A Y T, CHENG X, CHENG C K, NONG W, CHEUNG M K, CHAN R H, HUI J H L, KWAN H S. Discovery of microRNA-like RNAs during early fruiting body development in the model mushroom., 2018, 13(9): e0198234.

[38] ZENG W P, WANG J, WANG Y, LING J, FU Y P, XIE J T, JIANG D H, CHEN T, LIU H Q, CHENG J S. Dicer-like proteins regulate sexual development via the biogenesis of perithecium-specific microRNAs in a plant pathogenic fungus., 2018, 9: 818.

[39] LIU T, HU J, ZUO Y, JIN Y, HOU J. Identification of microRNA-like RNAs fromassociated with maize leaf spot by bioinformation analysis and deep sequencing., 2016, 291(2): 587-596.

[40] LAU S K, CHOW W N, WONG A Y, YEUNG J M, BAO J, ZHANG N, LOK S, WOO P C, YUEN K Y. Identification of microRNA-Like RNAs in mycelial and yeast phases of the thermal dimorphic fungus., 2013, 7(8): e2398.

[41] LI J, HULL J J, LIANG S, WANG Q, CHEN L, ZHANG Q, WANG M, MANSOOR S, ZHANG X, JIN S. Genome-wide analysis of cotton miRNAs during whitefly infestation offers new insights into plant-herbivore interaction.,2019, 20(21): e5357.

[42] 熊翠玲, 杜宇, 陳大福, 鄭燕珍, 付中民, 王海朋, 耿四海, 陳華枝, 周丁丁, 吳素珍, 石彩云, 郭睿. 意大利蜜蜂幼蟲腸道的miRNAs的生物信息學預測及分析. 應用昆蟲學報, 2018, 55(6): 1023-1033.

XIONG C L, DU Y, CHEN D F, ZHENG Y Z, FU Z M, WANG H P, GENG S H, CHEN H Z, ZHOU D D, WU S Z, SHI C Y, GUO R. Bioinformatic prediction and analysis of miRNAs inthelarval gut., 2018, 55(6): 1023-1033. (in Chinese)

[43] 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, KUMAR D. First identification of long non-coding RNAs in fungal parasite.,2018, 49(5): 660-670.

[44] 黃居敏, 張普照, 王芳, 李旸, 馮少俊, 楊明. 白僵菌的代謝產(chǎn)物及藥理活性研究進展. 中國生化藥物雜志, 2014, 34(9): 167-173.

HUANG J M, ZHANG P Z, WANG F, LI y, FENG S J, YANG MAdvanced studies on metabolites and pharmacological of.,2014, 34(9): 167-173. (in Chinese)

[45] CALVO A M, WILSON R A, BOK J W, KELLER N PRelationship between secondary metabolism and fungal development,2002, 66(3): 447-459.

[46] CORNMAN R S, BENNETT A K, MURRAY K D, EVANS J D, ELSIK C G, ARONSTEIN K. Transcriptome analysis of the honey bee fungal pathogen,: implications for host pathogenesis., 2012, 13: 285.

[47] WOLF J C, MIROCHA C J. Regulation of sexual reproduction in(“Graminearum”) by F-2 (Zearalenone)., 1973, 19(6): 725-734.

[48] 李瓊, 崔春來, 宋紅生, 王四寶. mro-miR-33在綠僵菌產(chǎn)孢中的作用. 菌物學報, 2017, 36(6): 671-678.

LI Q, CUI C L, SONG H S, WANG S B. The effects of mro-miR-33 on the conidial production in.,2017, 36(6): 671-678. (in Chinese)

[49] 閆思源, 姜學軍. 細胞自噬及真菌中自噬研究概述. 菌物學報, 2015, 34(5): 871-879.

YAN S Y, JIANG X J. Overview of autophagy and related research in fungi., 2015, 34(5): 871-879. (in Chinese)

[50] KIKUMA T, KITAMOTO K. Analysis of autophagy inby disruption of Aoatg13, Aoatg4, and Aoatg15 genes., 2011, 316(1): 61-69.

[51] RICHIE D L, FULLER K K, FORTWENDEL J, MILEY M D, MCCARTHY J W, FELDMESSER M, RHODES J C, ASKEW D S. Unexpected link between metal ion deficiency and autophagy in., 2007, 6(12):2437-2447.

[52] ZHAO X, MEHRABI R, XU J R. Mitogen-activated protein kinase pathways and fungal pathogenesis., 2007, 6(10): 1701-1714.

[53] IGBARIA A, LEV S, ROSE M S, LEE B N, HADAR R, DEGANI O, HORWITZ B A. Distinct and combined roles of the MAP kinases ofin virulence and stress responses., 2008, 21(6): 769-780.

[54] CHEN X X, XU C, QIAN Y, LIU R, ZHANG Q Q, ZENG G H, ZHANG X, ZHAO H, FANG W G. MAPK cascade-mediated regulation of pathogenicity, conidiation and tolerance to abiotic stresses in the entomopathogenic fungus., 2016, 18(3): 1048-1062.

[55] LECLERQUE A, WAN H, ABSCHüTZ A, CHEN S, MITINA G V, ZIMMERMANN G, SCHAIRER H.-mediated insertional mutagenesis (AIM) of the entomopathogenic fungus., 2004, 45(2): 111-119.

[56] FANG W, PEI Y, BIDOCHKA M J. Transformation ofmediated by., 2006, 52(7): 623-626.

[57] ZHANG Y J, ZHAO J J, XIE M, PENG D L.-mediated transformation in the entomopathogenic fungusand development of benzimidazole fungicide resistant strains.,2014, 105: 168-173.

Comparative Analysis of microRNAs and corresponding target mRNAs inmycelium and spore

CHEN HuaZhi, ZHU ZhiWei, JIANG HaiBin, WANG Jie, FAN YuanChan, FAN XiaoXue, WAN JieQi, LU JiaXuan, 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)

【】exclusively infects honeybee larvae, resulting in chalkbrood disease. The objective of this study is to clarify the differences of number, structure and expression pattern of miRNAs betweenmycelium and spore based on deep sequencing and comparative analysis of purified mycelia (AaM) and spores (AaS) using small RNA-seq (sRNA-seq) and bioinformatics, and reveal the potential relationship between common miRNAs, specific miRNAs, differentially expressed miRNAs (DEmiRNAs) and their target mRNAs and the growth and development of mycelium and spore as well as pathogenesis of【】The pure culture ofwas gained under lab condition. AaM and AaS were respectively sequenced using sRNA-seq technology. Clean tags were obtained after filtration and quality control of raw reads. Common miRNAs and specific miRNAs in AaM and AaS were screened out using Venn analysis. DEmiRNAs in the AaM vs AaS comparison group were filtered out following the criteria of≤0.05 and |log2fold change|≥1. Target mRNAs of common miRNAs, specific miRNAs and DEmiRNAs were predicted using related bioinformatic software. Target mRNAs mentioned above were respectively annotated to GO database and KEGG database. The regulatory network of DEmiRNAs and target mRNAs was constructed on basis of target binding relationship, followed by visualization with Cytoscape. RT-qPCR was conducted to verify the reliability of the sequencing data.【】In total, 12 982 320 and 12 708 832 raw reads were produced from AaM and AaS, and after strict quality control, 10 800 101 and 9 888 848 clean tags were gained, respectively. The length of specific miRNAs in AaM was distributed among 18-26 nt, while that in AaS was distributed among 18-24 nt. Additionally, most of the miRNAs were distributed in 18 nt. MiRNAs with the first base U in both AaM and AaS were the most abundant. miRNAs with the highest expression levels in both AaM and AaS were miR6478-x, miR10516-x and miR482-x. These common miRNAs could target 5 946 mRNAs, while specific miRNAs in AaM and AaS could bind to 6 141 and 6 346 mRNAs, respectively. Targets of common miRNAs were annotated to 42 functional terms such as metabolism process, cellular process and catalytic activity, and 120 pathways including translation, carbohydrate metabolism and energy metabolism. In addition, a total of 93 DEmiRNAs were identified in AaM vs AaS comparison group, targeting 6 090 mRNAs annotated to 38 functional terms and 120 pathways. Moreover, complicated regulatory networks were formed between DEmiRNAs and target mRNAs, with miR-4968-y located in the center and linked to as many as 118 mRNAs. RT-qPCR result demonstrated the expression trend of five DEmiRNAs was consistent with that in the sequencing result, confirming the reliability of our sequencing data.【】The structures of miRNAs inmycelium and spore were similar, whereas their expression patterns were obviously different; mycelium and spore may specifically and differentially express part of miRNAs to regulate their growth, development and reproduction.

; mycelium; spore; miRNA; mRNA; target binding

10.3864/j.issn.0578-1752.2020.17.017

2019-12-25；

2020-02-04

國家自然科學基金（31702190）、國家現(xiàn)代農(nóng)業(yè)產(chǎn)業(yè)技術體系建設專項（CARS-44-KXJ7）、福建省自然科學基金（2018J05042）、福建省教育廳中青年教師教育科研項目（JAT170158）、福建農(nóng)林大學杰出青年科研人才計劃（xjq201814）、福建農(nóng)林大學科技創(chuàng)新專項（CXZX2017342, CXZX2017343）、福建農(nóng)林大學優(yōu)秀碩士學位論文資助基金（陳華枝）

陳華枝，E-mail：CHZ0720@outlook.com。祝智威，E-mail：zzw15235470398@163.com。陳華枝和祝智威為同等貢獻作者。通信作者郭睿，E-mail：ruiguo@fafu.edu.cn

（責任編輯 岳梅）