付中民，陳華枝，劉思亞，祝智威，范小雪，范元嬋，萬潔琦，張璐，熊翠玲，徐國鈞，陳大福，郭睿

意大利蜜蜂響應東方蜜蜂微孢子蟲脅迫的免疫應答

付中民，陳華枝，劉思亞，祝智威，范小雪，范元嬋，萬潔琦，張璐，熊翠玲，徐國鈞，陳大福，郭睿

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

【】通過對意大利蜜蜂（，簡稱意蜂）工蜂中腸響應東方蜜蜂微孢子蟲（）脅迫的差異表達基因（differentially expressed gene，DEG）及免疫通路進行深入細致的分析，揭示宿主的細胞和體液免疫應答，為關(guān)鍵免疫應答基因的篩選及功能研究打下基礎?；谇捌讷@得的正常及東方蜜蜂微孢子蟲脅迫的意蜂7日齡和10日齡工蜂中腸（Am7CK、Am7T、Am10CK、Am10T）轉(zhuǎn)錄組數(shù)據(jù)，按照FDR≤1，≤0.05和|log2fold change|≥1的標準篩選出各組的DEG，利用相關(guān)生物信息學軟件對DEG進行Pearson相關(guān)性分析、Venn分析、GO分類和KEGG代謝通路富集分析，進而對免疫通路富集基因進行統(tǒng)計和分析，利用實時熒光定量PCR（real-time quantitative PCR，RT-qPCR）驗證轉(zhuǎn)錄組數(shù)據(jù)的可靠性。差異分析結(jié)果顯示，Am7CK vs Am7T比較組包含472個上調(diào)基因和385個下調(diào)基因；Am10CK vs Am10T比較組包含611個上調(diào)基因和360個下調(diào)基因。Venn分析結(jié)果顯示，兩個比較組特有的DEG分別為739和853個，共有的DEG為118個。GO分類結(jié)果顯示，Am7CK vs Am7T比較組中上調(diào)和下調(diào)基因分別涉及23和29個功能條目，其中富集上調(diào)基因數(shù)最多的前5位分別為結(jié)合、催化活性、代謝進程、細胞進程和單組織進程；富集下調(diào)基因數(shù)最多的前5位分別為代謝進程、單組織進程、催化活性、細胞進程和結(jié)合。Am10CK vs Am10T比較組中上調(diào)和下調(diào)基因分別涉及36和26個功能條目，其中富集上調(diào)基因數(shù)最多的前5位分別為單組織進程、結(jié)合、細胞進程、催化活性和代謝活性；富集下調(diào)基因數(shù)最多的前5位分別為結(jié)合、細胞進程、催化活性、代謝進程和單組織進程。KEGG代謝通路富集分析結(jié)果顯示，Am7CK vs Am7T比較組中上調(diào)和下調(diào)基因分別富集在38和33條代謝通路，富集上調(diào)基因數(shù)最多的前5條分別為膽汁分泌、內(nèi)質(zhì)網(wǎng)蛋白加工、泛素介導的蛋白水解、PI3K-Akt信號通路和神經(jīng)營養(yǎng)因子信號通路；富集下調(diào)基因數(shù)最多的前5條分別為胞質(zhì)DNA傳感通路、嘌呤代謝、嘧啶代謝、RNA聚合酶和核糖體；涉及泛素介導的蛋白水解等3條細胞免疫通路，以及PI3K-Akt信號通路等7條體液免疫通路。Am10CK vs Am10T比較組中上調(diào)和下調(diào)基因分別富集在54和43條代謝通路，富集上調(diào)基因數(shù)最多的前5條分別為Hippo信號通路、藥物代謝-細胞色素P450、P450對外源物質(zhì)的代謝、泛素介導的蛋白水解和鞘脂類代謝；富集下調(diào)基因數(shù)最多的前5條分別為mRNA監(jiān)測、鞘脂類信號通路、果糖和甘露糖代謝、半乳糖代謝和鞘脂類代謝；涉及泛素介導的蛋白水解等7條細胞免疫通路，以及NF-B信號通路等2條體液免疫通路。RT-qPCR驗證結(jié)果顯示6個隨機挑選的DEG的表達量變化趨勢與測序結(jié)果一致，證實了本研究中測序數(shù)據(jù)的可靠性。進一步分析發(fā)現(xiàn)意蜂7日齡和10日齡工蜂中腸的NF-B信號通路均被東方蜜蜂微孢子蟲激活，隨即啟動了3種抗菌肽apidaecin、defensin-1和hymenoptaecin的合成，體現(xiàn)了它們在宿主抵御東方蜜蜂微孢子蟲入侵中的重要性。在轉(zhuǎn)錄組水平解析了意蜂工蜂中腸對東方蜜蜂微孢子蟲的免疫應答，揭示宿主在脅迫前期同時作出細胞和體液免疫應答，前者可能在抵御病原入侵方面發(fā)揮主要作用；宿主的細胞免疫在脅迫后期持續(xù)增強，但體液免疫較大程度地減弱；泛素介導的蛋白水解通路及富集DEG，NF-B信號通路及富集DEG，以及抗菌肽編碼基因、和可能在意蜂工蜂對東方蜜蜂微孢子蟲的免疫應答中起到關(guān)鍵作用。

意大利蜜蜂；東方蜜蜂微孢子蟲；中腸；免疫應答；細胞免疫；體液免疫

0 引言

