劉一名, 徐新建, 周姝婧, 姚 丹, 李 寒, 朱晨煜, 李 想, 王茗琦, 周冰峰, 朱翔杰

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

溫度是影響昆蟲的外部形態(tài)[1-2]、行為[3]、繁殖能力[4]、生理功能[5-6]、能量代謝[7]等的環(huán)境因子之一.作為典型的社會性昆蟲，蜜蜂具有調(diào)節(jié)巢溫的能力，蜂群可以通過群體行為將蜂群子區(qū)的中心溫度嚴格維持在35 ℃,以使卵、蟲、蛹正常發(fā)育[8-10].當溫度偏離最適溫度35 ℃時，蜜蜂會出現(xiàn)死亡率增大、發(fā)育歷期延長、外部形態(tài)異常、記憶能力減弱等現(xiàn)象[11-15].可見，蜜蜂對溫度比較敏感[16].

轉(zhuǎn)錄組學(xué)是從轉(zhuǎn)錄層面上研究機體在某一狀態(tài)下的基因表達水平及轉(zhuǎn)錄調(diào)控規(guī)律[17]，是目前在昆蟲中廣泛應(yīng)用的組學(xué)研究手段[18].研究表明，低溫會誘導(dǎo)斜紋夜蛾(Spodopteralitura)能量代謝從碳水化合物代謝向脂質(zhì)代謝轉(zhuǎn)變，并且能降低游離氨基酸的水平[19].中華蜜蜂(Apisceranacerana)成蜂經(jīng)低溫脅迫后的差異表達基因(differentially expressed genes, DEGs)在代謝通路上富集最多，且多數(shù)是與糖、氨基酸的生物合成和代謝相關(guān)的通路，如甘氨酸-絲氨酸-蘇氨酸代謝、碳代謝、磷酸肌醇代謝、脂肪酸代謝等，同時在Wnt、Hippo、FoxO等與生長發(fā)育相關(guān)的信號通路上也有富集[20].西方蜜蜂(A.mellifera)預(yù)蛹受到低溫脅迫后的DEGs也集中在代謝和生長發(fā)育通路上[21].

當蜜蜂幼蟲發(fā)育到6日齡時，蜂群會將巢房用蜂蠟封閉，不再需要工蜂飼喂食物，所以該齡期幼蟲是研究低溫脅迫影響的理想蟲態(tài).本研究團隊前期研究顯示，與預(yù)蛹和蛹早期相比，蜜蜂6日齡幼蟲對低溫的敏感性降低，20 ℃低溫脅迫48 h對其死亡率沒有影響[11].Mucci et al[22]研究顯示，西方蜜蜂幼蟲受到低溫脅迫時，抗氧化能力提高.因此，推測西方蜜蜂6日齡幼蟲可能通過提高抗氧化能力調(diào)節(jié)自身響應(yīng)低溫的能力.但具體通過哪些通路上的關(guān)鍵基因表達變化提高抗氧化能力，以及蜜蜂幼蟲的生長發(fā)育在轉(zhuǎn)錄組水平上會受到怎樣的影響仍需進一步明確.本試驗將西方蜜蜂工蜂6日齡幼蟲進行20 ℃低溫處理，利用轉(zhuǎn)錄組測序技術(shù)分析受低溫脅迫后的DEGs，為后續(xù)深入研究該蟲響應(yīng)低溫脅迫的機制提供依據(jù).

1 材料與方法

1.1 蟲源

西方蜜蜂工蜂幼蟲取自福建省福州市福建農(nóng)林大學(xué)實驗蜂場(N 26°5′，E 119°13′)4—6月的蜂群.蜂群群勢12足框，蜂王產(chǎn)卵正常.將空脾分別置于3個產(chǎn)卵蜂群中限王產(chǎn)卵1 d，以獲取足夠多卵齡一致的卵，再將卵脾放入同一哺育群哺育幼蟲，幼蟲經(jīng)哺育6 d后封蓋.為保證樣本發(fā)育一致性，取2 h內(nèi)封蓋的6日齡幼蟲，分別在正常發(fā)育條件[(35±0.2) ℃，RH 75%, CK]和低溫條件[(20±0.2) ℃，RH 75%, T]下培養(yǎng)4 h后，取樣本立即用液氮冷凍，-80 ℃保存.從3個蜂群各取10頭幼蟲作為每組樣本的3個生物學(xué)重復(fù).

1.2 高通量測序及測序數(shù)據(jù)質(zhì)控

