姚 丹, 周姝婧, 周冰峰, 徐新建, 朱翔杰

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

蜜蜂不僅為人類提供了豐富多樣的蜂產(chǎn)品，而且是農(nóng)業(yè)中最重要的授粉昆蟲，在維護生態(tài)平衡發(fā)揮重要作用(Allen-Wardelletal., 1998)。近年來，全球氣候變化和病蟲害危害導致蜜蜂的種群數(shù)量大幅下降(Milbrathetal., 2015)，蜜蜂蜂群的減少引起了人們對蜜蜂健康的極大關注。

環(huán)境溫度是影響昆蟲分布、形態(tài)、生理、行為等的主要生態(tài)因素之一(Burnett, 1949; Ponsonby and Copland, 1996)。 獨居性昆蟲可通過休眠或滯育等方式耐受低溫等不良環(huán)境，生態(tài)幅一般較寬(Tauber and Tauber, 1976)。而蜜蜂的卵、蟲、蛹(統(tǒng)稱為蜂子)是在蜂群維持較穩(wěn)定的巢溫下發(fā)育的，因此形成了狹溫性特點，發(fā)育溫區(qū)為29～38℃，最適發(fā)育溫度為35℃(Tautzetal., 2003; Grohetal., 2004; 朱翔杰等, 2006)。偏離最適溫度，羽化后成蜂的外部形態(tài)(Medrzyckietal., 2010)、翅脈發(fā)育(周冰峰等, 2011; Zhuetal., 2018)、神經(jīng)突觸數(shù)量(Grohetal., 2004)、壽命(Wangetal., 2016)和行為(Jonesetal., 2005; Becheretal., 2009)等均受到影響。蜜蜂發(fā)育對溫度的敏感性，使得蜂子成為研究溫度對個體發(fā)育影響的理想材料，有助于認識溫度與生命的關系。

低溫對昆蟲的影響主要表現(xiàn)在糖代謝、脂類代謝、氨基酸代謝、嘌呤代謝和硫胺素代謝等代謝通路上富集大量差異表達基因(differentially expressed genes, DEGs)(庫爾班·吐松等, 2016; Tangetal., 2017; Wangetal., 2017; Xuetal., 2017; Qinetal., 2019)。在0℃下20日齡中華蜜蜂Apisceranacerana成蜂轉(zhuǎn)錄組測序檢測到與環(huán)境進程相關的Ca2+信號通路、FoxO信號通路、Hippo信號通路等，同時熱激蛋白、鋅指蛋白和絲氨酸/蘇氨酸蛋白激酶等差異表達(Xuetal., 2017)。

相對其他日齡的蜜蜂封蓋子，3日齡封蓋子處于幼蟲-蛹變態(tài)的過渡期，表現(xiàn)出對低溫脅迫最敏感。前期研究發(fā)現(xiàn)，不同日齡封蓋子經(jīng)20℃低溫脅迫后3日齡封蓋子最早出現(xiàn)大量封蓋子死亡，羽化蜜蜂的壽命也最短(Wangetal., 2016)。到目前為止，尚無研究對蜜蜂封蓋子發(fā)育期間低溫敏感的機制給出解釋。為了深入研究低溫對蜜蜂化蛹影響的分子機制，本研究采用轉(zhuǎn)錄組學方法對低溫處理不同時間的意大利蜜蜂Apismelliferaligustica3日齡封蓋子的共有差異表達基因進行趨勢分析，從分子水平探討蜜蜂化蛹時受低溫影響的關鍵基因和關鍵通路。本研究對豐富昆蟲生態(tài)學和昆蟲發(fā)育生物學有重要意義。

1 材料與方法

1.1 供試蜜蜂

樣本取自春季增長階段的意大利蜜蜂A.m.ligustica蜂群，通過限制蜂王產(chǎn)卵空間，獲得卵齡一致的蜜蜂卵。3 d后將卵脾放入同一個哺育群哺育。移入哺育群5 d后這些幼蟲即將封蓋，獲取4 h內(nèi)同時封蓋的封蓋子(封蓋巢房內(nèi)的幼蟲和蛹)，割取樣本巢脾，放入恒溫恒濕箱(35±0.2℃, RH 75%)[CTHI-250B, 施都凱儀器設備(上海)有限公司]中培育至3日齡(對照組，CK)。再將一部分3日齡封蓋子放入溫度20±0.2℃、RH 75%(Wangetal., 2016)的恒溫恒濕箱中分別處理18 h和36 h(分別標記為T18和T36)，作為低溫處理組。所取的封蓋子樣本立即用液氮冷凍，放入-80℃保存?zhèn)溆?。每個樣本用30只封蓋子(3個蜂群，10頭/群)混合建庫和高通量測序。

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