【研究意義】蜜蜂是一種專食植物蜜粉的昆蟲，可為70%農(nóng)作物和大部分野生植物授粉[1]，對農(nóng)業(yè)生產(chǎn)、生態(tài)環(huán)境和生物多樣性至關(guān)重要[2]。意大利蜜蜂（，簡稱意蜂）屬于西方蜜蜂（），為我國飼養(yǎng)規(guī)模最大的蜂種，經(jīng)濟效益顯著[3]。東方蜜蜂微孢子蟲（）特異性地侵染成年蜜蜂中腸上皮細胞，對蜂王、工蜂和雄蜂均有侵染性[4-5]，可對蜜蜂造成諸多方面的影響，如細胞破壞[6]、能量脅迫[7-8]和壽命縮短[9]等，還能與其他生物或非生物脅迫因子共同危害蜂群[10-11]。近年來，歐洲和美國等地的西方蜜蜂相繼暴發(fā)蜂群崩潰綜合癥，即蜂群中只剩下蜂王、卵及一些未成年的工蜂和大量蜜粉殘留于巢脾內(nèi)，大量的成年工蜂短期內(nèi)突然在巢外消失，該病已被證明與東方蜜蜂微孢子蟲密切相關(guān)[12-14]。目前，關(guān)于蜜蜂響應東方蜜蜂微孢子蟲的免疫應答機制尚不明確。利用最新的組學技術(shù)對意蜂工蜂中腸響應東方蜜蜂微孢子蟲脅迫的差異表達基因（differentially expressed gene，DEG）及其潛在功能進行系統(tǒng)分析，可在組學水平揭示意蜂對東方蜜蜂微孢子蟲的免疫應答，為在分子水平闡明宿主的免疫應答機制提供重要依據(jù)。【前人研究進展】前人對東方蜜蜂微孢子蟲的研究主要集中在流行病學、致病機理及診斷防治[15-17]等方面。Cornman等于2009年測序、組裝并注釋了東方蜜蜂微孢子蟲的參考基因組[18]，為其組學及分子生物學研究奠定了基礎。此后，在意蜂對東方蜜蜂微孢子蟲的脅迫應答方面進行了一些組學研究，但相關(guān)信息總體仍較為有限[19-20]。Badaoui等通過RNA-seq技術(shù)研究了東方蜜蜂微孢子蟲感染后5、10和15 d西方蜜蜂的基因表達情況，發(fā)現(xiàn)宿主的抗菌肽、角質(zhì)層蛋白和氣味結(jié)合蛋白編碼基因表達量均為下調(diào)，導致宿主的免疫功能下降，從而促進微孢子蟲在蜂群中的存活和增殖[19]；Aufauvre等利用二代測序和實時熒光定量PCR（real-time quantitative PCR，RT-qPCR）技術(shù)對東方蜜蜂微孢子蟲單獨脅迫、殺蟲劑單獨脅迫及二者共同脅迫1 d和7 d的西方蜜蜂中腸樣品進行研究，結(jié)果顯示東方蜜蜂微孢子蟲與殺蟲劑共同脅迫雖然導致西方蜜蜂的死亡率顯著上升，但二者不存在協(xié)同效應；此外，無論是東方蜜蜂微孢子蟲單獨脅迫還是與殺蟲劑共同脅迫，都能導致宿主中腸免疫相關(guān)基因、角質(zhì)層基因和海藻糖代謝相關(guān)基因表達量的顯著改變[20]。筆者團隊前期已利用RNA-seq技術(shù)對正常及東方蜜蜂微孢子蟲脅迫的意蜂7 d和10 d工蜂中腸進行了全轉(zhuǎn)錄組測序，在高質(zhì)量組學數(shù)據(jù)的基礎上對高表達基因（highly expressed gene，HEG）表達譜進行了初步分析[21]；全面解析了意蜂工蜂中腸發(fā)育過程中長鏈非編碼RNA（long non-coding RNA，lncRNA）和環(huán)狀RNA（circular RNA，circRNA）的差異表達譜、調(diào)控網(wǎng)絡及潛在作用機制[22-23]。上述研究結(jié)果為本研究提供了較好的數(shù)據(jù)支持和前期基礎?！颈狙芯壳腥朦c】前期研究中，筆者團隊在mRNA組學水平解析了意蜂10日齡工蜂中腸響應東方蜜蜂微孢子蟲脅迫的HEG表達譜及HEG的潛在作用，為揭示宿主的脅迫應答機制提供了一些有益信息。但僅從基因表達量角度進行分析具有一定的局限性。DEG與東方蜜蜂微孢子蟲脅迫存在直接的聯(lián)系，因此，本研究在前期分析的基礎上，進一步從基因表達差異變化倍數(shù)的角度對意蜂工蜂中腸響應東方蜜蜂微孢子蟲脅迫的DEG表達譜及免疫應答進行分析和探討?！緮M解決的關(guān)鍵問題】基于前期獲得的高質(zhì)量全轉(zhuǎn)錄組數(shù)據(jù)[21-23]，對意蜂工蜂中腸響應東方蜜蜂微孢子蟲脅迫的DEG表達譜及潛在功能進行細致的生物信息學分析和分子生物學驗證，進而對宿主的細胞和體液免疫通路及富集DEG進行深入分析，在轉(zhuǎn)錄組層面揭示意蜂對東方蜜蜂微孢子蟲脅迫的免疫應答。

1 材料與方法

供試意蜂工蜂取自福建農(nóng)林大學蜂學學院教學蜂場。東方蜜蜂微孢子蟲孢子由福建農(nóng)林大學蜂學學院蜜蜂保護實驗室純化和保存。試驗于2017年在福建農(nóng)林大學完成。

1.1 意大利蜜蜂工蜂中腸全轉(zhuǎn)錄組數(shù)據(jù)來源

筆者所在課題組前期已利用Illumina技術(shù)對正常及東方蜜蜂微孢子蟲脅迫的意蜂工蜂中腸進行全轉(zhuǎn)錄組測序，通過嚴格的質(zhì)控和過濾獲得了高質(zhì)量的測序數(shù)據(jù)[22-23]。

工蜂中腸的制備及測序簡述如下：將剛出房的工蜂放入滅菌的四周打孔的塑料盒（20只/盒），每個盒子上方插一支裝有50%（w/v）無菌糖水的飼喂器。將工蜂出房當天記為0 d，（34±0.5）℃培養(yǎng)24 h。饑餓處理2 h，單頭固定，處理組工蜂每頭飼喂含純化孢子（1×106個）的糖水5 μL，對照組工蜂每頭飼喂不含孢子的糖水5 μL，24 h后均飼喂不含孢子的糖水，每日更換糖水，及時清理未食盡糖水的工蜂及死去的工蜂。飼喂至7 d和10 d，分別拉取中腸，置于RNA-Free的EP管，放入液氮速凍。然后迅速轉(zhuǎn)移到-80℃超低溫冰箱保存?zhèn)溆?。每次取樣時間均控制在20 s內(nèi)。處理組和對照組均包含3個生物學重復。意蜂7日齡工蜂中腸對照組和處理組樣品為Am7CK：Am7CK1、Am7CK2、Am7CK3；Am7T：Am7T1、Am7T2、Am7T3；意蜂10日齡工蜂中腸對照組和處理組樣品為Am10CK：Am10CK1、Am10CK2、Am10CK3；Am10T：Am10T1、Am10T2、Am10T3。參照郭睿和陳大福等[22,24]的方法對上述12個中腸樣品進行RNA抽提和cDNA文庫構(gòu)建，委托廣州基迪奧生物技術(shù)有限公司進行雙端測序，測序平臺為Illumina HiSeq 4000。

