潘 霞 徐永健 寧 燕 沈錫權

(寧波大學海洋學院,浙江 寧波 315211)

海馬(Hippocampus)是中國水產(chǎn)經(jīng)濟中重要的高價值海水魚類,在中藥、海洋水族以及紀念品市場有很大的貿(mào)易成交量[1-2]。中醫(yī)認為,海馬具有補氣活血、補腎壯陽的功效,是治療陽痿、不育等疾病的佳品[3-4]。大海馬(Hippocampus kudaBleeker)廣泛分布在熱帶印度洋-太平洋地區(qū),通常生活在較淺的近海棲息地,包括紅樹林、海草床和河口等[5]。海馬的生長適溫范圍為12～33℃,最適溫度為19～28℃,但其對溫度急性變化的適應能力差,常由于水溫的急劇變化而患病甚至死亡[6]。水溫是影響魚類生長發(fā)育的重要環(huán)境因素之一,水溫驟變影響魚類的攝食、代謝、生長,還可能造成分解代謝加強,合成代謝減弱[7]。為了減輕和避免海馬養(yǎng)殖過程中溫度驟變所產(chǎn)生的有害影響,有必要研究高/低溫影響大海馬的分子生理機制。

分析溫度驟變前后的轉(zhuǎn)錄組學,可以較好地反映大海馬對溫度的應激調(diào)控。轉(zhuǎn)錄組測序(RNA-Seq)結合了轉(zhuǎn)錄組測序的實驗方法和數(shù)字基因表達譜的信息分析方法,可用于研究特定生物過程中基因的差異表達,具有通量高、重復性好、檢測范圍廣等優(yōu)點[8],已被證明是許多非模型魚對各種試驗條件轉(zhuǎn)錄反應的可靠評估方法[9],其中在硬骨魚研究中主要涉及發(fā)育、適應性進化、免疫應答和應激反應等[10]。此外,RNASeq 也被用于研究硬骨魚溫度適應和馴化所涉及的分子機制,如鯰魚[11]、虹鱒魚[12]和斑馬魚[13]。當前,對海馬的研究主要集中在溫度、鹽度、疾病等生理生化層面,鮮有研究采用RNA-Seq 技術分析溫度脅迫下大海馬的轉(zhuǎn)錄表達變化。本研究以大海馬幼體為研究對象,采用RNA-Seq 技術對溫度驟變前后大海馬肝臟組織進行轉(zhuǎn)錄組測序,分析不同試驗組之間基因轉(zhuǎn)錄表達的差異,以期挖掘溫度應激調(diào)節(jié)的候選標志基因,為闡明其分子機制提供理論基礎。

1 材料與方法

1.1 試驗材料

大海馬(H.kuda)幼體由寧波大學海馬養(yǎng)殖研究實驗室培育。挑選健康、活力良好的幼體,體重為0.301 2±0.027 6 g,體長為5.35±0.18 cm,兩者均無顯著性差異。在循環(huán)水養(yǎng)殖系統(tǒng)中暫養(yǎng)1 周,以避免其他環(huán)境因素的干擾。暫養(yǎng)溫度為25±0.5℃,鹽度25‰,充氧DO＞5 mg·L-1,試驗用海水為經(jīng)砂濾后的天然海水,每天定時投喂2 次,攝食2 h 后進行吸污。試驗水溫設置：對照組(CK)25℃,高溫脅迫組(H-test)31℃,低溫脅迫組(L-test)19℃,鹽度均為25‰。溫度脅迫12 h 后,每組隨機取5 尾大海馬的肝臟置于凍存管,液氮速凍,于-80℃冰箱中保存?zhèn)溆谩?/p>

1.2 RNA 提取、文庫構建和測序

