蜜蜂球囊菌菌絲和孢子中長鏈非編碼RNA的比較及潛在功能分析

陳華枝，王杰，祝智威，蔣海賓，范元嬋，范小雪，萬潔琦，盧家軒，鄭燕珍，付中民,2,3，徐國鈞，陳大福,2,3，郭睿,2,3

蜜蜂球囊菌菌絲和孢子中長鏈非編碼RNA的比較及潛在功能分析

陳華枝1，王杰1，祝智威1，蔣海賓1，范元嬋1，范小雪1，萬潔琦1，盧家軒1，鄭燕珍1，付中民1,2,3，徐國鈞1，陳大福1,2,3，郭睿1,2,3

1福建農(nóng)林大學(xué)動物科學(xué)學(xué)院（蜂學(xué)學(xué)院），福州 350002；2福建農(nóng)林大學(xué)蜂療研究所，福州 350002；3福建農(nóng)林大學(xué)蜂產(chǎn)品加工與應(yīng)用教育部工程研究中心，福州 350002

蜜蜂球囊菌（，簡稱球囊菌）是一種專性侵染蜜蜂幼蟲的真菌病原，可引起成年蜜蜂數(shù)量和蜂群群勢的急劇下降。長鏈非編碼RNA（long non-coding RNA，lncRNA）是一類新近發(fā)現(xiàn)的非編碼RNA，在表觀遺傳、細(xì)胞周期、劑量補償?shù)缺姸嗌顒又邪l(fā)揮重要生物學(xué)功能?！尽棵鞔_球囊菌菌絲和孢子中l(wèi)ncRNA的數(shù)量、種類和表達(dá)譜差異，并探究共有l(wèi)ncRNA、特有l(wèi)ncRNA和差異表達(dá)lncRNA（differentially expressed lncRNA，DElncRNA）在菌絲與孢子中的潛在功能。利用基于鏈特異性建庫的lncRNA-seq技術(shù)對球囊菌的純化菌絲（AaM）和純化孢子（AaS）分別進(jìn)行測序。根據(jù)FPKM（Fragment Per Kilobase of per Million mapped reads）法計算lncRNA在AaM和AaS中的表達(dá)水平。通過Venn分析篩選AaM與AaS的共有l(wèi)ncRNA和特有l(wèi)ncRNA。按照≤0.05且|log2fold change|≥1的標(biāo)準(zhǔn)篩選AaM vs AaS比較組中的DElncRNA。通過Blast工具將共有l(wèi)ncRNA、特有l(wèi)ncRNA和DElncRNA的上下游基因比對GO和KEGG數(shù)據(jù)庫，以進(jìn)行功能及通路注釋。根據(jù)靶向結(jié)合關(guān)系構(gòu)建共有l(wèi)ncRNA、特有l(wèi)ncRNA和DElncRNA的競爭性內(nèi)源RNA（competing endogenous RNA，ceRNA）調(diào)控網(wǎng)絡(luò)并利用Cytoscape軟件進(jìn)行可視化。利用RT-qPCR驗證測序數(shù)據(jù)的可靠性。AaM和AaS分別測得108 614 646和105 675 408條原始讀段（raw reads），經(jīng)嚴(yán)格過濾得到107 780 032和104 621 402條有效讀段（clean reads），20分別為98.76%和98.72%，30分別達(dá)到95.84%和95.78%。共鑒定到850個lncRNA。AaM和AaS的共有l(wèi)ncRNA為701個，二者的特有l(wèi)ncRNA分別為39和110個。上述共有l(wèi)ncRNA通過順式作用調(diào)控3 992個上下游基因，它們涉及細(xì)胞進(jìn)程、代謝進(jìn)程和催化活性等42個功能條目，以及代謝途徑、次生代謝物的生物合成和抗生素的生物合成等117條通路；AaM的特有l(wèi)ncRNA和AaS的特有l(wèi)ncRNA通過順式作用分別調(diào)控243和672個上下游基因。AaM vs AaS比較組包含的255個DElncRNA通過順式作用調(diào)控1 479個上下游基因，它們涉及代謝進(jìn)程、細(xì)胞進(jìn)程和催化活性等41個功能條目，以及代謝途徑、次生代謝產(chǎn)物的生物合成和抗生素的生物合成等107條通路。從共有l(wèi)ncRNA、孢子特有l(wèi)ncRNA和DElncRNA中分別預(yù)測到41、5和13個微小RNA（microRNA，miRNA）的前體序列。調(diào)控網(wǎng)絡(luò)分析結(jié)果顯示，菌絲lncRNA、孢子lncRNA形成較為復(fù)雜的ceRNA調(diào)控網(wǎng)絡(luò)；菌絲lncRNA可靶向結(jié)合8個miRNA，進(jìn)而調(diào)控77個mRNA；孢子lncRNA可靶向結(jié)合7個miRNA，進(jìn)而調(diào)控87個mRNA；2個DElncRNA（TCONS_00008630與TCONS_00009302）可靶向結(jié)合miR-4968-y，進(jìn)而調(diào)控10個mRNA。RT-qPCR驗證結(jié)果顯示4個DElncRNA的差異表達(dá)趨勢與測序結(jié)果一致，表明測序數(shù)據(jù)真實可靠。共有l(wèi)ncRNA、特有l(wèi)ncRNA和DElncRNA可能通過調(diào)控上下游基因的表達(dá)，作為miRNA的前體，以及充當(dāng)ceRNA影響菌絲和孢子的物質(zhì)和能量代謝、自噬、轉(zhuǎn)錄、MAPK信號通路、泛素介導(dǎo)的蛋白水解、蛋白酶體以及次生代謝產(chǎn)物的生物合成等生物學(xué)過程，從而調(diào)節(jié)球囊菌的生長、發(fā)育、生殖和致病性。

蜜蜂球囊菌；長鏈非編碼RNA；菌絲；孢子；競爭性內(nèi)源RNA

0 引言

