低溫脅迫下兩種不同形態群體擬穴青蟹轉錄組分析

陳 潔 董 夢 崔 悅 房翠蓮 陳旭堂 劉子明

(1. 麗水學院生態學院, 麗水 323000; 2. 麗水市生態經濟研究所, 麗水 323000)

擬穴青蟹(Scylla paramamosain)屬于甲殼綱、梭子蟹科、青蟹屬, 主要分布于我國東南沿海海域,是青蟹屬的優勢種[1—3]。“螯足腕節外緣中部發達刺的有無”是鑒定青蟹屬不同種的主要分類特征之一。對福建漳州擬穴青蟹進行形態學分析, 發現擬穴青蟹大多數個體(80%左右)一對螯足腕節外緣中部都不具有發達的刺, 少數個體具有發達的刺[4]。于是, 將“一對螯足腕節外緣中部都不具有發達的刺”和“至少一只螯足腕節外緣中部具有發達的刺”的擬穴青蟹分別定義為Sp1和Sp2群體。通過研究上述兩個形態群體擬穴青蟹鰓、肝胰腺和肌肉線粒體呼吸速率的季節變化, 發現Sp1和Sp2群體線粒體呼吸代謝的季節變化存在差異, 其中冬季差異最為顯著[5]。

至今, 有關溫度對擬穴青蟹的生物學效應已有不少的研究報道。Kong等[6]的研究結果表明, 冬季青蟹鰓超氧化物歧化酶(SOD)、過氧化氫酶(CAT)和谷胱甘肽過氧化物酶(GPX)活性顯著低于春、夏和秋季, 且低溫馴化會降低擬穴青蟹肌肉SOD、CAT和GPX活性[7], 由此說明在低溫條件下青蟹代謝水平降低, 活性氧生成減少進而引起抗氧化酶活性的下降。近年來, 熱休克蛋白(HSPs)成為研究生物受溫度脅迫的一大熱點[8—10]。目前擬穴青蟹hsp90[11]、hsp70[12]、hsp60[13]和hsp10[14]等熱休克蛋白基因已被克隆, 初步的研究結果表明擬穴青蟹hsp70、hsp60和hsp10表達對溫度變化十分靈敏。

上述這些研究主要從個體水平考察溫度對擬穴青蟹的影響, 在群體水平上探索溫度脅迫對擬穴青蟹生物學效應的報道很少。本研究以Sp1和Sp2兩種形態群體擬穴青蟹為研究對象, 通過轉錄組測序分析, 篩選低溫脅迫下差異表達基因, 將初步揭示擬穴青蟹不同形態群體對低溫脅迫的適應性差異。這項工作的開展將有利于擬穴青蟹耐寒優質品種的選育及其種質資源的保護。

1 材料與方法

1.1 擬穴青蟹樣品采集與處理

擬穴青蟹采樣于廣東省汕頭海域(緯度約23°22′N)。采集到的擬穴青蟹活體運回實驗室, 按照本課題組建立的方法[4]篩選出Sp1和Sp2個體。Sp1和Sp2個體的殼寬平均值分別為 (12.56±1.37)和(12.28±1.85) cm, 雌雄比例為1﹕1。選取大小相似、肢體完整、體表無菌斑的個體暫養于200 L的培養箱中, 內含150 L沙濾海水, 鹽度22, 溫度20℃, 暫養期間不投喂飼料, 暫養3d。擬穴青蟹Sp1和Sp2群體, 每個群體48只, 設置低溫脅迫組(8.0±0.5)℃和對照組(20.0±0.5)℃, 每組設置3個平行組(每組每個群體8只), 分別放于6個培養箱中。6h后取肝胰腺, 液氮速凍, 隨后存放于–80℃備用。

1.2 轉錄組測序

利用Trizol(生工, 中國)提取擬穴青蟹肝胰腺總RNA。1%瓊脂糖凝膠電泳和 NanoDrop 2000(Thermo,美國)檢測總RNA質量和濃度。將檢測合格的總RNA送杭州聯川生物技術股份有限公司測序。采用mRNASeq sample preparation kit(Illumina, 美國)構建cDNA文庫, 利用IlluminaHiseq4000高通量平臺(LC Sciences, 美國), 通過Paired-end 300 bp雙末端測序模式對12個樣本進行測序。

1.3 序列組裝和基因注釋