參考Liu 等[11]的方法,通過混樣以消除個體之間的差異。采用RNAiso plus 試劑盒(大連寶生物有限公司)從肝臟樣品中提取總RNA。采用K5500 分光光度計(北京凱奧科技發(fā)展有限公司)測量260、280 nm波長處的吸光度值,并用1%瓊脂糖電泳檢查RNA 樣品的純度和完整性。采用Bioanalyzer 2100 系統(tǒng)(Agilent Technologies,美國)評估RNA 的濃度和完整性。每個樣品各取2 μg RNA,采用NEBNext? UltraTMRNA 文庫制備試劑盒(New England Biolabs,美國)按照制造商的建議構建測序文庫。在Illumina 平臺上測序,測序策略為PE150[安諾優(yōu)達基因科技(北京)有限公司]。

1.3 轉(zhuǎn)錄組的組裝、注釋和差異表達基因分析

使用Perl 腳本處理原始數(shù)據(jù)以確保用于信息分析的數(shù)據(jù)質(zhì)量,包括去除接頭污染的Reads、低質(zhì)量的Reads 以及含N 比例大于5%的Reads。過濾后得到的Clean Data,統(tǒng)計其質(zhì)量、數(shù)據(jù)量等性狀,包括過濾后的總序列中質(zhì)量值大于30(錯誤率小于0.1%)的堿基數(shù)的比例(clean Q30 bases rate,Q30)、數(shù)據(jù)量、堿基含量等。采用Trinity 軟件進行組裝,對得到的全長轉(zhuǎn)錄本進行評估,然后進行開放閱讀框(open reading frame,ORF)預測和功能分析。采用RPKM 方法計算基因表達量[14],DEGseq v1.18.0 軟件[15]用于差異表達基因(differentially expressed genes,DEGs)分析。選取Q≤0.05 和｜log2 ratio ｜≥1 的基因為DEGs,對其進行GO基因本體(Gene Ontology)富集分析和京都基因與基因組百科全書KEGG(Kyoto Encyclopedia of Genes and Genomes)通路富集分析。

1.4 實時熒光定量PCR(RT－qPCR)分析

選取13 個DEGs 對RNA-Seq 結果進行驗證,根據(jù)轉(zhuǎn)錄組文庫中的Unigene 序列,用Primer Premier 5 軟件設計特異性引物(表1)。以β2-微球蛋白(beta-2microglobulin,B2M)為內(nèi)參基因,采用2-ΔΔCt法[16]計算目的基因的相對表達水平。采用TB Green Premix Ex TaqTMⅡ(Tli RNaseH Plus)(TaKaRa,大連)進行定量表達檢測,在Eppendorf realplex4中進行。每組共9 尾魚,3 尾大海馬肝臟為一個混樣,設3 個重復。

2 結果與分析

2.1 大海馬肝臟mRNA 測序結果及從頭組裝

高通量測序顯示(表2),CK、H-test 和L-test 組分別獲得49 556 386、49 404 042、48 043 100 個原始讀段(raw reads)。過濾后分別獲得47 880 348、47 987 534和46 520 778 個clean reads。Q30 值均在94%以上,說明該數(shù)據(jù)可靠。由表3可知,高質(zhì)量的測序數(shù)據(jù)用Trinity 軟件組裝后生成了86 712 個轉(zhuǎn)錄本,長度在200～2 000 bp 之間,平均長度為1 443.25 bp,N50 長度為2 485 bp。所得轉(zhuǎn)錄本序列進一步組裝處理后獲得42 264 個單基因簇(Unigene),平均長度為1 122.71 bp,N50 為2 223 bp。其中,在200～600 bp 范圍內(nèi)的Unigenes 占53.27%,600～1 000 bp 占12.41%,1 000～2 000 bp 占16.30%,2 000 bp 以上的Unigenes 為18.02%。

2.2 基因功能注釋