【研究意義】蜜蜂球囊菌（，簡稱球囊菌）是一種能夠侵染蜜蜂幼蟲而導(dǎo)致白堊病的真菌病原[1]。哺育蜂將被球囊菌孢子污染的花粉飼喂幼蟲是導(dǎo)致該病的直接原因[2]。白堊病一般發(fā)生于春季和初夏時期，患病蜂群中的幼蟲會出現(xiàn)成批死亡的癥狀，雖然該病不會導(dǎo)致蜂群的滅亡，但能導(dǎo)致成蜂數(shù)量和蜂群生產(chǎn)力的急劇降低，同時會使其他病原微生物入侵的概率增加，嚴(yán)重影響蜜蜂健康和養(yǎng)蜂生產(chǎn)[3]。然而，由于球囊菌菌絲生長、孢子形成和萌發(fā)的相關(guān)分子機理尚未闡明，蜜蜂白堊病的防治成效有限。利用基于鏈特異性建庫的lncRNA-seq技術(shù)對球囊菌的孢子和菌絲分別進(jìn)行測序和組學(xué)分析，能夠明確球囊菌孢子和菌絲長鏈非編碼RNA（long non-coding RNA，lncRNA）的數(shù)量、種類和表達(dá)譜差異，解析lncRNA介導(dǎo)的菌絲生長、孢子萌發(fā)、雜交產(chǎn)孢和毒力因子合成的生物學(xué)過程，為闡明相關(guān)分子機制提供重要信息和數(shù)據(jù)基礎(chǔ)?！厩叭搜芯窟M(jìn)展】lncRNA指的是一類長度＞200 nt、含有2個及以上外顯子且不具有蛋白編碼能力的非編碼RNA（non-coding RNA，ncRNA），多由RNA聚合酶II轉(zhuǎn)錄形成，具有類似于mRNA的結(jié)構(gòu)如5′帽子和3′ polyA尾巴[4]。LncRNA曾經(jīng)被認(rèn)為是基因轉(zhuǎn)錄過程中的“暗物質(zhì)”和“噪聲”，但目前已被較多的研究證實能夠作為信號[5]、誘餌[6]、引導(dǎo)[7]及支架分子[8]，在表觀遺傳[9]、細(xì)胞周期調(diào)控[10]以及劑量補償[5]等諸多生命活動中發(fā)揮關(guān)鍵的調(diào)控功能。根據(jù)lncRNA基因在基因組上與蛋白編碼基因的相對位置，可將其分為基因間區(qū)lncRNA（intergenic lncRNA）、雙向lncRNA（bidirectional lncRNA）、內(nèi)含子RNA（intronic lncRNA）、反義鏈RNA（antisense lncRNA）和正義鏈RNA（sense lncRNA）[11]。LncRNA具有高度的物種、組織表達(dá)和發(fā)育時期表達(dá)特異性[12]。LncRNA與相鄰蛋白編碼基因的順式（）作用元件或3′UTR區(qū)域結(jié)合，從而通過順式作用在轉(zhuǎn)錄或轉(zhuǎn)錄后水平調(diào)控基因表達(dá)[13-14]。Kim等[15]曾鑒定出547個與禾谷鐮孢（）子實體形成相關(guān)的lncRNA，并發(fā)現(xiàn)其中多數(shù)lncRNA以順式作用調(diào)控mRNA的表達(dá)。LncRNA在釀酒酵母（）和裂殖酵母（）中可通過順式作用調(diào)節(jié)的啟動子，從而影響酵母的生殖細(xì)胞分化[16]。此外，lncRNA可作為微小RNA（microRNA，miRNA）前體，通過一系列加工和剪切形成成熟的miRNA，從而調(diào)控下游靶mRNA的表達(dá)[17-18]。再者，含有miRNA反應(yīng)元件（miRNA response element，MRE）的lncRNA還能作為競爭性內(nèi)源RNA（competing endogenous RNA，ceRNA）靶向結(jié)合miRNA，從而對靶基因的表達(dá)進(jìn)行間接調(diào)控[13,19]。Qin等[20]雖然早在2006年組裝并公布了球囊菌ARSEF7405菌株的基因組序列，但并未同時公布基因功能注釋信息，導(dǎo)致球囊菌的分子生物學(xué)和組學(xué)研究舉步維艱。前人對白堊病的研究主要涉及病原鑒定[21]、流行病學(xué)[22]、形態(tài)學(xué)[23]和快速診斷[24]等方面。2016年，Shang等[25]從頭組裝和注釋了球囊菌ARSEF 7405菌株的參考基因組（AAP 1.0），為深入開展球囊菌的分子生物學(xué)和組學(xué)研究打下了基礎(chǔ)。隨著高通量測序技術(shù)特別是基于鏈特異性建庫的lncRNA-seq技術(shù)的發(fā)展，人們已在動物[26]、植物[27]、細(xì)菌[28]、真菌[29]、病毒[30]中鑒定到大量lncRNA。【本研究切入點】筆者團(tuán)隊前期已利用lncRNA-seq技術(shù)對球囊菌的孢子和菌絲混合樣品進(jìn)行測序，預(yù)測和分析了球囊菌的lncRNA和環(huán)狀RNA（circular RNA，circRNA）的數(shù)量、種類及結(jié)構(gòu)特征，并驗證了二者的真實表達(dá)[31-32]。然而，球囊菌菌絲和孢子lncRNA的數(shù)量、結(jié)構(gòu)、表達(dá)譜差異還沒有明確，lncRNA與菌絲和孢子生長、發(fā)育、生殖和致病性的關(guān)系仍不清楚?！緮M解決的關(guān)鍵問題】利用基于鏈特異性建庫的lncRNA-seq技術(shù)對球囊菌的孢子和菌絲分別進(jìn)行測序，通過生物信息學(xué)方法篩選并深入分析孢子和菌絲的共有l(wèi)ncRNA、特有l(wèi)ncRNA和差異表達(dá)lncRNA（differentially expressed lncRNA，DElncRNA），以期解析球囊菌孢子和菌絲lncRNA的數(shù)量、種類和表達(dá)譜差異，挖掘菌絲生長、孢子萌發(fā)、雜交產(chǎn)孢、毒力因子等相關(guān)的lncRNA并探討其潛在功能。

1 材料與方法

試驗于2019年6月至2020年1月在福建農(nóng)林大學(xué)動物科學(xué)學(xué)院（蜂學(xué)學(xué)院）蜜蜂保護(hù)實驗室完成。

1.1 供試球囊菌

供試球囊菌菌株由福建農(nóng)林大學(xué)動物科學(xué)學(xué)院（蜂學(xué)學(xué)院）蜜蜂保護(hù)實驗室分離和保存[33-34]。

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

參照本課題組前期建立的方法[35]，將實驗室保存的球囊菌菌株接種到馬鈴薯葡萄糖瓊脂（potato dextrose agar，PDA）培養(yǎng)基，移入恒溫生化培養(yǎng)箱33℃條件下培養(yǎng)，培養(yǎng)7 d后分別制備純化的球囊菌菌絲和孢子。將上述菌絲樣品和孢子樣品分別命名為AaM和AaS。

1.3 RNA提取、cDNA鏈特異性文庫構(gòu)建及深度測序

從超低溫冰箱中取出菌絲和孢子樣品，分別置于兩個干凈的研缽，加入液氮充分研磨。利用AxyPrepTMMultisource Total RNA Miniprep試劑盒（TaKaRa公司，日本）分別提取球囊菌菌絲和孢子的總RNA。參照Guo等[31]的方法構(gòu)建鏈特異性cDNA文庫。采用Illumina HiSeqTM4000平臺對上述構(gòu)建好的文庫進(jìn)行雙端測序。測序數(shù)據(jù)已上傳到NCBI SRA數(shù)據(jù)庫，BioProject號：PRJNA560452。

1.4 高通量測序數(shù)據(jù)過濾質(zhì)控與lncRNA預(yù)測

參照GUO等[31]的方法對下機的原始讀段（raw reads）進(jìn)行嚴(yán)格質(zhì)控以保證測序數(shù)據(jù)的質(zhì)量。利用reads比對工具Bowtie2[36]將質(zhì)控得到的有效讀段（clean reads）比對到核糖體數(shù)據(jù)庫（允許錯配率為0），去除比對上的clean reads。再利用Tophat2[37]將剩余的clean reads比對到球囊菌的參考基因組（AAP 1.0）上，得到比對上的有效讀段（mapped clean reads）。使用Cufflinks軟件[38]進(jìn)行轉(zhuǎn)錄本的重構(gòu)，根據(jù)組裝出來的轉(zhuǎn)錄本在參考基因組上的位置，篩選長度≥200 bp且外顯子數(shù)目≥2的轉(zhuǎn)錄本，作為候選轉(zhuǎn)錄本。聯(lián)用CPC[39]和CNCI[40]軟件預(yù)測候選轉(zhuǎn)錄本的編碼能力，同時將轉(zhuǎn)錄本對比到SwissPort蛋白數(shù)據(jù)庫（https://www.uniprot.org），將沒有編碼能力且沒有蛋白注釋信息的轉(zhuǎn)錄本預(yù)測為lncRNA，取二者預(yù)測結(jié)果的交集作為可信的lncRNA。

1.5 LncRNA的Venn分析和差異分析

利用Omicshare在線工具集合（https://www. omicshare.com）對AaM與AaS的lncRNA進(jìn)行Venn分析，篩選出兩組樣品的共有和特有l(wèi)ncRNA。使用edgeR軟件[41]對AaM和AaS的lncRNA進(jìn)行差異分析，根據(jù)≤0.05（FDR矯正）且|log2fold change (FC)|≥1標(biāo)準(zhǔn)篩選DElncRNA。

1.6 共有l(wèi)ncRNA、特有l(wèi)ncRNA和DElncRNA的順式作用分析

