熊翠玲， 黃枳腱， 梁 勤， 徐細(xì)建， 張曌楠， 駱 群， 劉 敏， 陳大福， 郭 睿

(1.福建農(nóng)林大學(xué)蜂學(xué)學(xué)院，福建福州 350002； 2.江西省養(yǎng)蜂研究所，江西南昌 330200)

蜜蜂白堊病是由蜜蜂球囊菌侵染蜜蜂幼蟲而導(dǎo)致的致死性真菌病[1]。該病主要分布在歐洲、北美和中國等地區(qū)，早在1913年，Maassen等就報道了白堊病的發(fā)生，我國大陸1961年也曾報道發(fā)生類似蜜蜂白堊病的病害[2]。蜜蜂白堊病主要是使老熟幼蟲或封蓋幼蟲死亡，雄蜂幼蟲最易感染，因為雄蜂幼蟲多在巢脾的邊緣，易受低溫的影響。白堊病會導(dǎo)致大批幼蟲死亡，造成群勢下降，嚴(yán)重病群失去產(chǎn)蜂蜜和蜂王漿的能力，嚴(yán)重影響蜂群的生產(chǎn)力。球囊菌屬于球囊菌門球囊菌屬，迄今已鑒定出球囊菌屬的22個種，它們以致病性真菌或腐生菌存在于蜜蜂體內(nèi)；在蜂群里，患病幼蟲的尸體以及被污染的飼料與巢脾是疾病傳播的主要來源[3]。國內(nèi)養(yǎng)蜂生產(chǎn)中使用的主要蜂種是意大利蜜蜂和中華蜜蜂，意蜂極易受到球囊菌的侵染而暴發(fā)白堊病，而中蜂幾乎不受影響。2006年，Evans等公布了西方蜜蜂的基因組，為在分子水平研究西方蜜蜂奠定了基礎(chǔ)[4]。對于白堊病的前期研究主要是集中在生化[5]、病理[6]及檢測[7]等方面，分子水平的研究相對較少。目前，高通量測序技術(shù)應(yīng)用于白堊病的研究報道極少，蜜蜂幼蟲響應(yīng)球囊菌脅迫的應(yīng)答研究仍不深入。

以RNA seq為代表的高通量測序技術(shù)發(fā)展迅猛，現(xiàn)已廣泛應(yīng)用于昆蟲包括蜜蜂的相關(guān)研究[8-11]，然而，尚無利用二代測序技術(shù)研究白堊病過程蜜蜂幼蟲應(yīng)答的相關(guān)研究，本研究利用RNA-seq技術(shù)對球囊菌脅迫和正常意蜂幼蟲腸道進行深度測序，對表達水平最高的100個基因進行GO分類和KEGG代謝通路富集分析，并比較分析球囊菌侵染幼蟲腸道與正常幼蟲腸道之間高表達基因的差異。

1 材料與方法

1.1 生物材料

本研究中使用的意大利蜜蜂幼蟲取自福建農(nóng)林大學(xué)蜂學(xué)學(xué)院教學(xué)蜂場，蜜蜂球囊菌菌株由福建農(nóng)林大學(xué)蜂學(xué)學(xué)院蜜蜂保護實驗室保存并活化。

1.2 主要試驗試劑及儀器

RNase-free水購自中國上海生工生物公司，DNaseⅠ和Oligotex mRNA Kits Midi試劑盒購自德國Qiagen公司，Dynal M280磁珠購自Invitrogen公司，高碘酸鈉購自美國Sigma公司，DNA ligase購自美國Thermo公司，RNA Reagent抽提試劑盒、ExTaqpolymerase及Superscript Ⅱ reverse transcriptase均購自日本TaKaRa公司。純化cDNA的Ampure beads為美國Agencourt產(chǎn)品，cDNA文庫構(gòu)建試劑盒TruSeqTMDNA Sample Prep Kit -Set A為美國Illumina公司產(chǎn)品。

