劉思嘉 田菲 張存芳 喬志剛 趙凱
(1. 中國科學院西北高原生物研究所,西寧 810008;2. 中國科學院大學,北京 100049;3. 河南師范大學水產(chǎn)學院,新鄉(xiāng) 453007)
溫度是影響魚類生理、行為及進化過程中最重要的環(huán)境因素。作為變溫動物,由于缺乏自身的體溫調(diào)控系統(tǒng),魚類的體溫與環(huán)境溫度基本一致,機體內(nèi)多種生理、生化過程直接受到環(huán)境溫度的調(diào)控,這使得魚類對環(huán)境溫度的依賴性極高。魚類的低溫適應(yīng)是一個長期進化的過程,在與低溫環(huán)境的長期博弈中魚類也進化形成了多種適應(yīng)機制以保證其能夠在晝夜、季節(jié)甚至常年低溫的環(huán)境下生存[1-3]。鯉(Cyprinus carpio)是一種分布于世界各地的古老經(jīng)濟魚類,通過長期的人工馴化,已形成上百個品種,具有豐富的遺傳多態(tài)性[4]。由于其適中的體型、與模式生物斑馬魚(Danio rerio)近緣、廣溫適應(yīng)等特點,鯉也是研究變溫動物低溫適應(yīng)的理想實驗材料[5]。
轉(zhuǎn)錄組測序技術(shù)(RNA-seq technology)能夠在整體水平上較為全面地挖掘低溫響應(yīng)基因集,并能夠動態(tài)監(jiān)測其轉(zhuǎn)錄變化規(guī)律。Gracey等[6]利用由13 440個EST探針組成的芯片對鯉響應(yīng)低溫脅迫時多組織的轉(zhuǎn)錄變化進行檢測,獲得3 400個低溫響應(yīng)基因,其中多數(shù)基因的表達調(diào)控具有器官、組織特異性。Liang等[7]通過RNA-seq對具有耐低溫性狀的黑龍江野鯉多組織的轉(zhuǎn)錄活動進行檢測,發(fā)現(xiàn)大量低溫響應(yīng)基因,其中參與蛋白定位及蛋白轉(zhuǎn)運等生物學過程的基因被顯著富集,參與生物節(jié)律調(diào)控、內(nèi)質(zhì)網(wǎng)蛋白加工、內(nèi)吞作用、胰島素信號通路及溶酶體等有關(guān)通路的基因被顯著富集。
能量代謝調(diào)控是保證魚類在長期低溫環(huán)境下維持基本生理活動的基礎(chǔ)。肝臟是脊椎動物能量代謝過程最重要的功能器官,也是衡量魚類生理狀態(tài)時重點監(jiān)測的指征器官[8]。由于魚類肝臟(鯉科魚類中則為肝胰腺)在糖原儲備上有一定的缺陷,肝臟中糖質(zhì)新生等內(nèi)源性生糖途徑在魚類維持機體葡萄糖平衡具有重要的作用[9]。然而,目前尚無針對低溫耐受差異鯉品種在低溫脅迫條件下肝胰腺代謝活動相關(guān)分子機制的研究。因此,本研究以低溫耐受品種松浦鏡鯉,低溫敏感品種荷包紅鯉為研究對象,通過RNA-seq技術(shù)縱向比較各品種肝胰腺中低溫響應(yīng)基因集,同時橫向比較品種間低溫響應(yīng)差異基因,并將二者相結(jié)合相對全面地進行分析,有利于鑒定鯉低溫耐受相關(guān)重要基因及其參與的功能通路,從而探索魚類低溫耐受的潛在分子調(diào)控機制。
荷包紅鯉20尾(無直接親緣關(guān)系),松浦鏡鯉20尾(無直接親緣關(guān)系)均來自新鄉(xiāng)金龍生態(tài)養(yǎng)殖基地,為人工同期授精得到的6月齡幼魚。
1.2.1 低溫試驗處理及樣品制備 將40尾參試魚從暫養(yǎng)池(水溫24℃)隨機分為兩組(10尾/品種):常溫對照組(24℃)、低溫組(4℃)。對于低溫組的降溫過程為:以4℃ /h的降溫過程進行,直到降至實驗設(shè)定溫度。實驗期間兩組除溫度差異外,飼喂、供氧、水循環(huán)及管理一致。在低溫馴化后的15d對各組進行采樣,每組隨機取6尾健康個體進行采樣。取樣時用丁香酚將鯉麻醉后立即解剖取肝胰腺,液氮暫存。提取各樣本總RNA(提取試劑盒為:RNAprep Pure Tissue Kit,TIANGEN ),對提取的總RNA進行純度和質(zhì)量檢測,后將每組的總RNA樣本進行等含量混合成3個生物學重復(fù)(每組中的2個RNA樣本等量混為1個樣本)。
1.2.2 cDNA文庫構(gòu)建及Illumina測序 通過帶有OligodT(25)的磁珠富集、純化mRNA。以mRNA為模板合成cDNA。按照標準的Illumina建庫流程對每個測序樣本的cDNA構(gòu)建小片段測序文庫。采用2×100 pair-end的模式在Illumina HiSeq2000測序平臺進行測序。
1.2.3 轉(zhuǎn)錄本重構(gòu)及基因注釋 利用比對軟件STAR1(2.5.0)將每個樣本的clean reads比對到參考基因組上(基因組來源:http://www.carpbase.org/index.php),使用cuffmerge與原有的基因組注釋信息進行比較合并,除去小于200 bp的轉(zhuǎn)錄本,獲得最終的轉(zhuǎn)錄本。取每個基因中最長的轉(zhuǎn)錄本作為Unigene,基于Unigene進行基因的功能注釋。通過Trinotate軟件對Unigene行功能注釋。Trinotate軟件整合了多個數(shù)據(jù)庫,分別是蛋白數(shù)據(jù)庫NCBI-nr、SwissProt、KEGG和COG/KOG數(shù)據(jù)庫。比對得到跟Unigene序列相似性最高的蛋白,從而得到該Unigene的蛋白功能注釋信息。
1.2.4 表達豐度估計 對重構(gòu)的轉(zhuǎn)錄本進行表達豐度估計。首先使用bowtie 2.0將每個樣品的reads分別比對到組裝出來的轉(zhuǎn)錄本上,然后使用RSEM計算比對到特定基因或轉(zhuǎn)錄本上的reads數(shù)目。為了便于不同基因之間進行表達量的比較,使用FPKM對表達量進行標準化。FPKM(Expected number of Fragments Per Kilobase of transcript sequence per Millions base pairs sequenced)是每百萬fragments中來自某一基因每千個堿基長度的fragments數(shù)目,其同時考慮了測序深度和基因長度對fragments計數(shù)的影響。
1.2.5 差異表達分析 通過Cuf fl ink軟件對兩兩實驗處理組之間的比較,鑒定差異表達的基因?;虿町惐磉_的計算是基于基因表達水平分析中得到的read count數(shù)據(jù)。使用DESeq程序進行差異分析,篩選閾值為 FDR<0.05且 |log2FoldChange|>1,從而獲得差異表達矩陣?;诓町惐磉_矩陣,對差異表達基因進行功能及代謝通路的富集分析(GO和KEGG富集分析)。
表1 qRT-PCR所用引物信息
表2 樣品測序數(shù)據(jù)質(zhì)量評估統(tǒng)計表
1.2.6 差異表達基因的qRT-PCR驗證 從糖酵解/糖質(zhì)新生代謝通路(Pathway ID:ko00010)中選取5個編碼酶蛋白的差異基因(cycg001314、cycg034145、cycg035940、cycg021976和 cycg041259)進行熒光定量PCR驗證分析,樣品取材與RNA-seq一致。根據(jù)相關(guān)文獻[10]選擇鯉β-actin基因為內(nèi)參,使用Primer Primers 6.0設(shè)計擴增引物,引物信息見表1。反應(yīng)體系如下:上游引物1 μL,下游引物1 μL,SYBR green 10 μL,cDNA 模板 1 μL 以及 ddH2O 7 μL。反應(yīng)程序設(shè)置為:95℃ 30 s;95℃ 5 s,60℃ 30 s,40個循環(huán)。相對表達量按照2-ΔΔCT方法計算,利用R語言基礎(chǔ)軟件包中的LSD法對樣品中各基因表達量變化做統(tǒng)計學分析,以P<0.05為具有顯著性差異。所有反應(yīng)設(shè)置3次技術(shù)重復(fù),所有樣本設(shè)置3次生物學重復(fù)。
采用Illumina HiSeq2000高通量測序平臺對2個鯉品種在低溫脅迫15 d 后肝胰腺進行轉(zhuǎn)錄組測序,如表2所示,各測序樣品GC含量在46.03%-47.49%,較為一致,Q30在85%以上,說明RNA-seq數(shù)據(jù)量和質(zhì)量都較高,為后續(xù)的數(shù)據(jù)拼接組裝提供了可靠的數(shù)據(jù)源。
對樣品reads在參考基因組的比對效率進行分析(表3),各樣本reads在參考基因組的比對率在64.55%-67.90%,且品種間無顯著差異;各樣本reads比對到注釋基因上的比對率在42.51%-46.27%,較為一致。對轉(zhuǎn)錄組測序所獲得的reads通過組裝和拼接,獲得轉(zhuǎn)錄本總長度達366100425bp(約349.14 Mb),平均長度為2 260 bp,N50為3 301 bp。對轉(zhuǎn)錄本提取Unigene后共獲得73 088個Unigene,其中包括17088個新基因。
表3 樣品比對結(jié)果統(tǒng)計表
將Unigene序列比對到蛋白數(shù)據(jù)庫NR、SwissProt、KEGG 和COG/KOG進行基因注釋,如圖1顯示。73 088個Unigene序列在NR數(shù)據(jù)庫中有57 237(占78.31%)個基因被成功注釋,在SwissProt數(shù)據(jù)庫中有38 985(53.34%)個基因得到注釋;在KEGG數(shù)據(jù)庫中有31 357(42.90%)個基因得到注釋;在COG/KOG 數(shù)據(jù)庫中有23 042(31.53%)個基因得到注釋;4個數(shù)據(jù)庫一共注釋了58 153(79.57%)個Unigene序列,未能夠得到注釋的Unigene序列有14 935(20.43%)個,有20 189(27.62%)個Unigene 序列序列被4大數(shù)據(jù)庫同時注釋。絕大多數(shù)基因可以通過NR 數(shù)據(jù)庫得到注釋,只注釋到KEGG、COG 數(shù)據(jù)庫的基因相對較少,分別為360 和995。多數(shù)基因可以同時至少在2個數(shù)據(jù)庫得到注釋。
對重構(gòu)的轉(zhuǎn)錄本進行表達定量,篩選FDR<0.05且至少在兩組之間表達差異2倍以上(|log2FoldChange|>1)的基因后,共獲得由10 521個基因組成差異表達矩陣。根據(jù)矩陣中基因的FPKM值進行PCA分析。在PC1-PC2二維坐標系中(圖2),相同品種、溫度處理的樣本能夠聚在一起,說明樣本重復(fù)性良好;而不同品種、溫度處理的樣本分布分散,表明低溫脅迫后能夠明顯改變鯉肝胰腺的整體轉(zhuǎn)錄活動,且低溫耐受鯉品種和低溫敏感鯉品種在響應(yīng)低溫脅迫時肝胰腺的轉(zhuǎn)錄應(yīng)答模式差異較大。
通過溫恩分析,對10 521個差異基因在品種間的分布進行統(tǒng)計。如圖3所示,與24℃ 對照組相比,荷包紅鯉中有2 367個特異基因在低溫下的表達量發(fā)生了明顯變化,其中911個基因表達上調(diào);松浦鏡鯉中共有5 426個特異基因在低溫下的表達量發(fā)生了明顯變化,其中3 961個基因表達上調(diào)。另外,共有2 728個基因是在松浦鏡鯉和荷包紅鯉中在低溫脅迫后都發(fā)生表達變化的基因,其中,1 998個基因在松浦鏡鯉表達上調(diào),1 303個基因在荷包紅鯉表達上調(diào),1 017個基因在兩個鯉品種中均表現(xiàn)為低溫誘導(dǎo)上調(diào)。
圖1 Unigene功能注釋韋恩統(tǒng)計圖
圖2 差異表達矩陣PCA分析
圖3 差異表達基因韋恩統(tǒng)計圖
基于差異表達矩陣,對品種內(nèi)低溫脅迫差異表達基因(與各品種24℃對照組比較)及品種間差異表達基因進行功能及代謝通路的富集分析(GO和KEGG富集分析)。
GO分類可分為生物學過程、細胞組分和分子功能3大類條目。分別對荷包紅鯉、松浦鏡鯉各自的差異表達基因及其共有的差異基因進行GO富集分析。對3大類條目中顯著富集基因最多的5個類別進行統(tǒng)計分析,如圖4所示,在生物學過程條目中,基因多富集在轉(zhuǎn)錄、轉(zhuǎn)錄調(diào)控、小分子代謝過程、基于RNA聚合酶II型啟動子的轉(zhuǎn)錄調(diào)控等方面;在細胞組分條目中,基因多富集在胞核、胞質(zhì)及胞包膜成分;在分子功能條目中,參與ATP結(jié)合、鋅離子結(jié)合、金屬離子結(jié)合的基因所占比例最高。在富集基因最多的GO分類條目上,松浦鏡鯉與和荷包紅鯉基本一致,但在松浦鏡鯉鯉中參與上述生物學過程的基因數(shù)目都顯著多于荷包紅鯉。
分別對荷包紅鯉、松浦鏡鯉的差異表達基因在KEGG數(shù)據(jù)庫中進行代謝通路富集分析。對富集基因最多的10個KEGG二級條目進行統(tǒng)計分析,如圖5所示,富集在多種信號傳導(dǎo)、翻譯、內(nèi)分泌系統(tǒng)等代謝通路中的基因所占比例最大。在松浦鏡鯉中,蛋白質(zhì)折疊、排序及降解過程、糖酵解/糖質(zhì)新生過程被顯著富集。在糖酵解/糖質(zhì)新生代謝途徑中(Pathway ID:ko00010),多種編碼關(guān)鍵酶或催化亞基基因的表達量發(fā)生了顯著的變化(圖6),其中包括乳酸脫氫酶、1,6-二磷酸果糖激酶、6-磷酸果糖激酶、葡萄糖-6-磷酸酶、己糖激酶、葡萄糖激酶、葡萄糖磷酸變位酶、葡萄糖-6-磷酸異構(gòu)酶、丙酮酸脫氫酶等多種關(guān)鍵酶。
進一步對表達上調(diào)/下調(diào)基因進行GO和KEGG富集分析(表 4),發(fā)現(xiàn)在富集基因最多的GO和KEGG條目中,松浦鏡鯉低溫誘導(dǎo)表達上調(diào)的基因數(shù)目顯著多于荷包紅鯉。其中,參與轉(zhuǎn)錄調(diào)控等生物學過程的基因在松浦鏡鯉中的上調(diào)/下調(diào)比顯著高于荷包紅鯉,富集在DNA模板介導(dǎo)轉(zhuǎn)錄過程的基因,在荷包紅鯉表現(xiàn)為明顯的低溫表達下調(diào)的趨勢,而在松浦鏡鯉表現(xiàn)為低溫表達上調(diào);參與細胞組分的基因在兩個鯉品種中的上調(diào)/下調(diào)比值較為接近,且近似于1∶1;而在分子功能條目中,荷包紅鯉多數(shù)基因的上調(diào)/下調(diào)比接近于1∶1,而松浦鏡鯉中的表達上調(diào)基因數(shù)顯著多于表達下調(diào)基因數(shù);在細胞凋亡通路中,荷包紅鯉比松浦鏡鯉有更多的基因發(fā)生低溫表達上調(diào)和下調(diào)。
圖4 GO富集分析
圖5 KEGG富集分析
表4 表達上調(diào)/下調(diào)基因的富集分析
從糖酵解/糖質(zhì)新生代謝通路(Pathway ID:ko00010)中選取5個編碼葡萄糖-6-磷酸酶、1,6-二磷酸果糖激酶、磷酸烯醇式丙酮酸激酶、檸檬酸合酶及乳酸脫氫酶等關(guān)鍵酶蛋白的差異基因:g6pc、fbp、pck、cisy和ldhb2( 基 因 ID 分 別 是:cycg001314、cycg034145、cycg035940、cycg021976 和cycg041259),運用qRT-PCR的方法,對低溫脅迫下荷包紅鯉及松浦鏡鯉幼魚肝胰腺中的表達模式進行分析。通過qRT-PCR驗證(圖7-A-E),低溫脅迫后5個基因在松浦鏡鯉中的表達量都顯著的高于荷包紅鯉,其表達趨勢與RNA-seq定量結(jié)果(圖7-F)一致。對兩種方法下基因表達變化倍數(shù)值的相關(guān)性分析,兩種方法所得結(jié)果的相關(guān)系數(shù)達到0.89,結(jié)果高度一致,進一步驗證了轉(zhuǎn)錄組測序數(shù)據(jù)的可靠性。
魚類作為變溫脊椎動物,其生長、繁殖等重要的生命活動更易受到環(huán)境溫度影響。長期低溫脅迫對魚類的攝食活動、應(yīng)激反應(yīng)、營養(yǎng)代謝等生理、生化活動會產(chǎn)生顯著地影響[11],同種魚類由于遺傳背景和地理分布的差異,也會表現(xiàn)出迥然不同的低溫適應(yīng)能力[12]。荷包紅鯉的育成地我國南方低緯度溫暖地區(qū),常年自然水溫都保持在10℃以上,而松浦鏡鯉主產(chǎn)區(qū)在我國東北高緯度的寒冷地區(qū),冰封期長達6個月之久,水溫穩(wěn)定維持在0-4℃。本研究通過RNA-seq技術(shù)對上述2個低溫耐受能力不同的鯉品種在低溫脅迫條件下肝胰腺的轉(zhuǎn)錄活動進行檢測,發(fā)現(xiàn)在轉(zhuǎn)錄模式上二者具有顯著的差異。通過差異表達分析,松浦鏡鯉在低溫脅迫條件下差異表達基因數(shù)目顯著多于荷包紅鯉,這個可能是松浦鏡鯉在長期低溫適應(yīng)過程中進化而來的一種復(fù)雜網(wǎng)絡(luò)調(diào)控機制,但也有可能是由于兩個參試品種與參考基因組來源品種[4]親緣關(guān)系不同,reads比對效率差異導(dǎo)致的結(jié)果(松浦鏡鯉與參考基因組品種具有較近的親緣關(guān)系)。
圖6 糖酵解/糖質(zhì)新生代謝通路富集結(jié)果
圖7 qRT-PCR定量結(jié)果及相關(guān)性驗證
為了排除reads比對效率對結(jié)果的影響,對兩種鯉RNA-seq reads在參考基因組上的比對率進行分析(表3),發(fā)現(xiàn)兩個鯉品種間的基因組比對率及基因比對率并無顯著的差異,故而可以排除品種間序列比對效率差異對差異表達基因挖掘造成的影響。因此推測松浦鏡鯉在低溫脅迫過程中調(diào)動更多的基因參與低溫應(yīng)答可能是保證其在長期低溫環(huán)境下生存的重要分子調(diào)控過程。值得注意的是,Xin等[13]在研究葡萄低溫適應(yīng)的分子調(diào)控機制時同樣發(fā)現(xiàn)在低溫脅迫下,耐低溫葡萄品種比非耐低溫葡萄品種轉(zhuǎn)錄譜的表達變化更為廣泛。
相似地結(jié)果表明,在長期的進化過程中,低溫適應(yīng)的物種會調(diào)動更多的功能基因參與抵抗低溫帶來的不良影響,而這種特殊的轉(zhuǎn)錄變化在不同的物種間都較為普遍。耐低溫鯉品種不但比低溫敏感鯉品種具有更為廣泛的低溫轉(zhuǎn)錄譜變化,其差異表達基因的表達上調(diào)/下調(diào)模式與后者也有顯著地差異。例如,多數(shù)富集在DNA模板介導(dǎo)轉(zhuǎn)錄過程的基因,在耐低溫的松浦鏡鯉中表現(xiàn)為低溫表達上調(diào),而在低溫敏感的荷包紅鯉中則是表達下調(diào)的趨勢。這說明在低溫脅迫下,低溫耐受品種表現(xiàn)出更為強烈的正向轉(zhuǎn)錄調(diào)控。
更多參與轉(zhuǎn)錄正調(diào)控、信號轉(zhuǎn)導(dǎo)、能量代謝、內(nèi)分泌等生理活動的基因表達激活、上調(diào)都有助于耐低溫對低溫環(huán)境的適應(yīng);而低溫敏感品種中大量參與細胞凋亡、轉(zhuǎn)錄負調(diào)控等生理活動的基因表達上調(diào)則可能是低溫對其細胞造成損害導(dǎo)致的。除此之外,大量參與細胞細胞膜、細胞質(zhì)膜及細胞基質(zhì)等組分的功能基因在低溫脅迫下發(fā)生轉(zhuǎn)錄變化,而低溫敏感品種與低溫耐受品種的上調(diào)/下調(diào)比較為相似,而以往的研究也證實了低溫會改變魚類細胞膜的流動性和通透性[14],以上說明低溫環(huán)境造成的細胞結(jié)構(gòu)組分的改變、重組是普遍存在于魚類中的。
當環(huán)境溫度低于適宜的生理溫度范圍時,魚類的攝食行為會減弱,而食物轉(zhuǎn)化率也會受到抑制[15]。同時,魚類為抵抗低溫環(huán)境而發(fā)生的一系列行為、生理生化變化本身也是一個耗能的過程。當外源性能源的獲取受到抑制時,機體需要動員更多的內(nèi)源性能源物質(zhì)或潛在產(chǎn)能底物來提供能量應(yīng)對低溫環(huán)境。前人的研究發(fā)現(xiàn),魚類在響應(yīng)低溫脅迫時,血糖會顯著的升高以應(yīng)對不良脅迫造成的生理壓力[16-17]。由于魚類對直接產(chǎn)能物質(zhì)(肝糖原、肌糖原等)的貯備能力有限[9],糖質(zhì)新生作用能夠利用機體內(nèi)源產(chǎn)能物質(zhì)保證機體重要生理活動的能源補給,在長期的低溫適應(yīng)過程中需要持續(xù)葡萄糖供給,這對于維持鯉機體血糖穩(wěn)定起關(guān)鍵作用。與Long等[18-19]在斑馬魚(Danio rerio)低溫適應(yīng)的研究結(jié)果類似,本研究發(fā)現(xiàn)糖質(zhì)新生過程的關(guān)鍵酶編碼基因g6pc、fbp和pck的表達量在鯉低溫適應(yīng)過程中都顯著地升高了,而其在松浦鏡鯉的表達量都顯著高于荷包紅鯉,也進一步說明糖質(zhì)新生途徑對鯉低溫適應(yīng)的重要性。當魚類應(yīng)對低溫等非生物脅迫時耗氧量會增加并且產(chǎn)生乳酸[20],乳酸脫氫酶將乳酸催化生成丙酮酸,而后者是糖質(zhì)新生作用的重要底物。利用乳酸作為產(chǎn)能物質(zhì)既是一種經(jīng)濟的能源利用策略也能防止機體的乳酸或丙酮酸中毒。Kuo等[16]發(fā)現(xiàn)具有耐低溫性狀的草魚(Ctenopharyngodon idella)在低溫脅迫下乳酸脫氫酶的活性顯著高于適溫性的虱目魚(Chanos chanos)。本研究發(fā)現(xiàn)在低溫脅迫后,編碼乳酸脫氫酶催化亞基的ldhb2在松浦鏡鯉中的表達量顯著高于荷包紅鯉,也說明乳酸脫氫酶在魚類的低溫適應(yīng)中起關(guān)鍵作用,且其活性及mRNA轉(zhuǎn)錄水平可以作為評價魚類低溫耐受能力的生化指標。
本研究采用RNA-Seq技術(shù)比較分析了低溫敏感鯉品種荷包紅鯉與低溫耐受鯉品種松浦鏡鯉在低溫脅迫15 d后的肝胰腺轉(zhuǎn)錄組數(shù)據(jù),通過差異表達分析篩選出10 521個差異基因,其中5 426個基因僅在松浦鏡鯉中表達變化顯著。低溫脅迫會顯著改變鯉肝胰腺的轉(zhuǎn)錄活動,編碼糖酵解/糖質(zhì)新生途徑的多個關(guān)鍵酶基因在低溫耐受品種的表達量顯著高于低溫敏感品種,推測這些基因編碼蛋白功能有助于其低溫適應(yīng)。通過對低溫耐受、低溫敏感性狀極端鯉品種的比較轉(zhuǎn)錄分析,有利于鑒定鯉低溫耐受相關(guān)重要基因及其參與的功能通路,從而探索魚類低溫耐受的潛在分子調(diào)控機制。