利用Trizol法提取全蟲總RNA后，委托廣州基迪奧生物科技有限公司使用Illumina HiSeqTM 4000平臺對樣品進行建庫和RNA-seq測序.下機數(shù)據(jù)使用fastp v0.18.0過濾后，去除含adapter、N比例大于10%、全部都是A堿基、質(zhì)量值Q≤20的堿基占50%以上的reads，得到clean reads[23]，然后與西方蜜蜂的基因組Amel_HAv3.1(NCBI Assembly: GCF_003254395.2)進行比對.

1.3 樣本相關(guān)性分析

比較CK組和T組3個生物學(xué)重復(fù)樣品測序得到的基因之間的pearson相關(guān)系數(shù)，使用Omicshare生信云工具平臺(http://www.omicshare.com/tools/)進行分析.

1.4 DEGs篩選及功能分析

DEGs的篩選使用DESeq2[24]軟件，篩選條件為|log2(FC)|≥1且錯誤發(fā)現(xiàn)率(false discovery rate, FDR)≤0.05.其中，F(xiàn)C為差異倍數(shù)(fold change)，是T組與CK組之間3個生物學(xué)重復(fù)的平均FPKM值的比值.使用Omicshare生信云工具平臺對DEGs進行GO(gene ontology)功能分類和KEGG(Kyoto Encyclopedia of Genes and Genomes)富集分析，P<0.05為顯著富集GO功能注釋和KEGG通路.

1.5 實時熒光定量PCR(qPCR)驗證

挑選5個在KEGG通路或GO功能注釋上富集的DEGs進行qPCR定量檢測.引物使用NCBI Primer BLAST設(shè)計(表1).將上述測序所用的總RNA反轉(zhuǎn)錄得到cDNA，以此作為模板進行目的基因的擴增，qPCR反應(yīng)體系參照PerfectStartTMGreen qPCR SuperMix試劑盒(北京全式金生物科技有限公司)說明書.qPCR程序：95 ℃ 30 s；95 ℃ 5 s，60 ℃ 30 s(40個循環(huán))；熔解曲線條件為95 ℃ 5 s，60 ℃ 30 s，95 ℃ 1 s.以Actin作為內(nèi)參基因，利用2-ΔΔCt法[25]計算目的基因相對表達量，每個目的基因進行3次生物學(xué)重復(fù)和3次技術(shù)重復(fù)后取平均值，再進行l(wèi)og2(平均值)的轉(zhuǎn)換后，與測序數(shù)據(jù)進行對比.

表1 西方蜜蜂5個DEGs的qPCR引物信息Table 1 qPCR primers for 5 selected DEGs in A.mellifera

2 結(jié)果與分析

2.1 測序數(shù)據(jù)質(zhì)量分析

通過轉(zhuǎn)錄組測序得到CK組和T組的原始reads(單位長度150 bp)分別為87 593 339±6 594 537和106 975 591±5 116 984，過濾后得到高質(zhì)量有效reads分別為87 434 485±6 586 202和106 799 168±5 107 747，高質(zhì)量reads占比分別為99.82%±0.04%和99.84%±0.04%.各樣本有效堿基數(shù)據(jù)量均大于11 Gb，CK組和T組各個樣品兩端的平均Q20、 Q30分別為97.42%、92.83%和92.40%、92.68%(表2)，表明測序數(shù)據(jù)質(zhì)量較好，可用于后續(xù)分析.CK組和T組的3個生物學(xué)重復(fù)樣本間的pearson相關(guān)系數(shù)均大于0.968 8，顯示各組內(nèi)樣本極強相關(guān)，樣本重復(fù)性好.

表2 堿基過濾前后質(zhì)量統(tǒng)計1)Table 2 Base quality before and after filtration

2.2 DEGs篩選

經(jīng)過篩選，CK組與T組之間共檢測到98個響應(yīng)低溫脅迫的DEGs，其中87個上調(diào)表達，11個下調(diào)表達(圖1).

FC(fold change)為基因表達差異倍數(shù)，F(xiàn)DR(false discovery rate)為錯誤發(fā)現(xiàn)率.圖1 低溫處理組與正常發(fā)育組6日齡幼蟲的差異基因火山圖Fig.1 Volcano map of DEGs in 6-day-old larva grew under low temperature stress and optimal temperature

2.3 DEGs的功能富集分析

2.3.1 GO功能分類 西方蜜蜂工蜂6日齡幼蟲響應(yīng)低溫脅迫的基因富集在21條二級GO功能注釋上.其中，上調(diào)DEGs富集在21條注釋上，下調(diào)DEGs富集在15條注釋上(圖2).

圖2 低溫脅迫后DEGs的GO功能分類Fig.2 GO functional classification of DEGs in response to low temperature stress