將原始序列中低質量、帶接頭和N比例大于10%的序列去除后, 利用Trinity 2.4.0軟件進行無參拼接, 將拼接得到的最長序列作為基因的參考序列(Unigene)[15]。使用BUSCO(v5.1.2)軟件對轉錄組進行完整性評估。隨后用DIAMOND軟件將Unigene分別與NCBI NR(http://www.ncbi.nlm.nih.gov/)、eggNOG(http://eggnogdb.embl.de/)、KEGG(http://www.genome.jp/kegg/)、Pfam(http://pfam.sanger.ac.uk/)、GO(http://www.geneontology.org/)和Swiss-prot(http://www.expasy.ch/sprot/)數據庫進行比對和注釋。

1.4 基因差異表達分析

基因表達量采用TPM(Transcripts Per Million)方法計算[16]。差異表達基因采用edgeR[17]軟件包篩選, 篩選閾值為log2(fold change)>1和FDR<0.05。差異表達基因的GO和KEGG富集分析分別利用GOseq R和KOBAS軟件包進行(FDR<0.05)[18,19]。

1.5 熒光定量PCR

選取上調和下調基因各5個。用TRIzol抽提樣本總RNA, 利用AMV反轉錄酶(TaKaRa, 日本)反轉錄獲得cDNA第一鏈。使用BioRad CFX96 Touch實時熒光定量PCR系統(BioRad, 美國)進行qPCR, 反應體系為: SYBR Premix ExTaq(2×)緩沖液(TaKaRa)12.5 μL, cDNA模板0.5 μL, 引物(10 μmol/L)各1 μL,以ddH2O補足至25 μ L。反應程序為: 94℃預變性5min; 95℃ 30s, 60℃ 30s, 72℃ 30s, 40個循環。qPCR檢測結果使用2–ΔΔCt方法分析[20]。以18S rRNA作為內參基因, 所用引物見表 1, 引物的擴增效率根據熒光定量PCR的MIQE指南, 使用 cDNA樣品的標準曲線確定, 引物的特異性根據熔解曲線分析確定。每組每個群體4個樣品, 每個樣品技術重復4次。

表1 引物序列Tab. 1 Primers used in this study

2 結果

2.1 測序和組裝結果

低溫脅迫組和對照組2個群體各取3個肝胰腺經過轉錄組測序、原始數據過濾后共計獲得66.00 Gb有效數據(表 2)。有效數據已上傳至NCBI Sequence Read Archive (SRA)數據庫(登錄號: PRJNA 743858)。Trinity組裝去冗余后獲得81853個Unigene。Unigene序列長度在201—17151 bp, 平均長度為420 bp,N50為1460 bp。BUSCO分析顯示,C值(Complete)為917占比90.5%, 其中S值(Complete and single-copy)544,D值(Complete and duplicated)373, 分別占比53.7%和36.8%, 說明拼接的轉錄組完整性較好。

表2 不同形態群體擬穴青蟹肝胰腺轉錄組數據Tab. 2 Hepatopancreatic transcriptome data of Scylla paramamosain with different morphological populations

2.2 轉錄組序列功能注釋

本次轉錄組測序共計獲得81853個Unigene, 將Unigene分別與NR、Pfam、Swissprot、eggNOG、GO和KEEG數據庫比對, 其中比對到NR數據庫的基因數最多為18278, 注釋率為22.33%; 比對到swissprot庫的基因數最少為12556, 注釋率為15.34%(表 3)。GO分析結果顯示, 共有13951個Unigene得到注釋。總共匹配到50個GO類, 歸為3大類, 分別是生物進程(25類, 50.00%)、細胞組分(15類, 30.00%)和分子功能(10類, 20.00%), 其中富集在細胞質和細胞核等細胞組分的基因數最多。有12848條Unigene注釋到KEGG數據庫, 分屬于40個代謝通路, 歸為6大類。其中注釋序列數量較多的代謝途徑分別為運輸和分解代謝; 翻譯、信號轉導; 折疊、分選和降解和糖代謝等。

表3 DIAMOND注釋結果Tab. 3 DIAMOND annotation results

2.3 差異表達基因篩選

采用edgeR分析發現Sp1群體在低溫脅迫下共有15773個差異表達基因, 其中上調145個, 下調15628個(圖 1A); Sp2群體在低溫脅迫下共有323個差異表達基因, 其中上調114個, 下調209個(圖 1B)。

圖1 低溫脅迫下不同形態群體擬穴青蟹肝胰腺基因表達差異火山圖Fig. 1 The ‘volcano plot’ picture of differentially expressed genes of Scylla paramamosain with different morphological populations under cold challenge

依據FDR值越小, 差異越顯著的原則, 從兩種形態擬穴青蟹在低溫脅迫下差異表達基因中篩選到差異表達最顯著的基因各10個(表 4)。

2.4 差異基因GO和KEGG通路富集分析

對所有差異表達基因進行GO富集和KEGG通路分析, 以獲得其潛在的功能注釋。Sp1群體在低溫脅迫下差異表達基因在細胞質基質、RNA結合、核糖體結構組分和翻譯等富集量較高(圖 2A);Sp2群體在低溫脅迫下差異表達基因在甾體羥化酶活性、DNA整合和碳酸氫鹽跨膜轉運體活性等富集較多(圖 2B)。

圖2 差異表達基因的GO富集分布散點圖Fig. 2 GO enrichment results of differentially expressed genes

對差異表達基因進行KEGG富集分析, 結果顯示, Sp1群體在低溫脅迫下差異表達基因在核糖體、剪接體和內質網蛋白質加工等信號通路富集較多(圖 3A), 而Sp2群體在低溫脅迫下差異表達基因在賴氨酸降解、甘氨酸、絲氨酸和蘇氨酸的代謝、精氨酸和脯氨酸的代謝等代謝通路富集較多(圖 3B)。

圖3 差異表達基因的KEEG富集分布散點圖Fig. 3 KEGG enrichment results of differentially expressed genes

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

為了驗證轉錄組測序結果的可靠性, 從表 4中挑選了5個上調基因和5個下調基因。如圖 4所示,轉錄組表達分析結果和熒光定量PCR結果基本一致, 表明本研究轉錄組測序分析結果是可靠的。

圖4 RT-qPCR驗證差異表達基因Fig. 4 Validation of differentially expressed genes by RT-qPCR

表4 低溫脅迫下兩種不同形態群體擬穴青蟹差異表達最顯著的前10個基因Tab. 4 The top 10 differentially expressed genes of Scylla paramamosain with different morphological populations under cold challenge

3 討論

肝胰腺作為甲殼動物重要的組成部分, 具有消化吸收、解毒、清除氧自由基和免疫防護等功能[21,22]。前期研究已明確溫度脅迫可影響擬穴青蟹肝胰腺的膜脂肪酸組成、線粒體呼吸速率、線粒體酶活性、抗氧化酶活性和熱休克蛋白表達等[4,23,24]。因此本研究以擬穴青蟹肝胰腺為材料, 通過轉錄組測序分析比較Sp1和Sp2群體在低溫脅迫下轉錄水平的差異。測序共獲得81853條Unigene, 其中22.33%為已知基因, 這在一定程度上豐富了擬穴青蟹轉錄組數據庫。

變溫動物在應對低溫時, 通常通過溫度補償機制來提高代謝速率, 以維持正常生命活動。但當環境溫度劇烈變化時, 一些變溫動物則會停止代謝進入休眠狀態, 以度過惡劣的環境[25―27]。水溫從20℃驟降至8℃, 在Sp2群體擬穴青蟹肝胰腺中檢測到

323個Unigene顯著差異表達, 其中上調114個, 下調209個。這些Unigene在甾體羥化酶活性、DNA整合和碳酸氫鹽跨膜轉運體活性等GO富集較多, 并進一步富集在賴氨酸降解、甘氨酸、絲氨酸和蘇氨酸的代謝、精氨酸和脯氨酸的代謝等代謝通路。水溫從24℃緩慢降至18℃(每天降溫2℃), 凡納濱對蝦(Litopenaeus vannamei)肝胰腺的轉錄組測序結果同樣表明低溫脅迫下差異表達基因富集在氨基酸的生物合成途徑以及精氨酸和脯氨酸的代謝中[28]。凡納濱對蝦3個家系的差異表達Unigene分別為1973、1677和687個[28], 多于經歷水溫驟降12℃的擬穴青蟹Sp2群體, 由此推測Sp2群體的擬穴青蟹對水溫快速下降的代謝響應相對遲緩。

對于Sp1群體, 差異表達基因篩選的結果顯示低溫脅迫下共有15773個Unigene顯著差異表達。廣溫性草魚(Ctenopharyngodon idellus)從27℃水溫緩慢降至12℃和4℃, 腦轉錄組測序分別檢測到1273和1987個差異表達Unigene[29]。對低溫敏感的熱帶魚虎皮魚(Puntius tetrazona)從27℃水溫驟降至13℃, 腦、鰓和肝臟轉錄組測序分別檢測到10272、7771和4632個差異表達Unigene[30]。由此可看出,面對水溫驟降, 擬穴青蟹Sp1群體的代謝響應遠比Sp2群體敏感。Sp1群體肝胰腺下調表達的Unigene在細胞質基質、RNA結合、核糖體結構組分和翻譯等GO富集量較高, 并進一步富集在核糖體、剪接體和內質網蛋白質加工等信號通路, 這說明Sp1群體擬穴青蟹蛋白質合成加工過程大幅度減弱。

面對不利環境, Sp1群體擬穴青蟹能快速地降低生理活動, 減緩低溫對機體的損傷, 而Sp2群體擬穴青蟹在應對低溫時, 僅調整了物質代謝, 這可能使得Sp2群體在應對低溫時機體更容易受到傷害,甚至導致死亡現象的發生。在早期的研究中, 我們調研了福建漳州海區的Sp1和Sp2群體, 其中Sp1群體的數量是Sp2群體數量的3倍多, 這可能與兩個群體對溫度變化存在適應差異有關。在本研究的基礎上, 我們將進一步挖掘轉錄組數據同時結合蛋白質組學分析尋找其他低溫脅迫適應性指標來進一步揭示Sp1和Sp2群體對溫度的適應差異。