中華蜜蜂工蜂中腸響應東方蜜蜂微孢子蟲脅迫的高表達基因分析

付中民， 周丁丁， 陳華枝 ， 耿四海 ， 陳大福 ， 鄭燕珍 ， 熊翠玲 ， 徐國鈞 ， 張 曦 ， 郭 睿

(1. 福建農林大學動物科學學院蜂學學院，福州 350002； 2. 棗莊職業學院，棗莊 277100)

1 引 言

蜜蜂是自然界中最重要的授粉昆蟲，能顯著提高作物、果蔬的產量和品質[1]，具有無可比擬的經濟和生態價值.蜜蜂也是一種高度進化的社會性昆蟲，在社會行為學和種群遺傳學等研究領域應用廣泛[2].中華蜜蜂(Apisceranacerana，簡稱中蜂)是東方蜜蜂(Apiscerana)的指名亞種，也是適應我國本土自然環境的特有蜂種，具有善于利用零星蜜源、采集力強及抗病性強等優勢[3].

微孢子蟲是專性胞內寄生的真菌，宿主范圍包括哺乳動物、魚類和昆蟲等脊椎動物和無脊椎動物[4].微孢子蟲對家蠶和蜜蜂等經濟昆蟲的影響尤為嚴重[5].東方蜜蜂微孢子蟲(Nosemaceranae)特異性侵染成年蜜蜂的中腸上皮細胞，可導致蜜蜂免疫力下降，工蜂離巢飛行提前，壽命縮減，以及產蜜量下降，嚴重危害蜜蜂健康[6-7].隨著東方蜜蜂微孢子蟲基因組(大小：7.86 Mb)的公布和東方蜜蜂基因組(大小：228.332 Mb)完善，二者互作的組學研究不斷深入[8-10].較多的研究結果表明，N.ceranae是導致蜂群崩潰綜合癥的關鍵因素之一[11-12].前人對N.ceranae與西方蜜蜂(Apismellifera)的相互作用進行了較多的組學研究，取得了一些重要進展.Aufauvre等[13]利用二代測序技術對N.ceranae、殺蟲劑脅迫及二者共同脅迫1 d和7 d的西方蜜蜂的中腸樣品進行轉錄組測序，分析發現7 d時宿主的中腸免疫調節系統受到了強烈抑制，表皮防御功能以及海藻糖代謝通路受到了顯著影響.Badaoui等[14]利用RNA-seq技術研究了N.ceranae感染5、10和15 d的西方蜜蜂的基因表達譜，發現差異表達基因可歸入3種表達模式，進一步分析發現N.ceranae感染對宿主的氨基酸代謝產生干擾及抗菌肽基因的表達產生抑制，表明N.ceranae感染抑制了西方蜜蜂的免疫系統.但有關N.ceranae-東方蜜蜂互作的組學研究較為滯后，相關信息匱乏.

筆者所在課題組前期對N.ceranae脅迫意大利蜜蜂(Apismelliferaligustica，簡稱意蜂)工蜂過程中宿主和病原的高表達基因(highly expressed gene, HEG)進行全面分析，揭示了二者的HEG的表達譜及潛在作用[15].本研究利用Illumina測序技術對正常及N.ceranae脅迫的中蜂7 d和10 d工蜂中腸進行測序，通過生物信息學方法對宿主的HEG進行深入分析，以期在轉錄組水平探究中蜂對N.ceranae的脅迫應答.研究結果可揭示HEG在中蜂脅迫應答中的作用，為探明背后的分子機制提供有益的信息和線索.

2 材料與方法

2.1 實驗材料

中蜂工蜂取自福建農林大學動物科學學院蜂學學院教學蜂場；N.ceranae孢子由福建農林大學動物科學學院蜂學學院蜜蜂保護學實驗室保存并純化[16].

2.2 方 法