利用Trizol法提取1.1節(jié)樣品全蟲總RNA，用帶有Oligo(dT)的磁珠富集mRNA后，加入fragmentation buffer使其片段成為短片段的mRNA為模板。用六堿基隨機引物(random hexamers)合成cDNA第1鏈，并加入緩沖液、dNTPs、RNase H和DNA polymerase I合成cDNA第2鏈，經(jīng)過QiaQuick PCR試劑盒純化并加EB緩沖液洗脫經(jīng)末端修復、加堿基A，加測序接頭，再經(jīng)瓊脂糖凝膠電泳回收目的大小片段，并進行PCR擴增，從而完成整個文庫制備工作，構建好cDNA文庫，高通量測序委托廣州基迪奧生物科技有限公司完成，測序平臺為Illumina HiSeqTM。對下機數(shù)據(jù)的raw reads進行過濾，去除含adapter的reads、含N比例大于10%的reads和低質(zhì)量的reads，最終得到clean reads，然后利用Tophat2(2.1.1)軟件將reads比對到西方蜜蜂Apismellifera的參考基因組(ncbi_GCF_000002195.4)上，并利用Cufflinks組裝轉(zhuǎn)錄本。 轉(zhuǎn)錄組原始數(shù)據(jù)已上傳NCBI的SRA數(shù)據(jù)庫(獲取號: SRR10912937-SRR10912945)。

1.3 差異表達基因趨勢分析

使用DESeq2軟件(V1.20.0)對有生物學重復的處理組與對照組間(T18vsCK, T36vsCK)基因表達量進行差異分析(Loveetal., 2014)， DEGs篩選條件為FDR(錯誤發(fā)現(xiàn)率)≤0.05, |log2Fold change| (Fold change, FC, 指基因表達水平的變化倍數(shù), FC=T/CK，其中T為處理組的基因表達水平， CK為對照組的基因表達水平)≥1。利用STEM(Short Time-series Expression Miner)，選擇參數(shù){-pro 20 -ratio 1.0[log2(2)=1, log2(1.5)=0.5849625, log2(1.2)=0.2630344]}(Ernst and Bar-Joseph, 2006),對2個比較組間共有的差異表達基因進行趨勢分析，在軟件中設置P≤0.05為差異顯著的表達模式。

1.4 差異表達基因GO及KEGG pathway富集分析

對各個趨勢中的基因進行GO(Gene Ontology)和KEGG(Kyoto Encyclopedia of Genes and Genomes)代謝通路進行功能注釋(Omicshare tools, http:∥www.omicshare. com/tools/home/index/index.html)，并通過假設檢驗計算得到P值。得到的P值通過FDR校正之后，以P值≤0.05為閾值，滿足此條件的GO條目(term)和通路(pathway)定義為在該趨勢中顯著富集的GO term和Pathway。

1.5 RT-qPCR基因驗證

隨機挑選5個趨勢顯著的基因進行RT-qPCR驗證,β-actin為內(nèi)參基因(引物見表1)。利用RNA抽提試劑盒提取樣本總RNA，利用反轉(zhuǎn)錄試劑盒合成cDNA第1鏈(北京全式金生物技術有限公司)。以上述的cDNA作為模板進行qPCR，反應體系(20 μL): H2O 6.6 μL, 2×TransStart SuperMix 10 μL, 上下游引物(5 μmol/L)各1 μL, cDNA模板1 μL,Dye II 0.4 μL。qPCR反應條件: 95℃預變性4 min; 95℃變性15 s, 60℃退火和延伸45 s，共40個循環(huán)；熔解曲線默認系統(tǒng)程序(QuantStudio5, 美國)。樣本取自1.1節(jié)，對照組和處理組的封蓋子均來自3個蜂群，每個蜂群取1頭進行混合提取總RNA。每個樣本設置3次目的基因的重復擴增。

表1 RT-qPCR引物信息Table 1 Primers for RT-qPCR

1.6 數(shù)據(jù)分析

目的基因的RT-qPCR表達水平按照內(nèi)參基因的2-ΔΔCt法進行計算，即ΔCt=Ct目的基因-Ct內(nèi)參基因(Livak and Schmittgen, 2001)，各基因的相對表達量經(jīng)log2(X)轉(zhuǎn)換后與高通量測序中的表達水平進行比較。應用SPSS V21.0軟件對每個基因處理組與對照組之間平均相對表達量的差異顯著性進行非配對T檢驗。