上調(diào)DEGs中，26個富集在生物過程，23個富集在分子功能，12個富集在細胞組分.上調(diào)DEGs富集數(shù)最多的前5條二級GO功能注釋分別為細胞進程(19個)、結(jié)合(18個)、代謝進程(16個)、催化活性(14個)、單有機體進程(12個).顯著富集的GO功能注釋共有43條，其中，生物過程上的GO功能注釋有38條，多數(shù)與級聯(lián)反應(yīng)、應(yīng)激反應(yīng)、激酶活性相關(guān)(圖3).

富集因子指DEGs中位于該功能注釋的基因數(shù)與所有基因中位于該功能注釋的基因總數(shù)的比值.圖3 低溫脅迫后上調(diào)DEGs在生物過程中顯著富集的GO功能注釋Fig.3 GO enrichment analysis of up-regulated DEGs in response to low temperature stress

下調(diào)DEGs中，3個富集在生物過程，5個富集在分子功能，3個富集在細胞組分(圖2).顯著富集的GO功能注釋有6條，其中4條屬于生物過程.

2.3.2 KEGG富集分析 KEGG富集分析顯示，西方蜜蜂工蜂6日齡幼蟲受低溫脅迫后，上調(diào)DEGs富集在14條B級通路(圖4)的25條KEGG通路上，其中，代謝相關(guān)通路富集的DEGs最多.代謝相關(guān)的KEGG通路有氨基酸代謝上的賴氨酸降解，碳水化合物代謝的磷酸肌醇代謝，輔酶因子和維生素代謝的核黃素代謝，萜類和聚酮肽代謝的萜類骨架生物合成、昆蟲激素生物合成，核苷酸代謝的嘧啶代謝、嘌呤代謝.其中，萜類骨架生物合成、核黃素代謝顯著富集(圖5).除代謝相關(guān)通路外，富集在信號轉(zhuǎn)導(dǎo)B級通路上的DEGs也較多(圖4)，信號轉(zhuǎn)導(dǎo)的KEGG通路有Hippo信號通路、Notch信號通路、磷脂酰肌醇信號系統(tǒng).

縱坐標中彩色字體表示KEGG A級通路，黑色字體表示相應(yīng)A級通路下的B級通路.圖4 低溫脅迫后DEGs的KEGG富集分析Fig.4 KEGG enrichment analysis for DEGs in response to low temperature stress

西方蜜蜂工蜂6日齡幼蟲受低溫脅迫后，下調(diào)DEGs富集在3條B級通路(圖4)的3條KEGG通路上，分別為代謝方面的代謝途徑、甘氨酸-絲氨酸-蘇氨酸代謝和信號轉(zhuǎn)導(dǎo)方面的MAPK信號通路.其中，屬于氨基酸代謝的甘氨酸-絲氨酸-蘇氨酸代謝顯著富集(圖5).葡萄糖脫氫酶(GCDH, LOC551044)和葡萄糖氧化酶(GOD, LOC406081)同時富集在代謝途徑和甘氨酸-絲氨酸-蘇氨酸代謝上.

富集因子指DEGs中位于該通路的基因數(shù)與所有基因中位于該通路的基因總數(shù)的比值.圖5 低溫脅迫后DEGs顯著富集的KEGG通路Fig.5 Significant enriched KEGG pathways of DEGs in response to low temperature stress

2.4 DEGs的qPCR驗證

挑選IP3K、Ex、CYP18A1、GOD、GCDH5個DEGs進行qPCR檢測，結(jié)果顯示，與對照組相比，低溫脅迫后基因表達差異顯著，基因IP3K、Ex、CYP18A1均顯著上調(diào)表達，基因GOD和GCDH均顯著下調(diào)表達.這與轉(zhuǎn)錄組測序結(jié)果的基因表達變化趨勢一致(圖6).

圖中qPCR數(shù)據(jù)為平均值±標準差.柱上星號代表處理組與對照組之間基因表達差異顯著(非配對T檢驗)：*P<0.05;**P<0.01;***P<0.001.圖6 5個西方蜜蜂DEGs的qPCR檢測結(jié)果Fig.6 qPCR verification of 5 selected DEGs in A.mellifera

3 討論

3.1 西方蜜蜂工蜂6日齡幼蟲可能通過提高抗氧化能力應(yīng)對低溫脅迫

