轉(zhuǎn)錄組分析揭示東方蜜蜂微孢子蟲侵染意大利蜜蜂的分子機(jī)制

耿四海, 周丁丁, 范小雪, 蔣海賓, 祝智威, 王 杰, 范元嬋,王心蕊, 熊翠玲, 鄭燕珍, 付中民, 陳大福, 郭 睿

(福建農(nóng)林大學(xué)動(dòng)物科學(xué)學(xué)院(蜂學(xué)學(xué)院), 福州 350002)

微孢子蟲是一種細(xì)胞內(nèi)寄生的單細(xì)胞真菌，廣泛寄生脊椎動(dòng)物和無脊椎動(dòng)物(Cuomoetal., 2012)。自1857年首次在患病家蠶Bombyxmori中發(fā)現(xiàn)家蠶微孢子蟲Nosemabombycis(Nageli, 1857)以來，至今已經(jīng)有1 400多種微孢子蟲被鑒定出來(吳杰等, 2017)。微孢子蟲的基因組緊湊，缺少典型的線粒體，取而代之是一種微粒體，因而高度依賴宿主細(xì)胞提供能量(Burrietal., 2006)。

東方蜜蜂微孢子蟲Nosemaceranae是專性寄生蜜蜂中腸上皮細(xì)胞的單細(xì)胞真菌病原，可引起蜜蜂微孢子蟲病(Forsgren and Fries, 2010)。Freis等(1996)首次在北京地區(qū)的東方蜜蜂Apiscerana蜂群中分離和鑒定出東方蜜蜂微孢子蟲，之后該孢子蟲逐漸擴(kuò)散到歐洲(Higesetal., 2006)和中國臺(tái)灣地區(qū)(Huangetal., 2007)的西方蜜蜂Apismellifera蜂群，現(xiàn)已蔓延至世界各養(yǎng)蜂國家和地區(qū)(Shumkovaetal., 2018)。前人的研究結(jié)果表明東方蜜蜂微孢子蟲感染可導(dǎo)致成年蜜蜂腸道破壞、消化系統(tǒng)紊亂,脂質(zhì)代謝加快和脂質(zhì)儲(chǔ)存減少,免疫基因和細(xì)胞凋亡被抑制,能量脅迫,采集日齡提前及采集能力下降等；東方蜜蜂微孢子蟲感染嚴(yán)重可造成蜂群的群勢降低、采蜜量減少以及越冬失敗(Higesetal., 2007; Antúnezetal., 2009; Mayack and Naug, 2009; Botíasetal., 2013; Chenetal., 2013; Mikeetal., 2013; Kurzeetal., 2015)。此外，有研究表明東方蜜蜂微孢子蟲與蜂群崩潰失調(diào)癥之間存在密切的關(guān)聯(lián)(Higesetal., 2009; Simeunovicetal., 2014)。東方蜜蜂微孢子蟲在體外以休眠態(tài)的孢子存在，孢子壁由幾丁質(zhì)和孢壁蛋白組成，能夠抵御外界的不良環(huán)境，孢子被蜜蜂攝取之后，在蜜蜂中腸的特殊理化環(huán)境的刺激下發(fā)芽，內(nèi)部高度壓縮盤旋的極絲彈射而出，刺穿中腸上皮細(xì)胞并將有侵染性的胞原質(zhì)注入，隨后進(jìn)入繁殖期，利用宿主的物質(zhì)與能量合成系統(tǒng)進(jìn)行增殖，孢子數(shù)量逐漸增多，最終導(dǎo)致宿主細(xì)胞裂解，釋放出的成熟孢子繼續(xù)感染臨近的上皮細(xì)胞，或者隨糞便排出體外，感染新的宿主(Simeunovicetal., 2014)。

隨著高通量測序技術(shù)和生物信息學(xué)算法的發(fā)展，RNA-Seq技術(shù)不僅廣泛應(yīng)用于人類(Violletetal., 2017)、動(dòng)物(李鵬飛等, 2018)和植物(Munetal., 2017)的研究，在微生物的相關(guān)研究中也得到了較好的應(yīng)用，包括綠僵菌Metarhiziumanisopliae(Wangetal., 2014)、里氏木霉Trichodermareesei(Riesetal., 2013)、家蠶微孢子蟲(Liuetal., 2016a)和蜜蜂球囊菌Ascosphaeraapis(郭睿等, 2017)等。Huang等(2016)對感染東方蜜蜂微孢子蟲后1-6 d的西方蜜蜂中腸進(jìn)行測序，從混合數(shù)據(jù)中篩選出東方蜜蜂微孢子蟲的數(shù)據(jù)，進(jìn)而通過統(tǒng)計(jì)東方蜜蜂微孢子蟲的基因表達(dá)量，發(fā)現(xiàn)1 122個(gè)基因可分為4種表達(dá)模式；隨后，該團(tuán)隊(duì)對東方蜜蜂微孢子蟲的Dicer進(jìn)行RNAi，然后通過轉(zhuǎn)錄組測序和分析發(fā)現(xiàn)Dicer調(diào)控東方蜜蜂微孢子蟲的基因表達(dá)(Huangetal., 2018)。筆者團(tuán)隊(duì)前期結(jié)合鏈特異性建庫方法和二代測序技術(shù)對東方蜜蜂微孢子蟲孢子進(jìn)行了深度測序，在全基因組水平預(yù)測、鑒定和分析了東方蜜蜂微孢子蟲的長鏈非編碼RNA(long non-coding RNA, lncRNA)和環(huán)狀RNA(circular RNA, circRNA)(Guoetal., 2018a, 2018b)。但對于東方蜜蜂微孢子蟲的組學(xué)研究較之蜜蜂仍然較為滯后，東方蜜蜂微孢子蟲的侵染和增殖機(jī)制還未闡明。