1.2 測序數(shù)據(jù)質(zhì)控、差異表達基因（DEG）篩選及生物信息學分析

利用短reads比對工具bowtie將前期測序得到的各中腸樣品的有效讀段（clean reads）映射到核糖體數(shù)據(jù)庫，去除比對上核糖體的reads后，將剩余的reads映射到西方蜜蜂參考基因組（Aml_4.5），比對上的clean reads用于后續(xù)的生物信息學分析。

利用公式FPKM=Total exon fragment/[Mapped fragment (millions)×exon length (kB)]計算所有表達基因的FPKM值。利用R軟件計算各樣品之間的相關(guān)性系數(shù)。DEG的篩選標準為FDR≤1，≤0.05和|log2fold change|≥1。利用在線分析工具集合OmicsShare（http://www.omicshare.com/）進行GO分類和KEGG代謝通路富集分析等相關(guān)生物信息學分析，采用默認參數(shù)。

1.3 DEG的實時熒光定量PCR（RT-qPCR）驗證

為了驗證全轉(zhuǎn)錄組測序數(shù)據(jù)的可靠性，隨機選取6個DEG（TCONS_00035591、TCONS_00021864、XM_006569665.2、TCONS_00030754、XM_003251768.3和XM_016912533.1）進行RT-qPCR驗證。參照郭睿等[21]的方法，利用DNAMAN軟件根據(jù)所選DEG的堿基序列設計相應的特異性引物，委托上海生工生物工程有限公司進行合成，相關(guān)引物信息詳見表1。利用RNA抽提試劑盒（TaKaRa公司，中國）分別提取Am7CK、Am7T、Am10CK和Am10T樣品的總RNA，以Oligo (dT)23作為引物進行反轉(zhuǎn)錄，得到對應的cDNA，作為DEG和內(nèi)參基因的qPCR模板。反應體系（20 μL）包含上下游向引物（10.0 μmol·L-1）各1 μL，cDNA模板DNA 1 μL，SYBR Green Dye 10 μL，DEPC水7 μL。反應在QuantStudio熒光定量PCR儀（ThermoFisher公司，美國）上進行，按照SYBR Green Dye試劑盒（Vazyme公司，中國）操作說明書上的方法，每個反應進行3次重復。反應條件：95℃預變性1 min，95℃變性15 s，60℃退火30 s，共40個循環(huán)，最后72℃延伸45 s。利用2-ΔΔCt法對上述DEG的相對表達量進行計算。通過GraphPad Prism軟件進行相關(guān)數(shù)據(jù)分析及繪圖。

表1 RT-qPCR驗證使用的引物信息

2 結(jié)果

2.1 比對西方蜜蜂參考基因組數(shù)據(jù)分析

前期研究結(jié)果表明Am7CK和Am10CK的組內(nèi)各樣品相關(guān)性均在0.95以上[23]。Am7T的組內(nèi)各樣品相關(guān)性介于0.976—1.000，Am10T的組內(nèi)各樣品相關(guān)性均介于0.977—1.000（圖1），表明處理組的樣品重復性較好。通過映射核糖體數(shù)據(jù)庫去除rRNA的reads，共得到1 668 076 502條clean reads，進而通過映射西方蜜蜂基因組得到505 845 865條能夠比對上的clean reads，Q30均在95.16%以上，說明測序數(shù)據(jù)質(zhì)量良好，可用于后續(xù)分析。

圖1 各東方蜜蜂微孢子蟲感染的意蜂工蜂中腸樣品組內(nèi)不同生物學重復間的Pearson相關(guān)性

2.2 差異表達基因分析

差異分析結(jié)果顯示，Am7CK vs Am7T比較組中共有857個基因差異表達，其中上調(diào)基因為472個，下調(diào)基因為385個（圖2-A）；Am10CK vs Am10T比較組包含971個DEG，其中上調(diào)和下調(diào)基因數(shù)分別為611和360個（圖2-B）。進一步對上述兩個比較組的DEG進行Venn分析，結(jié)果顯示二者特有DEG分別為739和853個，共有DEG為118個（圖2-C）。

2.3 DEG的GO分類

GO分類結(jié)果顯示，Am7CK vs Am7T中的DEG主要分為生物學進程、細胞組分和分子功能3類，其中上調(diào)基因涉及23個功能條目，富集基因數(shù)最多的前5位分別是結(jié)合（26）、催化活性（23）、代謝進程（17）、細胞進程（17）和單組織進程（14）；下調(diào)基因分布在29個功能條目，富集基因數(shù)最多的前5位分別是代謝進程（12）、單組織進程（11）、催化活性（11）、細胞進程（10）和結(jié)合（9）（圖3-A—C）。

Am10CK vs Am10T中的DEG同樣分為生物學進程、細胞組分和分子功能，其中上調(diào)基因涉及36個功能條目，富集基因數(shù)最多的前5位分別是單組織進程（27）、結(jié)合（27）、細胞進程（25）、催化活性（23）和代謝進程（20）；下調(diào)基因分布在26個功能條目，富集基因數(shù)最多的前5位分別是結(jié)合（23）、細胞進程（22）、催化活性（21）、代謝進程（21）和單組織進程（17）（圖3-D—E）。

2.4 DEG的KEGG代謝通路富集分析

代謝通路富集分析結(jié)果顯示，Am7CK vs Am7T中上調(diào)基因涉及38條代謝通路，富集基因數(shù)最多的分別是膽汁分泌（2）、內(nèi)質(zhì)網(wǎng)蛋白加工（2）、泛素介導的蛋白水解（2）、PI3K-Akt信號通路（2）、神經(jīng)營養(yǎng)因子信號通路（1）、趨化因子信號通路（1）、近端小管碳酸氫鹽再生（1）、集管酸分泌（1）、催乳激素分泌（1）和胃酸分泌（1）；下調(diào)基因涉及33條代謝通路，富集基因數(shù)最多的分別是胞質(zhì)DNA傳感通路（2）、嘌呤代謝（2）、嘧啶代謝（2）、RNA聚合酶（2）、核糖體（2）、脂肪的消化和吸收（1）、血管平滑肌收縮（1）、蛋白質(zhì)的消化和吸收（1）、胰液分泌（1）和氧化磷酸化（1）。

