張文馨 潘 霞 沈錫權(quán) 徐永健
(寧波大學(xué)海洋學(xué)院, 應(yīng)用海洋生物技術(shù)教育部重點實驗室, 寧波 315211)
海馬是一種深受東南亞國家喜愛的傳統(tǒng)中藥[1]。當(dāng)前, 因全球野生海馬種群的棲息地遭到破壞而導(dǎo)致其資源數(shù)量驟減, 但海馬的市場需求卻與日俱增。因此, 海馬養(yǎng)殖業(yè)極具潛力, 不僅可以滿足市場需求,還能夠保護其野生資源[2]。盡管多年來針對海馬的繁殖和培育開展了許多研究工作, 但仍然存在諸多養(yǎng)殖技術(shù)難題[3], 這些問題主要由養(yǎng)殖環(huán)境的脅迫所致。
鹽度是魚類生長發(fā)育的重要環(huán)境因素之一, 鹽度變化對海馬的生長發(fā)育及生理生化等方面存在很大影響[4]。已有一些研究涉及其他海洋魚類, 證明鹽度變化會影響魚類的存活和生長[5], 甚至影響魚類的性別決定和早期發(fā)育[6]。而鹽度脅迫會引起滲透脅迫, 導(dǎo)致魚類血漿應(yīng)激相關(guān)激素增多, 能量代謝受到影響, 電解質(zhì)平衡被破壞, 進而影響生物體內(nèi)的滲透壓穩(wěn)態(tài)和免疫穩(wěn)態(tài)[7]。Lu等[8]的研究發(fā)現(xiàn), 高鹽脅迫使黃鱔(Nibea albiflora)肝胰腺抗氧化酶和非特異性免疫酶活性發(fā)生顯著變化。有報道指出, 長期低鹽脅迫顯著降低了小龍蝦(Cherax quadricarinatus)免疫相關(guān)基因(Toll基因和ProPO基因)的表達水平[9]。這些工作對于海馬養(yǎng)殖條件的研究都有很好的借鑒。海馬屬于廣鹽性魚類, 可以存活于10‰—32‰鹽度的環(huán)境中[10], 但海馬幼體對鹽度變化的適應(yīng)能力差, 鹽度急劇變化會導(dǎo)致滲透脅迫,并增加滲透調(diào)節(jié)和其他幾種生物過程所需的能量需求[11,12], 因此鹽度波動對海馬幼體的影響較大。為減輕和避免鹽度變化帶來的不利影響, 本研究以養(yǎng)殖大海馬(Hippocampus kuda)幼體為研究對象,采用RNA-Seq技術(shù)對鹽度脅迫下大海馬幼體的肝臟組織進行了轉(zhuǎn)錄組測序, 分析其肝臟組織在高/低鹽度脅迫下的轉(zhuǎn)錄表達情況, 揭示了大海馬鹽度應(yīng)激的調(diào)控機制, 并從中挖掘響應(yīng)鹽度應(yīng)激的候選基因, 作為海馬養(yǎng)殖過程中的生物標(biāo)記指示物。
本實驗用生物大海馬幼體由寧波大學(xué)海馬養(yǎng)殖研究實驗室培育。挑選健康無病且活力良好的個體 [個體質(zhì)量: (0.46±0.07) g, 體長: (5.78±0.42) cm(P>0.05)], 于實驗前置于循環(huán)水養(yǎng)殖系統(tǒng)中暫養(yǎng)1周。暫養(yǎng)期間使用鹽度25‰、溫度(25±0.5)℃、含氧量DO>5 mg/L過經(jīng)砂濾后的天然海水, 每日定時投喂 2 次, 攝食 2h 后吸污并換水50%。實驗開始后, 各組鹽度設(shè)置如下: 對照組鹽度25‰(Control,CK)、高鹽脅迫組31‰(HS-test)和低鹽脅迫組17‰(LS-test), 各組水溫均保持25℃。急性鹽度脅迫 12h 后, 每組隨機取 3 尾大海馬的肝臟置于同一凍存管中制成一個混樣以消除個體誤差, 以此每組各設(shè)置3個生物學(xué)重復(fù)[13]。樣品組織以液氮速凍,置于-80℃ 冰箱中保存。
采用 RNAiso plus 試劑盒(大連寶生生物有限公司)提取肝臟總 RNA。采用Multiskan SkyHigh微孔板分光光度計(Thermo Fisher Scientific, 美國)檢驗RNA樣品的濃度(ng/μL)和純度(A260/A280值), 利用 1%的瓊脂糖凝膠電泳檢查 RNA 樣品的完整性,采用 Bioanalyzer 2100 系統(tǒng)(Agilent Technologies,美國)評估 RNA 的濃度和完整性。取每組2 μg RNA樣品, 采用 NEBNext UltraTMRNA 文庫制備試劑盒(New England Biolabs, 美國)構(gòu)建測序文庫并在 Illumina 平臺上測序, 測序策略為 PE150 [安諾優(yōu)達基因科技(北京)有限公司]。
使用 Perl 腳本處理原始數(shù)據(jù), 以確保后續(xù)用于信息分析的數(shù)據(jù)質(zhì)量, 包括去除接頭污染的 Reads、低質(zhì)量的 Reads 及含N 比例大于5%的Reads。對過濾后所得到的 Clean Data的質(zhì)量和數(shù)據(jù)量等性狀進行統(tǒng)計, 包括過濾后的總序列中質(zhì)量值大于30(錯誤率<0.1%) 的堿基數(shù)的比例(Clean Q30 bases rate,Q30)、數(shù)據(jù)量和堿基含量等。采用 Trinity 軟件對進行組裝, 對得到的全長轉(zhuǎn)錄本進行評估, 然后進行開放閱讀框(Open reading frame, ORF)預(yù)測和功能分析[14]。采用 RPKM法計算基因表達量[15], 采用DEGseq v1.18.0 軟件[16]對差異表達基因進行分析。對DEGs(Q≤0.05 且| log2 ratio |≥1)進行基因本體(Gene Ontology, GO)富集分析及京都基因與基因組百科全書(Kyoto Encyclopedia of Genes and Genomes, KEGG)數(shù)據(jù)庫通路富集分析。
本實驗以 β2-微球蛋白(beta-2 microglobulin,B2M)作為內(nèi)參基因[14], 對10 個差異表達基因的RNA-Seq 結(jié)果進行了驗證。根據(jù)轉(zhuǎn)錄組文庫中的Unigene 序列, 用 Primer Premier 5 軟件設(shè)計特異性引物(表 1)。使用Eppendorf Mastercycler ep realplex4(Eppendorf, 德國)和TB Green Premix ExTaqTMⅡ(Tli RNaseH Plus; TaKaRa, 大連)試劑盒對差異表達基因的表達量進行定量測定, 采用 2-ΔΔCt法[17]計算目的基因的相對表達水平。
表1 RT-qPCR所用基因及其引物序列Tab. 1 Primers used for RT-qPCR
高通量測序結(jié)果顯示(表 2), CK、HS-test、CK和 LS-test 組分別獲得55510236、58040146和55840746個原始讀段(Raw reads)。過濾后分別獲得54453888(CK)、56489512(HS-test)和54407308(LS-test)個Clean reads。Q30分別為97.80%(CK)、95.88%(HS-test)和95.84%(LS-test), 均大于95%, 證明數(shù)據(jù)可靠。Trinity 軟件組裝高質(zhì)量的測序數(shù)據(jù)后(表 3), 共生成154430個轉(zhuǎn)錄本, 最大長度為15866 bp, 最小長度為201 bp, 平均長度為1234.69 bp,N50 長度為2318 bp。所得轉(zhuǎn)錄本序列進一步組裝處理后獲得71794個單基因簇(Unigene), 平均長度為820.71 bp, N50為1780 bp。其中68.07%的Unigenes在200—600 bp, 9.96%的Unigenes在600—1000 bp, 1000—2000 bp占10.87%, 11.10%的Unigenes在2000 bp以上。
表2 大海馬測序數(shù)據(jù)的統(tǒng)計匯總Tab. 2 Statistic summary of the sequencing data of H. kuda
表3 組裝結(jié)果統(tǒng)計表Tab. 3 Statistics of assembly results
將獲得的 Unigene 數(shù)據(jù)信息分別在蛋白質(zhì)家族域數(shù)據(jù)庫(Pfam)、蛋白質(zhì)序列數(shù)據(jù)庫(Swiss-Prot)、非冗余蛋白數(shù)據(jù)(Nr)、真核生物蛋白相鄰類的聚簇(KOG)、基因本體論(GO)數(shù)據(jù)庫和東京基因與基因組百科全書(KEGG)等數(shù)據(jù)庫比對注釋。結(jié)果顯示(表 4), 37.61%的 Unigenes在Nr 數(shù)據(jù)庫中注釋, 占比最多;核酸序列數(shù)據(jù)庫(Nt)注釋的Unigenes占 31.19%; 4781 個Unigenes在KEGG數(shù)據(jù)庫中注釋, 僅占比6.66%; 24450個Unigenes 在GO數(shù)據(jù)庫中得到了注釋, 占比34.06%。GO數(shù)據(jù)庫涵蓋3個方面(圖 1): 生物過程(Biological process)、細胞組成(Cellular component)和分子功能(Molecular function)。利用Trinotate的注釋結(jié)果, 統(tǒng)計每個GO條目中注釋到的基因, 并根據(jù)二級GO條目得到統(tǒng)計結(jié)果。生物學(xué)過程中, 占比前5的主要類別依次為細胞過程(72.75%)、生物調(diào)節(jié)(50.84%)、代謝過程(51.46%)、生物過程調(diào)節(jié)(47.50%)和刺激反應(yīng)(36.05%); 在分子功能類別中代表性的類別分別是結(jié)合功能(63.10%)、催化活性(34.67%)和分子功能調(diào)節(jié)(7.62%); 細胞組成部分中最主要的類別為分別為細胞(79.42%)、細胞組分(79.26%)、細胞器(63.61%)和膜(46.63%)。
圖1 GO 統(tǒng)計柱狀圖Fig. 1 GO statistics histogram
表4 各個數(shù)據(jù)庫的注釋匯總表Tab. 4 Summary of comments in each database
在 KOG 數(shù)據(jù)庫中, 14428個基因注釋到了 25個直系同源組(圖 2), 其中信號轉(zhuǎn)導(dǎo)機制(T,23.61%)與一般功能預(yù)測(R, 13.54%)代表了最大的類群, 其次是轉(zhuǎn)錄(K, 11.73%)、翻譯后修飾/蛋白質(zhì)周轉(zhuǎn)/分子伴侶(O, 8.18%)和細胞內(nèi)運輸/分泌和囊泡運輸(U, 6.10%)。
圖2 KOG統(tǒng)計圖Fig. 2 KOG chart
由圖 3 可知,HS-testvs. CK 共獲得2740 個差異表達基因(P<0.05), 其中顯著下調(diào)495 個, 顯著上調(diào) 2245 個; LS-testvs. CK 共獲得 3715 個差異表達基因, 1854 個顯著上調(diào), 1861 個顯著下調(diào)。高溫脅迫組的上調(diào)基因數(shù)量低于下調(diào)基因; 而在低溫脅迫組中, 上下調(diào)基因數(shù)量占比相當(dāng)。
圖3 大海馬肝臟轉(zhuǎn)錄組中差異表達基因個數(shù)統(tǒng)計圖Fig. 3 Statistical diagram of differentially expressed genes in liver transcriptome of H. kuda
肝臟轉(zhuǎn)錄組中差異表達基因個數(shù)統(tǒng)計結(jié)果顯示(圖 4), 兩兩比對后三組均顯著差異表達的基因共計161個, HS-testvs. CK 與 LS-testvs. CK有1260個相同的差異表達基因, 而 1480 個差異基因僅在 HS-testvs. CK顯著表達, 共有2455 個差異基因僅在 LS-testvs. CK顯著表達。
圖4 大海馬肝臟轉(zhuǎn)錄組中差異基因韋恩圖Fig. 4 Venn diagram of differential gene in liver transcriptome of H. kuda
KEGG 是一個實用程序數(shù)據(jù)庫, 整合了系統(tǒng)功能信息、基因組信息和化學(xué)信息, 有助于了解生物系統(tǒng)的高級功能和實用性, 如在分子水平探究細胞、生物和生態(tài)系統(tǒng)的信息, 而差異表達基因的Pathway 注釋分析有助于在被整合的分子互動網(wǎng)絡(luò)中進一步解讀基因的功能[18], 是國際最常用的生物信息數(shù)據(jù)庫之一。將獲得的轉(zhuǎn)錄組差異表達基因在KEGG中進行富集比對時, 對 KEGG數(shù)據(jù)庫中的Pathway 信息應(yīng)用超幾何檢驗進行富集分析, 歸納差異表達基因中顯著富集的 Pathway(Q<0. 05)[14]。其中Q值為多重假設(shè)檢驗校正之后的P值,Q值越小, 表示差異表達基因在該通路中的富集顯著性越高[19]。在本研究中, HS-testvs. CK 的差異表達基因共富集到 227 條代謝調(diào)節(jié)通路, 其中顯著富集的通路是類固醇生物合成(Steroid biosynthesis)。LS-testvs. CK 的差異表達基因富集到252條代謝調(diào)節(jié)通路,其中顯著富集的通路有蛋白酶體(Proteasome)、氧化磷酸化(Oxidative phosphorylation)、類固醇生物合成(Steroid biosynthesis)及檸檬酸循環(huán)(Citrate cycle)等。
表 5和表 6 分別為HS-testvs.CK和LS-testvs.CK的差異表達基因富集得到的前 20 個 KEGG 代謝通路, 其中高鹽脅迫組中代謝途徑 (如類固醇生物合成、維生素消化吸收、甘氨酸、絲氨酸和蘇氨酸代謝等氨基酸代謝、甘油磷脂代謝、蛋白質(zhì)運輸、磷酸鹽和次磷酸鹽代謝、不飽和脂肪酸的生物合成、類固醇激素生物合成等)和免疫途徑(如PPAR信號通路、抗原加工提呈、NF-κB信號通路、HIF-1信號通路、IgA生成的腸道免疫網(wǎng)絡(luò)、細胞色素P450對異生素的代謝等)受到影響。低鹽脅迫組富集得到的代謝通路主要涉及氨基酸相關(guān)代謝有關(guān)通路 (如甘氨酸、絲氨酸、蘇氨酸、色氨酸、丙氨酸、天冬氨酸和谷氨酸代謝、β-丙氨酸代謝、賴氨酸降解、賴氨酸生物合成)、能量和脂肪酸代謝有關(guān)通路 [如類固醇生物合成、氧化磷酸化、檸檬酸循環(huán)、碳代謝、淀粉和蔗糖代謝、核糖體生物發(fā)生 (真核生物)]和免疫代謝相關(guān)通路(如蛋白酶體、抗原加工提呈、抗壞血酸和藻酸鹽代謝)等。
根據(jù)表 5和表 6富集通路中的有關(guān)基因, 結(jié)合所查閱的相關(guān)文獻, 最終篩選得到Gst、Hsp70、Hsp90、Sod、Bcl-2、Gadd45α、Tcrβ、Tap2和Traf3等免疫相關(guān)基因,Fadsd6、Fas、Sqle、Cyp51、Elovl6和Slc27a6等脂肪酸代謝相關(guān)基因,Vlcad、Pdha1、Mdh1、Idh3b、G6pd和Sdhd等能量代謝相關(guān)基因, 及Gldc、Atp6v1e1、Sms、Fadh、Asl、Ass1和Glud1等與氨基酸代謝以及滲透轉(zhuǎn)運有關(guān)的基因在鹽度應(yīng)激中發(fā)揮重要作用。
本研究共篩選出了Gst、Hsp70、Hsp90、Sod、Bcl-2、Gadd45α、Fadsd6、Fas、Sqle、Cyp51、Elovl6、Pdha1、Mdh1、Idh3b、G6pd、Sms、Asl和Glud1 等對鹽度脅迫敏感的候選基因。
為驗證RNA-Seq檢測結(jié)果的準(zhǔn)確性, 通過查閱文獻共選取 10個差異表達基因 (表 1)進行驗證。以B2M作為內(nèi)參基因, 采用 RT-qPCR 相對定量方法, 檢測這些差異表達基因在試驗組與處理組中的表達量差異。將 RT-qPCR 檢測結(jié)果與轉(zhuǎn)錄組分析結(jié)果進行比較, 結(jié)果表明, 10個差異表達基因的RNA-seq 相對表達水平與 RT-qPCR 相對表達趨勢基本一致 (圖 5和圖 6), 證明了轉(zhuǎn)錄組分析結(jié)果的可靠性。
圖5 高鹽脅迫組 (HS-test) RT-qPCR 及 RNA-Seq 的比較分析Fig. 5 Comparative analysis of RT-qPCR and RNA-Seq in HStest
圖6 低鹽脅迫組 (LS-test) RT-qPCR 及 RNA-Seq 的比較分析Fig. 6 Comparative analysis of RT-qPCR and RNA-Seq in LS-test
水環(huán)境中鹽度的變化能夠?qū)λa(chǎn)動物的行為和生理生化等方面產(chǎn)生十分復(fù)雜的影響, 從行為至分子水平的影響均十分嚴(yán)重, 且影響程度存在物種差異性[20]。鹽度變化一般會直接破壞水產(chǎn)動物的電解質(zhì)平衡, 導(dǎo)致滲透壓改變, 此外還會影響氨基酸代謝、能量代謝、免疫代謝和應(yīng)激激素水平等方面[21]。有研究表明, 魚類和甲殼類水產(chǎn)動物的鰓[22,23]、肌肉[24]和肝臟[25,26]等組織器官都會因鹽度變化而受到影響。本研究結(jié)果也表明, 高鹽和低鹽對大海馬幼體造成的影響是不同的 (表 5和表 6): 在暴露于低鹽時, 主要涉及與氨基酸代謝 (大部分基因上調(diào))、能量代謝 (大部分基因上調(diào))和免疫功能 (大部分基因上調(diào))相關(guān)的途徑發(fā)生改變; 而暴露于高鹽時, 免疫功能 (大部分基因上調(diào))和脂質(zhì)代謝 (大部分基因下調(diào))相關(guān)的途徑發(fā)生了明顯變化。
鹽度脅迫會造成水產(chǎn)動物體內(nèi)氧化損傷, 其中氧化應(yīng)激下的酶活性變化及其相關(guān)機制已有許多報道[22,27]。在本研究中的低鹽脅迫下, 蛋白酶體途徑顯著富集, 且富集于此途徑的基因均顯著上調(diào)(如26s蛋白酶體非ATP酶調(diào)節(jié)亞單位11A、B等)。蛋白酶體途徑可以響應(yīng)感染、熱休克或氧化損傷等多種細胞應(yīng)激, 此途徑具有重要且復(fù)雜的功能。熱休克蛋白 (HSP)可以識別錯誤折疊或未折疊的蛋白質(zhì), 通過蛋白酶體的水解作用降解蛋白[28]。本研究發(fā)現(xiàn), 低鹽脅迫會使大海馬幼體體內(nèi)產(chǎn)生氧化損傷, DEGs中抗氧化基因 (如Sod和Gst)和熱休克蛋白基因 (如Hsp70和Hsp90)的顯著上調(diào)證實了這一點。此外, 表 5和表 6的結(jié)果發(fā)現(xiàn)了多個免疫相關(guān)通路, 如抗原加工提呈、NF-κB信號通路和IgA生成的腸道免疫網(wǎng)絡(luò)等, 經(jīng)分析發(fā)現(xiàn)高/低鹽脅迫組中多數(shù)免疫相關(guān)基因顯著上調(diào), 如Tcrβ、Tap2、Traf3和Gadd45α等。在水體鹽度變化時, 大海馬免疫防御系統(tǒng)中的免疫分子轉(zhuǎn)錄增加, 推測是為環(huán)境應(yīng)激可能導(dǎo)致的疾病做準(zhǔn)備[29]。
表5 HS-test vs. CK差異基因富集的前20個KEGG代謝通路Tab. 5 Top 20 KEGG metabolic pathways enriched for HS-test vs. CK differential genes
游離氨基酸、離子通道和水通道蛋白是已知的3種鹽度脅迫效應(yīng)物[30]。據(jù)報道游離氨基酸在細胞體積變化和滲透調(diào)節(jié)過程中具有重要功能, 且大多數(shù)水生物種在環(huán)境鹽度變化時, 將體內(nèi)的游離氨基酸積累在細胞內(nèi)[31]。在本研究中, 表 6的前20個KEGG通路中低鹽脅迫組富集了多個氨基酸代謝通路, 如甘氨酸、絲氨酸、蘇氨酸、色氨酸、丙氨酸、天冬氨酸和谷氨酸代謝、賴氨酸降解和生物合成及β-丙氨酸代謝等代謝途徑。除表 6外, 本研究還富集到了半胱氨酸和蛋氨酸代謝、組氨酸代謝、精氨酸生物合成及精氨酸和脯氨酸代謝等通路。以上通路中富集的基因多數(shù)顯著上調(diào), 表明氨基酸代謝的增加可能是大海馬幼體應(yīng)對低鹽脅迫的一個重要調(diào)節(jié)途徑, 類似結(jié)果也在牡蠣 (Crassostrea gigas)[32]、尼羅羅非魚 (Oreochromis niloticus)[33]及凡納濱對蝦 (Litopenaeus vannamei)[26]等水產(chǎn)動物中被發(fā)現(xiàn)。但研究發(fā)現(xiàn)不同物種中重要的游離氨基酸種類并不一樣[32], 因此可以將氨基酸的定性及定量作為今后研究的重要內(nèi)容之一, 將氨基酸代謝相關(guān)的潛在標(biāo)志基因與氨基酸含量相聯(lián)系, 進一步了解鹽度脅迫下大海馬幼體滲透調(diào)解過程中氨基酸代謝的機制。
表6 LS-test vs. CK差異基因富集的前20個KEGG代謝通路Tab. 6 Top 20 KEGG metabolic pathways enriched for LS-test vs. CK differential genes
肝臟是魚類代謝的主要器官, 同時可以為滲透調(diào)節(jié)作用提供能量[34]。關(guān)于滲透調(diào)節(jié)過程中能量代謝的研究主要集中在魚類的生理生化反應(yīng)、氧氣消耗和氨排泄等方面, 而在分子水平上的研究較少。水產(chǎn)動物的能量代謝途徑在響應(yīng)鹽度應(yīng)激時起關(guān)鍵作用[35], 有研究發(fā)現(xiàn)在高鹽脅迫下斑點鱸(Lateolabrax maculatus)肝臟中能量代謝相關(guān)的許多基因顯著上調(diào)[25]。而本研究在低鹽脅迫下大部分能量代謝的關(guān)鍵基因 (如Vlcad、Pdha1、Mdh1、Idh3b和G6pd等)顯著上調(diào), 說明大海馬幼體肝臟在低鹽脅迫時的能量代謝過程可能承擔(dān)了重要的角色。推測大量能量用于合成關(guān)鍵蛋白、離子調(diào)節(jié)和滲透調(diào)節(jié)等過程, 以維持體內(nèi)的穩(wěn)態(tài)和滲透平衡[36]。
水產(chǎn)動物的脂肪酸代謝在鹽度脅迫下也承擔(dān)了多種功能, 一種是提供能量維持離子轉(zhuǎn)運并調(diào)節(jié)滲透壓, 另一種是調(diào)整細胞膜的結(jié)構(gòu)[26]。在本研究中高鹽脅迫下大海馬幼體肝臟中富集了多個脂肪酸合成和降解途徑 (如不飽和脂肪酸的生物合成、類固醇生物合成、類固醇激素生物合成、脂肪酸代謝、脂肪酸降解和甘油磷脂代謝等), 且通路中的大部分基因顯著下調(diào), 如Fadsd6、Cyp51、Sqle、Fabp和Slc27a6等脂肪代謝基因在高鹽脅迫下均顯著下調(diào), 由此推測大海馬幼體體內(nèi)可能有足夠的能量用于抵抗高鹽脅迫。