LncRNA的順式作用指的是lncRNA的功能與其上下游相鄰的蛋白編碼基因有關(guān)，處于上下游的lncRNA可以與共表達(dá)基因的順式作用元件或者3′UTR存在交集，從而起到對靶基因的調(diào)控作用[42]。利用Blast軟件將共有l(wèi)ncRNA、特有l(wèi)ncRNA和DElncRNA的上下游基因分別映射到GO和KEGG數(shù)據(jù)庫，對映射上的GO條目數(shù)、KEGG通路數(shù)及相應(yīng)的基因數(shù)分別進(jìn)行統(tǒng)計。

1.7 DElncRNA的ceRNA調(diào)控網(wǎng)絡(luò)構(gòu)建分析及miRNA前體預(yù)測

筆者團(tuán)隊前期同時利用small RNA-seq（sRNA- seq）技術(shù)對球囊菌純化菌絲、純化孢子樣品分別進(jìn)行測序，獲得了高質(zhì)量的miRNA組學(xué)數(shù)據(jù)[43]，可為本研究中mRNA、DEmRNA、lncRNA和DElncRNA的靶向預(yù)測及ceRNA調(diào)控網(wǎng)絡(luò)構(gòu)建提供數(shù)據(jù)支持。sRNA-seq的原始數(shù)據(jù)已上傳NCBI SRA數(shù)據(jù)庫，BioProject號：PRJNA560456。

參照本課題組前期建立的方法[44]，利用TargetFind軟件[45]預(yù)測lncRNA與miRNA、miRNA與mRNA的靶向結(jié)合關(guān)系，并構(gòu)建及可視化lncRNA- miRNA-mRNA和DElncRNA-DEmiRNA-DEmRNA調(diào)控網(wǎng)絡(luò)。將預(yù)測出的lncRNA對比到miRBase（version 21），尋找潛在的miRNA前體，選取對其中覆蓋度＞90%的lncRNA，同時利用miRPara（version 6.3）軟件預(yù)測miRNA及其前體。

1.8 DElncRNA實時熒光定量PCR（RT-qPCR）驗證

隨機選取4個DElncRNA（TCONS_00006988、TCONS_00003707、TCONS_00001814、TCONS_00007359）進(jìn)行RT-qPCR驗證。參照郭睿等[26]方法設(shè)計合成DElncRNA的特異性上游引物和下游引物。選擇（5417）作為內(nèi)參基因。相關(guān)的引物信息詳見表1。利用RNA提取試劑盒（TaKaRa公司，日本）分別提取球囊菌菌絲和孢子的總RNA，使用cDNA第1鏈合成試劑盒（上海翊圣，中國）對提取的RNA進(jìn)行反轉(zhuǎn)錄，得到相應(yīng)的cDNA作為qPCR的模板。qPCR反應(yīng)體系為20 μL：SYBR Green Dye 10 μL，上下游引物各1 μL（10 μmol·L-1），cDNA模板1 μL，DEPC水7 μL。配制的反應(yīng)體系在ABI QuanStudio 3熒光定量PCR儀（ABI公司，美國）中進(jìn)行反應(yīng)。qPCR程序：94℃預(yù)變性5 min；94℃變性30 s，53℃退火40 s，72℃延伸30 s，共40個循環(huán)。每組qPCR反應(yīng)設(shè)置3次技術(shù)重復(fù)。采用2-??Ct方法計算DElncRNA的相對表達(dá)量。利用GraphPad Prism 7軟件（GraphPad公司，美國）進(jìn)行數(shù)據(jù)的Student’s-test和相關(guān)繪圖。

表1 本研究所用的引物

2 結(jié)果

2.1 LncRNA-seq數(shù)據(jù)的質(zhì)控與評估

AaM和AaS通過測序得到的raw reads分別為108 614 646和105 675 408條，質(zhì)控后得到的clean reads分別為107 780 032與104 621 402條；20分別為98.76%和98.72%，30分別達(dá)到95.84%和95.78%（表2）。上述結(jié)果表明測序數(shù)據(jù)質(zhì)量良好，可用于進(jìn)一步分析。

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

2.2 球囊菌菌絲與孢子中l(wèi)ncRNA的數(shù)量和種類

基于AaM的測序數(shù)據(jù)預(yù)測出740個lncRNA，其中基因間區(qū)lncRNA、雙向lncRNA、反義鏈lncRNA和正義鏈lncRNA分別為302、199、66和31個；基于AaS的測序數(shù)據(jù)預(yù)測出811個lncRNA，包括324個基因間區(qū)lncRNA，230個雙向lncRNA，68個反義鏈lncRNA和36個正義鏈lncRNA。Venn分析結(jié)果顯示AaM和AaS包含701個共有l(wèi)ncRNA，特有l(wèi)ncRNA的數(shù)量分別為39和110個。

2.3 球囊菌菌絲與孢子共有和特有l(wèi)ncRNA的篩選及上下游基因的數(shù)據(jù)庫注釋

GO數(shù)據(jù)庫注釋結(jié)果顯示，701個共有l(wèi)ncRNA可以通過順式作用潛在調(diào)控3 992個上下游基因，這些上下游基因可注釋到細(xì)胞進(jìn)程（542）、代謝進(jìn)程（539）和催化活性（477）等42個GO功能條目（圖1-A）；AaM的特有l(wèi)ncRNA的上下游基因可注釋到細(xì)胞進(jìn)程（49）、代謝進(jìn)程（41）和催化活性（39）等32個功能條目（圖1-B）；AaS的特有l(wèi)ncRNA的上下游基因可注釋到代謝進(jìn)程（130）、細(xì)胞進(jìn)程（129）和催化活性（117）等36個功能條目（圖1-C）。括號內(nèi)的數(shù)字代表注釋到該條目（或通路）的上下游基因數(shù)。

KEGG數(shù)據(jù)庫注釋結(jié)果顯示，共有l(wèi)ncRNA的上下游基因還可注釋到代謝途徑（260）、次生代謝物的生物合成（122）及抗生素的生物合成（92）等117條通路（圖2-A）；AaM特有l(wèi)ncRNA的上下游基因能注釋到代謝途徑（21）、次生代謝產(chǎn)物的生物合成（8）及嘌呤代謝（6）等56條通路（圖2-B）；AaS特有l(wèi)ncRNA的上下游基因能注釋到代謝途徑（73）、次生代謝產(chǎn)物的生物合成（26）及抗生素的生物合成（20）等95條通路（圖2-C）。

2.4 AaM vs AaS比較組中DElncRNA的篩選及功能和通路注釋

差異分析結(jié)果顯示，AaM vs AaS比較組包含255個DElncRNA，其中上調(diào)和下調(diào)lncRNA的數(shù)量分別為181和74個；上調(diào)幅度最大的前3位分別是TCONS_00004900（log2fc=13.975）、TCONS_00006709（log2fc=13.785）和TCONS_00001694（log2fc= 13.258）；下調(diào)幅度最大的前3位分別是TCONS_ 00004890（log2fc=-13.374）、TCONS_00005341（log2fc=-12.481）和TCONS_00008852（log2fc=-12.388）。

255個DElncRNA共預(yù)測出1 479個上下游基因，包括773個上游基因和706個下游基因。這些上下游基因可注釋到41個功能條目，包括細(xì)胞（203）和細(xì)胞組件（203）等12個細(xì)胞組分相關(guān)條目；催化活性（248）和結(jié)合（171）等10個分子功能相關(guān)條目；代謝進(jìn)程（275）和細(xì)胞進(jìn)程（270）等19個生物學(xué)進(jìn)程相關(guān)條目；此外，還可注釋到代謝途徑（137）及次生代謝產(chǎn)物的生物合成（61）等107條通路。括號內(nèi)的數(shù)字代表注釋到該條目（或通路）的上下游基因數(shù)。

2.5 共有l(wèi)ncRNA、特有l(wèi)ncRNA和DElncRNA作為miRNA前體的分析