A：Am7CK vs Am7T中DEG的火山圖Volcano plot of DEGs in Am7CK vs Am7T；B：Am10CK vs Am10T中DEG的火山圖Volcano plot of DEGs in Am10CK vs Am10T。紅色圓點代表上調(diào)基因，綠色圓點代表下調(diào)基因Red dots indicate up-regulated genes, green dots indicate down-regulated genes。C：Am7CK vs Am7T和Am10CK vs Am10T中DEG的Venn圖 Venn diagram of DEGs in Am7CK vs Am7T and Am10CK vs Am10T

Am10CK vs Am10T中的上調(diào)基因涉及54條代謝通路，富集基因數(shù)最多的分別是Hippo信號通路（3）、藥物代謝-細胞色素P450（2）、P450對外源物質(zhì)的代謝（2）、泛素介導的蛋白水解（2）、鞘脂類代謝（2）、氨基糖和核苷酸代謝（2）、藥物代謝-其他酶（2）、鞘脂類信號通路（2）、Hedgehog信號通路（2）和血小板激活（1）；下調(diào)基因涉及43條代謝通路，富集基因數(shù)最多的分別是mRNA監(jiān)控（3）、鞘脂類信號通路（3）、果糖和甘露糖代謝（2）、半乳糖代謝（2）、鞘脂類代謝（2）、氨基糖和核苷酸糖代謝（2）、RNA轉(zhuǎn)運（2）、心肌細胞的腎上腺素能信號（1）、突出囊泡循環(huán)（1）和背腹軸形成（1）。

2.5 細胞和體液免疫通路分析

進一步對免疫通路及富集基因進行分析，發(fā)現(xiàn)Am7CK vs Am7T中的DEG可注釋到3條細胞免疫通路，包括泛素介導的蛋白水解、細胞凋亡和溶酶體（表2）；以及7條體液免疫通路，包括PI3K-Akt信號通路、趨化因子信號通路、NF-B信號通路、cAMP信號通路、Ras信號通路、MAPK信號通路-酵母和FoxO信號通路（表3）。對于細胞免疫通路，富集在泛素介導的蛋白水解信號通路上的2個基因均上調(diào)表達，富集在細胞凋亡上的1個基因也呈現(xiàn)出表達上調(diào)，但富集在溶酶體上的1個基因表達量下調(diào)。對于體液免疫通路，富集在趨化因子信號通路和NF-B信號通路上的基因均為上調(diào)表達，但富集在cAMP信號通路、Ras信號通路、MAPK信號通路-酵母和FoxO信號通路上的基因均表現(xiàn)為下調(diào)。

A—C：Am7CK vs Am7T中上調(diào)和下調(diào)基因分別在生物學進程、分子功能和細胞組分相關(guān)條目的分類Classifications of up- and down-regulated genes within Am7CK vs Am7T comparison group in GO terms associated with biological process, molecular function and cellular component；D—F：Am10CK vs Am10T中上調(diào)和下調(diào)基因分別在生物學進程、分子功能和細胞組分相關(guān)條目的分類Classifications of up- and down-regulated genes within Am10CK vs Am10T comparison group in GO terms associated with biological process, molecular function and cellular component

在Am10CK vs Am10T中，DEG分別注釋到7條細胞免疫通路和2條體液免疫通路。對于細胞免疫通路，富集在泛素介導的蛋白水解、P450對外源物質(zhì)代謝、藥物代謝-細胞色素P450通路、細胞凋亡、Fc--R介導的吞噬作用和溶酶體通路上的基因都呈現(xiàn)出上調(diào)表達，僅有1個富集在內(nèi)吞作用上的基因表現(xiàn)為下調(diào)（表4）。對于體液免疫通路，分別有1個上調(diào)基因和1個下調(diào)基因富集在NF-B信號通路和PI3K-Akt信號通路（表5）。

Am7CK vs Am7T和Am10CK vs Am10T比較組中DEG在泛素介導的蛋白水解上的富集詳情如圖4所示。

表2 Am7CK vs Am7T中DEG富集的細胞免疫相關(guān)通路

表3 Am7CK vs Am7T中DEG富集的體液免疫相關(guān)通路

表4 Am10CK vs Am10T中DEG富集的細胞免疫相關(guān)通路

表5 Am10CK vs Am10T中DEG富集的體液免疫相關(guān)通路

A：Am7CK vs Am7T中DEG在泛素介導蛋白水解上的富集情況Enrichment of DEGs within Am7CK vs Am7T comparison group in ubiquitin mediated proteolysis；B：Am10CK vs Am10T中DEG在泛素介導蛋白水解上的富集情況Enrichment of DEGs within Am10CK vs Am10T comparison group in ubiquitin mediated proteolysis。紅色方框代表上調(diào)基因Red boxes represent up-regulated genes

2.6 DEG的RT-qPCR驗證

RT-qPCR檢測結(jié)果顯示，6個DEG的表達水平變化趨勢與測序數(shù)據(jù)中的基因表達水平變化趨勢一致（圖5），證明本研究中的mRNA組學數(shù)據(jù)真實可靠。

A：乙酰轉(zhuǎn)移酶編碼基因（TCONS_00035591）predicted acetyltransferase encoded gene；B：牙本質(zhì)涎磷蛋白編碼基因（TCONS_00021864）dentin sialophosphoprotein encoded gene；C：富含甘氨酸的細胞壁結(jié)構(gòu)蛋白1.8-like型X2編碼基因（XM_006569665.2） glycine-rich cell wall structural protein 1.8-like isoform X2 encoded gene；D：金星激酶受體編碼基因（TCONS_00030754）venus kinase receptor encoded gene；E：G12蛋白編碼基因（XM_003251768.3） protein G12 encoded gene；F：內(nèi)切型2編碼基因（XM_016912533.1）endochitinase isoform 2encoded gene。RT-qPCR組中，*表示p＜0.05，**表示p＜0.01 In RT-qPCR group, * indicates p＜0.05, ** indicates p＜0.01

3 討論