筆者團(tuán)隊(duì)前期已對東方蜜蜂微孢子蟲感染7和10 d的意大利蜜蜂Apismelliferaligustica工蜂中腸進(jìn)行轉(zhuǎn)錄組深度測序，通過連續(xù)比對西方蜜蜂和東方蜜蜂微孢子蟲的基因組過濾得到東方蜜蜂微孢子蟲的純凈數(shù)據(jù)，結(jié)合東方蜜蜂微孢子蟲孢子的轉(zhuǎn)錄組數(shù)據(jù)(Guoetal., 2018a)開展了東方蜜蜂微孢子蟲的高表達(dá)基因(highly expressed gene, HEG)分析，解析了東方蜜蜂微孢子蟲侵染意大利蜜蜂過程的HEG表達(dá)譜及潛在作用(熊翠玲等, 2019)。本研究在前期研究基礎(chǔ)上對東方蜜蜂微孢子蟲進(jìn)行深入的轉(zhuǎn)錄組分析，進(jìn)一步通過分析病原的毒力因子和侵染相關(guān)因子，揭示東方蜜蜂微孢子蟲侵染意大利蜜蜂的分子機(jī)制。研究結(jié)果可為探明東方蜜蜂微孢子蟲的侵染機(jī)制提供基礎(chǔ)，為蜜蜂微孢子蟲病的治療提供潛在靶點(diǎn)。

1 材料與方法

1.1 轉(zhuǎn)錄組數(shù)據(jù)來源

前期研究中，筆者團(tuán)隊(duì)利用基于鏈特異性建庫的lncRNA-Seq技術(shù)對東方蜜蜂微孢子蟲的純化孢子進(jìn)行了深度測序(NcCK-1, NcCK-2和NcCK-3為3個(gè)生物學(xué)重復(fù))，獲得了高質(zhì)量的mRNA和lncRNA組學(xué)數(shù)據(jù)(Guoetal., 2018a)，其中的mRNA組學(xué)數(shù)據(jù)可作為本研究中的對照組數(shù)據(jù)；此外，還利用二代測序技術(shù)對東方蜜蜂微孢子蟲感染7 d(NcT1-1, NcT1-2和NcT1-3為3個(gè)生物學(xué)重復(fù))和10 d(NcT2-1, NcT2-2和NcT2-3為3個(gè)生物學(xué)重復(fù))的意大利蜜蜂工蜂中腸進(jìn)行了全轉(zhuǎn)錄組測序，獲得了高質(zhì)量的測序數(shù)據(jù)(熊翠玲等, 2019)，其中的mRNA組學(xué)數(shù)據(jù)包含宿主來源的數(shù)據(jù)和病原來源的數(shù)據(jù)，因此通過連續(xù)比對核糖體數(shù)據(jù)庫、西方蜜蜂基因組和東方蜜蜂微孢子蟲基因組對上述混合數(shù)據(jù)進(jìn)行過濾，篩濾得到的病原的mRNA組學(xué)數(shù)據(jù)(即侵染意大利蜜蜂工蜂的東方蜜蜂微孢子蟲的mRNA組學(xué)數(shù)據(jù))作為本研究中的處理組數(shù)據(jù)。

上述轉(zhuǎn)錄組測序的原始數(shù)據(jù)已上傳美國國家生物技術(shù)信息中心(NCBI)SRA數(shù)據(jù)庫，Bioproject號分別為PRJNA395264(NcCK)和PRJNA406998(NcT1和NcT2)。

1.2 差異表達(dá)基因(DEG)分析

采用FPKM算法計(jì)算和歸一化基因表達(dá)量。利用edgeR軟件篩選NcCKvsNcT1, NcCKvsNcT2和NcT1vsNcT2比較組的DEG，標(biāo)準(zhǔn)為P≤0.05且|log2(Fold change)|≥1。利用OmicShare在線平臺(tái)(www.omicshare.com)對各比較組中上調(diào)基因和下調(diào)基因進(jìn)行Venn分析，采用默認(rèn)參數(shù)。參照郭睿等(2017)的方法，利用BLASTall軟件將DEG比對到GO和KEGG數(shù)據(jù)庫，進(jìn)行功能分類和代謝通路注釋，采用默認(rèn)參數(shù)。

1.3 毒力因子和侵染相關(guān)因子分析

東方蜜蜂微孢子蟲是專性寄生成年蜜蜂中腸上皮細(xì)胞的真菌病原，其增殖所需的物質(zhì)和能量完全依賴宿主細(xì)胞提供(Burrietal., 2006)。根據(jù)DEG在Nr數(shù)據(jù)庫的注釋信息，同時(shí)參考前人關(guān)于東方蜜蜂微孢子蟲的基因組學(xué)、轉(zhuǎn)錄組學(xué)和分子生物學(xué)的研究報(bào)道(Cornmanetal., 2009; Paldietal., 2010; Huangetal., 2016, 2018; Pelinetal., 2016; Rodríguez-Garcíaetal., 2018)以及其他微孢子蟲的相關(guān)研究報(bào)道(Corradi and Slamovits, 2011; Katinkaetal., 2011; Nakjangetal., 2013; Simeunovicetal., 2014)，篩選出幾丁質(zhì)酶、孢壁蛋白、蓖麻毒素B凝集素、極管蛋白等毒力因子相關(guān)DEG，以及己糖激酶、丙酮酸激酶、磷酸果糖激酶、ATP/ADP移位酶、ABC轉(zhuǎn)運(yùn)蛋白等侵染相關(guān)因子相關(guān)DEG，并統(tǒng)計(jì)上述相關(guān)DEG在各組樣品中的表達(dá)量(FPKM)及在各比較組中的log2(Fold change)。

1.4 DEG的RT-qPCR驗(yàn)證

