鯉在低溫脅迫下肝胰腺轉錄組測序分析

劉思嘉 田菲 張存芳 喬志剛 趙凱

（1. 中國科學院西北高原生物研究所，西寧 810008；2. 中國科學院大學，北京 100049；3. 河南師范大學水產學院，新鄉 453007）

溫度是影響魚類生理、行為及進化過程中最重要的環境因素。作為變溫動物，由于缺乏自身的體溫調控系統，魚類的體溫與環境溫度基本一致，機體內多種生理、生化過程直接受到環境溫度的調控，這使得魚類對環境溫度的依賴性極高。魚類的低溫適應是一個長期進化的過程，在與低溫環境的長期博弈中魚類也進化形成了多種適應機制以保證其能夠在晝夜、季節甚至常年低溫的環境下生存［1-3］。鯉（Cyprinus carpio）是一種分布于世界各地的古老經濟魚類，通過長期的人工馴化，已形成上百個品種，具有豐富的遺傳多態性［4］。由于其適中的體型、與模式生物斑馬魚（Danio rerio）近緣、廣溫適應等特點，鯉也是研究變溫動物低溫適應的理想實驗材料［5］。

轉錄組測序技術（RNA-seq technology）能夠在整體水平上較為全面地挖掘低溫響應基因集，并能夠動態監測其轉錄變化規律。Gracey等［6］利用由13 440個EST探針組成的芯片對鯉響應低溫脅迫時多組織的轉錄變化進行檢測，獲得3 400個低溫響應基因，其中多數基因的表達調控具有器官、組織特異性。Liang等［7］通過RNA-seq對具有耐低溫性狀的黑龍江野鯉多組織的轉錄活動進行檢測，發現大量低溫響應基因，其中參與蛋白定位及蛋白轉運等生物學過程的基因被顯著富集，參與生物節律調控、內質網蛋白加工、內吞作用、胰島素信號通路及溶酶體等有關通路的基因被顯著富集。

能量代謝調控是保證魚類在長期低溫環境下維持基本生理活動的基礎。肝臟是脊椎動物能量代謝過程最重要的功能器官，也是衡量魚類生理狀態時重點監測的指征器官［8］。由于魚類肝臟（鯉科魚類中則為肝胰腺）在糖原儲備上有一定的缺陷，肝臟中糖質新生等內源性生糖途徑在魚類維持機體葡萄糖平衡具有重要的作用［9］。然而，目前尚無針對低溫耐受差異鯉品種在低溫脅迫條件下肝胰腺代謝活動相關分子機制的研究。因此，本研究以低溫耐受品種松浦鏡鯉，低溫敏感品種荷包紅鯉為研究對象，通過RNA-seq技術縱向比較各品種肝胰腺中低溫響應基因集，同時橫向比較品種間低溫響應差異基因，并將二者相結合相對全面地進行分析，有利于鑒定鯉低溫耐受相關重要基因及其參與的功能通路，從而探索魚類低溫耐受的潛在分子調控機制。

1 材料與方法

1.1 材料

荷包紅鯉20尾（無直接親緣關系），松浦鏡鯉20尾（無直接親緣關系）均來自新鄉金龍生態養殖基地，為人工同期授精得到的6月齡幼魚。

1.2 方法

1.2.1 低溫試驗處理及樣品制備 將40尾參試魚從暫養池（水溫24℃）隨機分為兩組（10尾/品種）：常溫對照組（24℃）、低溫組（4℃）。對于低溫組的降溫過程為：以4℃ /h的降溫過程進行，直到降至實驗設定溫度。實驗期間兩組除溫度差異外，飼喂、供氧、水循環及管理一致。在低溫馴化后的15d對各組進行采樣，每組隨機取6尾健康個體進行采樣。取樣時用丁香酚將鯉麻醉后立即解剖取肝胰腺，液氮暫存。提取各樣本總RNA（提取試劑盒為：RNAprep Pure Tissue Kit，TIANGEN ），對提取的總RNA進行純度和質量檢測，后將每組的總RNA樣本進行等含量混合成3個生物學重復（每組中的2個RNA樣本等量混為1個樣本）。