前人對于蜜蜂響應東方蜜蜂微孢子蟲的脅迫應答雖有一些研究[18-20]，但多集中在分子生物學水平，組學相關(guān)研究較為有限。前期研究中，筆者團隊利用RNA-seq技術(shù)對正常及東方蜜蜂微孢子蟲脅迫的意蜂工蜂中腸進行了全轉(zhuǎn)錄組數(shù)據(jù)，獲得了高質(zhì)量的mRNA、miRNA、lncRNA和circRNA組學數(shù)據(jù)[21-23]；通過生物信息學方法篩選出對照組和處理組的高表達基因（HEG），并通過深入分析揭示了HEG的表達譜及潛在作用[21]。HEG分析是從基因表達量角度對宿主的脅迫應答進行初步分析，雖然能得到一些有益信息，但具有不可避免的局限性，體現(xiàn)在有些基因?qū)τ谝夥涔し渲心c的生長發(fā)育十分重要，其表達水平自始至終都維持高量表達，這部分HEG會對相關(guān)分析造成一些影響；此外，HEG的分析只涉及被東方蜜蜂微孢子蟲激活的基因，但被抑制的基因同樣值得深入研究。免疫應答是一個復雜的動態(tài)過程，闡明其背后的分子機制需要從不同角度盡可能深入地進行分析，以獲取盡可能多的有效信息。與HEG相比，差異表達基因（DEG）與東方蜜蜂微孢子蟲脅迫存在直接的聯(lián)系。為進一步探討意蜂工蜂中腸響應東方蜜蜂微孢子蟲脅迫的免疫應答，本研究在前期研究基礎上對DEG及其潛在作用進行相關(guān)分析，發(fā)現(xiàn)Am7CK vs Am7T中的上調(diào)和下調(diào)基因數(shù)分別為472和385個；Am10CK vs Am10T中的上調(diào)和下調(diào)基因數(shù)分別為611和360個；表明隨著脅迫時間的延長，更多的宿主基因被激活，而受到抑制的基因數(shù)基本一致，脅迫應答水平呈增強趨勢。進一步分析發(fā)現(xiàn)，有118個DEG為上述兩個比較組所共有，特有DEG的數(shù)量分別為739和853個，推測這些共有DEG在宿主的脅迫應答過程發(fā)揮基礎性的作用，而特有DEG在宿主脅迫應答的不同階段具有特定功能。此外，將Am10CK vs Am10T中的DEG與前期篩選出的Am10T的HEG[21]進行Venn分析，結(jié)果顯示有5個基因既為高量表達又為差異表達，包括半胱氨酸調(diào)節(jié)子的轉(zhuǎn)錄調(diào)節(jié)劑（transcriptional regulator for cysteine regulon）編碼基因、蜂毒肽前體蛋白（melittin precursor）編碼基因、白細胞蛋白酶抑制劑（antileukoproteinase-like isoform X2）編碼基因、硫氧還原蛋白（thioredoxin-2 isoform 1）編碼基因和胰蛋白酶（trypsin-7）編碼基因，推測這些基因在意蜂工蜂響應東方蜜蜂微孢子蟲脅迫的過程中發(fā)揮關(guān)鍵作用，值得進一步深入研究。下一步將針對上述5個基因進行全長序列的分子克隆及功能研究（例如RNAi敲減），以在分子水平探究它們的功能。

富集在細胞殺傷和應激反應條目的DEG可能與免疫應答密切相關(guān)。本研究發(fā)現(xiàn)，Am7CK vs Am7T中分別有6個上調(diào)基因和4個下調(diào)基因富集在應激反應，2個下調(diào)基因富集在細胞殺傷；Am10CK vs Am10T中分別有3個上調(diào)基因和3個下調(diào)基因富集在應激反應，1個上調(diào)基因和1個下調(diào)基因富集在細胞殺傷；表明部分基因處于激活狀態(tài)、部分基因處于抑制狀態(tài)，暗示著意蜂工蜂在對東方蜜蜂微孢子蟲的脅迫響應過程，既能產(chǎn)生直接的免疫應答，也會受到病原一定程度的抑制，反映出宿主和病原之間存在復雜的互作關(guān)系。

蜜蜂的免疫防御分為群體免疫和個體免疫，后者又分為細胞免疫和體液免疫[25]。本研究中，在東方蜜蜂微孢子蟲脅迫的意蜂7日齡工蜂中腸中，分別有2個（log2fold change=1.766和1.685）和1個上調(diào)基因（log2fold change=1.685）富集在泛素介導的蛋白水解和細胞凋亡，僅有1個下調(diào)基因（log2fold change= -1.02）富集在溶酶體；分別有1個（log2fold change=12.805）和1個上調(diào)基因（log2fold change=1.685）富集在趨化因子和NF-B信號通路，1個（log2fold change=-9.103）、1個（log2fold change=-1.938）、1個（log2fold change= -13.391）和1個下調(diào)基因（log2fold change=-13.391）富集在cAMP、Ras、MAPK和FoxO信號通路；分別有2個上調(diào)基因（log2fold change=12.805和1.766）和1個下調(diào)基因（log2fold change=-2.023）富集在PI3K-Akt；表明意蜂工蜂中腸的細胞免疫很大程度被激活，而體液免疫部分程度被抑制、部分程度被激活。在東方蜜蜂微孢子蟲脅迫的意蜂10日齡工蜂中腸中，富集在細胞免疫通路的絕大多數(shù)DEG都為上調(diào)表達，僅有1個富集在內(nèi)吞作用的基因下調(diào)表達（log2fold change=-1.543）；分別有1個上調(diào)基因（log2fold change= 1.016）和1個下調(diào)基因（log2fold change= -2.569）富集在NF-B和PI3K-Akt信號通路；表明隨著東方蜜蜂微孢子蟲脅迫時間的延長，意蜂工蜂中腸的細胞免疫被持續(xù)激活且激活程度不斷增強，而體液免疫的種類減少、程度下降。上述結(jié)果說明意蜂工蜂中腸在東方蜜蜂微孢子蟲脅迫的前期，免疫系統(tǒng)識別到病原入侵，同時作出細胞和體液應答，此時宿主的細胞免疫很大程度被激活，在抵御病原入侵方面發(fā)揮主要作用，宿主的體液免疫既有部分被激活也有部分被抑制，體現(xiàn)了意蜂及東方蜜蜂微孢子蟲互作的復雜性；在東方蜜蜂微孢子蟲脅迫的后期，宿主的細胞免疫更大范圍、更大程度地被激活，在免疫防御中的作用繼續(xù)增強，而宿主的體液免疫通路較大程度地下降，其中NF-B信號通路呈現(xiàn)出激活狀態(tài)，PI3K-Akt信號通路呈現(xiàn)出抑制狀態(tài)。