2 結果

2.1 測序結果質(zhì)量

通過對CK, T18和T36測序，分別得到原始讀段平均為56 759 342, 45 977 301和54 196 127，過濾后得到高質(zhì)量有效讀段平均分別為55 578 624, 45 032 615和53 121 097。統(tǒng)計兩端的Q20和Q30分別平均在98.20%和94.85%以上。表明本研究測序數(shù)據(jù)質(zhì)量良好，可用于進一步分析。

2.2 差異表達基因趨勢

通過差異表達基因(DEGs)分析，T18vsCK有1 394個DEGs，T36vsCK有2 722個DEGs，說明隨著低溫處理的時間增長，受影響的基因數(shù)明顯增加。其中兩比較組共有DEGs 1 062個，T18vsCK特有DEGs有332個，T36vsCK特有DEGs有1 660個，上調(diào)趨勢的DEGs顯著多于下調(diào)趨勢的DEGs。

對兩比較組共有的1 062個DEGs進行趨勢分析，在8個基因表達模式中有3個為顯著基因表達模式(P<0.05)，包括2個上調(diào)表達模式，Profile 6模式為基因表達量先增加后保持不變，有539個基因；Profile 7模式為基因表達量隨低溫脅迫時間延長一直上調(diào)，有271個基因；1個下調(diào)表達模式Profile 1，為隨低溫脅迫時間延長基因表達量先下調(diào)表達，之后保持不變，有183個基因(圖1)。多數(shù)基因歸入上述3個顯著基因表達模式。

圖1 低溫脅迫不同時間3日齡意大利蜜蜂封蓋子差異表達基因表達趨勢Fig. 1 Trend analysis of differentially expressed genes (DEGs) in the 3 d-old post-capped brood of Apis mellifera ligustica subjected to low temperature stress for different timeA: DEGs的趨勢分析結果，趨勢類型方框內(nèi)數(shù)字為P值The result of trend analysis of DEGs, the number in the box denoting the P-value; B: Profile 6包含DEGs的表達模式Expression profile of DEGs in Profile 6; C: Profile 7包含DEGs的表達模式Expression profile of DEGs in Profile 7; D: Profile 1包含DEGs的表達模式Expression profile of DEGs in Profile 1. CK: 對照組(未經(jīng)低溫脅迫的3日齡封蓋子)Control group (3 d-old post-capped brood unexposed to low temperature); T18: 20℃低溫處理18 h的3日齡封蓋子 3 d-old post-capped brood exposed to low temperature (20℃) for 18 h； T36: 20℃低溫處理36 h的3日齡封蓋子 3 d-old post-capped brood exposed to low temperature (20℃) for 36 h. 下圖同The same below.

2.3 差異表達基因的GO分類

GO富集分析結果顯示，3日齡意蜂封蓋子在受到低溫脅迫后DEGs歸類到生物學進程的基因數(shù)最多。Profile 6趨勢上的DEGs分別富集在30個GO條目上，Profile 7趨勢上的DEGs富集在28個GO條目上。這兩種趨勢從生物學進程分析，主要富集在細胞進程、代謝進程、單有機體進程；從分子功能分析，主要富集在結合、催化活性、分子傳感器活性；從細胞組件分析，主要富集在細胞膜組分、細胞膜、細胞。除此之外Profile 6趨勢中的DEGs還富集在核酸結合轉(zhuǎn)錄因子活性、生物學進程負調(diào)控、發(fā)育進程、突觸和突觸組分上，而Profile 7趨勢中還富集在行為(1個)、轉(zhuǎn)錄因子活性、蛋白質(zhì)結合(1個)、細胞連接等條目上(圖2)。

Profile 1為隨低溫脅迫時間延長基因表達量先下調(diào)表達，之后保持不變，該趨勢上富集183個DEGs，這些基因富集在25個GO條目上，排在前10的GO條目分別為代謝進程(28個)、單有機體進程(26個)、結合(22個)、細胞進程(21個)、催化活性(18個)、細胞膜組分(12個)、細胞膜(12個)、細胞(12個)、細胞組分(12個)、應激(8個)條目上(圖2)。說明低溫脅迫主要影響3日齡封蓋子的代謝進程和細胞結構等。同時機體核酸結合轉(zhuǎn)錄因子活性、發(fā)育相關的基因不會隨著低溫脅迫時間延長表達量改變，與行為相關的基因表達量隨著低溫脅迫時間延長而升高。