從菌絲和孢子的22個共有l(wèi)ncRNA中預(yù)測出41個miRNA前體序列；從孢子的1個特有l(wèi)ncRNA中預(yù)測出5個miRNA前體序列；從菌絲的特有l(wèi)ncRNA中未能預(yù)測出miRNA前體序列；從5個DElncRNA中預(yù)測出13個miRNA前體序列；進(jìn)一步分析發(fā)現(xiàn)這些miRNA前體序列均含有典型的頸環(huán)結(jié)構(gòu)，并且包含成熟miRNA的序列。圖3展示了DElncRNA中預(yù)測出的5個miRNA前體序列及相應(yīng)的成熟miRNA序列。

2.6 菌絲lncRNA、孢子lncRNA和DElncRNA的調(diào)控網(wǎng)絡(luò)構(gòu)建及分析

利用Cytoscape軟件對菌絲全部lncRNA、孢子全部lncRNA的調(diào)控網(wǎng)絡(luò)進(jìn)行可視化，分析結(jié)果顯示菌絲的10個lncRNA靶向結(jié)合8個miRNA，上述miRNA靶向結(jié)合77個mRNA（圖4-A）；孢子的8個lncRNA靶向結(jié)合7個miRNA，這些miRNA靶向結(jié)合87個mRNA（圖4-B）；LncRNA和mRNA位于調(diào)控網(wǎng)絡(luò)的外周，而miRNA處于調(diào)控網(wǎng)絡(luò)的核心。

A：AaM與AaS共有l(wèi)ncRNA上下游基因的GO數(shù)據(jù)庫注釋GO database annotation of upstream and downstream genes of common lncRNAs in AaM and AaS；B：AaM特有l(wèi)ncRNA的GO數(shù)據(jù)庫注釋GO database annotation of upstream and downstream genes of specific lncRNAs in AaM；C：AaS特有l(wèi)ncRNA的GO數(shù)據(jù)庫注釋GO database annotation of upstream and downstream genes of specific lncRNAs in AaS

A：AaM與AaS共有l(wèi)ncRNA上下游基因的KEGG富集分析KEGG database annotation of upstream and downstream genes of common lncRNAs in AaM and AaS；B：AaM特有l(wèi)ncRNA的KEGG富集分析KEGG database annotation of upstream and downstream genes of specific lncRNAs in AaM；C：AaS特有l(wèi)ncRNA的KEGG富集分析KEGG database annotation of upstream and downstream genes of specific lncRNAs in AaS

進(jìn)一步對DElncRNA與DEmiRNA、DEmRNA的調(diào)控網(wǎng)絡(luò)進(jìn)行構(gòu)建和分析，結(jié)果顯示DElncRNA TCONS_00008630與TCONS_00009302共同靶向miR-4968-y，進(jìn)而調(diào)控10個DEmRNA（圖4-C）。上述10個DEmRNA可注釋到11個功能條目，包括細(xì)胞（2）、細(xì)胞組分（2）、細(xì)胞器（2）、細(xì)胞器組分（1）、腔上包膜（1）、結(jié)合（3）、催化活性（2）、代謝進(jìn)程（3）、細(xì)胞進(jìn)程（3）、單細(xì)胞進(jìn)程（1）、生物調(diào)節(jié)（1）；此外還能注釋到11條通路，包括代謝途徑（2）、次生代謝物的生物合成（1）、抗菌活性（1）、不同環(huán)境中的微生物代謝（1）、氨基酸生物合成（1）、嘌呤代謝（1）、嘧啶代謝（1）、丙氨酸、天門冬氨酸和谷氨酸代謝（1）、RNA聚合酶（1）、堿基切除修復(fù)（1）、氮代謝（1）。括號內(nèi)的數(shù)字代表注釋到該條目（或通路）的上下游基因數(shù)。

綠色堿基代表成熟miRNA的序列Green bases indicate sequences of mature miRNAs

2.7 DElncRNA 的RT-qPCR驗證

隨機挑選4個DElncRNA進(jìn)行RT-qPCR驗證，結(jié)果表明DElncRNA的變化趨勢均與測序結(jié)果一致（圖5），證明本研究的測序數(shù)據(jù)和lncRNA的差異表達(dá)真實可靠。

3 討論

近年來，較多的研究結(jié)果證實lncRNA在真核生物的生理、代謝、免疫和疾病等過程扮演重要角色[46-47]。本研究首先在實驗室條件下對球囊菌進(jìn)行純培養(yǎng)，然后分離得到純化菌絲樣品和純化孢子樣品，進(jìn)而利用基于鏈特異性建庫的lncRNA-seq技術(shù)對二者分別進(jìn)行測序，通過生物信息學(xué)方法從菌絲和孢子中分別預(yù)測出740和811個lncRNA，其中701個lncRNA為菌絲和孢子所共有，39和110個lncRNA在二者中特異性表達(dá)。鑒于lncRNA表達(dá)具有物種、組織、發(fā)育時期和脅迫階段特異性[48]，上述共有l(wèi)ncRNA可能在球囊菌的不同存在形態(tài)中發(fā)揮一般性的調(diào)控功能，而特有l(wèi)ncRNA在菌絲和孢子中分別發(fā)揮不同的調(diào)控功能。本研究中，球囊菌菌絲和孢子處于離體環(huán)境，推測仍有可能存在一些lncRNA在球囊菌侵染蜜蜂幼蟲的過程中特異性表達(dá)。菌絲和孢子中數(shù)量最多的lncRNA類型均為基因間區(qū)lncRNA，而內(nèi)含子lncRNA的數(shù)量很少。這與東方蜜蜂微孢子蟲（）、羅伯茨綠僵菌（）、釀酒酵母[17,49]等真菌的lncRNA主要類型相似。對于菌絲和孢子共有l(wèi)ncRNA的上下游基因，分別有542、539和477個涉及細(xì)胞進(jìn)程、代謝進(jìn)程和催化活性；74和35個涉及應(yīng)激反應(yīng)和信號；14與13個涉及生殖和生殖進(jìn)程。這表明共有l(wèi)ncRNA可能通過順式作用調(diào)控上下游基因的表達(dá)，從而參與球囊菌菌絲和孢子的細(xì)胞生命活動、新陳代謝、信號轉(zhuǎn)導(dǎo)和有性生殖等生物學(xué)過程。此外，共有l(wèi)ncRNA的上下游基因還能注釋到117條通路，其中有133個上下游基因注釋到酪氨酸代謝（9）等14條氨基酸代謝相關(guān)通路；151個注釋到糖酵解/糖異生通路（24）等15條碳水化合物代謝相關(guān)通路；78個注釋到甘油磷脂代謝（17）等13條脂質(zhì)代謝相關(guān)通路；62個注釋到嘌呤代謝（34）和嘧啶代謝（28）；還有43個上下游基因可注釋到4條能量代謝相關(guān)通路，包括氧化磷酸化（22）、甲烷代謝（10）、硫代謝（7）及氮代謝（4）。上述結(jié)果再次說明共有l(wèi)ncRNA通過順式作用廣泛參與球囊菌菌絲和孢子的物質(zhì)代謝和能量代謝方面的調(diào)控。

圖4 菌絲lncRNA（A）、孢子lncRNA（B）和DElncRNA（C）的調(diào)控網(wǎng)絡(luò)

*: p＜0.05; **: p＜0.01

3.1 球囊菌菌絲和孢子的共有l(wèi)ncRNA與自噬通路具有潛在的調(diào)控關(guān)系