從NcCKvsNcT1, NcCKvsNcT2和NcT1vsNcT2比較組中分別隨機(jī)選取7個(gè)DEG，其中NcCKvsNcT1中包含6個(gè)上調(diào)基因和1個(gè)下調(diào)基因，NcCKvsNcT2中包含6個(gè)上調(diào)基因和1個(gè)下調(diào)基因，NcT1vsNcT2中包含5個(gè)上調(diào)基因和2個(gè)下調(diào)基因，以actin(GenBank登錄號: 36320842)作為內(nèi)參基因，利用Primer Premier 5設(shè)計(jì)引物(表1)，委托福州擎科生物有限公司合成。利用RNA提取試劑盒(TaKaRa公司，大連)分別提取NcCK以及NcT1和NcT2中腸總RNA(熊翠玲等, 2019)，然后利用cDNA第1鏈合成試劑盒(Vazyme公司，南京)反轉(zhuǎn)錄成cDNA，作為模板進(jìn)行qPCR，反應(yīng)在QuantStudio 3(ThermoFisher公司，美國)上進(jìn)行。qPCR反應(yīng)按照SYBR Green Dye試劑盒(Vazyme公司，南京)操作說明書進(jìn)行，本實(shí)驗(yàn)進(jìn)行3次生物學(xué)重復(fù)，每個(gè)生物學(xué)重復(fù)包含6頭意大利蜜蜂工蜂中腸。反應(yīng)體系(20 μL)：正反向引物(2.5 μmol/L)各1 μL, cDNA模板1 μL, SYBR Green Dye 10 μL, DEPC水7 μL。反應(yīng)程序: 95℃預(yù)變性1 min; 95℃變性15 s, 60℃延伸30 s，共40個(gè)循環(huán)。通過2-△△Ct法對上述基因的相對表達(dá)量進(jìn)行計(jì)算。利用Graph Prism 7軟件進(jìn)行數(shù)據(jù)分析和繪圖。

2 結(jié)果

2.1 東方蜜蜂微孢子蟲轉(zhuǎn)錄組的DEG分析

差異分析結(jié)果顯示，NcCKvsNcT1, NcCKvsNcT2和NcT1vsNcT2比較組分別得到1 397, 1 497和52個(gè)DEG，其中上調(diào)和下調(diào)基因分別為497和900個(gè), 517和977個(gè)，38和13個(gè)。Venn分析結(jié)果顯示，有10個(gè)基因在上述3個(gè)比較組中共同上調(diào)表達(dá)(圖1: A)，它們涉及碳水化合物代謝、蓖麻毒素B凝集素、細(xì)胞表面糖基化蛋白、信號肽酶和假定蛋白(表2)；僅有1個(gè)下調(diào)表達(dá)的假定蛋白編碼基因?yàn)楦鞅容^組所共有(圖1: B; 表2)。

表1 本研究使用的引物Table 1 Primers used in this study

圖1 各比較組中東方蜜蜂微孢子蟲轉(zhuǎn)錄組差異表達(dá)基因(DEGs)Venn分析Fig. 1 Venn analysis of differentially expressed genes (DEGs) in Nosema ceranae transcriptome between various comparison groupsA: 上調(diào)基因Up-regulated genes; B: 下調(diào)基因Down-regulated genes. NcCK: 微孢子蟲純化的孢子Purified spores of N.ceranae; NcT1, NcT2: 分別為感染意大利蜜蜂工蜂7和10 d中腸內(nèi)的東方蜜蜂微孢子蟲N. ceranae in the midgut of Apis mellifera ligustica workers at 7 and 10 d post infection, respectively. 下同The same below.

表2 各比較組中東方蜜蜂微孢子蟲轉(zhuǎn)錄組共同上調(diào)與下調(diào)基因信息統(tǒng)計(jì)Table 2 Summary of commonly up- and down-regulated genes in Nosema ceranae transcriptome between various comparison groups

2.2 東方蜜蜂微孢子蟲轉(zhuǎn)錄組中DEG的GO分類及KEGG 通路富集分析

GO分類結(jié)果表明，NcCKvsNcT1中DEG涉及生物學(xué)進(jìn)程、細(xì)胞組分和分子功能相關(guān)的23個(gè)功能條目，其中代謝進(jìn)程(221)、催化活性(195)、結(jié)合(193)、細(xì)胞進(jìn)程(193)、細(xì)胞(116)、細(xì)胞組件(116)、單組織進(jìn)程(97)和細(xì)胞器(63)的條目富集基因數(shù)最多；NcCKvsNcT2中DEG也涉及23個(gè)功能條目，同樣有較多的DEG富集在代謝進(jìn)程(218)、催化活性(196)、細(xì)胞進(jìn)程(195)、結(jié)合(193)、細(xì)胞(119)、細(xì)胞組件(119)、單組織進(jìn)程(94)和細(xì)胞器(51)等；NcT1vsNcT2中DEG涉及生物學(xué)進(jìn)程(4)、細(xì)胞組分(5)和分子功能(2)相關(guān)的11個(gè)功能條目，其中催化活性、細(xì)胞進(jìn)程、代謝進(jìn)程、單組織進(jìn)程和結(jié)合的條目富集基因數(shù)最多，分別為11, 10, 10, 7和6個(gè)(圖2)。

進(jìn)一步對各比較組的DEG進(jìn)行KEGG通路富集分析，結(jié)果顯示NcCKvsNcT1中有355個(gè)DEG富集在80條通路，涉及代謝、遺傳信息加工、環(huán)境信息加工、細(xì)胞進(jìn)程和有機(jī)體系統(tǒng)等5個(gè)大類，其中富集基因數(shù)最多的是新陳代謝途徑(93)、次級代謝產(chǎn)物合成(37)、核糖體(33)、抗生素生物合成(29)和真核生物核糖體生物合成(27)(表3)；NcCKvsNcT2中有356個(gè)DEG富集在79條通路，其中富集基因數(shù)最多的是新陳代謝途徑、次級代謝產(chǎn)物合成、核糖體、抗生素合成和真核生物核糖體生物合成，分別為92, 39, 32, 31和26個(gè)(表3)。此外NcCKvsNcT1中有7個(gè)DEG富集在MAPK信號通路，其中有5個(gè)基因上調(diào)表達(dá)，GenBank登錄號分別為9423510, 9424462, 9422987, 9423963和9422306；另有2個(gè)基因下調(diào)表達(dá)，GenBank登錄號分別為9422598和9424007。 NcCKvsNcT2中同樣也有7個(gè)DEG富集在MAPK信號通路，包括5個(gè)上調(diào)基因，GenBank登錄號分別為9423510, 9424610, 9424462, 9422987和9422306；以及2個(gè)下調(diào)基因，GenBank登錄號分別為9411598和9424007。在NcCK, NcT1和NcT2中，GenBank登錄號為9424462的基因的FPKM值分別為966, 829和22，GenBank登錄號為9422306的基因的FPKM值分別為209, 205和5。