1.2.2 cDNA文庫構建及Illumina測序 通過帶有OligodT（25）的磁珠富集、純化mRNA。以mRNA為模板合成cDNA。按照標準的Illumina建庫流程對每個測序樣本的cDNA構建小片段測序文庫。采用2×100 pair-end的模式在Illumina HiSeq2000測序平臺進行測序。

1.2.3 轉錄本重構及基因注釋 利用比對軟件STAR1（2.5.0）將每個樣本的clean reads比對到參考基因組上（基因組來源：http://www.carpbase.org/index.php），使用cuffmerge與原有的基因組注釋信息進行比較合并，除去小于200 bp的轉錄本，獲得最終的轉錄本。取每個基因中最長的轉錄本作為Unigene，基于Unigene進行基因的功能注釋。通過Trinotate軟件對Unigene行功能注釋。Trinotate軟件整合了多個數據庫，分別是蛋白數據庫NCBI-nr、SwissProt、KEGG和COG/KOG數據庫。比對得到跟Unigene序列相似性最高的蛋白，從而得到該Unigene的蛋白功能注釋信息。

1.2.4 表達豐度估計 對重構的轉錄本進行表達豐度估計。首先使用bowtie 2.0將每個樣品的reads分別比對到組裝出來的轉錄本上，然后使用RSEM計算比對到特定基因或轉錄本上的reads數目。為了便于不同基因之間進行表達量的比較，使用FPKM對表達量進行標準化。FPKM（Expected number of Fragments Per Kilobase of transcript sequence per Millions base pairs sequenced）是每百萬fragments中來自某一基因每千個堿基長度的fragments數目，其同時考慮了測序深度和基因長度對fragments計數的影響。

1.2.5 差異表達分析 通過Cuf fl ink軟件對兩兩實驗處理組之間的比較，鑒定差異表達的基因。基因差異表達的計算是基于基因表達水平分析中得到的read count數據。使用DESeq程序進行差異分析，篩選閾值為 FDR＜0.05且 |log2FoldChange|＞1，從而獲得差異表達矩陣。基于差異表達矩陣，對差異表達基因進行功能及代謝通路的富集分析（GO和KEGG富集分析）。

表1 qRT-PCR所用引物信息

表2 樣品測序數據質量評估統計表

1.2.6 差異表達基因的qRT-PCR驗證 從糖酵解/糖質新生代謝通路（Pathway ID：ko00010）中選取5個編碼酶蛋白的差異基因（cycg001314、cycg034145、cycg035940、cycg021976和 cycg041259）進行熒光定量PCR驗證分析，樣品取材與RNA-seq一致。根據相關文獻［10］選擇鯉β-actin基因為內參，使用Primer Primers 6.0設計擴增引物，引物信息見表1。反應體系如下：上游引物1 μL，下游引物1 μL，SYBR green 10 μL，cDNA 模板 1 μL 以及 ddH2O 7 μL。反應程序設置為：95℃ 30 s；95℃ 5 s，60℃ 30 s，40個循環。相對表達量按照2-ΔΔCT方法計算，利用R語言基礎軟件包中的LSD法對樣品中各基因表達量變化做統計學分析，以P＜0.05為具有顯著性差異。所有反應設置3次技術重復，所有樣本設置3次生物學重復。

2 結果

2.1 轉錄組測序與組裝

采用Illumina HiSeq2000高通量測序平臺對2個鯉品種在低溫脅迫15 d 后肝胰腺進行轉錄組測序，如表2所示，各測序樣品GC含量在46.03%-47.49%，較為一致，Q30在85%以上，說明RNA-seq數據量和質量都較高，為后續的數據拼接組裝提供了可靠的數據源。

對樣品reads在參考基因組的比對效率進行分析（表3），各樣本reads在參考基因組的比對率在64.55%-67.90%，且品種間無顯著差異；各樣本reads比對到注釋基因上的比對率在42.51%-46.27%，較為一致。對轉錄組測序所獲得的reads通過組裝和拼接，獲得轉錄本總長度達366100425bp（約349.14 Mb），平均長度為2 260 bp，N50為3 301 bp。對轉錄本提取Unigene后共獲得73 088個Unigene，其中包括17088個新基因。

表3 樣品比對結果統計表

2.2 Unigene功能注釋