將獲得的Unigene 數(shù)據(jù)信息分別在蛋白質(zhì)家族域數(shù)據(jù)庫(Pfam)、蛋白質(zhì)序列數(shù)據(jù)庫(Swiss-Prot)、非冗余蛋白數(shù)據(jù)庫(Nr)、真核生物蛋白相鄰類的聚簇(KOG)、基因本體論數(shù)據(jù)庫(GO)和東京基因與基因組百科全書(KEGG)等數(shù)據(jù)庫比對注釋。由表4可知,Nr 數(shù)據(jù)庫中注釋的Unigenes 所占比例最多,達53.64%；核酸序列數(shù)據(jù)庫(Nt)注釋的Unigenes 占44.49%；19 506 個Unigenes 在GO 數(shù)據(jù)庫中得到了注釋,占46.15%；KEGG 數(shù)據(jù)庫注釋的Unigenes 最少,共12 046 個。GO 分為分子功能(molecular function)、生物過程(biological process)和細胞組成(cellular component)3 個部分。利用Trinotate 的注釋結果,統(tǒng)計每個GO 條目中注釋到的基因,并根據(jù)二級GO 條目得到統(tǒng)計結果(圖1)。生物學過程的主要類別是細胞過程,占74.74%,其次是生物調(diào)節(jié)和代謝過程,分別占56.7%和49.7%；在分子功能類別中最主要的類別是結合功能(72.88%)和催化活性(35.33%)；細胞組成部分中代表性的類別為細胞部分(87.25%)、細胞器(57.24%)和細胞器部分(45.36%)。在KOG 數(shù)據(jù)庫中15 454 個基因注釋到了25 個直系同源組(圖2),其中信號轉(zhuǎn)導機制(T,21.26%)與一般功能預測(R,21.06%)代表了最大的類群,其次是翻譯后修飾/蛋白質(zhì)周轉(zhuǎn)/分子伴侶(O,8.29%)。

表1 RT-qPCR 所用基因及其引物序列Table 1 Genes used in RT-qPCR and their primer sequences

表2 大海馬測序數(shù)據(jù)的統(tǒng)計匯總Table 2 Statistic summary of the sequencing data of H.kuda

由圖3可知,H-test vs CK 獲得14 009 個差異表達基因,其中顯著上調(diào)5 543 個,顯著下調(diào)8 466 個；Ltest vs CK 獲得20 030 個差異表達基因(14 016 個差異基因顯著上調(diào),6 014 個差異表達基因顯著下調(diào))。高溫脅迫組上調(diào)基因數(shù)量低于下調(diào)基因,而低溫脅迫組中上調(diào)基因數(shù)量占比更高。由圖4可知,兩兩比對后三組均顯著差異表達的基因有2 851 個,H-test vs CK 與L-test vs CK 之間有7 324 個差異表達基因相同,而6 685 個差異基因僅在H-test vs CK 顯著表達,12 706 個差異基因僅在L-test vs CK 顯著表達。

圖1 GO 統(tǒng)計柱狀圖Fig.1 GO statistics histogram

表3 組裝結果統(tǒng)計表Table 3 Assembly result statistics

表4 各個數(shù)據(jù)庫的注釋匯總表Table 4 Summary of comments in each database

2.3 差異基因的KEGG 富集分析

KEGG 是一個綜合數(shù)據(jù)庫,大致分為系統(tǒng)信息、基因組信息和化學信息三大類,有助于了解生物系統(tǒng)的高級功能和實用性,如在基因和分子水平了解細胞、生物和生態(tài)系統(tǒng)的信息,其對差異表達基因的Pathway注釋分析有助于進一步解讀基因的功能[17]。差異表達基因在KEGG 數(shù)據(jù)庫中富集,以Q＜0.05 為標準,對KEGG 中每個Pathway 應用超幾何檢驗進行富集分析,找出差異表達基因中顯著性富集的Pathway。其中Q 值為多重假設檢驗校正之后的P值,Q 值越小,表示差異表達基因在該通路中的富集顯著性越高[18]。本研究中,H-test vs CK 的差異表達基因富集到332 條通路,顯著富集的通路主要是蛋白質(zhì)消化吸收、類固醇生物合成以及細胞外基質(zhì)ECM-受體相互作用,L-test vs CK 的差異表達基因富集到332 條通路,顯著富集的通路為同源重組及DNA 復制。

圖2 KOG 功能分類統(tǒng)計圖Fig.2 KOG function classification chart