圖2 各比較組中東方蜜蜂微孢子蟲轉(zhuǎn)錄組差異表達(dá)基因(DEGs)的GO分類Fig. 2 GO categories of differentially expressed genes (DEGs) in Nosema ceranae transcriptome between various comparison groups

2.3 東方蜜蜂微孢子蟲的毒力因子分析

根據(jù)NcCKvsNcT1, NcCKvsNcT2和NcT1vsNcT2的DEG的Nr數(shù)據(jù)庫注釋情況對東方蜜蜂微孢子蟲的毒力因子進(jìn)行統(tǒng)計(jì)，分別發(fā)現(xiàn)5個(gè)孢壁蛋白基因、1個(gè)孢壁蛋白前體基因和1個(gè)孢壁和錨定盤復(fù)合蛋白基因，其中孢壁蛋白9基因(GenBank登錄號: 9422782)在NcCKvsNcT1和NcCKvsNcT2中下調(diào)表達(dá)[log2(Fold change)值分別為-2.10和-2.06]，孢壁蛋白12基因(GenBank登錄號: 9422437)也下調(diào)表達(dá)[log2(Fold change)值分別為-0.86和-0.67]；孢壁蛋白8基因(GenBank登錄號: 9422311)在NcCKvsNcT1中表達(dá)量下調(diào)[log2(Fold change)值為-2.05]，但在NcCKvsNcT2中其表達(dá)量不變；在NcCKvsNcT1和NcCKvsNcT2中均上調(diào)表達(dá)的基因包括2個(gè)孢壁蛋白基因[GenBank登錄號為9422653的基因的log2(Fold change)值分別為7.11和7.12; GenBank登錄號為9423533的基因的log2(Fold change)值分別為6.32和6.93]、1個(gè)孢壁蛋白前體基因[GenBank登錄號: 9422389, log2(Fold change)值分別為5.03和5.46]和1個(gè)孢壁和錨定盤復(fù)合蛋白基因[GenBank登錄號: 9424403, log2(Fold change)值分別為5.54和4.61]。此外，另有1個(gè)內(nèi)切幾丁質(zhì)酶基因[GenBank登錄號: 9422749, log2(Fold change)值分別為2.94和3.60]、1個(gè)幾丁質(zhì)合酶1基因[GenBank登錄號為9422288的基因的log2(Fold change)值分別為6.08和5.56]、3個(gè)極管蛋白基因[GenBank登錄號為9423198的基因的log2(Fold change)值分別為5.62和6.34; GenBank登錄號為9423201的基因的log2(Fold change)值分別為5.74和6.37; GenBank登錄號為9422340的基因的log2(Fold change)值分別為3.94和4.39]和7個(gè)蓖麻毒素B凝集素基因[GenBank登錄號為9422649的基因的log2(Fold change)值分別為4.76和5.98; GenBank登錄號為9422719的基因的log2(Fold change)值分別為4.89和6.32; GenBank登錄號為9422642的基因的log2(Fold change)值分別為4.77和4.98; GenBank登錄號為9423906的基因的log2(Fold change)值分別為2.32和3.23; GenBank登錄號為9422645的基因的log2(Fold change)值分別為3.43和4.74; GenBank登錄號為9422644的基因的log2(Fold change)值分別為2.72和3.21; GenBank登錄號為9422647的基因的log2(Fold change)值分別為5.62和3.38]在NcCKvsNcT1和NcCKvsNcT2中均上調(diào)表達(dá)(表4)。

表3 NcCK vs NcT1和NcCK vs NcT2中差異表達(dá)基因(DEGs)富集數(shù)前18的KEGG通路Table 3 Top 18 KEGG pathways enriched by differentially expressed genes (DEGs) in NcCK vs NcT1 and NcCK vs NcT2

2.4 東方蜜蜂微孢子蟲的侵染相關(guān)因子分析

進(jìn)一步對東方蜜蜂微孢子蟲的其他侵染相關(guān)因子進(jìn)行統(tǒng)計(jì)，發(fā)現(xiàn)在NcCKvsNcT1和NcCKvsNcT2中均上調(diào)表達(dá)的糖酵解關(guān)鍵酶基因包括己糖激酶基因[GenBank登錄號: 9423428, log2(Fold change)值分別為3.38和2.79]、丙酮酸激酶基因[GenBank登錄號: 9424331, log2(Fold change)值分別為1.81和2.69]和6-磷酸果糖激酶基因[GenBank登錄號: 9424162, log2(Fold change)值分別為1.91和2.79]等(表5)。在NcCKvsNcT1和NcCKvsNcT2中均上調(diào)表達(dá)的有3個(gè)ATP/ADP移位酶的編碼基因[GenBank登錄號為9424586的基因的log2(Fold change) 值分別為2.31和2.34; GenBank登錄號為9422880的基因的log2(Fold change)值分別為3.60和3.44; GenBank登錄號為9424363的基因的log2(Fold change)值分別為0.72和0.45]。此外，有1個(gè)ATP/ADP移位酶的編碼基因[GenBank登錄號:9422375, log2(Fold change)值分別為-3.65和-3.90]在NcCKvsNcT1和NcCKvsNcT2中均下調(diào)表達(dá)。另發(fā)現(xiàn)7個(gè)DEG涉及ABC轉(zhuǎn)運(yùn)蛋白，其中有2個(gè)[GenBank登錄號為9423041的基因的log2(Fold change)值分別為1.83和1.85; GenBank登錄號為9422518的基因的log2(Fold change)值分別為0.93和0.64]在NcCKvsNcT1和NcCKvsNcT2中均上調(diào)表達(dá)；還有1個(gè)DEG(GenBank登錄號: 9423588)在NcCKvsNcT1中下調(diào)表達(dá)[log2(Fold change)值為-0.79]，而在NcCKvsNcT2中表達(dá)量上調(diào)[log2(Fold change)值為0.67]；其余的4個(gè)ABC轉(zhuǎn)運(yùn)蛋白基因[GenBank登錄號為9423050的基因的log2(Fold change)值分別為-3.65和-3.34; GenBank登錄號為9423588的基因的log2(Fold change)值分別為-1.72和-0.13; GenBank登錄號為9422298的基因的log2(Fold change)值分別為-1.91和-1.56; GenBank登錄號為9422297的基因的log2(Fold change)值分別為-0.19和-0.35]在NcCKvsNcT1和NcCKvsNcT2中表現(xiàn)為表達(dá)量下調(diào)(表5)。

