東方蜜蜂微孢子蟲侵染中華蜜蜂過程中的高表達基因表達譜及其潛在作用

熊翠玲， 陳華枝， 耿四海， 周倪紅， 周丁丁， 祝智威， 陳大福，鄭燕珍， 徐國鈞， 張 曦， 郭 睿

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

1 引 言

微孢子蟲是一種真菌病原，廣泛感染人類、哺乳類、兩棲類和昆蟲等[1-3].東方蜜蜂微孢子蟲(Nosema ceranae)專性寄生成年蜜蜂的中腸上皮細胞，由Fries等人于1996年首次在中華蜜蜂(Apis cerana cerana，簡稱中蜂)中發現并報道[4]，隨后該病原迅速傳播蔓延至美國和歐洲地區，目前已擴散到世界各地的西方蜜蜂(Apis mellifera)蜂群[5].被N. ceranae感染的蜜蜂壽命減少1/3[6]，產蜜量降低40%，分泌蜂蠟量減少25%[7]，并伴隨有早熟行為[8].此外，N. ceranae還被認為是誘發蜂群崩潰綜合癥的主要病原之一[9].N. ceranae主要通過糞-口或口-口等途徑在蜜蜂個體間傳播，當孢子進入蜜蜂中腸，在堿性環境的刺激下，孢子發芽并彈射出高度壓縮的極絲，極絲穿入上皮細胞并將有感染性的孢質原漿注入其中；N. ceranae經過增殖形成孢子體，最終導致宿主細胞裂解，釋放出的成熟孢子繼而感染臨近的細胞，或者通過糞便排泄到外界環境中，侵染其他蜜蜂個體[10].

前人對于N. ceranae的研究集中在病原鑒定、流行病學、病理學、侵染機理等方面[10-12].Cornman等人于2009年完成N. ceranae的基因組組裝并公布其序列及注釋信息，為其組學及分子生物學研究提供了基礎.較之蜜蜂，N. ceranae的組學研究進展緩慢，相關信息較為缺乏.Huang[13]等通過對N. ceranae感染1-6d的A. mellifera中腸進行轉錄組測序和分析，發現在感染過程中宿主的細胞凋亡基因表達受到抑制，N. ceranae在一個繁殖周期中有1 122個基因差異表達，并分為4種表達模式.Evans等發現1 545個A. mellifera基因可能受到N. ceranae的5個miRNA的調控，這些miRNA聚集成6個共表達組，6個組中的4個(918個基因)與至少1個miRNA的表達顯著相關，N. ceranae的miRNA傾向于與病原基因正相關而與宿主基因負相關，表明N. ceranae的miRNA可調控自身基因表達也能跨界調控宿主的部分基因表達[14].目前，前人僅對侵染A. mellifera的N. ceranae進行了一些組學研究，而對于侵染東方蜜蜂(Apis cerana)的N. ceranae以及二者互作，還缺乏相關研究報道.近期，筆者團隊結合鏈特異性建庫方法和RNA-seq技術對N. ceranae孢子進行了全轉錄組測序，通過生物信息學方法分別預測出83個長鏈非編碼RNA(long non-coding RNA, lncRNA)和204個環狀RNA(circular RNA, circRNA)，并對它們的數量、種類、結構特征、表達譜和調控網絡進行了系統的分析和驗證[15-16].

A. cerana是N. ceranae的原始寄主，二者經過漫長的協調進化已相互適應.有研究表明N. ceranae對A. mellifera具有更強的毒力[5].因此，對侵染A. cerana的N. ceranae進行組學研究，可為解析N. ceranae的侵染機制和病原-宿主的適應性進化機制提供有價值的信息和線索.為了在轉錄組水平探究在N. ceranae侵染中蜂過程中的的作用，本研究利用基于鏈特異性建庫方法的RNA-seq技術對N. ceranae侵染7 d和10 d的中蜂工蜂中腸進行測序，篩選出病原的高表達基因(highly expressed gene, HEG)，并通過生物信息方法對N. ceranae的HEG進行深入分析，并通過RT-PCR技術對部分HEG進行驗證，以期揭示HEG在N. ceranae侵染中蜂過程的作用，為深入N. ceranae的侵染機制提供有益的信息和線索.