2.2.1 孢子純化、接種感染及人工飼養N.ceranae孢子的粗提和純化方法如下：(1)從N.ceranae感染嚴重的意蜂蜂群巢門口抓取外勤蜂若干只，拉取中腸放入研缽加水充分研磨，差速離心后棄上清液，純水懸浮沉淀，繼續多次差速離心以去除雜質，將得到的孢子粗提液，按照1∶10比例加入無菌水后分裝到1.5 mL滅菌的EP管內，4 ℃，3000 r/min離心5 min，吸去上清液后將5管渾濁液合為1管，4 ℃，1 200 r/min離心3 min；(2)分別配制25%、50%、75%、100%的percoll純化液.各濃度的percoll液取200 μL按照密度從高到低的濃度依次分裝到干凈的EP管中，最后加入離心后的渾濁液，4 ℃，14 000 r/min離心30 min；(3)將無菌的注射器小心吸取孢子，針頭插入乳白色的孢子帶中，注意不要吸到孢子帶上下方的液體.

參照本實驗室前期已建立的方法[17]進行中蜂工蜂的接種感染和人工飼養：(1)配制含孢子的50%(w/v)糖水，孢子濃度為2×108/mL；(2)選取群勢較強的成熟封蓋子脾迅速提至實驗室，放入34±0.5 ℃，(60%～70%)RH培養箱，將剛出房的工蜂放入滅菌的塑料盒(塑料盒四周打孔以通風，30只/盒)，將裝有50%(w/v)無菌糖水的飼喂器插入每個盒子上方；(3)將工蜂出房當天記為0 d，34±0.5 ℃培養24 h；將上述中蜂工蜂單頭固定后饑餓處理2 h，處理組工蜂單只飼喂含孢子的糖水5 μL，對照組工蜂單只飼喂不含孢子的糖水5 μL；(4)24 h后均飼喂不含孢子的糖水，每日更換糖水，及時清理未食盡糖水的工蜂及死去的工蜂，食盡糖水的工蜂用于后續實驗，用50%(w/v)的糖水飼喂工蜂，適宜條件飼喂至10 d.處理組和對照組均包含3個生物學重復，每組重復包含9只工蜂中腸.適宜條件飼喂至7 d時，處理組中的9只工蜂中腸，分別放入3個EP管，該混合樣品記為Ac7T(Ac7T-1、Ac7T-2、Ac7T-3)；對照組的9只工蜂中腸，分別放入3個EP管，該混合樣品記為Ac7CK(Ac7CK-1、Ac7CK-2、Ac7CK-3)；每次取樣在20 s內迅速放于液氮后，迅速轉移到-80℃超低溫冰箱保存備用.10 d的處理和對照組中腸樣品也按照上述方法進行取樣和保存，分別記為Ac10T(Ac10T-1、Ac10T-2、Ac10T-3)和Ac10CK(Ac10CK-1、Ac10CK-2、Ac10CK-3).

2.2.2 RNA提取、cDNA建庫及RNA-seq 首先利用RNAiso Reagent試劑盒分別抽提上述12個中腸組織樣品的總RNA，然后用RNase-Free DNase I去除殘留的基因組DNA.接著利用Oligotex mRNA Kits Midi對抽提的總RNA進行RNA純化，通過分光光度計對純化RNA進行定量檢測.RNA的質量通過瓊脂糖凝膠電泳和NanoDrop ND-1000( Nano Drop公司，美國)進行檢測.參照陳大福等[18]的方法進行cDNA文庫構建.委托廣州基迪奧生物技術有限公司，采用Illumina Hiseq 4000測序平臺進行雙端測序.

2.2.3 測序數據質控、HEG篩選及生物信息學分析 對于測序得到的原始讀段(raw reads)，去除含接頭(adapter)、含未知堿基(N)的比例大于10%、低質量的讀段，用Tophat軟件將過濾得到的有效讀段(clean reads)映射核糖體數據庫，進而將未比對上的clean reads映射東方蜜蜂基因組(assembly ACSNU-2.0)，將比對上的clean reads用于后續的數據分析.