表4 東方蜜蜂微孢子蟲的毒力因子統(tǒng)計(jì)Table 4 Summary of virulence factors of Nosema ceranae

表5 東方蜜蜂微孢子蟲的侵染相關(guān)因子統(tǒng)計(jì)Table 5 Summary of infection-associated factors of Nosema ceranae

2.5 東方蜜蜂微孢子蟲轉(zhuǎn)錄組DEG的RT-qPCR驗(yàn)證

為驗(yàn)證RNA-Seq結(jié)果，在NcCKvsNcT1, NcCKvsNcT2和NcT1vsNcT2比較組中分別隨機(jī)挑取7個(gè)DEG進(jìn)行RT-qPCR驗(yàn)證，實(shí)驗(yàn)結(jié)果與轉(zhuǎn)錄組數(shù)據(jù)的差異表達(dá)趨勢一致(圖3)，證明了測序數(shù)據(jù)和基因差異表達(dá)情況的可靠性。

圖3 各比較組中東方蜜蜂微孢子蟲轉(zhuǎn)錄組差異表達(dá)基因(DEGs)的RT-qPCR驗(yàn)證Fig. 3 RT-qPCR verification of differentially expressed genes (DEGs) in Nosema ceranae transcriptome between various comparison groupsA: NcCK vs NcT1; B: NcCK vs NcT2; C: NcT1 vs NcT2. 誤差棒表示各DEG的RT-qPCR結(jié)果的方差。Error bars represent the variance of RT-qPCR results of each DEG.

3 討論

前人研究報(bào)道東方蜜蜂微孢子蟲在感染西方蜜蜂后10-20 d內(nèi)都處于一個(gè)持續(xù)增殖的過程，在此過程病原的孢子數(shù)持續(xù)增多(Huang and Solter, 2013)。前期研究中，我們發(fā)現(xiàn)東方蜜蜂微孢子蟲感染的意大利蜜蜂工蜂的死亡率在感染后1-10 d范圍內(nèi)不斷上升，且與未感染工蜂的死亡率差異極顯著(Chenetal., 2019)，表明隨著宿主中腸內(nèi)東方蜜蜂微孢子蟲的不斷增殖，病原數(shù)量持續(xù)增多，對宿主的脅迫壓力持續(xù)上升。在另一項(xiàng)研究中，我們發(fā)現(xiàn)東方蜜蜂微孢子蟲的孢子數(shù)在4-13 d的范圍內(nèi)持續(xù)增加(未發(fā)表數(shù)據(jù))。本研究是屬于完整課題的一部分，完整課題的主要研究目的是探究意大利蜜蜂與東方蜜蜂微孢子蟲的相互作用機(jī)制，一方面是解析宿主的脅迫應(yīng)答和免疫應(yīng)答機(jī)制，另一方面是解析病原的侵染機(jī)制。鑒于東方蜜蜂微孢子蟲的增殖周期尚無定論，以及處于增殖過程的病原存在孢原質(zhì)、裂殖體和成熟孢子等各種狀態(tài)，綜合考慮后選取東方蜜蜂微孢子蟲感染后7和10 d作為本研究取樣、測序和數(shù)據(jù)分析的兩個(gè)時(shí)間點(diǎn)。目前仍缺乏不依賴于蜜蜂宿主的東方蜜蜂微孢子蟲體外培養(yǎng)體系。本研究中使用的東方蜜蜂微孢子蟲純化孢子是從患微孢子蟲病的意大利蜜蜂外勤蜂中腸制備而來，其中可能既含有經(jīng)增殖而致死宿主的孢子，也含有未侵入宿主細(xì)胞的孢子。鑒于目前沒有區(qū)分此二種類型孢子的技術(shù)手段，我們參照前人研究中采用的接種西方蜜蜂所采用的孢子濃度和接種方法(Antúnezetal., 2009; Huang and Solter, 2013; Huangetal., 2016, 2018; Rodríguez-Garcíaetal., 2018)，使每頭工蜂攝入等量的孢子。對于微孢子蟲，相對于代謝和生命活動(dòng)活躍的增殖階段(如裂殖體)，休眠態(tài)孢子中僅存在較低水平的代謝和生命活動(dòng)(Williamsetal., 2005; Corradietal., 2008)。本研究中，用于接種工蜂的孢子(其中已完成侵染的孢子仍然為休眠態(tài)孢子)來源于同一批制備的東方蜜蜂微孢子蟲純化孢子，通過對照組和處理組數(shù)據(jù)比較分析得到的基因差異表達(dá)信息應(yīng)由病原侵染所致。本研究中，在NcCKvsNcT1, NcCKvsNcT2和NcT1vsNcT2中分別鑒定出1 397, 1 497和52個(gè)顯著DEG，表明東方蜜蜂微孢子蟲在其侵染意大利蜜蜂工蜂的過程中伴隨著活躍的基因差異表達(dá)，暗示這些DEG中蘊(yùn)含著病原侵染相關(guān)的重要信息。