圖3 大海馬肝臟轉(zhuǎn)錄組中差異表達基因個數(shù)統(tǒng)計圖Fig.3 Statistical diagram of differentially expressed genes in liver transcriptome of H.kuda

圖4 大海馬肝臟轉(zhuǎn)錄組中差異基因韋恩圖Fig.4 Venn diagram of differential gene in liver transcriptome of H.kuda

表5、表6 為富集得到的前20 個KEGG 代謝通路,其中高溫脅迫組中免疫途徑(如ECM-受體相互作用、IgA 腸道免疫網(wǎng)絡、軍團菌病、PI3K-Akt 信號通路等)、代謝通路(如蛋白質(zhì)消化吸收、類固醇生物合成、脂肪消化吸收、甘油酯代謝等)及細胞凋亡(依賴于鐵的程序性細胞死亡)等受到影響。涉及的顯著表達基因主要有SOD、HSP70、HSP90、PIK3R、CYCS、CASP9、CASP3、P21、Akt、IL-10、TLR、CCR9、MHC2 等。低溫脅迫主要涉及DNA 損傷修復(如同源重組、DNA 復制、細胞周期、核苷酸切除修復、堿基切除修復等)及代謝(花生四烯酸代謝、谷胱甘肽代謝、糖酵解等)等過程。篩選發(fā)現(xiàn)SOD、HSP70、HSP90、P21、P53、BAX、HSP70、HSP90、StARD7、ApoA4、CYP51、Fadsd6、MSH、DDB2、XRCC2、RAD52、Ogg1、PMS2 等基因在低溫應激中發(fā)揮重要作用。

表5 H-test vs CK 差異基因富集的前20 個KEGG 代謝通路Table 5 Top 20 KEGG metabolic pathways enriched in H-test vs CK differential genes

2.4 實時熒光定量PCR 結果

為驗證基于RNA-Seq 方法檢測的結果是否準確,隨機選取13 個差異表達基因(表1)進行驗證。以B2M作為內(nèi)參基因,采用RT-qPCR 相對定量方法檢測這些基因在試驗組與處理組中的轉(zhuǎn)錄表達差異。將RT-qPCR 方法檢測結果與轉(zhuǎn)錄組分析結果進行比較,結果表明,13 個被測基因的RT-qPCR 相對表達與RNA-seq 相對表達水平趨勢基本一致(圖5、圖6),證明了轉(zhuǎn)錄組分析結果的可靠性。

表6 L-test vs CK 差異基因富集的前20 個KEGG 代謝通路Table 6 Top 20 KEGG metabolic pathways enriched in L-test vs CK differential genes

圖5 高溫脅迫組RT-qPCR 及RNA-Seq 的比較分析Fig.5 Comparative analysis of RT-qPCR and RNA-Seq in H-test

圖6 低溫脅迫組RT-qPCR 及RNA-Seq 的比較分析Fig.6 Comparative analysis of RT-qPCR and RNA-Seq in L-test

3 討論

