轉(zhuǎn)錄組學(xué)分析低溫脅迫對(duì)西方蜜蜂工蜂6日齡幼蟲(chóng)基因表達(dá)的影響

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

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

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

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

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

1 材料與方法

1.1 蟲(chóng)源

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

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

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

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

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

1.4 DEGs篩選及功能分析

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

1.5 實(shí)時(shí)熒光定量PCR(qPCR)驗(yàn)證

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

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

2 結(jié)果與分析

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

通過(guò)轉(zhuǎn)錄組測(cè)序得到CK組和T組的原始reads(單位長(zhǎng)度150 bp)分別為87 593 339±6 594 537和106 975 591±5 116 984，過(guò)濾后得到高質(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組各個(gè)樣品兩端的平均Q20、 Q30分別為97.42%、92.83%和92.40%、92.68%(表2)，表明測(cè)序數(shù)據(jù)質(zhì)量較好，可用于后續(xù)分析.CK組和T組的3個(gè)生物學(xué)重復(fù)樣本間的pearson相關(guān)系數(shù)均大于0.968 8，顯示各組內(nèi)樣本極強(qiáng)相關(guān)，樣本重復(fù)性好.

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

2.2 DEGs篩選

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

FC(fold change)為基因表達(dá)差異倍數(shù)，F(xiàn)DR(false discovery rate)為錯(cuò)誤發(fā)現(xiàn)率.圖1 低溫處理組與正常發(fā)育組6日齡幼蟲(chóng)的差異基因火山圖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日齡幼蟲(chóng)響應(yīng)低溫脅迫的基因富集在21條二級(jí)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個(gè)富集在生物過(guò)程，23個(gè)富集在分子功能，12個(gè)富集在細(xì)胞組分.上調(diào)DEGs富集數(shù)最多的前5條二級(jí)GO功能注釋分別為細(xì)胞進(jìn)程(19個(gè))、結(jié)合(18個(gè))、代謝進(jìn)程(16個(gè))、催化活性(14個(gè))、單有機(jī)體進(jìn)程(12個(gè)).顯著富集的GO功能注釋共有43條，其中，生物過(guò)程上的GO功能注釋有38條，多數(shù)與級(jí)聯(lián)反應(yīng)、應(yīng)激反應(yīng)、激酶活性相關(guān)(圖3).

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

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

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

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

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

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

2.4 DEGs的qPCR驗(yàn)證

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

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

3 討論

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

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

3.2 西方蜜蜂工蜂6日齡幼蟲(chóng)還可能通過(guò)抑制細(xì)胞生長(zhǎng)應(yīng)對(duì)低溫脅迫

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

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