Cornman等(2009)組裝及注釋了東方蜜蜂微孢子蟲的基因組，分析發(fā)現(xiàn)東方蜜蜂微孢子蟲的供能主要依賴碳水化合物代謝、糖酵解、氧化磷酸化和磷酸戊糖途徑。本研究發(fā)現(xiàn)，有10個(gè)基因在NcCKvsNcT1, NcCKvsNcT2和NcT1vsNcT2 3個(gè)比較組均上調(diào)表達(dá)，其中4個(gè)上調(diào)基因涉及碳水化合物代謝，包括己糖激酶基因(GenBank登錄號: 9423428)、甘油-3-磷酸脫氫酶基因(GenBank登錄號: 9423035)、磷酸甘油酸激酶基因(GenBank登錄號: 9423279)和丙酮酸脫氫酶e1-β亞基基因(GenBank登錄號: 9423243)(表2)，表明東方蜜蜂微孢子蟲在侵染過程中通過上調(diào)上述4個(gè)基因的表達(dá)量，滿足自身增殖的能量需要。蓖麻毒素是一種可抑制蛋白質(zhì)活性和改變宿主細(xì)胞激素編碼基因的表達(dá)水平的細(xì)胞毒素，具有結(jié)合并粘附細(xì)胞的蓖麻毒素B凝集素的結(jié)構(gòu)(Spooner and Lord, 2015)。有研究表明在東方蜜蜂微孢子蟲侵染家蠶Bombyxmori卵巢BmN細(xì)胞的過程中，蓖麻毒素B凝集素能夠增強(qiáng)家蠶微孢子蟲對宿主細(xì)胞的粘附能力(Liuetal., 2016b)。Huang等(2016)研究發(fā)現(xiàn)東方蜜蜂微孢子蟲感染西方蜜蜂的次日就能檢測到蓖麻毒素編碼基因的表達(dá)，進(jìn)一步分析發(fā)現(xiàn)蓖麻毒素含有信號肽，同時(shí)宿主的核糖體蛋白受到抑制。本研究中，3個(gè)比較組共有的上調(diào)基因包含3個(gè)蓖麻毒素B凝集素編碼基因(GenBank登錄號分別為9422649, 9422645和9422719)(表2)，暗示東方蜜蜂微孢子蟲在侵染意大利蜜蜂工蜂的過程中通過增強(qiáng)蓖麻毒素B凝集素的分泌以加強(qiáng)對中腸上皮細(xì)胞的粘附力，同時(shí)抑制宿主細(xì)胞蛋白質(zhì)的合成。細(xì)胞膜蛋白糖基化是細(xì)胞識(shí)別和應(yīng)答外界環(huán)境的重要方式之一(劉嘯塵等, 2018)。本研究中，細(xì)胞表面糖基化蛋白編碼基因(GenBank登錄號: 9423191)在NcCKvsNcT1, NcCKvsNcT2和NcT1vsNcT2 3個(gè)比較組中分別上調(diào)1.74, 3.06和1.32倍(表2)，推測東方蜜蜂微孢子蟲在侵染過程中通過上調(diào)細(xì)胞表面糖基化蛋白基因的表達(dá)對宿主細(xì)胞內(nèi)環(huán)境產(chǎn)生應(yīng)答。此外，共同上調(diào)基因中包含1個(gè)假定蛋白編碼基因(GenBank登錄號: 9422567)，共同下調(diào)的1個(gè)基因也為假定蛋白編碼基因(GenBank登錄號: 9422077)(表2)，說明東方蜜蜂微孢子蟲的基因組功能注釋信息尚不完全，有待通過分子克隆、基因和蛋白功能研究進(jìn)一步完善。到目前為止，包括東方蜜蜂微孢子蟲在內(nèi)的絕大多數(shù)微孢子蟲的基因功能未知，最大的制約因素是缺乏轉(zhuǎn)基因操作技術(shù)平臺(tái)。Guo等(2016)通過堿液刺激家蠶微孢子蟲孢子體外發(fā)芽，繼而感染家蠶BmN細(xì)胞系，3 d后利用脂質(zhì)體包裹piggyBac轉(zhuǎn)座子載體和pIZT/V5-His非轉(zhuǎn)座子載體轉(zhuǎn)染家蠶微孢子蟲感染的細(xì)胞，通過連續(xù)的熒光觀察和分子生物學(xué)檢測證實(shí)將2種載體攜帶的gfp基因成功導(dǎo)入家蠶微孢子蟲基因組；此外，作者還通過飼喂家蠶微孢子蟲孢子感染家蠶個(gè)體，繼而注射脂質(zhì)體包裹的上述兩種載體，進(jìn)而也通過連續(xù)的熒光觀察和分子生物學(xué)檢測證實(shí)家蠶微孢子蟲基因組插入了gfp基因。目前，筆者團(tuán)隊(duì)正在參考該方法嘗試建立東方蜜蜂微孢子蟲的替代細(xì)胞連續(xù)感染系和轉(zhuǎn)基因操作技術(shù)體系。

