陳大福，郭睿，熊翠玲，梁勤，鄭燕珍，徐細建，張曌楠，黃枳腱，張璐，王鴻權，解彥玲，童新宇

?

中華蜜蜂幼蟲腸道響應球囊菌早期脅迫的轉(zhuǎn)錄組學

陳大福，郭睿，熊翠玲，梁勤，鄭燕珍，徐細建，張曌楠，黃枳腱，張璐，王鴻權，解彥玲，童新宇

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

【目的】白堊病是困擾養(yǎng)蜂生產(chǎn)的頑疾。目前，尚無利用二代測序技術研究中華蜜蜂（,簡稱中蜂）幼蟲白堊病的報道。本研究利用RNA-seq技術對健康（AcCK）及球囊菌()脅迫的中蜂4日齡幼蟲腸道（AcT）進行深度測序，在轉(zhuǎn)錄組水平研究中蜂幼蟲在球囊菌脅迫早期的脅迫應答?！痉椒ā客ㄟ^Illumina HiSeq 2500平臺對AcCK和AcT進行雙端（PE125）測序，首先對測序數(shù)據(jù)進行質(zhì)控和評估，利用edgeR軟件進行差異表達基因（DEG）分析，進而對DEGs進行GO富集分析及KEGG代謝通路富集分析，最后，利用實時熒光定量PCR（qRT-PCR）驗證測序數(shù)據(jù)的可靠性?！窘Y(jié)果】 AcCK和AcT的轉(zhuǎn)錄組測序共得到188 457 338條原始讀段（raw reads）, 經(jīng)過濾得到182 088 448條有效讀段（clean reads），兩端Q20與Q30均在97.96%和94.97%以上，說明測序數(shù)據(jù)質(zhì)量良好；主成分分析（PCA）結(jié)果顯示第一與第二主成分可分別解釋基因表達總體差異的75.8%和10.7%；DEG分析結(jié)果顯示上調(diào)基因與下調(diào)基因的數(shù)量分別為344和239個；GO富集分析結(jié)果顯示DEGs共富集在36個GO分類（term）上，其中基因富集數(shù)最多的細胞（106 unigenes）、細胞組件（106 unigenes）和代謝進程（104 unigenes）；KEGG代謝通路富集分析結(jié)果顯示上調(diào)與下調(diào)基因分別富集在72個和45個代謝通路上，其中上調(diào)基因富集數(shù)最多的是核糖體（72 unigenes）、碳代謝（16 unigenes）和糖酵解（14 unigenes），而下調(diào)基因富集數(shù)最多的是碳代謝（9 unigenes）、二羧酸代謝（8 unigenes）和氨基酸生物合成（7 unigenes），進一步分析表明中蜂幼蟲腸道的部分細胞免疫對球囊菌的脅迫產(chǎn)生應答，而體液免疫不產(chǎn)生應答，宿主的代謝相關基因受到球囊菌的顯著抑制?！窘Y(jié)論】揭示了中蜂幼蟲在球囊菌入侵早期的脅迫應答，為深入解析中蜂幼蟲的脅迫應答機制提供了重要信息，也為在分子水平研究關鍵應答基因打下了基礎。

中華蜜蜂；幼蟲腸道；球囊菌；轉(zhuǎn)錄組；RNA-seq

0 引言