圖2 低溫脅迫不同時間3日齡意大利蜜蜂封蓋子差異表達基因GO富集分析Fig. 2 GO enrichment analysis of differentially expressed genes in the 3 d-old post-capped brood of Apis mellifera ligustica subjected to low temperature stress for different time

2.4 差異表達基因的KEGG pathway富集分析

通過KEGG pathway富集分析結果顯示，Profile 6趨勢上的DEGs富集在70條信號通路上，其中有48條與代謝相關的信號通路?；驍?shù)富集最多的為代謝途徑(22個)、神經(jīng)活性受體(11個)、次生代謝物的生物合成(11個)、不同環(huán)境中的微生物代謝(8個)、抗生素的生物合成(6個)、細胞外基質(zhì)-受體相互作用(5個)、氨基酸的生物合成(5個)、拼接體(5個)、胞吞作用(5個)等；找到基因ILP(insulin-like peptide receptor gene, ID411297)和基因FoxO(forkhead box protein O gene, ID 727091)同時富集在FoxO信號通路和壽命調(diào)節(jié)信號通路上(表2)。

表2 Profile 6中DEGs基因數(shù)富集前12的KEGG通路Table 2 Top 12 KEGG pathways enriched by DEGs of Profile 6

Profile 7趨勢中的271個DEGs富集在51條信號通路上，其中與代謝相關的信號通路有31條。基因數(shù)富集最多的為代謝進程(16個)、次生代謝物的生物合成(7個)、嘌呤代謝(4個)、賴氨酸降解(3個)、Jak-STAT信號通路(2個)、Notch信號通路(2個)、淀粉和蔗糖代謝(2個)、TGF-β信號通路(2個)、基礎轉(zhuǎn)錄因子(2個)、Hippo信號通路(2個)；其中基因CBP(CREB-binding protein gene, ID726332)同時富集在Jak-STAT, Notch, TGF-β, FoxO和Wnt 5條信號通路上(表3)。

表3 Profile 7中DEGs基因數(shù)富集前12的KEGG通路Table 3 Top 12 KEGG pathways enriched by DEGs in Profile 7

Profile 1模式中183個DEGs富集在43條信號通路上。DEGs主要富集在代謝進程模塊中，包括新陳代謝總覽中的代謝途徑(13個)、不同環(huán)境中微生物代謝(3個)、次生代謝物的生物合成(3個)；核苷酸代謝中有嘌呤代謝(3個)和嘧啶代謝(3個)；氨基酸代謝中有酪氨酸代謝(2個)、丙氨酸、天冬氨酸和谷氨酸代謝(2個)；有1個基因CYP450306a1 (phm, cytochrome P450 306a1 gene, ID08398)富集在昆蟲激素的生物合成信號通路上，還發(fā)現(xiàn)基因WNT1(protein Wnt-1 gene, ID為413502)同時富集在Hippo, Wnt和hedgehog(Hh)3條信號通路上(表4)。

表4 Profile 1中DEGs基因數(shù)富集前13的KEGG通路Table 4 Top 13 KEGG pathways enriched by DEGs in Profile 1

2.5 差異表達基因的RT-qPCR驗證

利用RT-qPCR技術檢測隨機挑選的5個DEGs (kibra,FoxO,CBP,Myo20和arm)的表達模式。結果顯示這些DEGs的表達水平與RNA-Seq測序數(shù)據(jù)的基因表達水平變化趨勢一致(圖3)，說明本次轉(zhuǎn)錄組的測序數(shù)據(jù)真實可靠。

圖3 低溫脅迫不同時間3日齡意大利蜜蜂封蓋子差異表達基因RNA-Seq測序數(shù)據(jù)的RT-qPCR驗證Fig. 3 Verification of RNA-Seq sequencing data of differentially expressed genes in the 3 d-old post-capped brood of Apis mellifera ligustica subjected to low temperature stress for different time by RT-qPCRA: T18 vs CK; B: T36 vs CK. 圖中數(shù)據(jù)為平均值±標準差，柱上雙星號代表同一基因處理組與對照組之間相對表達量差異極顯著(P≤0.01, 非配對T檢驗)。Data in theFigure are mean±SD, and the double asterisk above bars represents extremely significant difference between the treatment group and the control group of the same gene in the relative expression level (P≤0.01, unpaired T-test).

3 討論