NcCKvsNcT1和NcCKvsNcT2比較組包含的DEG富集在代謝進(jìn)程、細(xì)胞進(jìn)程、單組織進(jìn)程等生物學(xué)進(jìn)程相關(guān)功能條目，細(xì)胞、細(xì)胞組件和細(xì)胞器等細(xì)胞組分相關(guān)功能條目，結(jié)合、催化活性和轉(zhuǎn)運(yùn)活性等分子功能相關(guān)功能條目，表明東方蜜蜂微孢子蟲通過改變新陳代謝、細(xì)胞生命活動(dòng)相關(guān)基因的表達(dá)水平促進(jìn)自身增殖。NcT1vsNcT2比較組包含的DEG主要富集在催化活性、代謝進(jìn)程和細(xì)胞進(jìn)程等功能條目且多數(shù)基因上調(diào)表達(dá)，說明東方蜜蜂微孢子蟲在侵染意大利蜜蜂工蜂的過程中物質(zhì)和能量代謝活躍。進(jìn)一步對DEG進(jìn)行KEGG通路分析，發(fā)現(xiàn)NcCKvsNcT1和NcCKvsNcT2比較組包含的DEG大量富集在新陳代謝相關(guān)通路，包括新陳代謝途徑(分別為93和92個(gè))、次級代謝產(chǎn)物合成(分別為37和39個(gè))、抗生素合成(分別為29和31個(gè))、嘌呤代謝(分別為21和21個(gè))、嘧啶代謝(分別為20和19個(gè))、碳代謝(分別為16和18個(gè))、甘油磷脂代謝(分別為11和13個(gè))、糖酵解/糖異生(分別為11和12個(gè))等物質(zhì)代謝通路，以及氧化磷酸化(分別為10和10個(gè))和甲烷代謝(分別為4和6個(gè))等能量代謝通路(表3)。上述結(jié)果再次表明東方蜜蜂微孢子蟲侵染意大利蜜蜂工蜂的過程伴隨著復(fù)雜的物質(zhì)和能量代謝活動(dòng)，病原本身需要提升物質(zhì)和能量代謝滿足增殖需要，但同時(shí)也面臨宿主的抑制。真菌通過MAPK、cAMP和雙組份信號通路應(yīng)答外界環(huán)境，調(diào)控其生長發(fā)育和毒力 (侯彬彬等, 2015)。蜜蜂球囊菌是另一種常見的蜜蜂真菌病原，特異性侵染蜜蜂幼蟲而導(dǎo)致白堊病，筆者團(tuán)隊(duì)在前期研究中發(fā)現(xiàn)，對于侵染意大利蜜蜂幼蟲的蜜蜂球囊菌，有48個(gè)DEG富集在MAPK信號通路，它們的表達(dá)水平隨蜜蜂球囊菌脅迫時(shí)間的延長而逐漸增強(qiáng)(陳大福等, 2017)。本研究中，對于NcCKvsNcT1和NcCKvsNcT2比較組，分別有7和7個(gè)DEG富集在MAPK信號通路；筆者前期研究發(fā)現(xiàn)，對于NcCK, NcT1和NcT2，分別有5, 7和7個(gè)涉及MAPK信號通路的基因高量表達(dá)(熊翠玲等, 2019)；深入分析發(fā)現(xiàn)本研究中MAPK信號通路上有2個(gè)DEG與前期研究的高表達(dá)基因重疊，其中1個(gè)DEG(GenBank登錄號: 9424462)在NcCK, NcT1和NcT2均高量表達(dá)(FPKM值分別為966, 829和22)，另外1個(gè)DEG(GenBank登錄號: 9422306)在NcT1和NcT2高量表達(dá)(FPKM值分別為209和205)；上述結(jié)果表明MAPK信號通路及其富集DEG在東方蜜蜂微孢子蟲的侵染過程中發(fā)揮重要作用，下一步擬通過設(shè)計(jì)相應(yīng)的siRNA對GenBank登錄號分別為9424462和9422306的基因進(jìn)行RNAi，以探究其功能。

東方蜜蜂微孢子蟲侵染蜜蜂時(shí)通過中空的極管將具有侵染性的孢原質(zhì)注入宿主的中腸上皮細(xì)胞，因此極管蛋白在東方蜜蜂微孢子蟲的侵染過程扮演著關(guān)鍵角色(Rodríguez-Garcíaetal., 2018)。前人研究報(bào)道干擾極管蛋白3編碼基因的表達(dá)可顯著降低東方蜜蜂微孢子蟲孢子的數(shù)量(Rodríguez-Garcíaetal., 2018)。本研究中，在NcCK, NcT1和NcT2中極管蛋白1編碼基因(GenBank登錄號: 9423198)的表達(dá)量分別為171, 8 047和13 854；極管蛋白2編碼基因(GenBank登錄號: 9423201)的表達(dá)量分別為353, 17 919和27 591；極管蛋白編碼基因(GenBank登錄號: 9422340)的表達(dá)量分別為190, 2 923和3 994(表4)，表明東方蜜蜂微孢子蟲在侵染增殖過程中比孢子狀態(tài)上調(diào)表達(dá)，并且在侵染過程持續(xù)上調(diào)，印證了極管蛋白對東方蜜蜂微孢子蟲侵染的重要性(Rodríguez-Garcíaetal., 2018)。孢壁蛋白保護(hù)東方蜜蜂微孢子蟲在孢子狀態(tài)抵制外界不良環(huán)境，同時(shí)在應(yīng)答外界環(huán)境的變化傳遞信號和附著等作用(Southernetal., 2007)。孢壁蛋白還可以與極管蛋白互相作用，例如家蠶微孢子蟲的孢壁蛋白5與極管蛋白2和極管蛋白3相互作用，孢壁蛋白5對極管蛋白錨定在孢壁上起到非常重要的作用，孢壁蛋白9定位到極管上，并且可以與極管蛋白1和極管蛋白2互相作用(Yangetal., 2018)。本研究中，一共檢測到5個(gè)孢壁蛋白編碼基因(GenBank登錄號分別為9422782, 9422311, 9422437, 9422653和9423533)、1個(gè)孢壁蛋白前體編碼基因(GenBank登錄號: 9422389)和1個(gè)孢壁和錨定盤復(fù)合蛋白編碼基因(GenBank登錄號: 9424403)，其中孢壁蛋白編碼基因(GenBank登錄號分別為9422653和9423533)在NcCKvsNcT1和NcCKvsNcT2上調(diào)表達(dá)(表4)，推測它們主要在孢子侵染和增殖過程中起到重要的作用，而孢壁和錨定盤復(fù)合蛋白編碼基因(GenBank登錄號: 9424403)在NcCKvsNcT1和NcCKvsNcT2上調(diào)表達(dá)，推測其在東方蜜蜂微孢子蟲侵染增殖過程發(fā)揮重要的作用。此外，孢壁蛋白9編碼基因在NcCK中的表達(dá)量很高(FPKM值為12 059)，而在NcT1和NcT2的表達(dá)量相對較低(FPKM值分別為2 807和2 885)，暗示其在東方蜜蜂微孢子蟲孢子抵御外界不良環(huán)境中起特殊作用。本研究發(fā)現(xiàn)，另有4個(gè)蓖麻毒素B凝集素編碼基因(GenBank登錄號分別為9422642, 9423906, 9422644和9422647)分別在NcCKvsNcT1(分別為4.77, 2.32, 2.72和5.62倍), NcCKvsNcT2(分別為4.98, 3.23, 3.21和3.38倍)中上調(diào)表達(dá)(表4)，再次說明蓖麻毒素B凝集素在東方蜜蜂微孢子蟲侵染過程中的重要性。

