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

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

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

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

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

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

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

0 引言

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

1 材料與方法

試驗(yàn)于2019年在福建農(nóng)林大學(xué)動(dòng)物科學(xué)學(xué)院（蜂學(xué)學(xué)院）蜜蜂保護(hù)實(shí)驗(yàn)室完成。

1.1 供試真菌

球囊菌菌株由福建農(nóng)林大學(xué)動(dòng)物科學(xué)學(xué)院（蜂學(xué)學(xué)院）蜜蜂保護(hù)實(shí)驗(yàn)室分離和保存[22-24]。

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

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

1.3 cDNA文庫(kù)構(gòu)建及深度測(cè)序

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

1.4 測(cè)序數(shù)據(jù)質(zhì)控

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

1.5 菌絲和孢子miRNA的預(yù)測(cè)與Venn分析

將嚴(yán)格質(zhì)控后得到的small RNA有效標(biāo)簽序列（clean tags）分別比對(duì)到GenBank數(shù)據(jù)庫(kù)（http://www. ncbi.nlm.nih.gov/Web/Genbank/）、Rfam數(shù)據(jù)庫(kù)（http://rfam.xfam.org/）、核糖體數(shù)據(jù)庫(kù)（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，以及比對(duì)上的基因組外顯子區(qū)域、內(nèi)含子區(qū)域和重復(fù)序列區(qū)域的clean tags，將未比對(duì)上的clean tags與miRBase數(shù)據(jù)庫(kù)（www.mirbase.org/）中的miRNA前體序列進(jìn)行比對(duì)，從而鑒定已知miRNA序列。利用Mireap_V 0.2軟件將剩余的未比對(duì)上的clean tags比對(duì)球囊菌參考基因組（assembly AAP 1.0），得到可能的前體序列，根據(jù)clean tags在前體序列的分布信息和前體結(jié)構(gòu)能量信息，采用貝葉斯模型打分預(yù)測(cè)新miRNA。按照TPM=T×106/N（T表示miRNA的tags，N表示總miRNA的tags）算法對(duì)AaM和AaS中每個(gè)miRNA的表達(dá)量進(jìn)行歸一化處理。利用基迪奧在線工具集合（www.omicshare.com）對(duì)AaM和AaS的miRNA進(jìn)行Venn分析，采用默認(rèn)參數(shù)，分別篩選出AaM和AaS的共有miRNA及特有miRNA。

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

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

1.7 DEmiRNA的篩選和靶mRNA的預(yù)測(cè)及分析

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

1.8 DEmiRNA的Stem-loop RT-qPCR驗(yàn)證

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

表1 本研究使用的引物

2 結(jié)果

2.1 測(cè)序數(shù)據(jù)的質(zhì)控與評(píng)估

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

2.2 球囊菌菌絲和孢子miRNA的預(yù)測(cè)及比較分析

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

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

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

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

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

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

表3 AaM中表達(dá)量最高的前10位miRNA

表4 AaS中表達(dá)量最高的前10位miRNA

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

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

從AaM vs AaS中篩選出93個(gè)DEmiRNA，包括65個(gè)上調(diào)miRNA和28個(gè)下調(diào)miRNA；其中上調(diào)幅度最大的是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）；下調(diào)幅度最大的為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能靶向結(jié)合6 090個(gè)mRNA。

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

表5 AaM vs AaS比較組中前10位上調(diào)和下調(diào)miRNA

2.5 DEmiRNA的調(diào)控網(wǎng)絡(luò)構(gòu)建及分析

根據(jù)靶向結(jié)合關(guān)系構(gòu)建調(diào)控網(wǎng)絡(luò)，分析結(jié)果顯示13個(gè)DEmiRNA可靶向結(jié)合131個(gè)mRNA；同一個(gè)DEmiRNA可靶向結(jié)合多個(gè)mRNA，如miR-4968-y可靶向結(jié)合多達(dá)118個(gè)mRNA；同時(shí)，部分mRNA可靶向結(jié)合1—2個(gè)DEmiRNA，如KZZ97168.1分別靶向結(jié)合miR-4968-y和miR-6769-y，KZZ91616靶向結(jié)合miR5782-y（圖3）。

