基于RNA-seq技術分析意大利蜜蜂幼蟲腸道響應球囊菌脅迫的高表達基因

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

(1.福建農林大學蜂學學院，福建福州 350002； 2.江西省養蜂研究所，江西南昌 330200)

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

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

1 材料與方法

1.1 生物材料

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

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產品，cDNA文庫構建試劑盒TruSeqTMDNA Sample Prep Kit -Set A為美國Illumina公司產品。

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

1.3 試驗方法

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

1.3.2 測序樣品準備 腸道組織是球囊菌孢子萌發的主要場所，也是蜜蜂抵御病原入侵的重要免疫器官，故選擇幼蟲腸道作為測序材料。分別剖取處理組和對照組4、5、6日齡幼蟲腸道，為盡量減少腸道RNA的降解，將從解剖取樣到樣品放入液氮速凍的時間控制在30 s以內，每剖取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作為反轉錄引物，用1 000 U SuperscriptⅡ逆轉錄酶在42 ℃下孵育1 h合成第1鏈cDNA；隨后利用高碘酸鈉氧化mRNA的5′端帽子結構，并連接生物素；通過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公司，美國)構建文庫，并利用TruSeq PE Cluster Kit(Illumina公司，美國)進行擴增。最后Illumina Hiseq 2500測序儀進行測序反應。測序數據已上傳NCBI SRA數據庫，SRA Num：SRA456722。

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

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

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

2 結果與分析

2.1 RNA seq測序數據概述

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

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

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

表1 RNA-seq數據過濾信息統計

注：括號內為占比，表2同。

表2 Clean reads與參考基因組的比對信息統計

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

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

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

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

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

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

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

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

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

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

2.4 對照組與處理組間比較

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

2.5 處理組間比較

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

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

表3 AmCK組基因富集量前10位KEGG代謝通路統計

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

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

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

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

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

3 討論

本研究選取意蜂幼蟲腸道作為測序材料，因為腸道是球囊菌侵染的主要場所，其轉錄組變化能更準確地反映意蜂幼蟲對球囊菌入侵的應答。為了從全局水平了解意蜂幼蟲受球囊菌脅迫前后的基因表達譜，本研究利用RNA-seq技術對正常幼蟲腸道和球囊菌脅迫的幼蟲腸道進行深度測序，通過比對意蜂參考基因組，共發現有13 592個基因在幼蟲腸道中表達，其中已知基因有11 918個(87.68%)，另有197個新基因，這些新基因將為完善意蜂基因組信息提供一定補充。

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

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

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