夏季高溫或冬季寒潮[19]以及人工養(yǎng)殖過程中換水等水溫驟變對水產(chǎn)生物的生長、發(fā)育和攝食等方面均有較大影響[20]。本研究分析了高/低溫脅迫處理的幼體大海馬肝臟轉(zhuǎn)錄組變化,高/低溫脅迫組與對照組相比分別獲得差異表達基因14 009 個和20 030 個,表明溫度脅迫可能激活了多種細胞代謝,以應答溫度應激對大海馬幼體造成的傷害。高/低溫脅迫均造成了抗氧化以及細胞凋亡方面相關基因表達的顯著改變,如超氧化物歧化酶(SOD)、熱休克蛋白(HSP70、HSP90)、細胞色素C(CYCS)、半胱天冬酶- 9(CASP9)、半胱天冬酶-3(CASP3)、細胞周期蛋白依賴性激酶抑制劑(P21)、腫瘤蛋白P53、Bcl-2 相關的X蛋白(BAX)等。此外,高溫脅迫組中免疫途徑中磷脂酰肌醇3-激酶調(diào)節(jié)亞基(PIK3R)、P21、RAC-γ 絲氨酸/蘇氨酸-蛋白激酶(Akt)、白細胞介素10(IL-10)、Toll 樣受體(TLR)等基因表達發(fā)生顯著變化。而低溫脅迫組中脂肪酸合成和DNA 損傷修復的相關基因表達差異顯著,如StAR 相關脂質(zhì)轉(zhuǎn)移結構域蛋白7(StARD7)、載脂蛋白A-Ⅳ(ApoA4)、羊毛甾醇14α-脫甲基酶(CYP51)、Δ-6 脂肪?；ワ柡兔?Fadsd6)、DNA 錯配修復蛋白(MSH)、DNA 損傷結合蛋白DDB2、DNA 修復蛋白XRCC2、DNA 修復蛋白RAD52、8-氧代鳥嘌呤糖基化酶(Ogg1)、錯配修復核酸內(nèi)切酶PMS2 等。此外,本研究發(fā)現(xiàn)(圖3)低溫脅迫組上調(diào)基因數(shù)是高溫脅迫組的2.6 倍,說明低溫脅迫較高溫脅迫組誘導了更多基因表達,在斑馬魚幼體[21]和大黃魚幼體[19]中也發(fā)現(xiàn)了這樣的現(xiàn)象,可能與大海馬幼體階段的耐寒能力有關。

溫度脅迫會增加水生生物體內(nèi)的活性氧(reactive oxygen,ROS)含量,其積累到一定程度會激發(fā)體內(nèi)抗逆性和抗感染的重要生物分子HSP70、HSP90[22]和SOD等基因的表達量顯著上調(diào)[23]。本研究中,高/低溫脅迫處理均造成了HSP70、HSP90 和SOD等抗氧化相關基因顯著上調(diào),可能是為了消除體內(nèi)的ROS,而熱休克蛋白(heat shock proteins,HSPs)作為分子伴侶可防止有害蛋白質(zhì)聚集,促進蛋白質(zhì)正確折疊并幫助修復[24]。

細胞凋亡指為維持內(nèi)環(huán)境穩(wěn)定,由基因控制的細胞自主的有序死亡。高/低溫脅迫均會導致細胞凋亡相關基因表達差異顯著,但兩者表達差異基因不同。本研究發(fā)現(xiàn)高溫脅迫組中顯著上調(diào)的基因主要是P21、CYCS、CASP9 和CASP3。而低溫脅迫組差異表達的凋亡相關基因有P53、P21、BAX。細胞色素C 大量釋放會激活Caspase-9,進而激活Caspase-3 后誘導細胞凋亡[25]。P53 參與多種修復DNA 損傷途徑,如染色質(zhì)重塑、堿基切除修復等[26]。而P21 是位于P53 基因下游的細胞周期蛋白依賴性激酶抑制因子,主要作用是減少受損DNA 的復制和積累,抑制細胞凋亡[27]。此外,BAX表達Bcl-2 相關的X 蛋白,通過抗氧化機制影響細胞凋亡[28]。本研究結果表明,高溫脅迫激活了caspase 依賴的細胞凋亡途徑,而低溫脅迫則造成與抑制細胞凋亡的相關基因上調(diào),并激活大量DNA 損傷修復基因,以抵抗低溫傷害。