將Unigene序列比對到蛋白數據庫NR、SwissProt、KEGG 和COG/KOG進行基因注釋，如圖1顯示。73 088個Unigene序列在NR數據庫中有57 237（占78.31%）個基因被成功注釋，在SwissProt數據庫中有38 985（53.34%）個基因得到注釋；在KEGG數據庫中有31 357（42.90%）個基因得到注釋；在COG/KOG 數據庫中有23 042（31.53%）個基因得到注釋；4個數據庫一共注釋了58 153（79.57%）個Unigene序列，未能夠得到注釋的Unigene序列有14 935（20.43%）個，有20 189（27.62%）個Unigene 序列序列被4大數據庫同時注釋。絕大多數基因可以通過NR 數據庫得到注釋，只注釋到KEGG、COG 數據庫的基因相對較少，分別為360 和995。多數基因可以同時至少在2個數據庫得到注釋。

2.3 差異表達分析

對重構的轉錄本進行表達定量，篩選FDR＜0.05且至少在兩組之間表達差異2倍以上（|log2FoldChange|＞1）的基因后，共獲得由10 521個基因組成差異表達矩陣。根據矩陣中基因的FPKM值進行PCA分析。在PC1-PC2二維坐標系中（圖2），相同品種、溫度處理的樣本能夠聚在一起，說明樣本重復性良好；而不同品種、溫度處理的樣本分布分散，表明低溫脅迫后能夠明顯改變鯉肝胰腺的整體轉錄活動，且低溫耐受鯉品種和低溫敏感鯉品種在響應低溫脅迫時肝胰腺的轉錄應答模式差異較大。

通過溫恩分析，對10 521個差異基因在品種間的分布進行統計。如圖3所示，與24℃ 對照組相比，荷包紅鯉中有2 367個特異基因在低溫下的表達量發生了明顯變化，其中911個基因表達上調；松浦鏡鯉中共有5 426個特異基因在低溫下的表達量發生了明顯變化，其中3 961個基因表達上調。另外，共有2 728個基因是在松浦鏡鯉和荷包紅鯉中在低溫脅迫后都發生表達變化的基因，其中，1 998個基因在松浦鏡鯉表達上調，1 303個基因在荷包紅鯉表達上調，1 017個基因在兩個鯉品種中均表現為低溫誘導上調。

圖1 Unigene功能注釋韋恩統計圖

圖2 差異表達矩陣PCA分析

圖3 差異表達基因韋恩統計圖

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

基于差異表達矩陣，對品種內低溫脅迫差異表達基因（與各品種24℃對照組比較）及品種間差異表達基因進行功能及代謝通路的富集分析（GO和KEGG富集分析）。

GO分類可分為生物學過程、細胞組分和分子功能3大類條目。分別對荷包紅鯉、松浦鏡鯉各自的差異表達基因及其共有的差異基因進行GO富集分析。對3大類條目中顯著富集基因最多的5個類別進行統計分析，如圖4所示，在生物學過程條目中，基因多富集在轉錄、轉錄調控、小分子代謝過程、基于RNA聚合酶II型啟動子的轉錄調控等方面；在細胞組分條目中，基因多富集在胞核、胞質及胞包膜成分；在分子功能條目中，參與ATP結合、鋅離子結合、金屬離子結合的基因所占比例最高。在富集基因最多的GO分類條目上，松浦鏡鯉與和荷包紅鯉基本一致，但在松浦鏡鯉鯉中參與上述生物學過程的基因數目都顯著多于荷包紅鯉。

分別對荷包紅鯉、松浦鏡鯉的差異表達基因在KEGG數據庫中進行代謝通路富集分析。對富集基因最多的10個KEGG二級條目進行統計分析，如圖5所示，富集在多種信號傳導、翻譯、內分泌系統等代謝通路中的基因所占比例最大。在松浦鏡鯉中，蛋白質折疊、排序及降解過程、糖酵解/糖質新生過程被顯著富集。在糖酵解/糖質新生代謝途徑中（Pathway ID：ko00010），多種編碼關鍵酶或催化亞基基因的表達量發生了顯著的變化（圖6），其中包括乳酸脫氫酶、1，6-二磷酸果糖激酶、6-磷酸果糖激酶、葡萄糖-6-磷酸酶、己糖激酶、葡萄糖激酶、葡萄糖磷酸變位酶、葡萄糖-6-磷酸異構酶、丙酮酸脫氫酶等多種關鍵酶。