核黃素是黃素單核苷酸(flavin mononucleotide, FMN)和黃素腺嘌呤二核苷酸(flavin adenine dinucleotide, FAD)的組成部分和前體，其以FMN和FAD的形式在昆蟲體內(nèi)參與廣泛的氧化還原反應(yīng)[26]，對細胞功能、昆蟲生長和發(fā)育至關(guān)重要[27].有研究顯示，果蠅受到氧化損傷時，體內(nèi)核黃素會通過提高超氧化物歧化酶和過氧化氫酶的活性增強抗氧化能力，使果蠅壽命明顯延長[28].本研究表明，核黃素代謝通路在受低溫脅迫的西方蜜蜂工蜂6日齡幼蟲中顯著富集，黃素腺嘌呤二核苷酸合酶(FADS)基因上調(diào)表達.這說明西方蜜蜂可能以此提高體內(nèi)核黃素抗氧化能力，延長自身壽命.這與本團隊前期研究的蜜蜂6日齡幼蟲在所有蟲態(tài)低溫處理中壽命最長的結(jié)果[11]相吻合.GOD可將葡萄糖氧化為葡萄糖酸和過氧化氫[29]；有研究顯示，過氧化氫會抑制絲腺細胞的程序死亡[30].本研究中，西方蜜蜂工蜂6日齡幼蟲在受到低溫脅迫后，GOD基因下調(diào)表達，可能會減少葡萄糖的氧化，增強幼蟲抗氧化能力，且過氧化氫的減少可能會造成絲腺細胞程序性死亡，使幼蟲吐絲受到影響.還有研究顯示，幼蟲代謝旺盛時，GOD活性相對較高[31-32].低溫脅迫后GOD基因下調(diào)表達可能也是導(dǎo)致幼蟲取食量下降的原因.

3.2 西方蜜蜂工蜂6日齡幼蟲還可能通過抑制細胞生長應(yīng)對低溫脅迫

膨化蛋白(Ex)是調(diào)控Hippo信號通路的關(guān)鍵蛋白.Hippo信號通路是廣泛存在于多細胞生物中的信號通路[33]，其主要功能是抑制細胞增殖和促進細胞凋亡[34].在Hippo信號通路上游，Ex蛋白與Kibra和Mer蛋白結(jié)合構(gòu)成KEM復(fù)合體，起到激活Hippo信號通路、抑制細胞生長的作用[35].本研究中，西方蜜蜂工蜂6日齡幼蟲在受到低溫脅迫后，Ex基因顯著上調(diào)表達，激活Hippo信號通路，可能會抑制幼蟲的細胞增殖而使其發(fā)育變慢.Notch通路的功能是促進細胞的增殖和分化[36].Numb蛋白通過調(diào)控Notch delt樣配體4(Dll4)在細胞內(nèi)的位置和穩(wěn)定性抑制Notch通路表達[37].西方蜜蜂工蜂6日齡幼蟲在受到低溫脅迫后，富集在Notch通路上的Numb基因也顯著上調(diào)表達，可能以此抑制Notch通路而使幼蟲生長發(fā)育變慢.細胞色素P45018A1(CYP18A1)基因參與激素生物合成通路，編碼一種具有26-羥化酶活性的細胞色素P450酶；26-羥化酶對于20-羥基蛻皮激素(20E)分解代謝具有重要作用，其能使20E分解為20-羥基蛻皮激素酸，使其不能與下游基因結(jié)合而調(diào)控蛻皮；CYP18A1基因缺失會導(dǎo)致果蠅幼蟲期延長和蛹的死亡[38].西方蜜蜂工蜂6日齡幼蟲在受到低溫脅迫后，CYP18A1基因上調(diào)表達，可能會使幼蟲中20E的滴度少于正常發(fā)育的幼蟲，影響幼蟲蛻皮的過程.

細胞代謝與信號轉(zhuǎn)導(dǎo)緊密相關(guān)：信號傳遞網(wǎng)絡(luò)可以促進細胞的代謝活性[39]，代謝物也可以控制信號蛋白的活性[40].西方蜜蜂工蜂6日齡幼蟲通過代謝和信號轉(zhuǎn)導(dǎo)的共同作用對低溫做出響應(yīng).其響應(yīng)低溫脅迫的基因FADS、GOD、CYP18A1富集的通路分別參與輔酶及維生素代謝、氨基酸代謝、萜類和聚酮肽代謝；Ex、Numb分別富集的Hippo、Notch信號通路均屬于信號轉(zhuǎn)導(dǎo).因此，西方蜜蜂工蜂6日齡幼蟲對低溫的抵抗能力是由多種代謝反應(yīng)和信號轉(zhuǎn)導(dǎo)途徑共同調(diào)控的.后續(xù)可以從轉(zhuǎn)錄組學(xué)檢測到的代謝與信號轉(zhuǎn)導(dǎo)通路入手，研究二者在蜜蜂響應(yīng)低溫脅迫中的作用機制.