利用FPKM(Fragment Per Kilobase of per Million mapped reads)算法計算基因表達量，公式為FPKM =Total exon fragment/[Mapped fragment (millions) × exon length (KB)]，按照FPKM >15的標準從所有表達基因中篩選出HEG[19].利用在線工具平臺OmicsShare(http://www.omicshare.com/tools/index.php/Home/Index/index.html)對篩選出的共有及特有HEG進行GO分類及KEGG pathway富集分析等相關生物信息學分析.

表1 RT-PCR的引物信息

2.2.4 HEG的RT-PCR驗證 為驗證本研究RNA-seq的數據的可靠性，從4組中腸樣品的共有HEG中隨機選取8個HEG進行RT-PCR驗證.根據上述8個基因的序列，利用Primer Premier 5設計特異性引物，委托福州擎科生物技術有限公司進行引物合成，相關引物序列信息見表1.利用試劑盒AxyPrepTMMultisource Total RNA Miniprep(AXYGEN公司，中國)抽提上述4組中腸樣品的總RNA.以上述總RNA作為模板，利用試劑盒HiScript 1st Strand cDNA Synthsis Kit(Vazyme公司，中國)進行反轉錄，反應步驟按照說明書進行.PCR反應在T100熱循環儀(Bio-Rad公司，美國)上進行，反應體系包括：Mixture 10 μL，ddH2O 7 μL，上下游引物各1 μL，cDNA模板1 μL.反應條件：94 ℃預變性5 min；94 ℃變性50 s，60 ℃退火30 s，72 ℃延伸1 min，共34個循環；最后72 ℃延伸10 min，4 ℃保存.PCR產物經1.5%的瓊脂糖凝膠電泳檢測，核酸凝膠成像儀(上海培清，中國)觀察拍照.

3 結果與分析

3.1 高通量測序數據概述

本文中蜂工蜂中腸樣品共測得1 809 736 786條raw reads，經嚴格質控得到1 562 162 742條clean reads，將未比對上核糖體數據庫的clean reads映射到東方蜜蜂基因組，Ac7CK、Ac7T、Ac10CK及Ac10T各組的平均比對率為75.78%、55.01%、78.13%和44.19%；Q30均達93.34%及以上(表2).此外，Ac7CK、Ac7T、Ac10CK及Ac10T的Pearson相關系數都在0.872 8及以上(圖1).上述結果說明測序樣品的重復性較好，測序數據的質量可滿足進一步分析需要.

表2 RNA-seq數據信息統計

圖1 中蜂工蜂中腸各組樣品內不同生物學重復間的Pearson相關性

Fig.1 Pearson correlations among different biological repeats within eachA.c.ceranaworker’s midgut sample group

3.2 HEG的Venn分析

根據FPKM >15的標準，分別從Ac7CK、Ac7T、Ac10CK及Ac10T中篩選出3 162、3 311、3 305和2 463個HEG.Venn分析結果顯示四組樣品的共有HEG為2 074個，特有HEG分別為89、283、156和78個(圖2).

圖2 中蜂工蜂中腸各組樣品HEG的Venn分析

Fig.2 Venn investigation of HEGs in everyA.c.ceranaworker midgut sample group

3.3 共有及特有HEG的GO分類

GO富集分析結果顯示，共有HEG主要涉及生物學進程、細胞組分和分子功能三類，分布于45個GO條目.生物進程中，富集基因數最多的是代謝進程(409)、細胞進程(380)和單組織進程(314)(圖3A)；細胞組分中，富集基因數最多的是細胞(254)、細胞組件(254)和細胞膜(153)(圖3B)；分子功能中，富集基因數最多的是結合(385)、催化活性(329)和轉運器活性(60)(圖3C).

Ac7CK特有HEG涉及37個GO條目，基因富集數最多的是細胞進程(204)、代謝進程(196)、結合(191)、催化活性(181)和單組織進程(165)；Ac7T特有HEG涉及35個GO條目，基因富集數最多的是細胞進程(219)、代謝進程(210)、結合(193)、催化活性(189)和單組織進程(183)；Ac10CK特有HEG涉及37個GO條目，基因富集數最多的是細胞進程(230)、代謝進程(223)、結合(217)、催化活性(203)和單組織進程(183)；Ac10T特有HEG涉及28個GO條目，基因富集數最多的是細胞進程(71)、結合(63)、單組織進程(59)、代謝進程(58)和催化活性(51).

圖3 共有HEG的GO分類Fig.3 GO categorization of the shared HEGs

3.4 共有及特有HEG的KEGG代謝通路富集分析

KEGG富集分析結果顯示，共有HEG富集在39條通路，其中基因富集數最多的前3位分別是核糖體(90)、氧化磷酸化(63)和蛋白酶體(25)(圖4).

Ac7CK的特有HEG富集在39條通路，富集基因數最多的前4位分別是內質網蛋白加工(33)、氧化磷酸化(27)、核糖體(25)和內吞作用(25)；Ac7T的特有HEG富集在39條通路，富集基因數最多的前4位分別是核糖體(32)、氧化磷酸化(29)、內質網蛋白加工(26)和剪接體(23)；Ac10CK的特有HEG富集在39條通路，富集基因數最多的前4位分別是內質網蛋白加工(38)、氧化磷酸化(32)、核糖體(27)、和碳代謝(26)；Ac10T的特有HEG富集在37條通路，富集基因數最多的前4位分別是內質網蛋白加工(8)、PI3K-Akt信號通路(7)、Wnt信號通路(6)和TGF-β信號通路(5).

圖4 共有HEG的KEGG代謝通路富集分析

圓圈大小代表富集在某一通路的基因數多少，越大代表基因數越多；圓圈顏色代表某一通路的富集基因的顯著性高低，越紅代表顯著性越高

Fig.4 KEGG pathway enrichment analysis of the shared HEGs

The size of the circle indicates the number of enriched genes in a certain pathway, the bigger the more; the color of the circle indicates the significance of enriched genes in a certain pathway, and the redder the color, the higher the more significant

圖5 8個共有HEG的RT-PCR驗證

A:類易化海藻糖轉運體tret1編碼基因; B: 類胰蛋白酶-1編碼基因; C: 類磷脂酶A1編碼基因; D: 精氨酸激酶亞基X1編碼基因; E: 類α胰蛋白酶-3編碼基因; F: 40S核糖體蛋白S15編碼基因; G: 類鋅羧肽酶編碼基因; H: 類生命必需蛋白lethal(2)編碼基因

Fig.5 RT-PCR verification of eight shared HEGs

A:facilitated trehalose transporter Tret1-like encoded gene; B: trypsin-1-like encoded gene; C: phospholipase A1-like encoded gene; D: arginine kinase isoform X1 encoded gene; E: trypsin alpha-3-like encoded gene; F: 40S ribosomal protein S15 encoded gene; G: zinc carboxypeptidase-like encoded gene; H: protein lethal(2) essential for life-like encoded gene

從4個樣品共有的HEG中隨機挑選8個進行RT-PCR驗證，電泳結果顯示8對引物均能擴增出符合預期的目的條帶(圖5)，初步證實了本研究預測的HEG真實表達.

4 討 論

中蜂是我國養蜂生產的主要蜂種之一，經長期進化已高度適應我國的生態環境.中蜂是N.ceranae的原始寄主，二者經過長期的協同進化已相互適應，N.ceranae經跨種侵染已擴散至世界各地的西方蜜蜂蜂群.目前，關于N.ceranae與蜜蜂互作方面的研究主要集中在西方蜜蜂[20-22].前人也對中蜂與N.ceranae的互作進行了一些研究，但組學研究較為滯后，相關信息頗為匱乏.Chaimanee等[21]曾對N.ceranae侵染3、 6和12 d的西方蜜蜂(Apismellifera)工蜂進行研究，通過檢測4種抗菌肽基因(defensin、abaecin、apidaecin和hymenoptaecin)、eater基因和卵黃原蛋白編碼基因(vitellogenin)的表達水平，發現N.ceranae侵染對宿主的體液免疫造成了較強抑制.Sinpoo等[22]對西方蜜蜂和東方蜜蜂進行了蜜蜂微孢子蟲(Nosemaapis)和N.ceranae交叉感染試驗，發現二種蜜蜂的抗菌肽編碼基因和細胞免疫相關基因均顯著上調表達，但東方蜜蜂的上皮細胞表現出更強的體液免疫水平且含有更少的孢子數，作者推測東方蜜蜂對N.apis和N.ceranae可能具有更強的抵抗力.本研究利用RNA-seq技術對正常及N.ceranae感染的中蜂工蜂中腸進行測序，從正常的7和10 d工蜂中腸分別篩選出3 162及3 305個HEG，從N.ceranae脅迫的7和10 d工蜂中腸分別篩選出3 311及2 463個HEG.Ac7CK、Ac7T、Ac10CK和Ac10T的共有HEG數為2 074個，推測這些共有HEG對中蜂工蜂中腸的生長和發育中具有重要作用.此外，Ac7T和Ac10T的特有HEG數分別為283及78個，說明隨著脅迫時間的延長，宿主應答的特有HEG數呈下降趨勢，推測這些HEG在宿主響應N.ceranae脅迫的不同階段發揮特殊作用.本研究中，Ac7CK、Ac7T、Ac10CK及Ac10T各組的平均比對率為75.78%、55.01%、78.13%和44.19%，處理組的平均比對率明顯低于對照組，究其原因，處理組中腸樣品包含大量的病原孢子，其測序數據為宿主和病原的混合數據，因病原的clean reads無法比對上東方蜜蜂的參考基因組，勢必影響宿主clean reads的比對率；而對照組中腸樣品的測序數據中絕大部分都為宿主本身的數據，因而clean reads的比對率較高.進一步分析發現，處理組和對照組測序得到的raw reads，以及質控的clean reads，基本處于同一數量級，因此對后續生物信息學分析的影響很小.

本研究中，Ac7T和Ac10T的特有HEG中分別有219和71個注釋到細胞進程，210和58個注釋到代謝進程，189和51個注釋到催化活性，顯示伴隨N.ceranae的入侵，中蜂工蜂中腸的新陳代謝和細胞生命活動處于不同幅度的激活狀態，但激活程度隨著病原增殖而呈下降趨勢.病原微生物入侵昆蟲時會遭到宿主細胞免疫和體液免疫的聯合抵抗，包括內吞作用、吞噬作用、酶促水解作用、黑化作用以及抗菌肽的合成與釋放等[15].本研究發現，Ac7T和Ac10T的特有HEG中分別有95和68個富集在應激反應，1和1個富集在免疫系統進程，表明中蜂工蜂對N.ceranae產生了免疫應答且應答水平隨脅迫時間的延長有所下降.

本研究中，Ac7T和Ac10T的210和58個特有HEG可注釋到77和39條物質代謝通路，包含碳水化合物代謝、氨基酸代謝、輔助因子與維生素代謝、核苷酸代謝、多糖生物合成與代謝等，進一步說明中蜂工蜂中腸的物質代謝因N.ceranae脅迫而激活，但激活程度隨脅迫時間延長而降低.前期研究中，筆者所在課題組對脅迫中蜂7和10日齡工蜂中腸的N.ceranae的共有HEG進行了代謝通路分析，發現與物質代謝相關的信號通路多達43條，包含碳代謝(15)、糖酵解/糖異生(11)、氧化磷酸化(8)和磷酸戊糖途徑(5)等，表明N.ceranae在侵染過程中的物質和能量代謝較為旺盛(數據未發表).結合本研究中的分析結果，推測N.ceranae在其增殖過程中一方面增強自身物質和能量代謝以滿足核酸和蛋白合成需要，另一方面通過某種操縱策略提高宿主的代謝水平，促使其攝入更多的物質和能量滿足中蜂和N.ceranae的需要，從而利于病原的增殖.Li等[23]對新羽化的成年意蜂工蜂接種感染N.ceranae，通過檢測宿主體內的脂類含量發現宿主的脂質含量隨感染時間的延長而顯著降低，作者推測脂質可能是用于維持宿主生長和病原增殖的首選營養物質.本研究發現，Ac7T和Ac10T中分別有31和8個特有HEG富集到12和6條脂質代謝通路，包括脂肪酸降解、磷脂酰甘油代謝、鞘脂代謝等，表明宿主通過提高脂質的合成來滿足自身和N.ceranae的需要，與前人的研究結果一致.前人研究發現在N.ceranae感染的西方蜜蜂對低濃度的蔗糖溶液可產生伸吻反應，取食更多的蔗糖溶液，氧化磷酸化富集基因的數量增加，表明N.ceranae對西方蜜蜂造成了能量脅迫[24-26].本研究中，Ac7T和Ac10T中分別有29和1個特有HEG富集在氧化磷酸化，表明N.ceranae對宿主造成能量脅迫，與西方蜜蜂的相關研究結果一致.

Toll樣受體是一類富含亮氨酸的重復序列的模式識別受體蛋白質家族，廣泛分布于免疫細胞及某些體細胞的表面，作為宿主防御反應中的重要組成部分以及先天性免疫和獲得性免疫的連接點[27].MAPK屬于絲氨酸-蘇氨酸蛋白激酶，MAPK信號通路是一種包括MAPK、MAPK激酶激酶、MAPK激酶的保守三級激酶模式，共同調節著細胞的生長、分化、增殖及對環境的應激反應、免疫應答、炎癥反應等多種細胞生理及病理的過程[28].NF-κB信號通路通過調節某些抗死亡基因和分子的表達，介導調節巨噬細胞中的炎癥反應及免疫應答[29].本研究中，Ac7T中分別有3、6和1個特有HEG富集在Toll樣受體信號通路、MAPK信號通路和NF-κB信號通路；Ac10T中有1個特有HEG富集在NF-κB信號通路，但未發現有HEG富集在其余兩條體液免疫通路.上述結果表明Toll樣受體等3條信號通路在中蜂工蜂響應N.ceranae脅迫的前期表現活躍，可能發揮重要的免疫防御作用，隨著脅迫時間的延長，宿主的NF-κB信號通路在免疫防御中的角色漸趨關鍵.細胞因子介導的Jak-STAT信號通路參與細胞的增殖、分化、凋亡以及免疫調節等諸多生物學過程.本研究中，Ac10T中有2個特有HEG富集在Jak-STAT信號通路，表明該信號通路被N.ceranae激活.筆者團隊在前期研究中發現對于N.ceranae脅迫10 d的意蜂工蜂中腸，沒有HEG富集Jak-STAT和NF-κB通路[15]，表明中蜂和意蜂響應N.ceranae脅迫的免疫應答存在差異，值得進一步深入探討.N.ceranae感染可抑制西方蜜蜂中腸內穩態和細胞凋亡相關基因的表達水平[15].本研究發現，Ac7T中有2個特有HEG富集在細胞凋亡，Ac10T中無HEG富集在此通路，推測中蜂工蜂中腸細胞可通過提高相關基因的表達促進細胞凋亡，以限制N.ceranae的增殖.

綜上所述，本研究對正常及N.ceranae脅迫的中蜂工蜂中腸HEG進行了細致深入的生物信息學分析，提供了宿主的HEG表達譜，揭示了HEG在宿主響應N.ceranae脅迫中的潛在作用.研究結果為在分子水平深入解析中蜂響應N.ceranae脅迫的應答機制、宿主-病原互作機制提供了有益的信息和線索.