恒溫恒濕氣候箱購自中國寧波江南儀器廠，超純水儀購自中國四川沃特爾水處理設(shè)備有限公司，高速冷凍離心機購自德國Eppendorf公司，倒置顯微鏡為中國上海光學(xué)儀器五廠產(chǎn)品，超凈工作臺為中國蘇州安泰空氣技術(shù)有限公司產(chǎn)品，PCR儀為美國Bio Rad公司產(chǎn)品，凝膠成像系統(tǒng)為中國上海培清科技有限公司產(chǎn)品，超低溫冰箱為中科美菱低溫科技股份有限公司產(chǎn)品。

1.3 試驗方法

1.3.1 意蜂幼蟲飼養(yǎng) 中蜂幼蟲的飼料配方參照王倩等的方法[12]，其中將D-果糖和D-葡萄糖換為新鮮蜂蜜。預(yù)試驗中空白對照組中蜂幼蟲7日齡成活率達到90%以上，說明本實驗室建立的中華蜜蜂幼蟲人工飼養(yǎng)模型可為本研究提供符合要求的中蜂幼蟲。從福建農(nóng)林大學(xué)蜂學(xué)學(xué)院蜂場選擇20群健康蜂群(無白堊病癥狀且蜜蜂球囊菌特異性引物常規(guī)PCR檢測為陰性)。用移蟲針挑取2日齡幼蟲，放入無菌的24孔細(xì)胞培養(yǎng)板(每孔對應(yīng)1頭幼蟲，孔內(nèi)加有35 ℃預(yù)溫的幼蟲飼料)，將24孔板放入恒溫恒濕培養(yǎng)箱，35 ℃、70% RH條件下培養(yǎng)。每隔24 h吸去舊飼料，加入新飼料。3日齡時，處理組飼喂含有球囊菌孢子的人工飼料(終濃度為1×107孢子/mL)，空白對照組飼喂正常人工飼料。適宜條件下飼喂幼蟲至7日齡，處理組和對照組均設(shè)置3個生物學(xué)重復(fù)。

1.3.2 測序樣品準(zhǔn)備 腸道組織是球囊菌孢子萌發(fā)的主要場所，也是蜜蜂抵御病原入侵的重要免疫器官，故選擇幼蟲腸道作為測序材料。分別剖取處理組和對照組4、5、6日齡幼蟲腸道，為盡量減少腸道RNA的降解，將從解剖取樣到樣品放入液氮速凍的時間控制在30 s以內(nèi)，每剖取1組樣品，立即液氮速凍并迅速放入-80 ℃超低溫冰箱保存，通過Illumina Hiseq 2500平臺對處理組和對照組的12個幼蟲腸道樣品進行深度測序，其中處理組樣品包括4日齡幼蟲腸道(AmT1-1、AmT1-2、AmT1-3)、5日齡幼蟲腸道(AmT2-1、AmT2-2、AmT2-3)和6日齡幼蟲腸道(AmT3-1、AmT3-2、AmT3-3)，對照組為4日齡幼蟲腸道(AmCK-1、AmCK-2、AmCK-3)。

1.3.3 RNA抽提和cDNA合成 液氮冷卻研磨，按照RNA抽提試劑盒的說明書進行RNA抽提，最后將RNA溶解于200 μL RNase-free水中。取1 μL溶解的RNA，用Evolution 600分光光度計(Thermo Scientific公司，美國)定量確定濃度。

1.3.4 mRNA純化與cDNA合成 抽提的總RNA首先利用10 U DNase Ⅰ在37 ℃條件下消化1 h，然后利用Oligotex mRNA Kits Midi進行mRNA純化，純化后的mRNA用分光光度計進行定量。以10 μg mRNA作為模板，GsuI-oligo dT作為反轉(zhuǎn)錄引物，用1 000 U SuperscriptⅡ逆轉(zhuǎn)錄酶在42 ℃下孵育1 h合成第1鏈cDNA；隨后利用高碘酸鈉氧化mRNA的5′端帽子結(jié)構(gòu)，并連接生物素；通過Dynal M280磁珠篩選連接了生物素的mRNA/cDNA，并通過堿裂解釋放第1鏈cDNA；然后通過DNA連接酶在第1鏈cDNA的5′末端加上接頭，利用ExTaq聚合酶合成第2鏈cDNA。最后通過GsuⅠ酶切去除polyA和5′端接頭。