【研究意義】蜜蜂作為一種社會學模式昆蟲，對發(fā)育學、神經(jīng)生物學、行為學和病原-宿主互作等研究有重要價值[1-5]。蜜蜂也是最重要的授粉昆蟲，在農(nóng)業(yè)生產(chǎn)、糧食安全和生態(tài)維持中也發(fā)揮著巨大作用[6]。白堊病是球囊菌（）特異性侵染蜜蜂幼蟲的致死性真菌病，長期困擾養(yǎng)蜂生產(chǎn)，每年給養(yǎng)蜂業(yè)造成巨大損失。近年來，隨著蜂產(chǎn)品全球貿(mào)易的快速發(fā)展，白堊病呈逐年上升趨勢[7]。有關白堊病的研究主要集中于西方蜜蜂（）幼蟲，但關于西方蜜蜂幼蟲響應球囊菌脅迫的應答研究報道極少，目前尚無球囊菌侵染東方蜜蜂（）幼蟲過程中的應答機制研究，也無防治白堊病的有效方法。【前人研究進展】國內(nèi)養(yǎng)蜂生產(chǎn)中使用的主要蜂種是意大利蜜蜂（，簡稱意蜂）和中華蜜蜂（，簡稱中蜂），二者分別屬于西方蜜蜂和東方蜜蜂。意蜂極易暴發(fā)白堊病，但中蜂極少受到該病的影響。近20年來，國內(nèi)外學者對白堊病開展了一系列研究，主要集中在病原分類鑒定[8-9]、形態(tài)學[10-11]、病理學[12-13]、流行病學[14]、侵染過程[15-16]、蜜蜂防御[17-18]以及疾病防治[19]等方面。筆者課題組也在球囊菌的病理和檢測方面開展了較為系統(tǒng)的研究[20-24]。2006年，基因組信息的公布[25]為其分子生物學研究奠定了重要基礎。Evans等[26]在蜜蜂細胞內(nèi)鑒定出了大多數(shù)NF-B信號通路中的基因，如的2個同系物，但都不與直系同源，表明是僅存在于果蠅的特化分支而不存在于其他昆蟲；Aronstein等[27]利用cDNA-AFLP對健康及球囊菌感染西方蜜蜂幼蟲進行了測序比較，發(fā)現(xiàn)差異表達基因（DEGs）參與了能量代謝和蛋白轉(zhuǎn)運，其中的類幾丁質(zhì)編碼基因很可能參與了宿主對球囊菌的抵抗；Cornman等[28]利用RNA-seq技術對幼蟲芽孢桿菌（）感染后12 h和72 h的蜜蜂幼蟲進行測序，鑒定出75個顯著上調(diào)和6個顯著下調(diào)基因，其中若干編碼抗菌肽基因和2個編碼圍食膜基質(zhì)基因顯著上調(diào)；Julie等[29]通過RNA-seq研究了東方蜜蜂微孢子蟲（）和殺蟲劑單獨或共同飼喂意蜂后中腸的轉(zhuǎn)錄組變化，發(fā)現(xiàn)和殺蟲劑并無協(xié)同效應，長期暴露于殺蟲劑環(huán)境抑制了意蜂的免疫基因的表達。2015年，韓國研究人員測序并公布了的基因組[30]，但當時并未公布基因位置和功能注釋信息。為深入開展中蜂的轉(zhuǎn)錄組學研究，筆者課題組前期已組裝并注釋了中蜂幼蟲腸道的參考轉(zhuǎn)錄組[31]?！颈狙芯壳腥朦c】養(yǎng)蜂生產(chǎn)中，偶爾可見中蜂幼蟲罹患白堊病的現(xiàn)象，前期研究發(fā)現(xiàn)中蜂幼蟲可在實驗室條件下被球囊菌侵染。白堊病研究基本都集中于球囊菌侵染西方蜜蜂幼蟲，有關球囊菌侵染中蜂幼蟲的相關研究幾乎沒有。本研究利用RNA-seq技術對健康及球囊菌脅迫的中蜂幼蟲腸道進行深度測序，進而在轉(zhuǎn)錄組水平研究宿主球囊菌入侵早期的脅迫應答?！緮M解決的關鍵問題】以健康及球囊菌脅迫的中蜂4日齡幼蟲腸道作為研究對象，通過Illumina測序技術研究中蜂幼蟲腸道在球囊菌入侵早期的脅迫應答，為深入解析中蜂幼蟲響應球囊菌脅迫的應答機制及關鍵應答基因的功能研究提供重要信息。

1 材料與方法

試驗于2015年12月至2016年8月在福建農(nóng)林大學蜂學學院蜜蜂保護學實驗室進行。

1.1 材料

中蜂幼蟲取自福建農(nóng)林大學蜂學學院教學蜂場，球囊菌菌株由福建農(nóng)林大學蜂學學院蜜蜂保護學實驗室保存并活化。

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

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

1.2 方法

1.2.1 球囊菌活化、孢子純化及計數(shù) 將4℃保存的球囊菌培養(yǎng)皿置于已紫外滅菌的超凈臺，用酒精燈上灼燒片刻的鑷子夾取少量菌絲在平板中央2 cm左右圓心區(qū)域內(nèi)劃線操作，蓋上培養(yǎng)皿，置于37℃生化箱恒溫培養(yǎng)，3 d后觀察球囊菌生長情況。接種8—10 d以后，待培養(yǎng)皿上黑色孢子較多時，刮取孢子至干凈的EP管中，用研磨棒充分研磨，按照Jensen等[32]的方法離心純化球囊菌孢子，梯度稀釋后用血球計數(shù)板進行計數(shù)。

1.2.2 中蜂幼蟲的人工飼養(yǎng)及腸道樣品準備 按照已報道的飼料配方[33]配制中蜂幼蟲飼料并進行改良，預試驗結(jié)果顯示中蜂幼蟲的7日齡成活率可達70%以上。從學院教學蜂場群勢較強的健康中蜂蜂群（無白堊病癥狀且PCR鑒定為陰性）移取2日齡幼蟲至無菌的24孔細胞培養(yǎng)板（每孔對應1只幼蟲，孔內(nèi)加有35℃預溫的幼蟲飼料），35℃，相對濕度70%條件下飼養(yǎng)。每隔24 h更換飼料。處理組3日齡幼蟲飼喂含球囊菌孢子的人工飼料（孢子終濃度1×107個/mL），對照組3日齡幼蟲飼喂正常飼料。飼喂孢子后24 h，分別剖取處理組（AcT）和對照組（AcCK）4日齡幼蟲腸道組織，每剖取一只幼蟲腸道，迅速將腸道移至RNA Free的EP管，再投入液氮速凍，待一組腸道樣品（7只幼蟲腸道）集齊后，迅速轉(zhuǎn)移保存于-80℃。本試驗進行3次生物學重復。AcCK與AcT的3個生物學重復分別為AcCK-1、AcCK-2、AcCK-3和AcT-1、AcT-2、AcT-3。因幼蟲腸道是球囊菌的寄生場所，而本文旨在對中蜂幼蟲在球囊菌脅迫早期的應答進行研究，故選取4日齡幼蟲腸道作為測序?qū)ο蟆?/p>