自噬是真核生物非常重要且高度保守的蛋白質(zhì)降解過程，通過自噬體將細(xì)胞中的細(xì)胞器、蛋白質(zhì)及其他生物大分子包裹后運送到具有降解作用的細(xì)胞器進(jìn)行降解并重新利用，從而使得細(xì)胞在營養(yǎng)缺乏的調(diào)節(jié)下維持基本生存[50]。自噬已被證明與真菌的孢子產(chǎn)生存在直接關(guān)聯(lián)[51]。稻瘟病菌（）是嚴(yán)重危害水稻的真菌病原，Liu等[52]在稻瘟病菌中鑒定到與自噬相關(guān)的基因，并發(fā)現(xiàn)該基因與稻瘟病菌的分生孢子形成直接相關(guān)；Deng等[53]發(fā)現(xiàn)缺失的稻瘟病菌菌株產(chǎn)孢能力顯著降低。球囊菌通過異宗配合進(jìn)行有性生殖，PDA培養(yǎng)基上生長早期的球囊菌菌絲為白色，當(dāng)異宗菌絲接觸發(fā)生配合后才會形成內(nèi)含大量孢子的黑色孢子囊。本研究中，30個共有l(wèi)ncRNA（TCONS_00008713、TCONS_00008558和TCONS_00003019等）的25個上下游基因（KZZ86589.1、KZZ86859.1和KZZ86861.1等）可注釋到自噬通路，說明這些共有l(wèi)ncRNA可能通過順式作用對自噬相關(guān)的上下游基因進(jìn)行表達(dá)調(diào)控，從而影響球囊菌的雜交產(chǎn)孢。

3.2 球囊菌菌絲特有l(wèi)ncRNA可能參與病原對蜜蜂幼蟲的侵染過程

昆蟲中腸內(nèi)側(cè)有一層致密的圍食膜，主要由幾丁質(zhì)、多糖和蛋白質(zhì)等成分組成，能夠?qū)⑸掀ぜ?xì)胞與食物隔開，構(gòu)成昆蟲抵御病原微生物入侵的一道物理屏障[54]。在被感染的蜜蜂幼蟲腸道內(nèi)，球囊菌菌絲大量生長的過程一方面通過機械作用刺穿圍食膜和腸壁，同時合成和分泌一些幾丁質(zhì)酶、水解酶、脂酶和蛋白酶等毒力因子共同作用分解宿主中腸圍食膜，協(xié)同促進(jìn)病原的增殖和侵染[3]。泛素是一類存在于多數(shù)真核細(xì)胞中的高度保守的蛋白質(zhì)，對于蛋白合成與降解至關(guān)重要[55]。本研究中，菌絲的特有l(wèi)ncRNA的2個上下游基因（KZZ91868.1和KZZ97813.1）可注釋到泛素介導(dǎo)的蛋白水解通路，另有1個上下游基因（KZZ95233.1）可注釋到蛋白酶體通路，表明調(diào)控KZZ91868.1、KZZ97813.1和KZZ95233.1的菌絲特有l(wèi)ncRNA（TCONS_00004660、TCONS_00000056和TCONS_00002376）參與調(diào)節(jié)病原侵染宿主過程對圍食膜的分解。真菌侵染昆蟲的過程中會通過產(chǎn)生一些次生代謝產(chǎn)物（例如毒素）促進(jìn)自身侵染。卵孢白僵菌（）[56]在侵染血黑蝗（）過程中分泌的草酸能夠溶解宿主體壁；金龜子綠僵菌（）[57]在侵染大蠟螟（）過程中產(chǎn)生的細(xì)胞松弛素可降低血淋巴的吞噬能力，有助于其自身在宿主體內(nèi)的存活。本研究中，4個菌絲特有l(wèi)ncRNA調(diào)控的8個上下游基因可注釋到次生代謝產(chǎn)物的生物合成，包括KZZ86645.1（核糖磷酸焦磷酸激酶1編碼基因）和KZZ87421.1（2-異丙基蘋果酸合成酶編碼基因）等，說明這些特有l(wèi)ncRNA可能參與球囊菌侵染過程毒力因子的合成與分泌。但需要強調(diào)的是，本研究測序材料來源于實驗室條件下獲得的球囊菌純培養(yǎng)，其表達(dá)的lncRNA必然與處于侵染過程的球囊菌所表達(dá)的lncRNA存在差異，因此上述結(jié)果顯示出的是部分特有l(wèi)ncRNA與球囊菌致病力的潛在關(guān)聯(lián)，若要明確二者的直接關(guān)系，需要通過對處于侵染過程的球囊菌lncRNA組學(xué)數(shù)據(jù)、純化孢子的lncRNA組學(xué)數(shù)據(jù)進(jìn)行比較分析，從而篩選和挖掘DElncRNA的相關(guān)信息。目前，筆者團(tuán)隊已利用基于鏈特異性建庫的lncRNA-seq技術(shù)對正常及球囊菌脅迫的意蜂幼蟲腸道、中蜂幼蟲腸道進(jìn)行深度測序，獲得了高質(zhì)量的lncRNA組學(xué)數(shù)據(jù)（未發(fā)表數(shù)據(jù)），下一步將通過比對宿主參考基因組過濾掉來源于宿主的lncRNA組學(xué)數(shù)據(jù)，再將剩余數(shù)據(jù)比對球囊菌參考基因組，從而獲得來源于球囊菌本身的lncRNA組學(xué)數(shù)據(jù)，進(jìn)而結(jié)合本研究中球囊菌純化孢子的lncRNA組學(xué)數(shù)據(jù)進(jìn)行細(xì)致深入的比較分析。

3.3 球囊菌在孢子狀態(tài)仍具有l(wèi)ncRNA轉(zhuǎn)錄及較低水平的代謝活動

真菌孢子是一種休眠態(tài)的存在形式，相對于營養(yǎng)狀態(tài)的菌絲，人們對于孢子中核酸的轉(zhuǎn)錄、翻譯和代謝方面的認(rèn)識較為有限。球囊菌的孢子對不良環(huán)境有很強的抵抗力，可以存活并保持感染力超過15年以上[58]。本研究發(fā)現(xiàn)，僅在孢子中特異性表達(dá)的110個lncRNA可通過順式作用調(diào)控672個上下游基因，其中有300個上下游基因可注釋到物質(zhì)和能量代謝相關(guān)通路，包括絡(luò)氨酸代謝（6）等12條氨基酸代謝通路，氨基酸糖和核苷酸糖代謝（7）等15條碳代謝通路，氧化磷酸化（4）等4條能量代謝通路，以及脂肪酸降解（2）等9條脂質(zhì)代謝通路。上述結(jié)果暗示球囊菌在其孢子狀態(tài)仍進(jìn)行著較低水平的物質(zhì)和能量代謝活動，部分lncRNA參與了上述代謝活動的調(diào)控。東方蜜蜂微孢子蟲是另一種廣泛感染世界各地蜂群的蜜蜂真菌病原。筆者團(tuán)隊前期研究發(fā)現(xiàn)東方蜜蜂微孢子蟲的孢子中同樣存在一定水平的轉(zhuǎn)錄活動[29]。這與本研究的結(jié)果相似，表明蜜蜂真菌病原的孢子內(nèi)部同樣存在低水平的轉(zhuǎn)錄及新陳代謝活動，推測這對于維持孢子的完整性和活力是必要的。

3.4 球囊菌菌絲和孢子的DElncRNA具有通過調(diào)控MAPK信號通路影響病原致病性的潛在功能

本研究中，共鑒定到181個上調(diào)lncRNA和74個下調(diào)lncRNA。這些DElncRNA的上下游基因可注釋到41個功能條目，其中分別有3、10、21和36個上下游基因注釋到生長、發(fā)育進(jìn)程、生殖和應(yīng)激反應(yīng)，說明DElncRNA參與上述生物學(xué)過程的調(diào)控。真菌通過信號通路感受外界環(huán)境變化進(jìn)而調(diào)整其生理狀態(tài)[59]；此外，信號通路與真菌的繁殖、生長和毒力密切相關(guān)[60]。真菌中，MAPK信號通路是一條具有較高保守性的關(guān)鍵信號通路，已被證實在釀酒酵母的交配、侵染、增殖、細(xì)胞壁完整性、滲透調(diào)節(jié)和孢子形成等方面起到重要作用[61]。新月彎孢（）是一種侵染玉米的真菌，Ni等[62]通過對該病原的2個MAPK信號通路基因和進(jìn)行研究，發(fā)現(xiàn)與病原孢子形成和致病性有關(guān)，而與細(xì)胞壁形成及病原侵染過程相關(guān)。中華蜜蜂（）和意大利蜜蜂（）是我國養(yǎng)蜂生產(chǎn)中使用的主要蜂種，前者對球囊菌抗性較強，后者則較為易感。筆者團(tuán)隊前期研究發(fā)現(xiàn)，對于侵染意蜂幼蟲的球囊菌，有48個基因富集在MAPK通路且呈上調(diào)表達(dá)趨勢，而侵染中蜂的球囊菌只有11個下調(diào)基因富集于MAPK通路[13]。本研究發(fā)現(xiàn)，11個DElncRNA（TCONS_00000156、TCONS_00006998和TCONS_00008039等）通過順式作用調(diào)控的9個上下游基因注釋到MAPK信號通路，暗示這些DElncRNA可能在球囊菌侵染蜜蜂幼蟲的過程中參與對MAPK信號通路的調(diào)控，從而影響病原的致病性。