進一步對表達上調/下調基因進行GO和KEGG富集分析（表 4），發現在富集基因最多的GO和KEGG條目中，松浦鏡鯉低溫誘導表達上調的基因數目顯著多于荷包紅鯉。其中，參與轉錄調控等生物學過程的基因在松浦鏡鯉中的上調/下調比顯著高于荷包紅鯉，富集在DNA模板介導轉錄過程的基因，在荷包紅鯉表現為明顯的低溫表達下調的趨勢，而在松浦鏡鯉表現為低溫表達上調；參與細胞組分的基因在兩個鯉品種中的上調/下調比值較為接近，且近似于1∶1；而在分子功能條目中，荷包紅鯉多數基因的上調/下調比接近于1∶1，而松浦鏡鯉中的表達上調基因數顯著多于表達下調基因數；在細胞凋亡通路中，荷包紅鯉比松浦鏡鯉有更多的基因發生低溫表達上調和下調。

圖4 GO富集分析

圖5 KEGG富集分析

表4 表達上調/下調基因的富集分析

2.5 qRT-PCR驗證

從糖酵解/糖質新生代謝通路（Pathway ID：ko00010）中選取5個編碼葡萄糖-6-磷酸酶、1，6-二磷酸果糖激酶、磷酸烯醇式丙酮酸激酶、檸檬酸合酶及乳酸脫氫酶等關鍵酶蛋白的差異基因：g6pc、fbp、pck、cisy和ldhb2（ 基 因 ID 分 別 是：cycg001314、cycg034145、cycg035940、cycg021976 和cycg041259），運用qRT-PCR的方法，對低溫脅迫下荷包紅鯉及松浦鏡鯉幼魚肝胰腺中的表達模式進行分析。通過qRT-PCR驗證（圖7-A-E），低溫脅迫后5個基因在松浦鏡鯉中的表達量都顯著的高于荷包紅鯉，其表達趨勢與RNA-seq定量結果（圖7-F）一致。對兩種方法下基因表達變化倍數值的相關性分析，兩種方法所得結果的相關系數達到0.89，結果高度一致，進一步驗證了轉錄組測序數據的可靠性。

3 討論

魚類作為變溫脊椎動物，其生長、繁殖等重要的生命活動更易受到環境溫度影響。長期低溫脅迫對魚類的攝食活動、應激反應、營養代謝等生理、生化活動會產生顯著地影響［11］，同種魚類由于遺傳背景和地理分布的差異，也會表現出迥然不同的低溫適應能力［12］。荷包紅鯉的育成地我國南方低緯度溫暖地區，常年自然水溫都保持在10℃以上，而松浦鏡鯉主產區在我國東北高緯度的寒冷地區，冰封期長達6個月之久，水溫穩定維持在0-4℃。本研究通過RNA-seq技術對上述2個低溫耐受能力不同的鯉品種在低溫脅迫條件下肝胰腺的轉錄活動進行檢測，發現在轉錄模式上二者具有顯著的差異。通過差異表達分析，松浦鏡鯉在低溫脅迫條件下差異表達基因數目顯著多于荷包紅鯉，這個可能是松浦鏡鯉在長期低溫適應過程中進化而來的一種復雜網絡調控機制，但也有可能是由于兩個參試品種與參考基因組來源品種［4］親緣關系不同，reads比對效率差異導致的結果（松浦鏡鯉與參考基因組品種具有較近的親緣關系）。

圖6 糖酵解/糖質新生代謝通路富集結果

圖7 qRT-PCR定量結果及相關性驗證

為了排除reads比對效率對結果的影響，對兩種鯉RNA-seq reads在參考基因組上的比對率進行分析（表3），發現兩個鯉品種間的基因組比對率及基因比對率并無顯著的差異，故而可以排除品種間序列比對效率差異對差異表達基因挖掘造成的影響。因此推測松浦鏡鯉在低溫脅迫過程中調動更多的基因參與低溫應答可能是保證其在長期低溫環境下生存的重要分子調控過程。值得注意的是，Xin等［13］在研究葡萄低溫適應的分子調控機制時同樣發現在低溫脅迫下，耐低溫葡萄品種比非耐低溫葡萄品種轉錄譜的表達變化更為廣泛。