1.2.3 cDNA文庫構(gòu)建及Illumina測序 利用RNAiso Reagent試劑盒抽提處理組和對照組中蜂幼蟲腸道組織的總RNA，然后用RNase-free DNaseI去除基因組DNA殘留。RNA的質(zhì)量通過瓊脂糖凝膠電泳和NanoDrop ND-1000（NanoDrop，Wilmington，DE，USA）進行檢測。利用Oligotex mRNA Kits Midi試劑盒說明書，純化各樣品總RNA中的mRNA。以10 μg mRNA作為模板，GsuI-oligo dT作為反轉(zhuǎn)錄引物，用1000 U Superscript II reverse transcriptase在42℃下孵育1 h合成第1鏈cDNA；隨后利用高碘酸鈉氧化mRNA的5′端帽子結(jié)構(gòu)，并連接生物素；通過Dynal M280磁珠篩選連接了生物素的mRNA/cDNA，并通過堿裂解釋放第1鏈cDNA；然后通過DNA ligase在第1鏈cDNA的5′末端加上接頭，利用Ex Taq polymerase 合成第2鏈cDNA。最后，通過I酶切去除polyA和5′端接頭。利用Ampure beads對上述cDNA進行純化，cDNA文庫通過TruSeqTMDNA Sample Prep Kit-Set A進行構(gòu)建和TruSeq PE Cluster Kit進行擴增。上述6個幼蟲腸道樣品委托廣州基迪奧生物科技有限公司進行測序，測序平臺為Illumina HiSeq2500。本研究測得的轉(zhuǎn)錄組數(shù)據(jù)已上傳美國國家生物技術信息中心（NCBI）SRA數(shù)據(jù)庫，SRA號：SRA456721。

1.2.4 數(shù)據(jù)分析 對于下機數(shù)據(jù)，利用Perl腳本去除含有adaptor、未知核苷酸比例大于5%和低質(zhì)量reads，獲得有效讀段（clean reads）。利用R軟件（version 2.16.2）進行測序飽和度分析。使用短reads比對工具bowtie[34]將clean reads比對（mapping）到核糖體數(shù)據(jù)庫（最多允許5個錯配），去除mapping上核糖體的reads，將保留下來的數(shù)據(jù)用于轉(zhuǎn)錄組的組裝及分析，進而利用SOAP aligner/soap2軟件[35]將未mapping上核糖體的reads mapping到筆者課題組組裝的中蜂幼蟲腸道參考轉(zhuǎn)錄組[31]。

利用FPKM（Fragments Per Kilobase of Transcript Per Million Mapped Reads）法計算基因表達量。利用R軟件（version 2.16.2）計算各樣品之間的相關性系數(shù)。利用edgeR軟件[36]進行DEG分析。DEG的篩選標準為FDR≤0.05且|log2fold change≥1。將DEGs向GO數(shù)據(jù)庫（http://www.geneontology.org/）的各條目（term）映射，并計算每個term的基因數(shù)，從而得到具有某個GO功能的基因列表及基因數(shù)目統(tǒng)計，然后應用超幾何檢驗，找出與整個基因組背景相比，在DEGs中顯著富集的GO term。KEGG代謝通路顯著性富集分析以KEGG pathway為單位，應用超幾何檢驗，找出與整個基因組背景相比在差異表達基因中顯著性富集的代謝通路。

1.2.5 實時熒光定量PCR（qRT-PCR）驗證 為了驗證RNA-seq數(shù)據(jù)，隨機選取5個DEGs進行qRT-PCR。qRT-PCR反應按照SYBR Green Dye試劑盒（Vazyme公司，中國）操作說明書進行，每個反應進行3次重復。20 μL的反應體系中含SYBR Green Dye 10 μL，10.0 μmol·L-1正、反向引物各1 μL，cDNA模板DNA 1 μL，DEPC水8 μL。qRT-PCR反應在ABI 7500熒光定量PCR儀（ABI公司，美國）上進行，反應條件：95℃預變性1 min，95℃變性15 s，60℃延伸30 s，共40個循環(huán)，最后72℃延伸45 s。所選基因的相對表達量采用2-??Ct法[37]計算。

2 結(jié)果

2.1 Illumina測序數(shù)據(jù)質(zhì)控與評估

中蜂幼蟲腸道樣品的Illumina PE125測序共測得188 457 338條reads，經(jīng)過濾得到182 088 448條clean reads，各樣品clean reads數(shù)均在95.4%以上，兩端Q20均在97.96%以上，兩端Q30均在94.97%以上（表1）。飽和度分析結(jié)果顯示隨著測序量的增多，檢測到的基因數(shù)也隨之上升、增長速度趨于平緩，說明檢測到的基因數(shù)趨于飽和（附圖1）。上述結(jié)果表明本研究中的測序數(shù)據(jù)質(zhì)量良好，可用于進一步分析。

Pearson相關性分析結(jié)果顯示AcCK與AcT的組內(nèi)各生物學重復之間的相關性均接近1，說明樣本的重復性較高（附圖2）。進一步對AcCK與AcT的各樣品進行主成分分析（PCA），結(jié)果顯示第一主成分（PC1）與第二主成分（PC2）共能解釋樣品基因表達總體差異的86.5%，AcCK與AcT的組內(nèi)各生物學重復聚類良好（圖1），表明AcCK與AcT的總體基因表達模式差別顯著。