3.5 球囊菌菌絲和孢子的DElncRNA可能作為ceRNA調(diào)控物質(zhì)和能量代謝以及遺傳信息傳遞

Salmena等[63]提出了ceRNA假說，該假說認(rèn)為具有miRNA反應(yīng)元件的轉(zhuǎn)錄本，例如lncRNA、circRNA、假基因等均可以通過競爭結(jié)合miRNA減弱miRNA對mRNA的抑制或降解作用。此后，該假說已被較多的研究[19,64]所證實。本研究中構(gòu)建的DElncRNA-DEmiRNA-DEmRNA調(diào)控網(wǎng)絡(luò)并不復(fù)雜，其中僅有2個DElncRNA TCONS_00008630（log2fc= -1.67）和TCONS_00009302（log2fc=-11.48）靶向結(jié)合miR-4968-y（log2fc=3.37），進(jìn)而調(diào)控10個靶mRNA，這些靶mRNA涉及能量代謝（1）、復(fù)制與修復(fù)（1）、轉(zhuǎn)錄（1）、氨基酸代謝（1）等活動。上述結(jié)果表明球囊菌菌絲和孢子的DElncRNA充當(dāng)ceRNA調(diào)控球囊菌的物質(zhì)和能量代謝以及遺傳信息傳遞。

4 結(jié)論

通過對球囊菌菌絲和孢子共有l(wèi)ncRNA、特有l(wèi)ncRNA和DElncRNA的全面分析和深入探討，解析了菌絲和孢子中l(wèi)ncRNA的數(shù)量、種類和表達(dá)譜差異，并揭示了部分lncRNA可能通過順式作用和ceRNA機制參與調(diào)控球囊菌菌絲和孢子的生長、發(fā)育和生殖等生物學(xué)過程。研究結(jié)果為深入理解球囊菌的基礎(chǔ)生物學(xué)提供了重要的參考信息，也為深入研究lncRNA及其調(diào)控網(wǎng)絡(luò)介導(dǎo)球囊菌對蜜蜂幼蟲的侵染過程提供了必要的數(shù)據(jù)基礎(chǔ)。

[1] Calderón RA, Rivera G, Sánchez LA, Zamora LG. Chalkbrood () and some other fungi associated with Africanized honey bees () in Costa Rica. Journal of Apicultural Research, 2004, 43(4): 187-188.

[2] WOOD M. Microbes help bees battle chalkbrood. Agricultural Research, 1998, 46(8): 16-17.

[3] ARONSTEINK K A, MURRAY D. Chalkbrood disease in honey bees.Journal of Invertebrate Pathology, 2010, 103(Suppl.1): 20-29.

[4] KAPRANOV P, CHENG J, DIKE S, NIX D A, Duttagupta R, Willingham A T, Stadler P F, HERTEL J, HACKERMüLLER J, HOFACKER I L,. RNA maps reveal new RNA classes and a possible function for pervasive transcription. Science, 2007, 316(5830): 1484-1488.

[5] PONTIER D B, GRIBNAU J. Xist regulation and function explored. Human Genetics, 2011, 130(2): 223-236.

[6] KINO T, HURT D E, ICHIJO T, NADER N, CHROUSOS G P. Noncoding RNA gas5 is a growth arrest- and starvation-associated repressor of the glucocorticoid receptor. Science Signaling, 2010, 3(107): ra8.

[7] HEO J B, SUNG S. Vernalization-mediated epigenetic silencing by a long intronic noncoding RNA.Science, 2011, 331(6013): 76-79.

[8] TSAI M C, MANOR O, WAN Y, MOSAMMAPARAST N, WANG J K, LAN F, SHI Y, SEGAL E, CHANG H Y. Long noncoding RNA as modular scaffold of histone modification complexes. Science, 2010, 329(5992): 689-693.

[9] MERCER T R, MATTICK J S. Structure and function of long noncoding RNAs in epigenetic regulation.Nature Structural and Molecular Biology, 2013, 20(3): 300-307.

[10] LI M, SUN X, CAI H, SUN Y, PLATH M, LI C, LAN X, LEI C, LIN F, BAI Y, CHEN H. Long non-coding RNA adncr suppresses adipogenic differentiation by targeting miR-204. Biochimica et Biophysica Acta, 2016, 1859(7): 871-882.

[11] ST LAURENT G, WAHLESTEDT C, KAPRANOV P. The landscape of long noncoding RNA classification. Trends in Genetics, 2015, 31(5): 239-251.

[12] LI Z X, HAN K W, ZHANG D F, CHEN J G, XU Z, HOU L J. The role of long noncoding RNA in traumatic brain injury. Neuropsychiatric Disease and Treatment, 2019, 15: 1671-1677.

[13] 郭睿, 王海朋, 陳華枝, 熊翠玲, 鄭燕珍, 付中民, 趙紅霞, 陳大福. 蜜蜂球囊菌的microRNA鑒定及其調(diào)控網(wǎng)絡(luò)分析. 微生物學(xué)報, 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. Acta Microbiologica Sinica, 2018, 58(6): 1077-1089. (in chinese)

[14] KOPP F, MENDELL J T. Functional classification and experimental dissection of long noncoding RNAs. Cell, 2018, 172(3): 393-407.

[15] KIM W, Miguel-Rojas C, Wang J, TownsendJ P, Trail F. Developmental dynamics of long noncoding RNA expression during sexual fruiting body formation in. mbio, 2018, 9(4): e01292-18.

[16] HiriartE, VerdelA. Long noncoding RNA-based chromatin control of germ cell differentiation: a yeast perspective. Chromosome Research, 2013, 21(6/7): 653-663.

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

[18] WILUSZ J E, SUNWOO H, SPECTOR D L. Long noncoding RNAs: functional surprises from the RNA world. Genes and Development, 2009, 23(13): 1494-1504.

[19] FAN G Q, WANG Z, ZHAI X Q, CAO Y B. ceRNA cross-talk in paulownia witches’ broom disease. International Journal of Molecular Sciences, 2018, 19(8): 2463.

[20] QIN X, EVANS J D, ARONSTEIN K A, MURRAY K D, WEINSTOCK G M. Genome sequences of the honey bee pathogensand. Insect Molecular Biology, 2006, 15(5): 715-718.

[21] SPILTOIR C F. Life cycle of(). American Journal of Botany, 1955, 42(6): 501-508.

[22] FLORES J M, SPIVAK M, GUTIéRREZ I. Spores ofcontained in wax foundation can infect honeybee brood. Veterinary Microbiology, 2005, 108(1/2): 141-144.

[23] SKOU J P. More details in support of the classAscosphaeromycetes. Mycotaxon, 1988, 31(1): 191-198.

[24] Anderson D, Giacon H Search articles by 'Giacon H' Giacon H, Gibson N Search articles by 'Gibson N' Gibson N. Detection and thermal destruction of the chalkbrood fungus () in honey. Journal of Apicultural Research, 1997, 36(3/4): 163-168.

[25] SHANG Y F, XIAO G H, ZHENG P, CEN K, ZHAN S, WANG C S. Divergent and convergent evolution of fungal pathogenicity. Genome Biology and Evolution, 2016, 8(5): 1374-1387.