相似地結果表明，在長期的進化過程中，低溫適應的物種會調動更多的功能基因參與抵抗低溫帶來的不良影響，而這種特殊的轉錄變化在不同的物種間都較為普遍。耐低溫鯉品種不但比低溫敏感鯉品種具有更為廣泛的低溫轉錄譜變化，其差異表達基因的表達上調/下調模式與后者也有顯著地差異。例如，多數富集在DNA模板介導轉錄過程的基因，在耐低溫的松浦鏡鯉中表現為低溫表達上調，而在低溫敏感的荷包紅鯉中則是表達下調的趨勢。這說明在低溫脅迫下，低溫耐受品種表現出更為強烈的正向轉錄調控。

更多參與轉錄正調控、信號轉導、能量代謝、內分泌等生理活動的基因表達激活、上調都有助于耐低溫對低溫環境的適應；而低溫敏感品種中大量參與細胞凋亡、轉錄負調控等生理活動的基因表達上調則可能是低溫對其細胞造成損害導致的。除此之外，大量參與細胞細胞膜、細胞質膜及細胞基質等組分的功能基因在低溫脅迫下發生轉錄變化，而低溫敏感品種與低溫耐受品種的上調/下調比較為相似，而以往的研究也證實了低溫會改變魚類細胞膜的流動性和通透性［14］，以上說明低溫環境造成的細胞結構組分的改變、重組是普遍存在于魚類中的。

當環境溫度低于適宜的生理溫度范圍時，魚類的攝食行為會減弱，而食物轉化率也會受到抑制［15］。同時，魚類為抵抗低溫環境而發生的一系列行為、生理生化變化本身也是一個耗能的過程。當外源性能源的獲取受到抑制時，機體需要動員更多的內源性能源物質或潛在產能底物來提供能量應對低溫環境。前人的研究發現，魚類在響應低溫脅迫時，血糖會顯著的升高以應對不良脅迫造成的生理壓力［16-17］。由于魚類對直接產能物質（肝糖原、肌糖原等）的貯備能力有限［9］，糖質新生作用能夠利用機體內源產能物質保證機體重要生理活動的能源補給，在長期的低溫適應過程中需要持續葡萄糖供給，這對于維持鯉機體血糖穩定起關鍵作用。與Long等［18-19］在斑馬魚（Danio rerio）低溫適應的研究結果類似，本研究發現糖質新生過程的關鍵酶編碼基因g6pc、fbp和pck的表達量在鯉低溫適應過程中都顯著地升高了，而其在松浦鏡鯉的表達量都顯著高于荷包紅鯉，也進一步說明糖質新生途徑對鯉低溫適應的重要性。當魚類應對低溫等非生物脅迫時耗氧量會增加并且產生乳酸［20］，乳酸脫氫酶將乳酸催化生成丙酮酸，而后者是糖質新生作用的重要底物。利用乳酸作為產能物質既是一種經濟的能源利用策略也能防止機體的乳酸或丙酮酸中毒。Kuo等［16］發現具有耐低溫性狀的草魚（Ctenopharyngodon idella）在低溫脅迫下乳酸脫氫酶的活性顯著高于適溫性的虱目魚（Chanos chanos）。本研究發現在低溫脅迫后，編碼乳酸脫氫酶催化亞基的ldhb2在松浦鏡鯉中的表達量顯著高于荷包紅鯉，也說明乳酸脫氫酶在魚類的低溫適應中起關鍵作用，且其活性及mRNA轉錄水平可以作為評價魚類低溫耐受能力的生化指標。

4 結論

本研究采用RNA-Seq技術比較分析了低溫敏感鯉品種荷包紅鯉與低溫耐受鯉品種松浦鏡鯉在低溫脅迫15 d后的肝胰腺轉錄組數據，通過差異表達分析篩選出10 521個差異基因，其中5 426個基因僅在松浦鏡鯉中表達變化顯著。低溫脅迫會顯著改變鯉肝胰腺的轉錄活動，編碼糖酵解/糖質新生途徑的多個關鍵酶基因在低溫耐受品種的表達量顯著高于低溫敏感品種，推測這些基因編碼蛋白功能有助于其低溫適應。通過對低溫耐受、低溫敏感性狀極端鯉品種的比較轉錄分析，有利于鑒定鯉低溫耐受相關重要基因及其參與的功能通路，從而探索魚類低溫耐受的潛在分子調控機制。