RNA-seq數(shù)據(jù)比對中蜂幼蟲腸道參考轉(zhuǎn)錄組情況統(tǒng)計顯示，各樣品clean reads比對上參考轉(zhuǎn)錄組的比例均在88.91%以上（表2），AcCK與AcT的表達基因數(shù)目為36 604和37 108，分別占中蜂幼蟲腸道參考轉(zhuǎn)錄組的84.05%和85.21 %。

圖1 各中蜂幼蟲腸道樣本的PCA分析

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

表2 RNA-seq數(shù)據(jù)比對參考基因組情況統(tǒng)計

2.2 差異表達基因（DEG）分析

利用edgeR進行DEG分析，結(jié)果顯示在球囊菌脅迫后，4日齡幼蟲腸道共有583個DEGs，其中上調(diào)基因與下調(diào)基因的數(shù)量分別為344和239個（圖2），說明有相當一部分基因受球囊菌脅迫而被激活，也有部分基因受到球囊菌的抑制。

2.3 差異表達基因的GO富集分析

DEGs的GO分類結(jié)果顯示，這些DEGs分為3類：生物學進程、細胞組分和分子功能，分布于36個GO term上，分別為代謝進程、生長、發(fā)育進程、生殖、生殖進程、運動、細胞成分組織或生物起源、免疫系統(tǒng)進程、多組織進程、定位、單組織進程、細胞進程、應激反應、信號、生物調(diào)控、結(jié)構(gòu)分子活性、電子載體活性、抗氧化活性、脒基-核苷酸交換因子活性、蛋白結(jié)合轉(zhuǎn)錄因子活性、催化活性、分子功能調(diào)節(jié)器、核酸結(jié)合轉(zhuǎn)錄因子活性、轉(zhuǎn)運器活性、結(jié)合、分子傳感器活性、細胞、細胞組件、大分子復合物、細胞器組件、細胞器、突觸、細胞膜內(nèi)腔、細胞膜和細胞膜組件（圖3）。GO富集分析結(jié)果顯示，基因富集數(shù)最多的GO term為細胞（106 unigenes），其次為細胞組件（106 unigenes）和代謝進程（104 unigenes）。

2.4 差異表達基因的KEGG代謝通路富集分析

分別對上調(diào)基因和下調(diào)基因進行KEGG代謝通路富集分析，結(jié)果表明上調(diào)基因富集在72個代謝通路上，其中基因富集數(shù)最多的是核糖體（ribosome）（72 unigenes），其次為碳代謝（carbon metabolism）（16 unigenes）和糖酵解（glycolysis/gluconeogenesis）（14 unigenes）（圖4-A），而下調(diào)基因富集在45個代謝通路上，其中基因富集數(shù)最多的是碳代謝（9 unigenes），其次為二羧酸代謝（glyoxylate and dicarboxylate metabolism）（8 unigenes）和氨基酸生物合成（biosynthesis of amino acids）（7 unigenes）（圖4-B）。值得關注的是，上調(diào)基因中分別有3、2、1和1個基因富集在內(nèi)吞作用（endocytosis）、吞噬體（phagosome）、MAPK信號通路（MAPK signaling pathway）和泛素介導的蛋白水解（ubiquitin mediated proteolysis），說明意蜂幼蟲的免疫相關通路被球囊菌脅迫所激活；多達32個下調(diào)基因富集在新陳代謝相關通路，如碳代謝（9 unigenes）、氧化磷酸化（oxidative phosphorylation）（5 unigenes）、精氨酸生物合成（arginine biosynthesis）（3 unigenes）、氮代謝（nitrogen metabolism）（2 unigenes）和丙酮酸代謝（pyruvate metabolism）（2 unigenes），說明對中蜂幼蟲的新陳代謝相關通路受到球囊菌的抑制，球囊菌入侵對宿主的代謝系統(tǒng)產(chǎn)生廣泛的影響。

圖2 球囊菌脅迫中蜂幼蟲腸道的DEGs

1：生物調(diào)控 Biological regulation；2：細胞成分組織或生物合成 Cellular component organization or biosynthesis；3：細胞進程 Cellular process；4：發(fā)育進程 Developmental process；5：生長 Growth；6：免疫系統(tǒng)進程 Immune system process；7：定位 Localization；8：運動 Locomotion；9：代謝進程 Metabolic process；10：多組織進程 Multi-organism process；11：多細胞生物進程 Multicellular organismal process；12：生殖 Reproduction；13：生殖進程 Reproductive process；14：應激反應 Response to stimulus；15：信號Signaling；16：單組織進程 Single-organism process；17：細胞 Cell；18：細胞組件 Cell part；19：大分子復合物 Macromolecular complex；20：細胞膜 Membrane；21：細胞膜組件 Membrane part；22：細胞膜內(nèi)腔 Membrane-enclosed lumen；23：細胞器Organelle；24：細胞器組件 Organelle part；25：突觸 Synapse；26：抗氧化活性 Antioxidant activity；27：結(jié)合 Binding；28：催化活性 Catalytic activity；29：電子載體活性 Electron carrier activity；30：脒基核苷酸交換因子活性 Guanyl-nucleotide exchange factor activity；31：分子功能調(diào)節(jié)器 Molecular function regulator；32：分子轉(zhuǎn)換器活性 Molecular transducer activity；33：核酸結(jié)合轉(zhuǎn)錄因子活性 Nucleic acid binding transcription factor activity；34：蛋白結(jié)合轉(zhuǎn)錄因子活性 Protein binding transcription factor activity；35：結(jié)構(gòu)分子活性 Structural molecule activity；36：轉(zhuǎn)運器活性 Transporter activity