生物體受到溫度、pH 值等刺激或細胞毒性損傷時,PI3K-Akt 信號傳導途徑會調(diào)節(jié)細胞凋亡、蛋白質(zhì)合成和細胞周期等重要過程[29]。本研究結果表明,在高溫脅迫下,PI3K-Akt 信號途徑相關基因的表達受到顯著影響(表5),共有62 個基因上調(diào),81 個基因顯著下調(diào)。此途徑中的2 個關鍵基因PIK3R1 以及AKT3的表達均顯著下調(diào)。而低溫脅迫下,PI3K-Akt 信號傳導途徑總體未受到較大影響。Sun 等[30]在雜色鮑(Haliotis diversicolor)中也發(fā)現(xiàn)在熱應激或缺氧及病原體攻擊時,PI3K、AKT基因表達顯著下調(diào)。此外,此途徑中涉及免疫的差異基因有IL-10、TLR等顯著上調(diào),顯著下調(diào)的差異基因有CC 趨化因子受體9 型(CCR9)、主要組織相容性復合物(MHC2)等。IL-10可下調(diào)T 輔助細胞因子(Th1 和Th2)、趨化因子和MHC Ⅱ類抗原的表達,是免疫和炎癥的關鍵調(diào)節(jié)因子[31]。TLR 信號通路也在先天免疫防御機制中起關鍵作用[32]。以上結果表明,高溫脅迫相較低溫脅迫更可能導致大海馬體內(nèi)免疫系統(tǒng)紊亂,最終導致患病率上升。

本研究中,低溫脅迫組中與脂肪酸代謝等有關的基因,如Fadsd6、ApoA4、StARD7 和CYP51 等基因表達顯著上調(diào),而高溫脅迫組中僅ApoA4 顯著上調(diào)。魚的脂肪酸代謝對溫度變化很敏感[33],膜脂肪酸去飽和被認為是魚適應低溫的一種重要機制,對維持膜的流動性、酶活力和細胞的正常功能至關重要[34]。研究發(fā)現(xiàn),低水溫可以增加硬骨魚類Δ6-去飽和酶基因表達,而水溫升高導致幾種淡水魚的Δ6-去飽和酶活性降低[35]。這可能是因為低溫脅迫導致魚體細胞中脂肪酸組成發(fā)生改變,增加膜脂的不飽和脂肪酸比率,以維持細胞膜流動性。載脂蛋白可能用于細胞增殖和組織修復期間支持細胞生長材料的運輸。而StAR 蛋白是介導膽固醇生物合成中線粒體轉(zhuǎn)運的限速蛋白[36],CYP51 是CYP 超家族中進化最保守的成員,而且是甾醇生物合成的關鍵酶[37]。本研究中低溫脅迫導致類固醇及不飽和脂肪酸等相關基因表達上調(diào),而在高溫脅迫中相反,可能是因為脂肪酸對于大海馬幼體抵抗低溫應激有重要作用。

環(huán)境刺激會造成氧化和抗氧化系統(tǒng)之間的平衡被打破,導致DNA 受到損傷。細胞周期調(diào)節(jié)、蛋白質(zhì)伴侶蛋白、DNA 修復等過程有降低DNA 損傷的作用[38]。本研究中,低溫脅迫組富集了多條DNA 損傷修復途徑(表6),如同源重組、DNA 復制、細胞周期、嘧啶代謝、核苷酸切除修復、堿基切除修復、錯配修復、嘌呤代謝等。篩選得到的上調(diào)基因有MSH、DDB2、XRCC2、RAD52、Ogg1、PMS2 等。上述基因中,高溫脅迫中顯著上調(diào)的僅有PMS2,說明低溫脅迫可能導致魚體內(nèi)產(chǎn)生DNA 損傷,誘導修復基因顯著表達以抵抗低溫應激。

4 結論

本研究采用RNA-Seq 技術,在新一代Illumina 測序平臺,采用PE150 測序策略,構建了幼體大海馬肝臟轉(zhuǎn)錄組數(shù)據(jù)庫,獲得了大量幼體大海馬的轉(zhuǎn)錄組數(shù)據(jù)信息,并分別構建高溫和低溫應激后的基因表達數(shù)據(jù)庫,對3 個數(shù)據(jù)庫進行從頭拼接、組裝和分析,挖掘了溫度脅迫相關的基因。通過轉(zhuǎn)錄組研究手段了解幼體大海馬的溫度應激反應,并初步了解其所涉及的重要過程和途徑,為進一步研究溫度應激的確切機制提供了基礎,而且有助于預防極端溫度對大海馬造成損傷。