2 材料與方法

2.1 生物材料

本研究中使用的中蜂取自福建農林大學動物科學學院(蜂學學院)教學蜂場；N. ceranae孢子[15]由福建農林大學蜂學學院蜜蜂保護課題組純化和保存.

2.2 方 法

2.2.1 蜂群選擇 挑選外觀健康、群勢較強的10群中蜂，在每個蜂箱的巢房門口隨機抓取6只工蜂，小心拉取中腸，加適量無菌水后充分研磨后取少許研磨液進行顯微制片和鏡檢，利用前人報道的N. ceranae特異性引物(F: CGGCGACGATGTGATATGAAAATATTAA, R: CCCGGTCATTCTCAAACAAAAAACCG)[17]和N. apis特異性引物(F: GGGGGCATGTCTTTGACGTACTATGTA, R: GGGGGGCGTTTAAAATGTGAAACAACTATG)[17]對鏡檢不到孢子的樣品進行PCR檢測，將鏡檢和PCR檢查都為陰性的中蜂蜂群作為本研究的實驗蜂群.

2.2.2 中蜂的飼喂感染及人工飼養 參照郭睿[17]等的方法進行中蜂的飼喂接種及人工飼養，簡述如下：將實驗蜂群內封蓋子較多的巢脾迅速提至實驗室，放入34±0.5 ℃恒溫培養箱，將剛出房的工蜂(當天記為0 d)放進干凈的塑料盒(塑料盒四周打孔通風，30只/盒)，每個盒子上方中心插入一支裝有50%(w/v)無菌糖水的飼喂器；(3)34±0.5 ℃培養24 h，將上述工蜂饑餓處理2 h，然后給每只工蜂飼喂5 μL糖水(含1×106個N. ceranae孢子)，舍棄未食盡糖水的工蜂，將食盡糖水的工蜂用于后續實驗；(4)每日檢查工蜂的存活情況，及時清理死亡的工蜂，每隔24 h更換糖水.本實驗進行三次生物學重復.

2.2.3 測序樣品準備及二代測序 因蜜蜂中腸是N. ceranae寄生和二者互作的主要部位，故本研究選取N. ceranae侵染的中蜂工蜂中腸作為測序材料.分別拉取N. ceranae感染7 d(即8日齡)和10 d(即11日齡)的工蜂中腸，迅速移至RNA free的EP管中，三只中腸并入一管，液氮速凍后轉移保存于-80 ℃超低溫冰箱中備用.將侵染中蜂8日齡和11日齡工蜂的N. ceranae分別設為Nc1和Nc2.Nc1組的3個生物學重復為：Nc1-1、Nc1-2和Nc1-3；Nc2組的3個生物學重復為：Nc2-1、Nc2-2、Nc2-3.利用RNAiso Reagent試劑盒(TaKaRa公司，中國)抽提上述樣品的總RNA，再用RNase-free DNase I去除上述樣品中基因組DNA殘留，然后用瓊脂糖凝膠電泳進行片段大小選擇.參照郭睿等[18]的方法構建cDNA文庫.建好的文庫委托廣州基迪奧生物技術有限公司進行雙端測序，測序平臺為Illumina HiSeqTM4000.原始數據已上傳NCBI的SRA數據庫，BioProject號：PRJNA562784.

2.2.4 數據質控及相關生物信息學分析 對于獲得的原始數據讀段(raw reads)，利用Perl腳本去除有adaptors、未知核苷酸比例大于5%和低質量的reads，獲得有效讀段(clean reads).按照下述步驟獲得純凈的N. ceranae數據：(1)用短reads比對工具bowtie將clean reads比對到核糖體數據庫(最多允許錯配5個)；(2)利用TopHat2軟件將未比對上的數據比對A. cerana的參考基因組(assembly ACSNU-2.0)以去除宿主本身的clean reads；(3)進而利用TopHat2軟件將未比對上剩余的clean reads繼續比對到N. ceranae的參考基因組(assembly ASM98816v1)，比對上的clean reads即為病原的數據.進一步再利用Cufflinks軟件對比對上的clean reads進行組裝.