A：上調(diào)基因Up-regulated genes；B：下調(diào)基因Down-regulated genes

進而對上述免疫相關基因和代謝相關基因的表達模式進行分析，結(jié)果顯示球囊菌脅迫后，中蜂幼蟲腸道的免疫相關基因表達量呈不同程度的上調(diào)（圖5-A），而代謝相關基因表達量不同程度地下調(diào)（圖5-B），進一步表明中蜂幼蟲腸道的免疫相關基因受球囊菌脅迫而激活、代謝相關基因受到球囊菌的抑制。

2.5 RNA-seq數(shù)據(jù)的qRT-PCR驗證

利用qRT-PCR檢測隨機選取的5個DEGs（4個上調(diào)基因，1個下調(diào)基因），結(jié)果顯示這些DEGs的表達水平變化趨勢與RNA-seq數(shù)據(jù)中的基因表達水平變化趨勢一致（圖6），證明了本研究轉(zhuǎn)錄組數(shù)據(jù)的可靠性。

A：中蜂幼蟲腸道的免疫相關基因A. c. cerana larval gut’s immune-related DEGs；B：中蜂幼蟲腸道的代謝相關基因A. c. cerana larval gut’s metabolism-related DEGs

A：類圍食膜因子3-C前體表皮蛋白編碼基因cuticular protein analogous to peritrophins 3-C precursor encoding gene；B：表皮蛋白1前體編碼基因cuticular protein 1 precursor encoding gene；C：類表皮蛋白同種型X2編碼基因cuticle protein-like isoform X2 encoding gene；D：類表皮蛋白38編碼基因cuticle protein 38-like encoding gene；E：G蛋白耦聯(lián)受體Mth-like 3編碼基因probable G-protein coupled receptor Mth-like 3 encoding gene

3 討論

白堊病對蜂群的危害嚴重，每年給養(yǎng)蜂業(yè)造成巨大損失。目前，白堊病的組學研究十分有限。前期研究中，筆者課題組組裝了中蜂幼蟲腸道的參考轉(zhuǎn)錄組，并對其進行了功能及代謝通路注釋[31]。在此基礎上，本研究利用RNA-seq技術對健康及球囊菌脅迫的中蜂4日齡幼蟲腸道進行深度測序，進而研究其響應球囊菌脅迫的應答。腸道是昆蟲的重要器官，在消化、發(fā)育和免疫等諸多方面中發(fā)揮重要作用。本研究針對中蜂幼蟲腸道進行測序，其轉(zhuǎn)錄組變化能更為精確地反映宿主響應病原脅迫的應答。在AcCK和AcT中分別檢測到已知基因36 604和37 108個，約占中蜂幼蟲腸道轉(zhuǎn)錄組基因總數(shù)的84.05%和85.21%，推測部分未檢測到的基因可能在中蜂4日齡幼蟲腸道不表達。

養(yǎng)蜂生產(chǎn)中，意蜂幼蟲極易被球囊菌侵染而暴發(fā)白堊病，而中蜂幼蟲極少受其影響。前期研究發(fā)現(xiàn)在球囊菌脅迫后的意蜂幼蟲腸道中共檢測到24個DEGs，其中上調(diào)基因與下調(diào)基因分別為4和20個（未發(fā)表數(shù)據(jù)）。本研究發(fā)現(xiàn)球囊菌脅迫后的中蜂幼蟲腸道有583個基因差異表達（上調(diào)基因344個，下調(diào)基因239個），遠多于球囊菌脅迫后的意蜂幼蟲腸道，這些DEGs或許與中蜂的球囊菌抗性密切相關。此前，大多數(shù)研究都集中于蜜蜂的群體防御，通常認為清理行為在蜜蜂抵御球囊菌入侵過程中至關重要。前期試驗結(jié)果表明，實驗室條件下對中蜂幼蟲和意蜂幼蟲接種球囊菌，兩種幼蟲皆能被侵染而發(fā)生白堊病，但前者的發(fā)病率遠低于后者，說明在中蜂幼蟲與意蜂幼蟲在個體水平上對球囊菌也存在較大的抗性差異。角質(zhì)層和圍食膜是昆蟲免疫防御的第一道防線[38]，當病原微生物突破這道防線后，將遭遇昆蟲體內(nèi)細胞免疫與體液免疫的抵抗，包括內(nèi)吞作用、黑化作用、吞噬作用、蛋白的酶促水解及抗菌肽等[39-40]。內(nèi)吞作用和吞噬作用在蜜蜂抵御真菌侵染的過程中發(fā)揮著重要作用[40]。本研究發(fā)現(xiàn)分別有3個和2個上調(diào)基因富集在內(nèi)吞作用和吞噬體，說明這兩個代謝通路受球囊菌脅迫而激活。泛素介導的蛋白水解對清除衰老和損傷細胞意義重大，同時也在調(diào)節(jié)細胞進程的諸多方面發(fā)揮著關鍵作用，包括細胞周期、分化、發(fā)育和免疫等[41]。本研究發(fā)現(xiàn)1個上調(diào)基因富集在泛素介導的蛋白水解，說明該通路同樣被球囊菌所激活。上述結(jié)果表明中蜂幼蟲的部分細胞免疫對球囊菌的脅迫產(chǎn)生應答，推測它們在球囊菌入侵的早期發(fā)揮重要作用。昆蟲的體液免疫也會因病原微生物的入侵而被激活，從而合成分泌抗菌肽殺滅病原，但在本研究中，未發(fā)現(xiàn)有DEGs富集在體液免疫如NF-B和Jak-STAT信號通路，也未發(fā)現(xiàn)蜜蜂4種抗菌肽（abaecin、apidaecin、defensin和hymenoptaecin）編碼基因的差異表達，表明中蜂幼蟲的體液免疫在球囊菌脅迫早期尚未被球囊菌激活。碳代謝對生物體內(nèi)氨基酸的合成與轉(zhuǎn)化至關重要。本研究發(fā)現(xiàn)在球囊菌脅迫的中蜂4日齡幼蟲腸道中，富集在碳代謝上的DEGs包括16個上調(diào)基因和9個下調(diào)基因，推測多數(shù)碳代謝相關基因上調(diào)表達可能與宿主的細胞免疫水平增強有關，而部分數(shù)碳代謝相關基因表達水平下調(diào)反映出球囊菌-中蜂幼蟲間互作的復雜性。本研究還發(fā)現(xiàn)多達72個上調(diào)基因富集在核糖體，說明在球囊菌脅迫早期，宿主的蛋白合成活躍，推測宿主通過提升蛋白合成滿足抵御球囊菌入侵及自我修復的需要。