1.3.5 cDNA測序 將3、4、5、6日齡意蜂幼蟲中腸組織的cDNA利用705型Fisher Scientific超聲波破碎儀打斷至 300～500 bp，利用Ampure beads(Illumina公司，美國)進行純化。隨后純化的cDNA利用TruSeqTMDNA Sample Prep Kit -Set A(Illumina公司，美國)構(gòu)建文庫，并利用TruSeq PE Cluster Kit(Illumina公司，美國)進行擴增。最后Illumina Hiseq 2500測序儀進行測序反應(yīng)。測序數(shù)據(jù)已上傳NCBI SRA數(shù)據(jù)庫，SRA Num：SRA456722。

1.3.6 基因表達豐度統(tǒng)計 通過Illumina PE125測序，獲得大量原始讀段，去除含接頭、含N比例大于10%的和低質(zhì)量的讀段，得到清晰讀段，進而利用Tophat軟件將清晰讀段比對至意蜂的參考基因組(Amel_4.5)，表達量的計算使用FPKM法[10]，其計算公式為：FPKM=(106C×103)/(NL)。

1.4 高表達基因的GO分類及KEGG代謝通路富集分析

根據(jù)基因的FPKM值，選取各腸道樣品表達量水平最高的前100個基因，利用在線分析工具OmicsShare(http：//www.omicshare.com/tools/index.php/Home/Index/index.html)進行GO分類及KEGG代謝通路富集分析。

2 結(jié)果與分析

2.1 RNA seq測序數(shù)據(jù)概述

對12個腸道樣品進行深度測序，獲得的讀段數(shù)都在 26 405 020 條以上，去除低質(zhì)量讀段后的清晰讀段數(shù)達到 26 405 020 條以上，Q20達到97%以上，說明RNA-seq數(shù)據(jù)質(zhì)量較好。首先將清晰讀段比對核糖體數(shù)據(jù)庫，最多允許5個錯配，去除比對上核糖體的讀段，剩余讀段用于比對意蜂參考基因組(表1、表2)。各樣品比對上的比例均在86%以上，說明測序數(shù)據(jù)真實可靠。并對各腸道樣品進行主成分分析(PCA)(圖1)，結(jié)果顯示處理組和對照組樣品的聚類情況良好，說明測序樣品具有較好的重復(fù)性。

2.2 各腸道樣品前100個高表達基因的GO分類

一般認(rèn)為FPKM 值在0.10～3.75之間的是低表達基因，在>3.75～15之間為中度表達基因，大于15的為高表達基因，根據(jù)FPKM值對各樣品基因進行排序，選取前100位的高表達基因進行下一步分析。

表1 RNA-seq數(shù)據(jù)過濾信息統(tǒng)計

注：括號內(nèi)為占比，表2同。

表2 Clean reads與參考基因組的比對信息統(tǒng)計

將意蜂幼蟲腸道表達水平最高的100個基因進行GO注釋(圖2)，這100個高表達基因主要集中在細(xì)胞組分、生物學(xué)進程和分子功能上發(fā)揮作用，AmCK、AmT1和AmT2的高表達基因GO分類數(shù)均為29個，AmT3的分類數(shù)有27個，在細(xì)胞組分方面，基因主要富集在細(xì)胞整體、細(xì)胞部分、胞外區(qū)、大分子復(fù)合物、隔膜、隔膜部分、細(xì)胞器和細(xì)胞器部分；在生物學(xué)進程方面，基因主要富集在細(xì)胞進程、行為、生物調(diào)節(jié)、細(xì)胞組分起源、發(fā)育進程、代謝進程、定位、多組織進程、多細(xì)胞組織進程、繁殖、繁殖進程、刺激反應(yīng)、單組織進程；在分子功能方面，基因主要富集在抗氧化活性、電子載體活性、結(jié)合(Binding)、催化活性、結(jié)構(gòu)分子活性。此外，富集在細(xì)胞組分中的基因比率高于在分子功能、生物學(xué)進程中的基因比率。

AmT1與AmCK比較，同是4日齡，高表達基因富集的GO分類數(shù)量相同，且分別富集在細(xì)胞組分、生物學(xué)進程和分子功能上的基因數(shù)量基本一致，說明在脅迫的早期階段，球囊菌尚未喚起意蜂幼蟲的應(yīng)答，AmT1、AmT2、AmT3之間比較，高表達基因富集的GO分類數(shù)量逐漸下降，且分別富集在細(xì)胞組分、生物學(xué)進程和分子功能上的基因數(shù)量也呈逐漸下降的趨勢，說明隨著時間延長，球囊菌脅迫增強，意蜂幼蟲的應(yīng)答也逐漸增強。

同時，GO富集分析結(jié)果顯示，隨著脅迫時間的延長，富集在代謝過程和催化活性上的基因數(shù)量呈下降趨勢，表明隨著球囊菌的增殖破壞了幼蟲腸道的正常代謝功能。

2.3 各腸道樣品前100個高表達基因的KEGG代謝通路富集分析

對樣品AmCK、AmT1、AmT2、AmT3表達水平最高的前100個基因進行KEGG在線分析，結(jié)果顯示這些基因分別富集在94、93、93、83個代謝通路上(圖3)。統(tǒng)計分析發(fā)現(xiàn)這些高表達基因主要集中在新陳代謝、有機體系統(tǒng)、遺傳信息加工、環(huán)境信息加工等6個途徑上。

2.3.1 新陳代謝 根據(jù)KEGG直系同源號進行分類，AmCK新陳代謝相關(guān)基因涉及脂代謝、碳水化合物代謝、其他氨基酸代謝、氨基酸代謝、能量代謝、萜類和酮類化合物的代謝、外來物質(zhì)的降解和代謝等8個代謝途徑。分析結(jié)果顯示，細(xì)胞色素C氧化酶Ⅰ、線粒體ATP合酶脂質(zhì)結(jié)合蛋白、細(xì)胞色素b、細(xì)胞色素C氧化酶Ⅱ(線粒體)4個編碼蛋白參與能量代謝，酮脂酰CoA硫解酶、谷胱甘肽S-轉(zhuǎn)移酶S1、類視黃醛脫氫酶參與脂類代謝過程、α-葡萄糖苷酶前體、山梨糖醇脫氫酶、果糖二磷酸醛縮酶參與碳代謝過程。精氨酸激酶同工酶、酮脂酰CoA硫解酶、類視黃醛脫氫酶參與氨基酸代謝過程，在細(xì)胞色素P450通路上的谷胱甘肽S-轉(zhuǎn)移酶S1基因、谷胱甘肽S-轉(zhuǎn)移酶D1型X2基因共同作用，細(xì)胞色素P450是一種含高鐵血紅素的蛋白，作為一種末端氧化酶，細(xì)胞色素P450主要催化機體內(nèi)源和外源性物質(zhì)在體內(nèi)的氧化反應(yīng)，參與大部分藥物的氧化代謝[13-14]。可以得到多個基因可以參與多個通路。

2.3.2 遺傳信息加工 遺傳信息加工通路中主要涉及轉(zhuǎn)錄、翻譯、折疊、組裝和降解，其中參與翻譯通路高表達基因占據(jù)一個很大的比例，其中主要有編碼40S核糖體蛋白基因和60S核糖體蛋白，總共有37個基因編碼這2類蛋白，翻譯通路中的RNA轉(zhuǎn)運的延伸因子1-α同工酶基因，該基因參與許多重要的細(xì)胞過程和疾病，包括信號傳導(dǎo)、翻譯控制、凋亡、細(xì)胞骨架組成、病毒復(fù)制等[15]。熱休克蛋白在涉及轉(zhuǎn)錄、轉(zhuǎn)移、折疊、組裝和降解過程，該基因參與新生蛋白折疊、組裝加工成有生物學(xué)活性的蛋白質(zhì)。

2.3.3 環(huán)境信息加工 在KEGG分析中，環(huán)境信息加工過程具有2個方面，分別是信號轉(zhuǎn)導(dǎo)和信號分子及其相互作用。40S核糖體蛋白S6、 肌動蛋白相關(guān)蛋白等參與相關(guān)信號轉(zhuǎn)導(dǎo)表達。

2.3.4 細(xì)胞進程 細(xì)胞進程中高表達基因主要富集在轉(zhuǎn)運和分解代謝、細(xì)胞運動、細(xì)胞通訊等3個方面。天冬氨酸蛋白酶、GL12416、NPC2蛋白質(zhì)、肌動蛋白、熱休克蛋白等參與轉(zhuǎn)運和分解代謝。通路中GL12416基因產(chǎn)物參與吞噬體的形成以及細(xì)胞通訊過程中的間隙連接，相關(guān)基因還參與溶酶體的形成和作用。肌動蛋白參與吞噬體通路過程，是細(xì)胞的一種重要骨架蛋白。同時肌動蛋白細(xì)胞分泌、吞噬、移動、胞質(zhì)流動和胞質(zhì)分離等細(xì)胞通訊過程中起重要作用。

2.3.5 有機體系統(tǒng) 有機體系統(tǒng)主要包括免疫系統(tǒng)、循環(huán)系統(tǒng)、消化系統(tǒng)等9個系統(tǒng)。超氧化物歧化酶、熱休克蛋白等蛋白參與有機體系統(tǒng)。熱休克蛋白參與MPKM通路，絲裂原活化蛋白激酶(MAPK)是生物體內(nèi)重要的信號轉(zhuǎn)導(dǎo)系統(tǒng)之一，參與介導(dǎo)生長、發(fā)育、分裂、分化、死亡以及細(xì)胞間的功能同步等多種細(xì)胞過程，MAPK能被多種炎性刺激激活，并對炎癥的發(fā)生、發(fā)展起重要調(diào)控作用[16-17]。

2.4 對照組與處理組間比較

意蜂幼蟲腸道AmCK、AmT1、AmT2之間的差異性不大，高表達基因富集在新陳代謝、有機體系統(tǒng)、遺傳信息加工、環(huán)境信息加工和細(xì)胞進程上的代謝通路數(shù)基本一致。AmT3相對于AmCK具有一些特異性基因，特異性表達出一些抗菌肽、防御素前體蛋白物質(zhì)等，抗菌肽指昆蟲體內(nèi)經(jīng)誘導(dǎo)而產(chǎn)生的一類具有抗菌活性的堿性多肽物質(zhì)，抗菌肽對部分真菌、原蟲、病毒及癌細(xì)胞等均具有強有力的殺傷作用[18-19]，說明隨著球囊菌脅迫時間的增加，幼蟲腸道誘導(dǎo)產(chǎn)生抗菌肽。

2.5 處理組間比較

在球囊菌脅迫前期(AmT1、AmT2)，高表達基因富集在新陳代謝、有機體系統(tǒng)、遺傳信息加工、環(huán)境信息加工和細(xì)胞進程上的代謝通路數(shù)基本一致，而在AmT3脅迫后期，高表達基因富集在新陳代謝上的代謝通路數(shù)大幅下降，并且免疫性物質(zhì)增加，再次說明球囊菌對意蜂幼蟲的生長代謝產(chǎn)生了嚴(yán)重影響，防御機制開始發(fā)揮作用。

KEGG代謝通路富集分析結(jié)果(表3至表6)顯示，AmCK、AmT1、AmT2和AmT3中基因富集量最多的都為核糖體通路，但是富集的基因數(shù)量呈逐漸下降的趨勢，說明球囊菌破壞了宿主的蛋白合成系統(tǒng)，推測球囊菌能通過互作減少宿主的營養(yǎng)供給，從而促進自身增殖。

表3 AmCK組基因富集量前10位KEGG代謝通路統(tǒng)計

表4 AmT1樣品基因富集量最多的前10個KEGG代謝通路

表5 AmT2樣品基因富集量最多的前10個KEGG代謝通路

2.6 各腸道樣品前100個高表達基因的Venn分析

對處理組和對照各幼蟲腸道樣品的Venn分析結(jié)果顯示， 72個高表達基因為AmCK、AmT1、AmT、AmT3所共有，各樣品中的特有高表達基因分別為1、2、2、13個，特有的高表達基因大多是編碼抗菌肽、防御素前體、鐵蛋白、抑制劑等免疫方面的基因，推測這些共有基因在球囊菌侵染的全過程都發(fā)揮著重要功能，而少數(shù)特有基因則在病程的各個階段分別發(fā)揮作用，對這些共有與特有高表達基因的深入研究將有助于揭示意蜂幼蟲的應(yīng)答機制(圖4)。

表6 AmT3樣品基因富集量最多的前10個KEGG代謝通路

3 討論

本研究選取意蜂幼蟲腸道作為測序材料，因為腸道是球囊菌侵染的主要場所，其轉(zhuǎn)錄組變化能更準(zhǔn)確地反映意蜂幼蟲對球囊菌入侵的應(yīng)答。為了從全局水平了解意蜂幼蟲受球囊菌脅迫前后的基因表達譜，本研究利用RNA-seq技術(shù)對正常幼蟲腸道和球囊菌脅迫的幼蟲腸道進行深度測序，通過比對意蜂參考基因組，共發(fā)現(xiàn)有13 592個基因在幼蟲腸道中表達，其中已知基因有11 918個(87.68%)，另有197個新基因，這些新基因?qū)橥晟埔夥浠蚪M信息提供一定補充。