泛素介導的蛋白水解是一種經(jīng)典的細胞免疫通路，該通路由E1酶（泛素活化酶）、E2酶（泛素結(jié)合酶）和E3酶（泛素連接酶）組成，其中E3可特異性識別不同底物，然后與E2酶結(jié)合，E2酶集合E1酶活化的泛素，泛素與底物形成多聚泛素鏈，最終被蛋白質(zhì)酶體識別降解[26]。本研究中，對于東方蜜蜂微孢子蟲脅迫的意蜂7日齡和10日齡工蜂中腸，共有4個上調(diào)基因編碼3種E3酶，推測宿主通過提高E3酶的合成以識別東方蜜蜂微孢子蟲的關(guān)鍵蛋白質(zhì)，從而對其產(chǎn)生抑制。另外，病原也能夠利用泛素介導的蛋白水解的不同環(huán)節(jié)以逃避宿主細胞免疫系統(tǒng)的監(jiān)控[27]。蜜蜂白堊病是另一種常見的蜜蜂真菌病，該病的病原蜜蜂球囊菌（，簡稱球囊菌）特異性侵染蜜蜂幼蟲，早期的增殖場所為幼蟲腸道。前期研究中，筆者團隊發(fā)現(xiàn)在球囊菌脅迫的中華蜜蜂（）幼蟲腸道中，靶向結(jié)合泛素介導的蛋白水解通路上的3個基因的miRNA-21-x表達量顯著上調(diào)，作者推測球囊菌通過調(diào)控宿主泛素介導的蛋白水解通路基因的表達以協(xié)助球囊菌逃避宿主免疫系統(tǒng)的監(jiān)控[28]。東方蜜蜂微孢子蟲特異性寄生蜜蜂中腸上皮細胞，其能否利用miRNA介導的方法調(diào)控泛素介導的蛋白水解通路基因從而逃脫宿主的免疫系統(tǒng)監(jiān)控，值得進一步深入研究。陳陽研究發(fā)現(xiàn)泛素介導的蛋白水解與Toll樣受體（TLRs）信號通路共同參與鴨的基因介導的抗病毒天然免疫信號的傳遞途徑[29]。果蠅的Toll免疫信號通路主要識別革蘭氏陽性菌和真菌等病原微生物，通過激活Toll受體等一系列信號分子激活NF-B同源蛋白Dif和Dorsal并釋放到核內(nèi)，最終啟動下游抗菌肽分子的表達[30]。本研究發(fā)現(xiàn)，意蜂7日齡和10日齡工蜂中腸的NF-B信號通路都被激活，且啟動了抗菌肽編碼基因、和（NM_001011613.1、NM_001011615.1和NM_001011616.2）的表達，表明NF-B信號通路在宿主免疫應答中表現(xiàn)活躍。綜上，推測泛素介導的蛋白水解通路及富集DEG，NF-B信號通路及富集DEG，以及、和在意蜂工蜂對東方蜜蜂微孢子蟲的免疫應答中發(fā)揮關(guān)鍵作用。

關(guān)于mRNA與非編碼RNA（non-coding RNA，ncRNA）之間的相互作用，Salmena等[31]提出了“競爭性內(nèi)源RNA”假說，即任何包含miRNA反應元件（miRNA response element）的RNA如mRNA、假基因、lncRNA和circRNA，均可競爭性地結(jié)合miRNA，從而影響基因的表達水平。前期研究中，筆者團隊基于高質(zhì)量的全轉(zhuǎn)錄組數(shù)據(jù)對意蜂工蜂中腸發(fā)育過程的ncRNA進行了系統(tǒng)的差異表達譜及調(diào)控網(wǎng)絡分析，篩選出若干與中腸發(fā)育密切相關(guān)的關(guān)鍵lncRNA和circRNA[22-23]。本研究對意蜂工蜂中腸響應東方蜜蜂微孢子蟲脅迫的DEG進行了全面分析，得到了宿主響應東方蜜蜂微孢子蟲脅迫的基因表達譜信息，但相關(guān)基因集仍較為復雜。目前，筆者團隊正在對意蜂工蜂響應東方蜜蜂微孢子蟲脅迫的差異表達miRNA、差異表達lncRNA和差異表達circRNA進行深入分析，擬通過構(gòu)建和分析DEG與差異表達ncRNA之間的ceRNA調(diào)控網(wǎng)絡，篩選出居于核心位置的關(guān)鍵應答基因，作為后續(xù)功能研究的分子靶標。

4 結(jié)論

在mRNA組學水平對意蜂工蜂中腸響應東方蜜蜂微孢子蟲脅迫的免疫應答進行了深入解析，發(fā)現(xiàn)意蜂工蜂中腸在脅迫前期同時作出細胞和體液免疫應答，前者可能在抵御病原入侵方面發(fā)揮主要作用；宿主的細胞免疫在脅迫后期的作用持續(xù)增強，但體液免疫的作用較大程度的減弱；泛素介導的蛋白水解通路及富集DEG，NF-B信號通路及富集DEG，以及、和在意蜂工蜂對東方蜜蜂微孢子蟲的免疫應答中起到關(guān)鍵作用。研究結(jié)果為在分子水平闡明意蜂響應東方蜜蜂微孢子蟲的免疫應答機制提供了數(shù)據(jù)支持。

[1] KLEIN A M, VAISSIèRE B E, CANE J H, STEFFAN-DEWENTER I, CUNNINGHAM S A, KREMEN C, TSCHARNTKE T. Importance of pollinators in changing landscapes for world crops., 2007, 274(1608): 303-313.

[2] POTTS S G, BIESMEIJER J C, KREMEN C, NEUMANN P, SCHWEIGER O, KUNIN W E. Global pollinator declines: trends, impacts and drivers., 2010, 25(6): 345-353.

[3] 曾志將, 顏偉玉, 饒波, 謝憲斌. 蜜蜂性比研究進展. 上海交通大學學報(農(nóng)業(yè)科學版), 2003, 21(2): 164-167.

ZENG Z J, YAN W Y, RAO B, XIE X B. Advances on the honeybee sex ratio., 2003, 21(2): 164-167. (in Chinese)

[4] TRAVER B E, FELL R D. Low natural levels ofinqueens., 2012, 110(3): 408-410.

[5] TRAVER B E, FELL R D.in drone honey bees ()., 2011, 107(3): 234-236.

[6] GARCIA-PALENCIA P, MARTIN-HERNANDEZ R, GONZALEZ- PORTO A V, MARIN P, MEANA A, HIGES M. Natural infection bycauses similar lesions as in experimentally infected caged-worker honey bees ()., 2010, 49(3): 278-283.