中蜂幼蟲在球囊菌入侵早期并不表現(xiàn)出明顯的白堊病癥狀，但宿主和病原之間發(fā)生著復雜的互作，本研究發(fā)現(xiàn)有多達76個DEGs富集在新陳代謝（32個新陳代謝相關通路），本研究的結(jié)果表明此時中蜂幼蟲的代謝系統(tǒng)已受到球囊菌一定程度地干擾和破壞，推測病原通過抑制宿主的新陳代謝減少營養(yǎng)物質(zhì)的消耗、降低宿主的免疫應答水平，從而促進病原自身的侵染。在球囊菌入侵早期，這種干擾和破壞作用的累積，可能會導致入侵晚期白堊病的發(fā)生。此前的白堊病研究集中于意蜂且聚焦在病程晚期，忽略了病程早期的研究。本研究利用RNA-seq技術開展中蜂幼蟲腸道響應球囊菌脅迫的轉(zhuǎn)錄組學研究，發(fā)現(xiàn)在球囊菌脅迫的早期階段，中蜂幼蟲腸道也伴隨著復雜的應答。目前尚無一種殺真菌劑被批準應用于養(yǎng)蜂生產(chǎn)[1]，養(yǎng)蜂生產(chǎn)實踐中主要通過選育抗病品系、改善養(yǎng)蜂管理和保持清潔衛(wèi)生來防治白堊病[41]，但是效果并不理想，亟需有效的防治策略。未來的研究方向是利用RNA-seq技術對球囊菌脅迫后的不同日齡中蜂幼蟲腸道進行深度測序，全局研究中蜂幼蟲在轉(zhuǎn)錄組水平的應答，并應用趨勢分析或WGCNA分析篩選出關鍵應答基因，進而利用RNAi等技術手段驗證其功能。

4 結(jié)論

利用RNA-seq技術對健康及球囊菌脅迫的中蜂4日齡幼蟲腸道進行深度測序，通過對DEGs的GO和KEGG富集分析，發(fā)現(xiàn)在球囊菌脅迫早期，中蜂幼蟲的部分細胞免疫被球囊菌激活，而體液免疫未被激活，宿主的新陳代謝系統(tǒng)受到球囊菌的顯著抑制。研究結(jié)果為深入解析中蜂幼蟲響應球囊菌脅迫的應答機制提供了重要信息，也為白堊病的有效防治打下了基礎。

References

[1] GALIZIA C G, EISENHARDT D, GIURFA M, MENZEL R.. Dordrecht Netherlands. New York: Springer, 2012.

[2] BEGNA D, Han B, FENG M, FANG Y, LI J. Differential expression of nuclear proteomes between honeybee (L.) queen and worker larvae: a deep insight into caste pathway decisions., 2012, 11(2): 1317-1329.

[3] ZAYED A, ROBINSON G. Understanding the relationship between brain gene expression and social behavior: lessons from the honey bee., 2012, 46(6): 591-615.

[4] FORET S, KUCHARSKI R, PELLEGRINI M, FENG S H, JACOBSEN S E, ROBINSON G E, MALESZKA R. DNA methylation dynamics, metabolic fluxes, gene splicing, and alternative phenotypes in honey bees., 2012, 109(13): 4968-4973.

[5] KURZE C, ROUTTU J, MORITZ R F. Parasite resistance and tolerance in honeybees at the individual and social level., 2016, 119(4): 290-297.

[6] Committee on the Status of Pollinators in North America.. National Academies Press, 2007.

[7] AIZEN M A, GARIBALDI L A, CUNNINGHAM S A, KLEIN A M. How much does agriculture depend on pollinators? Lessons from long-term trends in crop production., 2009, 103(9): 1579-1588.

[8] SPILTOIR C F. Life cycle of., 1955, 42(6): 501-518.