真社會性昆蟲蜜蜂蜂子具有典型的狹溫性發(fā)育特點，蜂群維持穩(wěn)定的合適巢溫對蜂子的健康發(fā)育具有十分重要意義。3日齡封蓋子是蜜蜂進行化蛹起始階段(預蛹期)，其低溫敏感性強。本研究采用趨勢分析對兩個梯度時間的低溫脅迫預蛹轉(zhuǎn)錄本表達模式進行聚類分析，得到了梯度低溫脅迫條件下的DEGs及歸類，篩選出具有生物學意義的基因集和蜜蜂響應低溫脅迫的關鍵基因。在此基礎上，我們重點討論了FoxO信號通路中的FoxO上調(diào)表達在低溫脅迫后蜜蜂預蛹無法蛻皮化蛹中的可能機制。

低溫脅迫可能導致蛻皮激素(20-hydroxyecdysone, 20E)的合成減少。20E是甾醇類化合物，是在昆蟲前胸腺(prothoracic glands, PG)中以膽固醇或植物固醇類物質(zhì)為底物，經(jīng)一系列酶，特別是碳22-羥化酶、碳2-羥化酶、碳25-羥化酶的催化作用合成的(Gilbertetetal., 2003; Gilbertet, 2004)。在對下調(diào)模式Profile 1中的DEGs進行KEGG富集分析顯示，在激素合成信號通路上存在一個下調(diào)基因CYP450306a1 (phm)，它通過調(diào)節(jié)碳25-羥化酶的活性來參與20E合成，phm表達量下調(diào)將會導致碳25-羥化酶活性降低，進而影響合成20E(Niwaetal., 2004; Warrenetal., 2004)。因此推測低溫會影響20E合成，導致其滴度降低。

低溫脅迫后，3日齡封蓋子無法化蛹，可能是由于ILP調(diào)控FoxO磷酸化留在細胞質(zhì)內(nèi)，無法行使其功能而導致無法蛻皮。昆蟲化蛹時ILP調(diào)控FoxO完成蛻皮的過程如下：在昆蟲中存在脊椎動物的胰島素類似物ILP，它的重要靶標蛋白為FoxO，在ILP上調(diào)時，F(xiàn)oxO發(fā)生磷酸化，無法進入細胞核而留在細胞質(zhì)內(nèi)，促進細胞生長，延遲化蛹(Panetal., 2018)；當ILP下調(diào)時，F(xiàn)oxO去磷酸化入核與20E和EcR-Usp聚合體激活后，誘導Br-C,E75A和E75B等初級及DHR3,DHR4和ftz-f1等次級應答基因表達，來激活羧肽酶-A(CPA)表達完成蛻皮發(fā)生(Gonzyetal., 2002; Zhou and Riddiford, 2002; Wu and Brown, 2006; Minakuchietal., 2008; Jinderetal., 2013)。在研究中上調(diào)模式Profile 6的差異表達基因的富集分析結果顯示，在蛻皮相關的FoxO信號通路(3個)、與壽命相關的壽命調(diào)節(jié)信號通路(2個)同時富集到了ILP和FoxO基因。ILP的表達量上調(diào)會導致FoxO磷酸化后無法入核留在細胞質(zhì)內(nèi)，無法被20E和EcR-Usp聚合體激活，導致20E的初級應答基因Br-C無法接受20E信號，蛻皮信號傳遞受阻。另外，20E具有拮抗胰島素的類似物ILP的作用，同時與ILP共同調(diào)控FoxO信號通路，從前一段論述低溫脅迫可能導致20E滴度降低，也正好與該結果相符。因此推測在3日齡封蓋子在受到低溫脅迫后，F(xiàn)oxO發(fā)生磷酸化留在細胞質(zhì)內(nèi)不能入核，是阻礙蜜蜂化蛹的一個關鍵基因。化蛹是一個極其復雜的過程，蛻皮只是其中的一個環(huán)節(jié)，除了蛻皮受到影響之外，低溫脅迫帶來的影響應該還有很多，仍然需要繼續(xù)關注。

我們前期研究結果顯示3日齡意蜂封蓋子在受到低溫脅迫后壽命顯著縮短(Wangetal., 2016)，F(xiàn)oxO在果蠅的壽命調(diào)節(jié)起到重要作用(Tataretal., 2004)。轉(zhuǎn)錄組數(shù)據(jù)顯示FoxO基因在受到低溫脅迫時表達量發(fā)生改變，這也可能是受到低溫脅迫后，蜜蜂壽命縮短的原因之一。

致謝福建農(nóng)林大學科研助手劉也齊參與樣本獲取工作，碩士研究生劉一名、馮睿蓉和李寒參與部分文獻研究工作，本科生馮海成、蔣心怡、赫昱暢、常斐然參與部分實驗，在此表示感謝！