[7] MAYACK C, NAUG D. Energetic stress in the honeybeefrominfection., 2009, 100(3): 185-188.

[8] MARTIN-HERNANDEZ R, BOTIAS C, BARRIOS L, MARTINEZ- SALVADOR A, MEANA A, MAYACK C, HIGES M. Comparison of the energetic stress associated with experimentalandinfection of honeybees ()., 2011, 109(3): 605-612.

[9] LECOCQ A, JENSEN A B, KRYGER P, NIEH J C. Parasite infection accelerates age polyethism in young honey bees., 2016, 6: 22042.

[10] SHUTLER D, HEAD K, BURGHER-MACLELLAN K L, COLWELL M J, LEVITT A L, OSTIGUY N, WILLIAMS G R. Honey beeparasites in the absence offungi andmites., 2014, 9(6): e98599.

[11] WU J Y, SMART M D, ANELLI C M, SHEPPARD W S. Honey bees () reared in brood combs containing high levels of pesticide residues exhibit increased susceptibility to(Microsporidia) infection., 2012, 109(3): 326-329.

[12] COX-FOSTER D L, CONLAN S, HOLMES E C, PALACIOS G, EVANS J D, MORAN N A, QUAN P L, BRIESE T, HORNIG M, GEISER D M, MARTINSON V, VANENGELSDORP D, KALKSTEIN A L, DRYSDALE A, HUI J, ZHAI J, CUI L, HUTCHISON S K, SIMONS J F, EGHOLM M, PETTIS J S, LIPKIN W I. A metagenomic survey of microbes in honey bee colony collapse disorder., 2007, 318(5848): 283-287.

[13] CHEJANOVSKY N, OPHIR R, SCHWAGER M S, SLABEZKI Y, GROSSMAN S, COX-FOSTER D. Characterization of viral siRNA populations in honey bee colony collapse disorder., 2014, 454/455: 176-183.

[14] KRIBS-ZALETA C M, MITCHELL C. Modeling colony collapse disorder in honeybees as a contagion., 2014, 11(6): 1275-1294.

[15] MAYACK C, NATSOPOULOU M E, MCMAHON D P.alters a highly conserved hormonal stress pathway in honeybees., 2015, 24(6): 662-670.

[16] NATSOPOULOU M E, DOUBLET V, PAXTON R J. European isolates of the microsporidiaandhave similar virulence in laboratory tests on European worker honey bees., 2016, 47(1): 57-65.

[17] BRAVO J, CARBONELL V, SEPúLVEDA B, DELPORTE C, VALDOVINOS C E, MARTíN-HERNáNDEZ R, HIGES M. Antifungal activity of the essential oil obtained fromagainst infection in honey bees by., 2017, 149: 141-147.

[18] CORNMAN R S, CHEN Y P, SCHATZ M C, STREET C, ZHAO Y, DESANY B, EGHOLM M, HUTCHISON S, PETTIS J S, LIPKIN W I, EVANS J D. Genomic analyses of the microsporidian, an emergent pathogen of honey bees., 2009, 5(6): e1000466.

[19] BADAOUI B, FOUGEROUX A, PETIT F, ANSELMO A, GORNI C, CUCURACHI M, CERSINI A, GRANATO A, CARDETI G, FORMATO G, MUTINELLI F, GIUFFRA E, WILLIAMS J L, BOTTI S. RNA-sequence analysis of gene expression from honeybees ()infected with., 2017, 12(3): e0173438.

[20] AUFAUVRE J, MISME-AUCOUTURIER B, Viques B, TEXIER C, DELBAC F, BLOT N. Transcriptome analyses of the honeybee response toand insecticides., 2014, 9(3): e91686.

[21] 郭睿, 刀晨, 熊翠玲, 鄭燕珍, 付中民, 耿四海, 陳大福. 意大利蜜蜂工蜂中腸響應脅迫的高表達基因分析. 環(huán)境昆蟲學報, 2018, 40(5): 1106-1112.

GUO R, DAO C, XIONG C L, ZHENG Y Z, FU Z M, GENG S H, CHEN D F. Analysis of the highly expressed genes in the midgut ofworker under the stress of., 2018, 40(5): 1106-1112. (in Chinese)

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

[23] 郭睿, 陳華枝, 熊翠玲, 鄭燕珍, 付中民, 徐國鈞, 杜宇, 王海朋, 耿四海, 周丁丁, 劉思亞, 陳大福. 意大利蜜蜂工蜂中腸發(fā)育過程中的差異表達環(huán)狀RNA及其調(diào)控網(wǎng)絡分析. 中國農(nóng)業(yè)科學, 2018, 51(23): 4575-4590.

GUO R, CHEN H Z, XIONG C L, ZHENG Y Z, FU Z M, XU G J, DU Y, WANG H P, GENG S H, ZHOU D D, LIU S Y, CHEN D FAnalysis of differentially expressed circular RNAs and their regulation networks during the developmental process ofworker’s midgut, 2018, 51(23): 4575-4590. (in Chinese)

[24] 陳大福, 郭睿, 熊翠玲, 梁勤, 鄭燕珍, 徐細建, 張曌楠, 黃枳腱, 張璐, 王鴻權(quán), 解彥玲, 童新宇. 中華蜜蜂幼蟲腸道響應球囊菌早期脅迫的轉(zhuǎn)錄組學. 中國農(nóng)業(yè)科學, 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 L, TONG X Y. Transcriptome oflarval gut under the stress of., 2017, 50(13): 2614-2623. (in Chinese)

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

[26] CIECHANOVER A, 葛亮. 泛素介導的蛋白質(zhì)降解系統(tǒng)——從基礎研究到臨床應用. 生命科學, 2010, 22(3): 212-215.

CIECHANOVER A, GE L. The ubiquitin proteolytic system—from bench to the bedside., 2010, 22(3): 212-215. (in Chinese)

[27] 李朝飛, 龐義. 泛素-蛋白水解酶復合體通路與病毒侵染. 生物工程學報, 2004, 20(2): 151-156.

LI C F, PANG Y. Ubiquitin-proteasome pathway and virus infection., 2004, 20(2): 151-156. (in Chinese)

[28] 杜宇, 童新宇, 周丁丁, 陳大福, 熊翠玲, 鄭燕珍, 徐國鈞, 王海朋, 陳華枝, 郭意龍, 隆琦, 郭睿. 中華蜜蜂幼蟲腸道響應球囊菌脅迫的microRNA應答分析. 微生物學報, 2019. DOI: 10.13343/j.cnki. wsxb.20180401.