[9] LUMBSCH H T, Huhndorf S M. Outline of Ascomycota., 2007, 13: 1-58.

[10] BISSETT J. Contribution toward a monograph of the genus., 1988, 66(12): 2541-2560.

[11] CHORBI?SKI P. Enzymatic activity of strains of., 2003, 59(11): 1019-1022.

[12] BAMFORD S, HEATH L A F. The effects of temperature and pH on the germination of spores of the chalkbrood fungus,., 1989, 28(1): 36-40.

[13] WINSTON M L. The biology of the honey bee//. Harvard University Press, Cambridge, USA, 1991.

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

[15] BAILEY L.. Academic Press, London, UK, 1991.

[16] THEANTANA T, CHANTAWANNAKUL P. Protease and beta-N- acetylglucosaminidase of honey bee chalkbrood pathogen., 2008, 47(1): 68-76.

[17] EVANS J D, SPIVAK M. Socialized medicine: individual and communal disease barriers in honey bees.,2010, 103(Suppl. 1): S62-S72.

[18] TANJI T, HU X, WEBER A N, IP Y. Toll and IMD pathways synergistically activate an innate immune response in., 2007, 27(12): 4578-4588.

[19] HORNITZKY M. Literature review of chalkbrood. A report for the RIRDC. Publication No. 01/150, Kingston, ACT, AU. 2001.

[20] 梁勤, 陳大福, 王建鼎. 營養(yǎng)生態(tài)條件對蜜蜂球囊菌生長及產(chǎn)孢的影響. 中國生態(tài)農(nóng)業(yè)學報,2001, 9(4): 31-34.

LIANG Q, CHEN D F, WANG J D. Effects on the mycelia growth and spore-forming ofunder ecological condition of nutrients.,2001, 9(4): 31-34. (in Chinese)

[21] 李江紅, 鄭志陽, 陳大福, 梁勤. 影響蜜蜂球囊菌侵染蜜蜂幼蟲的因素及侵染過程觀察. 昆蟲學報, 2012, 55(7): 790-797.

LI J H, ZHENG Z Y, CHEN D F, LIANG Q. Factors influencinginfection on honeybee larvae and observation on the infection process., 2012, 55(7): 790-797. (in Chinese)

[22] 熊翠玲, 陳大福, 付中民, 馬曉云, 梁勤. 蜜蜂球囊菌檢測分子標記的靈敏性測定. 中國蜂業(yè), 2010, 61(11): 16-18.

XIONG C L, CHEN D F, FU Z M, MA X Y, LIANG Q. Determination of the sensitivity of DNA molecular marker for detectings., 2010, 61(11): 16-18. (in Chinese)

[23] 席偉軍, 李江紅, 陳大福, 梁勤. 環(huán)介導等溫擴增(LAMP)技術檢測蜜蜂球囊菌. 中國農(nóng)業(yè)科學, 2016, 49(4): 765-774.

XI W J, LI J H, CHEN D F, LIANG Q. Diagnosis of theby the loop-mediated isothermal amplification., 2016, 49(4): 765-774. (in Chinese)

[24] 鄭志陽, 李江紅, 梁勤, 陳大福. 蜜蜂球囊菌分泌多種胞外酶侵染蜜蜂幼蟲. 福建農(nóng)林大學學報 (自然科學版), 2011, 40(3): 280-284.

ZHENG Z Y, LI J H, LIANG Q, CHEN D F.secretes multiple extracellular enzymes to infect honeybee larvae., 2011, 40(3): 280-284. (in Chinese)

[25] Honeybee Genome Sequencing Consortium. Insights into social insects from the genome of the honeybee., 2006, 443(7114): 931-949.

[26] EVANS J D, ARONSTEIN K A, CHEN Y P, HETRU C, IMLER J L, JIANG H, KANOST M, THOMPSON G J, ZOU Z, HULTMARK D. Immune pathways and defense mechanisms in honey bees,., 2006, 15(5): 645-656.

[27] ARONSTEIN K A, MURRAY K D, SALDIVAR E. Transcriptional responses in honey bee larvae infected with chalkbrood fungus., 2010, 11: 391.

[28] CORNMAN R S, LOPEZ D, EVANS J D. Transcriptional response of honey bee larvae infected with the bacterial pathogen., 2013, 8(6): e65424.

[29] JULIE A, BARBARA M A, BERNARD V, CATHERINE T, FRéDéRIC D, NICOLAS B. Transcriptome analyses of the honeybee response toand insecticides., 2014, 9(3): e91686.

[30] PARK D, JUNG J W, CHOI B S, JAYAKODI M, LEE J, LIM J, YU Y, CHOI Y S, LEE M L, PARK Y, CHOI I Y, YANG T J, EDWARDS O R, NAH G, KWON H W. Uncovering the novel characteristics of Asian honey bee,, by whole genome sequencing., 2015, 16: 1.

[31] 徐細建, 郭睿, 駱群, 熊翠玲, 梁勤, 張串聯(lián), 鄭燕珍, 張曌楠, 黃枳腱, 張璐, 李汶東, 陳大福. 中華蜜蜂幼蟲腸道參考轉(zhuǎn)錄組的組裝及SSR分子標記鑒定. 中國農(nóng)業(yè)科學, 2017, 50(6): 1157-1166.