[26] 郭睿, 耿四海, 熊翠玲, 鄭燕珍, 付中民, 王海朋, 杜宇, 童新宇, 趙紅霞, 陳大福. 意大利蜜蜂工蜂中腸發(fā)育過程中長鏈非編碼RNA的差異表達(dá)分析. 中國農(nóng)業(yè)科學(xué), 2018, 51(18): 3600-3613.

GUO R, GENG S H, XIONG C L, ZHENG Y Z, FU Z M, WANG H P, DU Y, TONG X Y, ZHAO H X, CHEN D F. Differential expression analysis of long non-coding RNAs during the developmental process ofworker’s midgut. Scientia Agricultura Sinica, 2018, 51(18): 3600-3613. (in chinese)

[27] Flórez-ZapataN M V, Reyes-ValdésM H, MartínezO. Long non-coding RNAs are major contributors to transcriptome changes in sunflower meiocytes with different recombination rates. BMC Genomics, 2016, 17: 490.

[28] YANG L, LONG Y P, LI C, CAO L, GAN H Y, HUANG K L, JIA Y J. Genome-wide analysis of long noncoding RNA profile in human gastric epithelial cell response to. Japanese Journal of Infectious Diseases, 2015, 68(1): 63-66.

[29] 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. Apidologie, 2018, 49(5): 660-670.

[30] WANG Z Q, ZHAO Y W, ZHANG Y. Viral lncRNA: A regulatory molecule for controlling virus life cycle. Noncoding RNA Research, 2017, 2(1): 38-44.

[31] 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, DHIRAJ K, LIANG Q. Identification of long non-coding RNAs in the chalkbrood disease pathogen. Journal of Invertebrate Pathology, 2018, 156: 1-5.

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

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

[34] 張曌楠, 熊翠玲, 徐細(xì)建, 黃枳腱, 鄭燕珍, 駱群, 劉敏, 李汶東, 童新宇, 張琦, 梁勤, 郭睿, 陳大福. 蜜蜂球囊菌的參考轉(zhuǎn)錄組組裝及SSR分子標(biāo)記開發(fā). 昆蟲學(xué)報, 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. Acta Entomologica Sinica, 2017, 60(1): 34-44. (in chinese)

[35] 陳大福, 郭睿, 熊翠玲, 梁勤, 鄭燕珍, 徐細(xì)建, 張曌楠, 黃枳腱, 張璐, 王鴻權(quán), 解彥玲, 童新宇.中華蜜蜂幼蟲腸道響應(yīng)球囊菌早期脅迫的轉(zhuǎn)錄組學(xué). 中國農(nóng)業(yè)科學(xué), 2017, 50(13): 2614-2623.

CHEN D F, GUO R, XIONG C L, LIANG Q, ZHENG Y Z, XU X J, ZHANG Z N, HUANG Z J, ZHANG L, WANG H Q, XIE Y N, TONG X Y. Transcriptome oflarval gut under the stress of. Scientia Agricultura Sinica, 2017, 50(13): 2614-2623. (in chinese)

[36] LANGMEAD B, TRAPNELL C, POP M, SALZBERG S L. Ultrafast and memory-efficient alignment of short DNA sequences to the human genome. Genome Biology, 2009, 10(3): R25.

[37] KIM D, PERTEA G, TRAPNELL C, PIMENTEL H, KELLEY R, SALZBERG S L. Tophat2: accurate alignment of transcriptomes in the presence of insertions, deletions and gene fusions. Genome Biology, 2013, 14(4): R36.

[38] TRAPNELL C, ROBERTS A, GOFF L, PERTEA G, KIM D, KELLEY D R, PIMENTEL H, SALZBERG S L, RINN J L, PACHTER L. Differential gene and transcript expression analysis of RNA-seq experiments with TopHat and Cufflinks. Nature Protocols, 2012, 7(3): 562-578.

[39] KONG L, ZHANG Y, YE Z Q, LIU X Q, ZHAO S Q, WEI L P, GAO G. CPC: assess the protein-coding potential of transcripts using sequence features and support vector machine. Nucleic Acids Research, 2007, 35: W345-349.

[40] SUN L, LUO H T, BU D C, ZHAO G G, YU K T, ZHANG C H, LIU Y N, CHEN R S, ZHAO Y. Utilizing sequence intrinsic composition to classify protein-coding and long non-coding transcripts. Nucleic Acids Research, 2013, 41(17): e166.

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

[42] SHAO J, CHEN H, YANG D, JIANG M, ZHANG H, WU B, LI J, YUAN L, LIU C. Genome-wide identification and characterization of natural antisense transcripts by strand-specific RNA sequencing in. Scientific Reports, 2017, 7(1): 5711.

[43] 陳華枝, 祝智威, 蔣海賓, 王杰, 范元嬋, 范小雪, 萬潔琦, 盧家軒, 熊翠玲, 鄭燕珍, 付中民, 陳大福, 郭睿. 蜜蜂球囊菌菌絲和孢子中微小RNA及其靶mRNA的比較分析. 中國農(nóng)業(yè)科學(xué), 2020, 53(17): 3606-3619.

CHEN H Z, ZHU Z W, JIANG H B, WANG J, FAN Y C, FAN X X, WAN J Q, LU J X, XIONG C L, ZHENG Y Z, FU Z M, CHEN D F, GUO R. Comparative analysis of microRNAs and corresponding target mRNAs inmycelium and spore. Scientia Agricultura Sinica, 2020, 53(17): 3606-3619. (in Chinese)

[44] 熊翠玲, 陳華枝, 陳大福, 鄭燕珍, 付中民, 徐國鈞, 杜宇, 王海朋, 耿四海, 周丁丁, 劉思亞, 郭睿. 意大利蜜蜂工蜂中腸的環(huán)狀RNA及其調(diào)控網(wǎng)絡(luò)分析. 昆蟲學(xué)報, 2018, 61(12): 1363-1375.

XIONG C L, CHEN H Z, CHEN D F, ZHENG Y Z, FU Z M, XU G J, DU Y, WANG H P, GENG S H, ZHOU D D, LIU S Y, GUO R. Analysis of circular RNAs and their regulatory networks in the midgut ofworkers. Acta Entomologica Sinica, 2018, 61(12): 1363-1375. (in chinese)

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

[46] ATIANAND M K, FITZGERALD K A. Long non-coding RNAs and control of gene expression in the immune system. Trends in Molecular Medicine, 2014, 20(11): 623-631.

[47] QURESHI I A, MEHLER M F. Emerging roles of non-coding RNAs in brain evolution, development, plasticity and disease. Nature Reviews. Neuroscience, 2012, 13(8): 528-541.

[48] Pauli A, Valen E, Lin M F, Garber M, Vastenhouw N L, Levin J Z, Fan L, Sandelin A, Rinn J L, Regev A, Schier A F. Systematic identification of long noncoding RNAs expressed during zebrafish embryogenesis. Genome Research, 2012, 22(3): 577-591.

[49] WANG Z, JIANG Y, WU H, XIE X, HUANG B. Genome-wide identification and functional prediction of long non-coding RNAs involved in the heat stress response in. Frontiers in Microbiology, 2019, 10: 2336.

[50] YORIMITSU T, KLIONSKY D J. Autophagy: molecular machinery for self-eating. Cell Death and Differentiation, 2005, 12(Suppl. 2): 1542-1552.

[51] Fan C L, Chang A N, Liu T B. Role of autophagy in the reproduction of pathogenic fungi. Acta Microbiologica Sinica, 2019, 59(2): 224-234.

[52] LIU X H, ZHAO Y H, ZHU X M, ZENG X Q, HUANG L Y, DONG B, SU Z Z, WANG Y, LU J P, LIN F C. Autophagy-related proteinmoatg14 is involved in differentiation, development and pathogenicity in the rice blast fungus. Scientific Reports, 2017, 7: 40018.