Venn分析結(jié)果顯示，4個腸道樣品中的共有高表達基因有72個，比如部分參與核糖體通路的基因，合成加工成供幼蟲自身生長發(fā)育需要的物質(zhì)，也有參與核糖體通路的U-box結(jié)構(gòu)域蛋白基因，可作為分子伴侶或輔分子伴侶，還可能作為剪接體的部分因子，識別細(xì)胞內(nèi)錯誤折疊或異常蛋白質(zhì)，從而使其正確折疊或降解。另外，U-BOX蛋白質(zhì)還與癌病變以及自身免疫性疾病等病理過程有關(guān)，該類蛋白質(zhì)通常還在將底物轉(zhuǎn)運給蛋白體的過程中起作用，調(diào)節(jié)著體內(nèi)蛋白質(zhì)的降解[20]。各樣品中的特有高表達基因分別有1、2、2、13個，對比處理組樣品和對照組樣品中的特有高表達基因，將有可能是意蜂幼蟲響應(yīng)球囊菌脅迫的關(guān)鍵應(yīng)答基因，基因編碼的核苷酸結(jié)合蛋白β亞基1在AmT3上特異性表達，該蛋白特異性結(jié)合抗原，可以起到抗體的作用，在抵御侵染上起到一定作用。對比處理組各樣品中的特有高表達基因，或許能為明確球囊菌脅迫過程不同階段的關(guān)鍵應(yīng)答基因提供重要線索，這些特有高表達基因值得進一步研究。

由于經(jīng)費限制，本研究只設(shè)4日齡正常幼蟲腸道作為對照組，會導(dǎo)致一些基因表達信息的缺失，因此，對4、5、6日齡正常幼蟲腸道分別取樣并同時進行深度測序，便于更全面地獲得意蜂幼蟲響應(yīng)球囊菌脅迫的基因表達譜信息。

本研究中，在進行RNA-seq數(shù)據(jù)分析時，發(fā)現(xiàn)有部分reads不能定位于已有的意蜂參考基因組，它們有可能是尚未克隆的新基因，這些新基因也許在意蜂幼蟲腸道發(fā)育中發(fā)揮重要作用。本研究利用RNA-seq技術(shù)對正常和球囊菌脅迫意蜂幼蟲腸道進行深度測序，通過對各樣品表達量水平最高的前100個基因進行GO分類和KEGG代謝通路富集分析，初步解析了意蜂幼蟲響應(yīng)球囊菌脅迫的應(yīng)答，推測在6日齡之后免疫機制迅速增強，應(yīng)答機制的深入闡釋有賴于差異表達基因的進一步研究。