DU Y, TONG X Y, ZHOU D D, CHEN D F, XIONG C L, ZHENG Y Z, XU G J, WANG H P, CHEN H Z, GUO Y L, LONG Q, GUO R. MicroRNA responses in the larval gut oftostress., 2019. DOI: 10. 13343/j.cnki.wsxb.20180401. (in Chinese)

[29] 陳陽. 鴨RIG-I在免疫調(diào)節(jié)中的作用及分子調(diào)控機制[D]. 揚州: 揚州大學, 2015.

CHEN Y. Function and molecular mechanism of duck RIG-I gene in immune regulation[D]. Yangzhou: Yangzhou University, 2015. (in Chinese)

[30] 徐嬌. miR-964在果蠅Toll免疫響應中的調(diào)控作用研究[D]. 南京: 南京師范大學, 2018.

XU J. Study on the regulation of miR-964 in the immune response of drosophila Toll[D]. Nanjing: Nanjing Normal University, 2018. (in Chinese)

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

Immune responses oftostress

FU zhongMin, CHEN HuaZhi, LIU SiYa, ZHU ZhiWei, FAN XiaoXue, FAN YuanChan, WAN JieQi, ZHANG Lu, XIONG CuiLing, XU GuoJun, CHEN DaFu, GUO Rui

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

【】The objective of this study is to reveal the cellular and humoral immune responses ofworker’s midgut tostress, and to provide a basis for screening and functional study of key immune response genes by deep investigation of host differentially expressed genes (DEGs) and immune pathways.【】Base on the previously obtained transcriptome datasets of normal and-stressed midguts of7- and 10-day-old workers’ midgut samples (Am7CK, Am7T, Am10CK and Am10T), DEGs in each comparison group were filtered out following the standard FDR≤1,≤0.05 and |log2fold change|≥1. Additionally, Pearson correlations, Venn analysis, GO classification and KEGG pathway enrichment analysis were conducted by using related bioinformatic softwares. Further, DEGs enriched in immune pathways were summarized and analyzed. Finally, the reliability of transcriptome data was verified via real-time quantitative PCR (RT-qPCR).【】Differential expression analysis showed that there were 472 up-regulated and 385 down-regulated genes in Am7CK vs Am7T comparison group, and 611 up-regulated and 360 down-regulated genes in Am10CK vs Am10T comparison group. Venn analysis suggested that the number of specific DEGs of the aforementioned two comparison groups was 739 and 853, respectively, and 118 DEGs were shared. GO classification indicated that the up- and down-regulated genes in Am7CK vs Am7T were respectively annotated to 23 and 29 functional terms; among them, the top five terms of up-regulated genes were binding, catalytic activity, metabolic process, cellular process and single tissue process; while the top five terms of down-regulated genes were metabolic process, single tissue process, catalytic activity, cellular process and binding. In addition, the up- and down-regulated genes in Am10CK vs Am10T were respectively involved in 36 and 26 functional terms; among them, the top five terms of up-regulated genes were single process, binding, cellular process, catalytic activity and metabolic activity; while the top five terms of down-regulated genes were binding, cellular process, catalytic activity, metabolic process and single process. KEGG metabolic pathway enrichment analysis showed that the up- and down-regulated genes in Am7CK vs Am7T were respectively enriched in 38 and 33 pathways, and among them the top five enriched pathways of up-regulated genes were bile secretion, endoplasmic reticulum protein processing, ubiquitin-mediated proteolysis, PI3K-Akt signaling pathway and neurotrophic factor signaling pathway; while the top five enriched pathways of down-regulated genes were cytoplasmic DNA-sensing pathway, purine metabolism, pyrimidine metabolism, RNA polymerase and ribosome; three cellular immune pathways including ubiquitin-mediated proteolysis, and seven humoral pathways including PI3K-Akt signaling pathway were detected. Additionally, the up- and down-regulated genes in Am10CK vs Am10T were respectively associated with 54 and 43 pathways, and among them the top five enriched pathways of up-regulated genes were Hippo signaling pathway, drug metabolism-cytochrome P450, metabolism of xenobiotics by cytochrome P450, ubiquitin-mediated proteolysis, and sphingolipid metabolism; while the top five enriched pathways of down-regulated genes were mRNA surveillance pathway, sphingolipid signaling pathway, fructose and mannose metabolism, galactose metabolism, and sphingolipid metabolism; seven cellular immune pathways including ubiquitin-mediated proteolysis, and two humoral pathways including NF-B signaling pathway were found. The result of RT-qPCR verification demonstrated the expression trend of six randomly selected DEGs was consistent with that of the sequencing data, which confirmed the reliability of the transcriptome data. Further analysis noted that NF-B signaling pathway was activated byin both 7- and 10-day-old workers’ midguts of, followed by immediate initiation of three antimicrobial peptides such as apidaecin, defensin-1 and hymenoptaecin, indicating their key importance in host defense againstinvasion.【】The immune responses ofworker’s midgut towere uncovered at transcriptomic level, and it was revealed that the host made both cellular and humoral immune responses at the early stage ofstress. the host cellular immune responses may play a major role in resisting pathogen invasion. The host cellular immune continued to increase at the later stage of fungal stress, while the humoral immune drastically weakened. The ubiquitin mediated proteolysis and enriched DEGs, NF-B signaling pathway and enriched DEGs, and antibacterial peptide-encoded genes such as,andwere likely to play pivotal roles in host immune response tostress.

;; midgut; immune response; cellular immune; humoral immune

2019-04-17；

2019-05-24

國家現(xiàn)代農(nóng)業(yè)產(chǎn)業(yè)技術(shù)體系建設專項資金（CARS-44-KXJ7）、福建省自然科學基金（2018J05042）、福建省教育廳中青年教師教育科研項目（JAT170158）、福建農(nóng)林大學杰出青年科研人才計劃（xjq201814）、福建農(nóng)林大學科技創(chuàng)新專項基金（CXZX2017342，CXZX2017343）、福建省大學生創(chuàng)新創(chuàng)業(yè)訓練計劃（3165602032，3155006018）

付中民，E-mail：369699776@qq.com。陳華枝，E-mail：18965015689@163.com。付中民和陳華枝為同等貢獻作者。通信作者郭睿，E-mail：ruiguo@fafu.edu.cn

（責任編輯 岳梅）