各組基因的表達量用公式FPKM =Total exon fragment/[Mapped fragment (millions) × exon length (KB)]計算.利用R包計算個樣品之間的相關性系數.利用OmicShare在線分析工具(http://www.omicshare.com/tools/index.php/Home/Index/index.html)對各個樣品的HEG(FPKM >15)進行Venn分析，以及對共有HEG進行GO分類和KEGG代謝通路(pathway)富集分析等生物信息學相關分析.

2.2.5 HEG的RT-PCR驗證 為證明測序數據的可靠性，從共有HEG中隨機選擇8個進行RT-PCR驗證.根據所選HEG的序列信息，使用Primer Premier 5軟件設計特異性引物，委托上海生工生物工程有限公司進行引物合成，相關引物序列信息詳見表1.利用AxyPrepTMMultisource Total RNA Miniprep試劑盒(AXYGEN公司，中國)分別抽提上述受N. ceranae感染的中蜂工蜂中腸樣品的總RNA，作為模板進行反轉錄，得到的cDNA作為模板進行PCR擴增.PCR反應體系為20 μL，包括PCR Mixture 10 μL，正、反向引物各1 μL，cDNA模板1 μL，無菌水7 μL.PCR反應條件為：94 ℃預變性5 min；94 ℃50 s，58 ℃退火30 s，72 ℃延伸1 min，共35個循環；72 ℃延伸10 min.PCR擴增產物經1.5%瓊脂糖凝膠電泳檢測.

表1 RT-PCR的引物信息

3 結果與分析

3.1 測序數據概覽

比對上A. cerana基因組的clean reads數分別為180 237 114和36 054 033條，Q20和Q30分別在98.41%和95.28%及以上(表2).此外，各組內各生物學重復之間的Pearson相關性均在0.969 1及以上(圖1).以上結果說明本研究中的測序樣品重復性較好，高通量測序數據質量好，可滿足本研究中N. ceranae的HEG分析.

表2 測序數據的信息統計

Tab.2 Summary of sequencing data information

樣本Samples比對上N. ceranae基因組的有效讀段Clean reads mapped to the N. ceranae genome99%堿基正確率/%Q20/%99.9%堿基正確率/%Q30/%Nc1-138 984 05898.7596.30Nc1-235 848 94498.7496.23Nc1-3105 404 11298.6095.85Nc2-17 568 77498.4195.28Nc2-213 016 30398.8096.42Nc2-315 468 95698.6696.02

3.2 HEG的篩選及Venn分析

從Nc1和Nc2中分別篩選出1 482和901個HEG.共有HEG為890(59.6%)個，特有HEG分別為592(39.7%)和11(0.7%)個(圖2).Nc1和Nc2的部分特有HEG信息詳見表3和表4.

圖1 Nc1和Nc2中不同生物學重復間的相關性Fig.1 Pearson correlations between every two biological repeats within Nc1 and Nc2 groups

圖2 Nc1和Nc2中HEG的Venn分析

Fig.2 Venn analysis of HEGs in Nc1 and Nc2 groups

3.3 共有及特有HEG的GO分類

GO分類結果顯示分別有298、302和175個共有HEG富集在生物學進程、細胞組分和分子功能，共涉及24個功能條目，富集前5位分別是代謝進程(230個HEG)、催化活性(213個HEG)、結合(200個HEG)、細胞進程(198個HEG)、細胞(139個HEG)(圖3).

Nc1的特有HEG涉及21個功能條目，富集前10位分別是代謝進程(67個HEG)、細胞進程(64個HEG)、催化活性(63個HEG)、結合(61個HEG)、單組織進程(32個HEG).Nc2的特有HEG涉及4個功能條目，分別為代謝進程(1個HEG)、細胞進程(1個HEG)、單組織進程(1個HEG)和催化活性(1個HEG).

表3 Nc1的前11位特有HEGs

表4 Nc2的特有HEGs

圖3 Nc1和Nc2的共有HEGs的GO分類Fig.3 GO categorization of shared HEGs in Nc1and Nc2 groups

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

KEGG代謝通路富集分析結果顯示，共有HEG富集在新陳代謝、遺傳信息加工、環境信息加工、細胞組分、有機體系統等5大類的72條代謝通路，其中富集前10位分別是核糖體(53個HEG)、內質網蛋白加工(33個HEG)、蛋白酶體(27個HEG)、真核生物的核糖體生物合成(24個HEG)、嘧啶代謝(23個HEG)、RNA轉運(21個HEG)、嘌呤代謝(21個HEG)(圖4).進一步分析發現2個共有HEG富集在MAPK信號通路(圖5).

Nc1的特有HEG富集在51條代謝通路，富集前10位分別是泛素介導的蛋白水解(12個HEG)、RNA轉運(11個HEG)、細胞周期-酵母(11個HEG)、嘌呤代謝(11個HEG)、真核生物的核糖體生物合成(11個HEG).Nc2的特有HEG僅富集在1條代謝通路，即淀粉和蔗糖代謝.

圖4 Nc1和Nc2的共有HEG的KEGG代謝通路分析Fig.4 KEGG pathway enrichment analysis of shared HEGs in Nc1 and Nc2 groups

圖5 N. ceranae的MAPK信號通路概貌紅框表示Nc1和Nc2的共有HEF.Fig.5 Overview map of MAPK signaling pathway in N. ceranaeRed boxes indicate the shared HEF between Nc1 and Nc2.

圖6 8個共有HEG的RT-PCR驗證 A:磷酸甘露糖酶編碼基因; B: 假設蛋白NCER_102326編碼基因; C: RNA聚合酶II CTD磷酸酶編碼基因; D: 水通道樣蛋白編碼基因; E: 26s蛋白酶體調節亞基4編碼基因; F: 蛋白酶體亞基-1型編碼基因; G: 蓖麻毒素b凝集素編碼基因; H: 鳥苷酸結合蛋白sar1編碼基因Fig. 6 RT-PCR verification of eight shared HEGsA:phosphomannomutase encoded gene; B: hypothetical protein NCER_102326 encoded gene; C: RNA polymerase ii ctd phosphatase encoded gene; D: aquaporin-like protein encoded gene; E: 26s proteasome regulatory subunit 4 encoded gene; F: proteasome subunit beta type-1 encoded gene; G: ricin b lectin encoded gene; H: gtp-binding protein sar1 encoded gene

3.5 HEG的RT-PCR驗證

電泳結果顯示8個隨機挑選的共有HEG在Nc1和Nc2中均擴增出目的條帶(圖6)，初步證實了本研究預測的HEG真實存在.

4 討論與結論

N.ceranae廣泛感染世界各地的蜂群，除能對蜜蜂造成能量脅迫和免疫抑制等影響外，還能與其他生物或非生物脅迫因子共同危害蜜蜂健康[6-8].此前的組學相關研究主要集中在N. ceranae與A. mellifera的互作方面，有關N. ceranae與A. cerana互作的組學研究鮮有報道.為探究N. ceranae在侵染中蜂過程的HEG表達譜及功能，本研究利用RNA-seq技術對N. ceranae侵染7 d和10 d的中蜂工蜂中腸進行測序，得到高質量的測序數據，并通過連續比對A. cerana和N. ceranae的基因組獲得病原本身的轉錄組數據.從Nc1和Nc2中分別篩選出1 482和901個HEG，其中有890個在二者中皆高量表達(FPKM分別介于72 073.050-15.76和75 951.357-15.737)，分別有592個HEG(FPKM介于2 962.700-15.003)和11個HEG(FPKM介于51.873-16.387)特異性高量表達.推測上述共有HEG在N. ceranae的侵染過程發揮基礎性作用，而特有HEG在侵染的不同階段具有特定功能，但具體的功能有待進一步研究.Cornman等曾利用羅氏454平臺對N. ceranae孢子進行測序，通過相關生物信息學軟件對基因組進行了組裝和注釋[19].然而，由于缺失基因功能研究技術體系，包括N. ceranae在內的微孢子蟲的轉基因操作技術尚未完全建立，僅在家蠶微孢子蟲中有過一例轉基因研究報道[20]，因而，包括N. ceranae在內的微孢子蟲的基因功能研究近乎停滯，加上N. ceranae的組學數據相對有限，導致其基因功能注釋信息較不完善.本研究發現，Nc1和Nc2的HEG共包含362個假定蛋白編碼基因，這些HEG在各大主流蛋白功能數據庫尚缺乏功能注釋信息，其完善需要更多的組學數據的支持以及轉基因操作技術體系的建立.

本研究中，Nc1和Nc2的共有HEG涉及24個功能條目；分別有139、139、64、34、30和17個HEG富集在細胞、細胞組件、細胞器、細胞膜、細胞膜組件和細胞器組件，另有198和27個HEG富集在細胞進程和細胞組分組織或生物合成，表明N. ceranae在侵染中蜂的過程中伴隨著活躍的細胞生命活動，以促進自身在宿主細胞內的增殖.此外，上述共有HEG參與了72條代謝通路，包含碳代謝(15)和糖酵解/糖異生(11)等41條物質代謝通路，氧化磷酸化(8)和磷酸戊糖途徑(5)2條能量代謝通路.上述結果表明N. ceranae在侵染過程中的物質和能量代謝較為旺盛.Ponton等[21]的發現N. ceranae侵染A. mellifera后即利用宿主的物質和能量進行自身遺傳物質的復制、轉錄和翻譯.本研究發現分別有23、21和21個N. ceranae的HEG富集在嘧啶代謝、RNA轉運和嘌呤代謝；此外，分別有53、33和19個HEG富集在核糖體、內質網蛋白質加工和氨?；?tRNA生物合成.這些結果表明N. ceranae在中蜂工蜂中腸上皮細胞內的增殖過程中同時進行著較為活躍的DNA復制、RNA轉錄和蛋白質翻譯等活動.

信號通路指能將細胞外的信號經細胞膜傳入細胞內的一系列酶促反應通路，真菌通過MAPK等信號通路轉導系統對外界環境的變化作出相應的應答，進而對生長和毒力等相關基因的表達進行調控，從而影響孢子形成、營養感知和形態建成等生物學過程[22].陳大福等[23]發現對于侵染中蜂6日齡幼蟲的蜜蜂球囊菌(Ascopshaera apis，簡稱球囊菌)，富集在MAPK信號通路的HEG遠多于體外培養的球囊菌，推測MAPK信號通路在球囊菌的侵染過程作用關鍵.本研究中，2個共有HEG富集在MAPK信號通路，暗示此信號通路在N. ceranae的侵染過程中發揮重要作用.Jaroenlak等[24]在蝦肝腸微孢子蟲(Enterocytozoon hepatopenaei)中發現一種分布在內生孢子和外孢子層的新孢壁蛋白(EhSWP1)，作為病原在侵染早期的重要毒力因子，可能通過結合肝素附著到宿主細胞表面.極管蛋白與N. ceranae毒力大小有關，有研究表明干擾N. ceranae的極性管蛋白3基因(ptp3)可抑制N. ceranae增殖[25].Paldi等發現宿主攝入靶向N. ceranae的ADP/ATP轉移蛋白基因的siRNA能導致特定的基因沉默，進而會影響N. ceranae感染水平和宿主的生理狀況[26].本研究中，孢壁蛋白基因(XM_002996303.1, Nc1中FPKM=24 434.550, Nc2中FPKM=33 275.467)、極管蛋白基因(XM_002995446.1, Nc1中FPKM=17 490.823, Nc2中FPKM=29 944.273)、蓖麻毒素B鏈基因(XM_002995387.1, Nc1中FPKM=698.757, Nc2中FPKM=1 547.050)、ABC轉運體基因(XM_002995069.1，Nc1中FPKM=133.697, Nc2中FPKM=265.673)、ATP/ADP轉移酶基因(XM_002995682.1，Nc1中FPKM=10 832.593, Nc2中FPKM=7 752.717)和絲氨酸蘇氨酸蛋白激酶基因(XM_002995071.1, Nc1中FPKM=337.033, Nc2中FPKM=266.083)等毒力基因均為高量表達，表明N. ceranae在其侵染中蜂工蜂的過程中通過維持上述毒力基因的高表達水平以提升毒力蛋白的合成代謝，進而促進自身在宿主細胞內的增殖.海藻糖廣泛存在于酵母等真菌的孢子和子實體中，在脫水、干旱、高溫、冷凍、高滲透壓和有毒試劑等不良環境的脅迫下起到保護作用[27].海藻糖與海藻糖酶結合會產生大量葡萄糖，導致細胞內滲透壓升高促進微孢子蟲對宿主的侵染，在微孢子蟲孢子的發芽過程中海藻糖酶催化海藻糖水解為葡萄糖為孢子正常發育提供能量[28].本研究發現，海藻糖磷酸合酶1基因(XM_002996572.1, Nc1中FPKM=239.363, Nc2中FPKM=301.7)和海藻糖-6-磷酸合酶基因(XM_002996577.1, Nc1中FPKM=73.413, Nc2中FPKM=200.397)在Nc1和Nc2中皆高量表達，推測N. ceranae通過調控海藻糖磷酸合酶合成速率為自身增殖提供能量.筆者團隊在前期研究中發現對于N. ceranae侵染7d和10d的中蜂工蜂中腸，隨著侵染時間的延長，與宿主物質和能量代謝相關的HEG數量持續降低[29]，表明宿主的物質和能量代謝受到N. ceranae的抑制，體現出N. ceranae通過與中蜂互作而對后者的部分生命活動進行操縱.前人研究發現通過RNA干擾(RNAi)對N. ceranae的ATP/ADP轉移酶和裸角質層(Naked cuticle)等毒力因子基因進行沉默的效果良好[26, 30].未來可針對N. ceranae的海藻糖磷酸合酶1基因和海藻糖-6-磷酸合酶基因設計特異性的小干擾RNA(siRNA)，通過注射或飼喂感染蜜蜂對該上述病原基因進行沉默，有望用于蜜蜂微孢子蟲病的治療.值得一提的是，本研究發現若干表達水平相對特別高的基因，如部分SR22蛋白基因(XM_002996758.1, Nc1中FPKM=61 410.717, Nc2中FPKM=75 951.357)、假定蛋白NCER_100570基因(XM_002996307.1, Nc1中FPKM=69 904.740, Nc2中FPKM=75 053.357)和延伸因子1基因(XM_002995284.1, Nc1中FPKM=33 358.027, Nc2中FPKM=34 716.303)，下一步擬通過RT-PCR和RACE技術得到上述HEG的cDNA全長序列，進而在核酸和蛋白水平開展相關的功能研究.

綜上所述，N. ceranae可能通過MAPK信號通路對中蜂中腸的堿性環境作出應答，促進孢壁蛋白、極管蛋白、蓖麻毒素B鏈、ABC轉運體、ADP/ATP轉移蛋白和絲氨酸蘇氨酸蛋白激酶等毒力因子基因的高量表達，刺激孢子萌發和菌絲彈射，促進N. ceranae對中蜂中腸上皮細胞的侵染；在中蜂中腸細胞內，激活ADP/ATP轉移蛋白和海藻糖磷酸合酶基因的表達，對宿主造成能量脅迫，為自身增殖提供能量.本研究通過深入分析侵染中蜂工蜂的N. ceranae的HEG表達譜，揭示了HEG在病原侵染過程中的潛在作用，揭示N. ceranae通過維持孢壁蛋白、極管蛋白、蓖麻毒素B鏈、ABC轉運體、ADP/ATP轉移蛋白和絲氨酸蘇氨酸蛋白激酶等毒力因子基因，以及ATP/ADP轉移酶、海藻糖磷酸合酶1和海藻糖-6-磷酸合酶等能量代謝相關蛋白編碼基因的高量表達，促進自身在中蜂工蜂中腸的侵染過程.研究結果為明確N. ceranae對中蜂的侵染機制提供了有益信息和線索.