2.6 DEmiRNA的RT-qPCR驗(yàn)證

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

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

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

3 討論

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

圖3 DEmiRNA-mRNA的調(diào)控網(wǎng)絡(luò)

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

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

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

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

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

4 結(jié)論

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

[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] 陳大福, 郭睿, 熊翠玲, 梁勤, 鄭燕珍, 徐細(xì)建, 黃枳腱, 張曌楠, 張璐, 李汶東, 童新宇, 席偉軍. 脅迫意大利蜜蜂幼蟲(chóng)腸道的球囊菌的轉(zhuǎn)錄組分析. 昆蟲(chóng)學(xué)報(bào), 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] 郭睿, 陳大福, 黃枳腱, 梁勤, 熊翠玲, 徐細(xì)建, 鄭燕珍, 張曌楠, 解彥玲, 童新宇, 侯志賢, 江亮亮, 刀晨. 球囊菌脅迫中華蜜蜂幼蟲(chóng)腸道過(guò)程中病原的轉(zhuǎn)錄組學(xué)研究. 微生物學(xué)報(bào), 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] 張曌楠, 熊翠玲, 徐細(xì)建, 黃枳腱, 鄭燕珍, 駱群, 劉敏, 李汶東, 童新宇, 張琦, 梁勤, 郭睿, 陳大福. 蜜蜂球囊菌的參考轉(zhuǎn)錄組組裝及SSR分子標(biāo)記開(kāi)發(fā). 昆蟲(chóng)學(xué)報(bào), 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鑒定及其調(diào)控網(wǎng)絡(luò)分析. 微生物學(xué)報(bào), 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] 郭睿, 杜宇, 熊翠玲, 鄭燕珍, 付中民, 徐國(guó)鈞, 王海朋, 陳華枝, 耿四海, 周丁丁, 石彩云, 趙紅霞, 陳大福. 意大利蜜蜂幼蟲(chóng)腸道發(fā)育過(guò)程中的差異表達(dá)microRNA及其調(diào)控網(wǎng)絡(luò). 中國(guó)農(nóng)業(yè)科學(xué), 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] 杜宇, 范小雪, 蔣海賓, 王杰, 范元嬋, 祝智威, 周丁丁, 萬(wàn)潔琦, 盧家軒, 熊翠玲, 鄭燕珍, 陳大福, 郭睿. 微小RNA及其介導(dǎo)的競(jìng)爭(zhēng)性內(nèi)源RNA調(diào)控網(wǎng)絡(luò)在意大利蜜蜂工蜂中腸發(fā)育過(guò)程中的潛在作用. 中國(guó)農(nóng)業(yè)科學(xué), 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] 郭睿, 杜宇, 周倪紅, 劉思亞, 熊翠玲, 鄭燕珍, 付中民, 徐國(guó)鈞, 王海朋, 耿四海, 周丁丁, 陳大福. 意大利蜜蜂幼蟲(chóng)腸道在球囊菌脅迫后期的差異表達(dá)微小RNA及其靶基因分析. 昆蟲(chóng)學(xué)報(bào), 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] 熊翠玲, 杜宇, 陳大福, 鄭燕珍, 付中民, 王海朋, 耿四海, 陳華枝, 周丁丁, 吳素珍, 石彩云, 郭睿. 意大利蜜蜂幼蟲(chóng)腸道的miRNAs的生物信息學(xué)預(yù)測(cè)及分析. 應(yīng)用昆蟲(chóng)學(xué)報(bào), 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)物及藥理活性研究進(jìn)展. 中國(guó)生化藥物雜志, 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] 李瓊, 崔春來(lái), 宋紅生, 王四寶. mro-miR-33在綠僵菌產(chǎn)孢中的作用. 菌物學(xué)報(bào), 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] 閆思源, 姜學(xué)軍. 細(xì)胞自噬及真菌中自噬研究概述. 菌物學(xué)報(bào), 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

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

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

（責(zé)任編輯 岳梅）