[53] DENG Y Z, RAMOS-PAMPLONA M, NAQVI N I. Autophagy- ssisted glycogen catabolism regulates asexual differentiation in. Autophagy, 2009, 5(1): 33-43.

[54] WANG P, GRANADOS R R. Observations on the presence of the peritrophic membrane in larvaland its role in limiting baculovirus infection. Journal of Invertebrate Pathology, 1998, 72(1): 57-62.

[55] SURESH B, LEE J, KIM K S, RAMAKRISHNA S. The importance of ubiquitination and deubiquitination in cellular reprogramming. Stem Cells International, 2016, 2016:6705927.

[56] Bidochka M J, Khachatourians G G. The implication of metabolic acids produced byin pathogenesis of the migratory grasshopper,. Journal of Invertebrate Pathology, 1991, 58(1): 106-117.

[57] GoTZ P, MATHA V, VILCINSKAS A. Effects of the entomopathogenic fungusand its secondary metabolites on morphology and cytoskeleton of plasmatocytes isolated from the greater wax moth,. Journal of Invertebrate Pathology, 1997, 43(12): 1149-1159.

[58] Jensen A B, Aronstein K, Flores J M, Vojvodic S, Palacio M A, Spivak M. Standard methods for fungal brood disease research. Journal of Apicultural Research, 2013, 52(1): 10.3896/IBRA.1.52.1.13.

[59] Brown A J P, Budge S, Kaloriti D, Tillmann A, Jacobsen M D, Yin Z, Ene I V, Bohovych I, Sandai D, Kastora S, Potrykus J, Ballou E R, Childers D S, Shahana S, Leach M D. Stress adaptation in a pathogenic fungus. Journal of Experimental Biology, 2014, 217(1): 144-155.

[60] So K K, Kim D H. Role of MAPK signaling pathways in regulating the hydrophobin cryparin in the chestnut blight fungusMycobiology, 2017, 45(4): 362-369.

[61] JIANG C, ZHANG X, LIU H Q, XU J R. Mitogen-activated protein kinase signaling in plant pathogenic fungi. PLoS Pathogens, 2018, 14(3): e1006875.

[62] NI X, GAO J X, YU C J, WANG M, SUN J N, LI Y Q, CHEN J. MAPKs and acetyl-CoA are associated withpathogenicity and toxin production in maize.Journal of Integrative Agriculture, 2018, 17(1): 139-148.

[63] SALMENA L, POLISENO L, TAY Y, KATS L, PANDOLFI P P. A ceRNA hypothesis: the rosetta stone of a hidden RNA language? Cell, 2011, 146(3): 353-358.

[64] HUANG M J, ZHAO J Y, XU J J, LI Y, ZHUANG Y F, ZHANG X L. LncRNA ADAMTS9-AS2 controls human mesenchymal stem cell chondrogenic differentiation and functions as a ceRNA. Molecular Therapy-Nucleic Acids, 2019, 18: 533-545.

Comparison and potential functional analysis of long non-coding RNAs betweenmycelium and spore

Chen HuaZhi1, Wang Jie1, Zhu ZhiWei1, Jiang HaiBin1, Fan YuanChan1, Fan XiaoXue1, Wan JieQi1, Lu JiaXuan1, Zheng YanZhen1, Fu ZhongMin1,2,3, Xu GuoJun1, Chen DaFu1,2,3, Guo Rui1,2,3

1College of Animal Sciences (College of Bee Science), Fujian Agriculture and Forestry University, Fuzhou 350002;2Apitherapy Research Institution, Fujian Agriculture and Forestry University, Fuzhou 350002;3Engineering Research Center of Processing and Application of Bee Products of Ministry of Education, Fujian Agriculture and Forestry University, Fuzhou 350002

【】is a fungal pathogen that exclusively infects honeybee larvae, causing a sharp decrease of the population of adult honeybees and colony population. Long non-coding RNA (lncRNA), a recently discovered non-coding RNA (ncRNA), plays a vital biological role in various activities such as epigenetics, cell cycle and dose compensation.【】This study aimed to clarify the differences of number, type and expression profile of lncRNAs betweenmycelium and spore, and investigate the potential role of the common lncRNAs, specific lncRNAs and differentially expressed lncRNAs (DElncRNAs). 【】The purified mycelia (AaM) and purified spores (AaS) ofwere respectively sequenced using strand specific library-based lncRNA-seq technology. The expression levels of lncRNAs in AaM and AaS were calculated using FPKM (Fragment Per Kilobase of per Million mapped reads) method. Common lncRNAs and specific lncRNAs were filtered out following Venn analysis. DElncRNAs within AaM vs AaS comparison group were screened out following the standard of≤0.05 and |log2fold change|≥1. Upstream and downstream genes of common lncRNAs, specific lncRNAs and DElncRNAs were aligned against GO and KEGG databases to obtain function and pathway annotations. The competing endogenous RNA (ceRNA) regulation networks of common lncRNAs, specific lncRNAs and DElncRNAs were constructed following target binding relationships, followed by visualization using Cytoscape software. RT-qPCR was performed to verify the reliability of the sequencing data.【】In total, 108 614 646 and 105 675 408 raw reads were gained from AaM and AaS, and after strict filtering, 107 780 032 and 104 621 402 clean reads were obtained, respectively, with20 of 98.76% and 98.72%, and30 of 95.84% and 95.78%. A total of 850 lncRNAs were identified. Seven hundred and one lncRNAs were shared by AaM and AaS, and there were 39 and 110 specific lncRNAs. Viafunction, these shared lncRNAs could regulate 3 992 upstream and downstream genes involving in 42 functional terms such as cellular process, metabolism process and catalytic activity; and 117 pathways such as metabolism pathway, biosynthesis of secondary metabolites and biosynthesis of antibiotics. Specific lncRNAs for AaM and AaS could respectively regulate 243 and 672 upstream and downstream genes. InAaM vs AaS comparison group, 255 DElncRNAs were identified and found to regulate 1 479 upstream and downstream genes, which were associated with 41 functional terms including cellular process, metabolism process and catalytic activity; and 107 pathways including metabolism pathway, biosynthesis of secondary metabolites and biosynthesis of antibiotics. Forty one, five and 13 miRNA precursors were predicted from common lncRNAs, specific lncRNAs and DElncRNAs. The result of regulatory network analysis showed the formation of ceRNA networks of mycelium lncRNAs and spore lncRNAs; ten lncRNAs in mycelium could bind to eight miRNAs, further targeting 77 mRNAs; while eight lncRNAs in spore could link to seven miRNAs, further targeting 87 mRNAs; two DElncRNAs including TCONS_00008630 and TCONS_00009302 could simultaneously target miR-4968-y, further regulating ten mRNAs. The result of RT-qPCR suggested the differential expression trend of four DElncRNAs were in accordance of that in sequencing result, indicating the reliability of our sequencing data.【】Common lncRNAs, specific lncRNAs and DElncRNAs are likely to affect material and energy metabolisms, autophagy, transcription, MAPK signaling pathway, ubiquitin-mediated proteolysis, proteasome and biosynthesis of secondary metabolites, by regulating the expression of upstream and downstream genes, serving as miRNA precursors or ceRNAs, thus regulating the growth, development, reproduction and pathogenicity of

; long non-coding RNA (lncRNA); mycelium; spore; competing endogenous RNA (ceRNA)

10.3864/j.issn.0578-1752.2021.02.018

2020-02-10；

2020-02-25

國家自然科學(xué)基金（31702190）、國家現(xiàn)代農(nóng)業(yè)產(chǎn)業(yè)技術(shù)體系建設(shè)專項資金（CARS-44-KXJ7）、福建省自然科學(xué)基金（2018J05042）、福建農(nóng)林大學(xué)碩士生導(dǎo)師團(tuán)隊（郭睿）、福建農(nóng)林大學(xué)科技創(chuàng)新專項基金（CXZX2017342，CXZX2017343）、福建省大學(xué)生創(chuàng)新創(chuàng)業(yè)項目（202010389016，202010389162）

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

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