Xu X J, Guo R, Luo Q, Xiong C L, Liang Q, Zhang C L, Zheng Y Z, Zhang Z N, Huang Z J, Zhang L, Li W D, Chen D F.transcriptome assembly forlarval gut and identification of SSR molecular markers., 2017, 50(6): 1157-1166. (in Chinese)

[32] 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): 79-88.

[33] 王倩, 孫亮先, 肖培新, 劉鋒, 康明江, 胥保華. 室內(nèi)人工培育中華蜜蜂幼蟲技術研究. 山東農(nóng)業(yè)科學, 2009(11): 113-116.

WANG Q, SUN L X, XIAO P X, LIU F, KANG M J, XU B H. Study on technology for indoor artificial feeding oflarvae., 2009(11): 113-116. (in Chinese)

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

[35] HURGOBIN B. Short read alignment using SOAP2//. New York, NY , United States: Humana Press, 2016: 241-252.

[36] 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.

[37] LIVAK K J, SCHMITTGEN T D. Analysis of relative gene expression data using real-time quantitative PCR and the 2-??Ctmethod., 2001, 25: 402-408.

[38] ORIHEL T C. The peritrophic membrane: its role as a barrier to infection of the arthropod host//Maramorosch K, Shope R E.. Academic Press, New York, USA, 1975: 67-73.

[39] GLI?SKI Z, JAROSZ J. Infection and immunity in the honey bee., 2001, 36: 12-24.

[40] GLINSKI Z, BUCZEK K. Response of theto fungal infections., 2003, 38: 183-189.

[41] McBride W H, Iwamoto K S, Syljuasen R, Pervan M, Pajonk F. The role of the ubiquitin/proteasome system in cellular responses to radiation., 2003, 22: 5755-5773.

（責任編輯 岳梅）

A: AcCK-1; B: AcCK-2; C: AcCK-3; D: AcT-1; E: AcT-2; F: AcT-3

附圖1 各中蜂幼蟲腸道樣品的測序飽和度

Fig. S1 Sequencing saturation of everylarval gut sample

附圖2 各中蜂幼蟲腸道樣品不同生物學重復間的相關性

Fig. S2 Pearson correlation between every two biological repeats within eachlarval gut sample

Transcriptome oflarval gut Underthe Stressof

CHEN DaFu, GUO Rui, XIONG CuiLing, LIANG qin, ZHENG YanZhen, XU XiJian, ZHANG ZhaoNan, HUANG ZhiJian, ZHANG Lu, WANG HongQuan, XIE YanLing, TONG XinYu

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

【Objective】So far, no study on application of next-generation sequencing technology for the research of chalkbrood disease was reported. In the present research, RNA-seq technology was utilized to deep sequence of normal and-treated 4th instarlarval gut to mine larvae’s responses tochallenge.【Method】AcCK (un-treated group) and AcT (-treated group) were sequenced on Illumina HiSeq 2500 platform. After evaluation and filtration of raw data from RNA-seq, differentially expressed gene (DEG) analysis was performed using edgeR software, further, Gene Ontology (GO) and KEGG pathway enrichment analyses were carried out, and finally, real-time quantitative PCR (qRT-PCR) was conducted to validate the RNA-seq data.【Result】 In total, RNA-seq produced 188 457 338 raw reads, and after filtration, 182 088 448 clean reads were obtained, Q20 and Q30 of each sample were above 97.96% and 94.97%, respectively, indicating that RNA-seq data in this research were with high quality. principle component analysis (PCA) was performed on all genes level and the result showed PC1 and PC2 were able to account for 75.8% and 10.7% of the expressed genes’ overall differences, respectively. DEG analysis result displayed that there were 344 up-regulated genes and 239 down-regulated genes in AcCK VS AcT. GO enrichment analysis result showed that the DEGs were enriched in 36 GO terms, among them, the mostly enriched ones were cell (106 unigenes), cell part (106 unigenes) and metabolic process (104 unigenes). KEGG pathway enrichment analysis result suggested that up- and down-regulated genes were enriched in 72 and 45 pathways, respectively, and the mostly enriched pathways for up-regulated genes were ribosome (72 unigenes), carbon metabolism (16 unigenes) and glycolysis/gluconeogenesis (14 unigenes), while the mostly enriched pathways for down-regulated genes were carbon metabolism (9 unigenes), glyoxylate and dicarboxylate metabolism (8 unigenes) and amino acids biosynthesis (7 unigenes). further analysis demonstrated that the immune-related genes inlarval gut were activated, while the metabolism-related genes were greatly inhibited.【Conclusion】The findings of the study not only uncovered thelarval gut’s responses toduring the early stage of invasion, but also provided key information for clarifying the mechanism underlying the host’s responses to, thus laying a foundation for functional investigation of key responding genes.

; larval gut;; transcriptome; RNA-seq

2016-12-26；接受日期：2017-03-06

國家自然科學基金（30800806）、國家現(xiàn)代農(nóng)業(yè)產(chǎn)業(yè)技術體系（蜜蜂）建設專項（CARS-45-KXJ7）、福建農(nóng)林大學科技發(fā)展資金（KF2015123）

陳大福，Tel：0591-83726835；E-mail：dfchen826@fafu.edu.cn。通信作者郭睿，Tel：0591-87640197；E-mail：ruiguo@fafu.edu.cn