東方蜜蜂微孢子蟲寄生蜜蜂中腸上皮細(xì)胞，由于缺少典型的線粒體，所以自身只能通過糖酵解途徑產(chǎn)生少量ATP，絕大多數(shù)的ATP必需依賴宿主細(xì)胞提供(Cornmanetal., 2009)。Huang等(2016)通過對東方蜜蜂微孢子蟲脅迫的西方蜜蜂工蜂1-6 d的腸道進(jìn)行測序，發(fā)現(xiàn)東方蜜蜂微孢子蟲催化糖酵解的關(guān)鍵酶的基因在感染后1 d就開始表達(dá)，己糖激酶在感染后2 d被檢測到并且在2-6 d的表達(dá)量上調(diào)，糖酵解途徑的其他核心酶在感染后3 d也被檢測出來，同時(shí)發(fā)現(xiàn)ATP/ADP移位酶基因和ABC轉(zhuǎn)運(yùn)蛋白基因在東方蜜蜂微孢子蟲侵染過程中也維持高量表達(dá)。本研究中，涉及東方蜜蜂微孢子蟲糖酵解的幾種關(guān)鍵酶基因如己糖激酶基因(GenBank登錄號: 9423428)、丙酮酸激酶基因(GenBank登錄號: 9424331)和6-磷酸果糖激酶基因(GenBank登錄號: 9424162)在侵染過程較孢子狀態(tài)表達(dá)量上調(diào)(表5)，與前人研究結(jié)果(Huangetal., 2016)一致，再次說明東方蜜蜂微孢子蟲在其增殖過程中通過糖酵解途徑提供部分能量。有研究表明家蠶微孢子蟲和蝗蟲微孢子蟲Nosemalocustae的己糖激酶都包含信號肽，均為分泌蛋白(Timofeevetal., 2017; 陳紅麗等, 2018)。本研究發(fā)現(xiàn)，東方蜜蜂微孢子蟲的己糖激酶同樣也包含信號肽，推測東方蜜蜂微孢子蟲在侵染過程中將己糖激酶分泌到宿主細(xì)胞質(zhì)以調(diào)控宿主的糖酵解。ATP/ADP移位酶是微孢子蟲竊取宿主能量的關(guān)鍵酶之一，Paldi等(2010)通過RNAi對東方蜜蜂微孢子蟲的ATP/ADP移位酶基因進(jìn)行敲減，發(fā)現(xiàn)東方蜜蜂微孢子蟲感染的西方蜜蜂體內(nèi)的病原增殖被顯著抑制。本研究中，有3個(gè)ATP/ADP移位酶的編碼基因(GenBank登錄號分別為9424363, 9422880和9424586)在NcCKvsNcT1 [log2(Fold change)值分別為0.72, 3.60和2.31]和NcCKvsNcT2[log2(Fold change)值分別為0.45, 3.44和2.34]中的表達(dá)水平上調(diào)(表5)，表明ATP/ADP移位酶對于東方蜜蜂微孢子蟲竊取宿主細(xì)胞能量、促進(jìn)自身增殖的重要性。ABC轉(zhuǎn)運(yùn)蛋白的主要功能是結(jié)合ATP，利用ATP水解產(chǎn)生能量將物質(zhì)逆濃度梯度運(yùn)輸(馮振月等, 2018)。本研究中，GenBank登錄號分別為9422518和9423041的ABC轉(zhuǎn)運(yùn)蛋白編碼基因在NcCKvsNcT1 [log2(Fold change)值分別為0.93和1.83]和NcCKvsNcT2[log2(Fold change)值分別為0.64和1.85]中上調(diào)表達(dá)，體現(xiàn)出對東方蜜蜂微孢子蟲侵染過程中物質(zhì)運(yùn)輸?shù)闹匾裕坏硗?個(gè)GenBank登錄號分別為9423050, 9423588, 9422298和9422297的ABC轉(zhuǎn)運(yùn)蛋白編碼基因在NcCKvsNcT1 [log2(Fold change)值分別為-3.65, -1.72, -1.91和-0.19] 和NcCKvsNcT2[log2(Fold change)值分別為-3.34, -0.13, -1.56和-0.35]中的表達(dá)水平也為下調(diào)(表5)，鑒于東方蜜蜂微孢子蟲與蜜蜂之間存在密切的相互作用，上述3個(gè)基因的下調(diào)表達(dá)可能是由于東方蜜蜂微孢子蟲受到宿主的抑制，具體的機(jī)制有待于進(jìn)一步深入研究。

本研究解析了東方蜜蜂微孢子蟲在侵染意大利蜜蜂工蜂過程的差異基因表達(dá)譜，在mRNA組學(xué)層面揭示了東方蜜蜂微孢子蟲的侵染機(jī)制，篩選出病原的毒力因子和侵染相關(guān)因子基因可為后續(xù)功能研究提供候選靶標(biāo)，并為治療蜜蜂微孢子蟲病提供潛在的分子靶點(diǎn)。