不同抗流能力大黃魚(Larimichthys crocea)肌肉轉錄組學差異分析*

張靜靜 王亞冰 王 倩 韓多彩 彭士明①

(1.中國水產科學研究院東海水產研究所 上海 200090; 2.中國農業科學院研究生院 北京 100081)

大黃魚是我國養殖規模最大的海水魚類, 被列為我國八大優勢出口養殖水產品之一(農業農村部漁業漁政管理局等, 2021; 韓承義等, 2022)。傳統的大黃魚養殖主要集中在近海海域, 近海網箱養殖過程中,除陸源輸入外, 排泄物及殘餌以及養殖區自身污染物的大量排放可引起局部海域氮磷營養鹽、有機質含量過高, 造成水質惡化。隨著大黃魚養殖規模的擴大,病害頻發、品質退化、成活率低等問題日益突出, 嚴重影響了大黃魚的產業效益(Caoet al, 2015; Maet al,2020; 徐鵬等, 2022)。因此, 推動大黃魚深遠海養殖產業發展是實現大黃魚養殖提質增效、產業轉型升級的有效途徑, 也是踐行大食物觀, 向江河湖海要食物,推進我國深藍漁業發展的重要戰略舉措(徐皓等,2021)。深遠海養殖海域營養豐富、水體交換率高、污染物含量少, 配合現代化養殖裝備及先進的養殖技術, 可高效獲得高品質的養殖大黃魚(鄒國華等,2021)。然而, 深遠海養殖水域浪高流急、海況惡劣,且易受臺風大浪的侵襲, 這對于大黃魚的養殖是一個巨大的挑戰, 要求養殖大黃魚需要具備較強的抗流能力。目前, 通過國審的4 個大黃魚新品種其選育的性狀主要聚焦在生長、條形及耐低溫等方面, 這些品種無法應對復雜的海域環境, 特別是較快的水流條件。深遠海養殖魚類新品種現已被列入《2022 農業農村產業發展重大技術需求》。因此, 為滿足深遠海養殖大黃魚種業發展需求, 需培育抗流能力強適宜深遠海養殖的大黃魚新品種, 以推動我國大黃魚深遠海養殖產業的高質量發展。

魚類性狀的遺傳解析可為新品種創制及其育種技術開發提供重要的理論基礎。目前已有的針對大黃魚性狀遺傳基礎解析的研究主要聚焦在疾病(Zhaoet al, 2021; Baiet al, 2022; Zhanget al, 2022)、低氧(Muet al, 2020)、低溫(曾霖等, 2022)、生長(Zhouet al,2019)、性別控制(Luoet al, 2019)等方面, 針對大黃魚抗流性狀的遺傳基礎解析尚未見有相關報道。本研究基于轉錄組技術(RNA-seq)對不同抗流能力大黃魚肌肉組織進行轉錄組分析, 以期從轉錄組學的層面解析大黃魚抗流性狀的分子基礎, 研究結果可為培育適宜深遠海養殖大黃魚抗流新品種奠定理論基礎。

1 材料與方法

1.1 實驗魚抗流分組及樣本采集

2021 年9 月從福建省福鼎市沙埕灣主養區收集1 000 余尾大黃魚, 購買于同一家養殖戶, 這批魚來自同一批受精卵(來源于寧德市富發水產有限公司),養殖區位于沙埕灣后港村的養殖漁排, 體重(274.58±14.09) g, 體長(23.9±0.96) cm, 為1 齡大黃魚, 其養殖周期、養殖模式及所處水域的水文特征均為一致。大黃魚暫養于2 口室內水泥池(直徑6 m, 水深1.5 m,500 尾/池)中, 進行一周的。暫養條件為: 水溫(27.0±1.0) °C, 鹽度25～26, 溶氧量 DO>6 mg/L。每天早晚兩次投喂商業飼料至飽食狀態, 日換水量1/2。

在進行抗流篩選實驗前大黃魚禁食24 h, 通過抗流實驗裝置(圖1)進行大黃魚抗流能力強弱的分篩,抗流裝置內水深0.5 m, 每次取50 條魚放入水槽中(E區), 用充氣泵保持足夠的溶氧水平, 水流控制采用無極變速, 以實驗魚體長(約20 cm/s)的流速作為初始流速, 使實驗魚適應20～30 min, 隨后在1 min 內使流速增加到1.0 m/s, 觀察并記錄大黃魚的抗流時間,將抗流時間>30 min 的大黃魚歸為HM 組, 抗流時間<5 min 且以抵觸隔網30 s 作為篩選條件確定為SM組。用于進行樣品采集及統計48 h 累計死亡率的實驗魚其具體操作流程如下: 在分篩實驗前, 根據大黃魚的游動能力, 將所有大黃魚初步分為了游動能力強組和游動能力弱組。操作的目的除了提高分篩的效率外, 主要是考慮到抗流裝置每次分篩只能一次性放置50 尾魚, 實驗魚被篩選出來的先后順序是不同的。為了準確統計抗流組和非抗流組的48 h 累計死亡率, 避免相同實驗組內實驗魚在時間上的不同步性, 在初步分組的基礎上, 分別針對游動能力強組和游動能力弱組各進行3 批次的抗流分篩。針對游動能力強組的3 批次抗流分篩, 分別篩選到35、38、39尾HM 組實驗魚, 每次分篩結束后將篩選到的HM 組實驗魚放置于單獨的一個水泥池中(直徑4 m, 水深1 m), 然后分別統計3 個水泥池中實驗魚的48 h 累計死亡率。針對游動能力弱組的3 批次抗流分篩, 分別篩選到37、40、36 尾SM 組實驗魚, 采用上述相同的方法統計48 h 累計死亡率。48 h 后, 分別從HM 組、SM 組隨機收集18 尾大黃魚, 以20 mg/L 丁香酚溶液進行麻醉, 統一收集背鰭中間下方1 cm 左右的肌肉組織0.2 g, 每6 尾大黃魚肌肉組織等量混合放置于2 mL 的無菌凍存管中, 每個處理組共收集3 管, 經液氮速凍后置于–80 °C 冰箱保存。

圖1 大黃魚抗流實驗裝置示意圖Fig.1 The schematic diagram of anti-flow experiment of large yellow croaker

1.2 轉錄組測序

1.2.1 總RNA 提取、建庫與測序 根據Trizol 試劑盒(Invitrogen, 美國)提取每管肌肉組織的總RNA,并分別使用 Nanodrop?分光光度計、Agilent 2100 生物分析儀進行RNA 的純度、濃度及完整性檢測。用帶有Oligo(dT)的磁珠富集mRNA, 以mRNA 為模板,合成cDNA, 利用AM Pure XP beads 純化cDNA, 純化的雙鏈cDNA 再進行末端修復、加poly(A)尾并連接測序接頭, 然后用AMPure XP beads 進行片段大小選擇;最后通過PCR 富集得到cDNA 文庫。檢測合格的cDNA 文庫在Illumina HiseqTM2500 平臺上測序。

1.2.2 轉錄組分析 利用Fast QC 軟件對測序得到的原始數據進行質量評估, 去除原始數據中含有帶接頭和低質量的片段, 保證獲得的數據具有品質較高的有效測序數據(clean reads)。利HISAT2 軟件, 使用BWT 算法將clean reads 與大黃魚參考基組(ftp://ftp.ncbi.nlm.nih.gov/genomes/all/GCF/000/972/845/GCF_000972845.2_L_crocea_2.0)進行比對。使用 DESeq2軟件采用 FPKM 計算方法將 clean data 標準化, 以差異倍數|log2fold change|≥1 且P<0.05 為標準篩選差異表達基因(DEGs)。對組裝的序列進行遺傳相似性比較, 以進行蛋白質功能注釋, 將單基因與GO、COG(clusters of orthologous groups) 、KOG (eukaryotic orthologous groups)、KEGG、Nr (non- redundant)、Swiss-Pro 和Pfam 數據庫進行比對并注釋。根據注釋結果, 分別對DEGs 的功能和生物學途徑進行分類。將DEGs 與 GO 數據庫和 KEGG 數據庫進行比對,從而獲得 GO 功能注釋和KEGG 通路富集分析。

1.2.3 DEGs RT-qPCR 驗證 隨機選取DEGs 中的8 個基因, 以β-actin為內參基因, 進行RT-qPCR 檢驗,特異性引物使用軟件NCBI 在線設計(表1)。使用Quant Studio 6 Flex 進行RT-qPCR, 反應體系按 Ultra SYBR Mixture (康為, 中國)試劑盒說明書進行。反應體系為: 2×Ultra SYBR Mixture 12.5 μL; Forward primer (10 μmol/L) 0.5 μL; Reverse primer (10 μmol/L)0.5 μL; Template cDNA 1 μL; ddH2O 10.5 μL, 總體積25 μL。采用兩步法PCR 擴增標準程序, 循環參數為反應程序: 95 °C 預變性 10 min, 95 °C 變性15 s,60 °C 退火延伸 1 min, 共40 個循環。對所得數據進行溶解曲線分析, 用2–??Ct法計算基因的相對表達量。

表1 實時熒光定量PCR 引物Tab.1 Primers used for real-time PCR analysis

1.3 數據分析

2 結果與分析

2.1 轉錄組測序質量及注釋分析

利用RNA-seq 對HM組、SM組進行測序分析, 對測得的原始數據進一步質控后得到 21 522 509～24 869 415 條有效測序量, 與大黃魚參考基因組比對后得到41 412 309～47 651 205 條匹配的序列, 比對率為 95.27%～96.37%, 拼接共獲得 33 191 個單基因(unigene)。本試驗中GC 含量區間為51.51%～52.15%,Q30>92.9% (表2), 表明測序數據質量可用于對后續數據分析。

表2 轉錄組測序數據統計Tab.2 Summary statistics of transcriptome sequencing date

將unigene 與常用數據庫進行比對并注釋, 發現共有24 252 條unigene 被注釋, 其中有18 494 個unigene在GO 中注釋(76.26%), 長度大于1 000 bp 有12 895個 unigene; 在 KEGG 注釋到 24 177 個 unigene(99.69%), 長度大于1 000 bp 有15 871個單基因(表3)。

表3 Unigene 注釋統計Tab.3 Summary of unigene annotation

2.2 DEGs 分析

將HM 組與SM 組肌肉組織DEGs 進行對比分析。將基因表達量倍數設置為二倍以上,P<0.05 的基因作為DEGs。由HM 組與SM 大黃魚肌肉組織的RNA 測序和基因表達譜分析可知: DEGs 的總數一共有1 806, 其中有1 090 個基因顯著上調(P<0.05; 圖2, 紅色), 716 個基因顯著下調(P<0.05; 圖2, 綠色)。聚類分析熱圖顯示, HM 與SM 組間差異基因表達模式相差較大, HM組個體間表達模式相似, 但是SM 組個體間基因表達模式存在明顯差異(圖3)。

圖2 差異表達基因火山圖Fig.2 Volcano diagram of differentially expressed genes

圖3 組間差異表達基因分層聚類熱圖Fig.3 Heatmap of differentially expressed genes between groups

2.3 DEGs 的GO 富集分析

為進一步探究HM 組與SM 組DEGs 的功能, 對DEGs 進行GO 功能注釋。DEGs 在細胞組分(cellular component, CC)、分子功能(molecular function, MF)和生物學過程(biological process, BP)中均有富集。在CC 中顯著富集在肌鈣蛋白復合物、蛋白酶體輔助復合物、肌球蛋白、線粒體內膜等條目; 在MF 中顯著富集在MAP 激酶酪氨酸/絲氨酸/蘇氨酸磷酸酶活性條目; 在BP 中顯著富集在Rho 蛋白信號轉導的調節條目(圖4)。

圖4 轉錄組差異表達基因的GO 分析Fig.4 Gene ontology (GO) item analyses of differentially expressed genes in transcriptome

2.4 DEGs 的KEGG 富集分析

本研究中共有24 177 個基因在KEGG 數據庫中匹配, 并將富集程度最顯著的20 條通路進行散點圖繪制。HM 組與SM 組相比, 上調的DEGs 富集到的信號通路主要與能量代謝和信號轉導相關, 包括心肌收縮、氧化磷酸化信號通路、黏著斑、ECM-受體相互作用、AGE-RAGE、MAPK、心肌細胞中的腎上腺素能等信號通路, 而下調的DEGs 富集的信號通路主要為真核生物核糖體的發生、蛋白酶體、內質網中的蛋白加工、RNA 轉運以及泛素介導的蛋白水解等(圖5)。

圖5 肌肉組織差異表達基因KEGG 功能富集分布Fig.5 KEGG enrichment analysis of differentially expressed genes in muscle tissue

2.5 RT-qPCR 驗證

隨機選取轉錄組數據中8 個DEGs 做RT-qPCR驗證。RT-qPCR 所得到的基因表達情況與轉錄組獲得的基因表達趨勢一致, 說明轉錄組結果可靠(圖6)。

圖6 RNA-Seq 與RT-qPCR 的基因表達比較分析Fig.6 Comparative analysis of differentially expressed genes of RNA-Seq and RT-qPCR

2.6 抗流分組后48 h 累計死亡率

抗流分組后, HM 組和SM 組大黃魚均出現死亡情況, 但SM 組的死亡率在各個時間點均極顯著高于HM 組(P<0.01)。其中SM 組在抗流分組后12 h 內就出現死亡現象, 48 h 內累計死亡率達到了18.00%, 而HM 組在抗流試驗24 h 內才出現死亡情況, 48 h 內累計死亡率為5.33% (表4)。

表4 大黃魚累計死亡率Tab.4 Cumulative mortality rate of large yellow croaker

3 討論

3.1 DEGs 聚類分析

本研究采用RNA-Seq 技術對比揭示了不同抗流能力大黃魚群體間肌肉組織的基因表達差異, 為進一步挖掘定位大黃魚抗流性狀關鍵控制基因奠定了重要基礎。通過DEGs 聚類分析的結果發現, 抗流弱的SM 組DEGs 并沒有全部聚類在一起, 個體間差異較為明顯, 其原因可能與SM 組個體間的抗流能力存在較大差異有關。在本實驗條件下, SM 組大黃魚其抗流時間的跨度從0～5 min, 換言之, 有的個體抗流的時間可以維持將近5 min, 但有的個體其抗流的時間可能不足 1 min, 甚至個別個體的抗流時間僅有20～30 s, 由此導致了SM 組個體間的抗流能力差異較大; 而抗流能力強的HM 組中, 所有個體的抗流時間均在30 min 以上, 意味著所有個體均具有較強的抗流能力, 個體間的抗流能力差異不大, 所以HM 組DEGs 聚類分析的結果顯示其一致性較好。眾所周知,魚類表型性狀受遺傳與環境的雙重控制(Jonssonet al,2019; Raffiniet al, 2020; Downieet al, 2021), 本研究中大黃魚的來源相同, 因此可以斷定不同大黃魚個體間抗流能力的差異主要受遺傳因素的影響, SM 組個體間的抗流能力存在較大差異進而導致了DEGs 的表達模式未能呈現出較好的一致性。

3.2 DEGs 的GO 富集分析

肌球蛋白、肌鈣蛋白復合物是構成肌肉結構的關鍵組成部分, 肌節是肌肉收縮的基本單位, 參與各種類型的細胞運動和細胞骨架的維持(Drazicet al, 2018;Ochiaiet al, 2020) 。本研究發現, 不同抗流能力大黃魚群體間的DEGs 主要富集在肌球蛋白、肌鈣蛋白復合物、肌節和Rho 蛋白信號轉導的調控等與肌肉收縮功能有關的條目, 這表明了大黃魚肌肉組織中肌球蛋白、肌鈣蛋白復合物與肌節等重要組分的生理功能與其抗流能力的強弱之間存在密切的關聯。這一點在斑馬魚(Danio rerio)和金頭鯛(Sparus aurata)的研究中也得到了證實, 其游泳能力與骨骼肌肌球蛋白和肌鈣蛋白基因的表達密切相關(Van Der Meulenet al,2006; Palstraet al, 2020)。當然, 魚類肌肉組織中重要組分的生理功能與其相關基因的表達不僅與遺傳因素相關, 與環境因素同樣密切相關。Moya 等(2019)指出持續性的游泳會使金頭鯛(Sparus aurata)幼魚肌纖維特征發生明顯變化, 使虹鱒(Oncorhynchus mykiss)肌肉收縮發育相關的多個基因表達發生改變(Palstraet al, 2013), 此外, Totland 等(1987)研究表明較強的水流條件會使大西洋鮭魚(Salmo salar)白肌中單位面積的纖維數量減少, 單個肌肉纖維的平均面積增加。本研究中, 實驗所用大黃魚的養殖環境因素是一致的, 因此可以確定, 遺傳因素的差異影響了大黃魚個體間肌肉組織中重要組分的生理功能進而決定了其抗流能力的強弱。

3.3 DEGs 的KEGG 通路分析

心肌收縮對魚類維持正常的生理活動具有重要作用(Singhet al, 2016; Kublyet al, 2019; Nabeebaccus,2022), 周盈穎(2001)發現心肌收縮受心肌細胞中腎上腺素能信號通路的調控。本研究發現, 在抗流能力強的大黃魚肌肉組織中心肌收縮相關基因的表達顯著上調, 表明大黃魚心肌收縮功能與其抗流能力間存在正相關關系。心輸出量和血氧含量是魚類游泳能力的重要指標, 持續游泳中魚類通過心輸出量和血氧含量的改變影響能量代謝和氧氣利用, 來維持細胞穩態。Castro 等(2013)發現大西洋鮭魚的游泳訓練可以通過激活特定基因來誘導心臟的肥大和增生參與對心臟生長的調節。反過來, 心臟生長可以潛在地增加心輸出量, 促進氧氣和營養物質的輸送。Li 等(2023)利用臨界游泳速度對篩選得到的快速和慢速游泳的大黃魚幼魚進行轉錄組分析發現差異基因主要富集在氧化磷酸化、心肌收縮等途徑, 說明游泳過程中大黃魚通過增加心輸出量影響能量代謝和氧氣利用, 來維持細胞穩態。同樣地, 具有較強游泳能力的藍鰭金槍魚(Thunnus orientalis), 強大的心肌收縮能力是其在快速游泳狀態下保持心室充盈、維持心輸出量的重要生理基礎(Ciezareket al, 2020)。除此之外,較強的水流條件也會引起魚類新陳代謝的顯著變化,進而引起機體能量需求的增加(袁喜等, 2012; Schrecket al, 2016; Yuet al, 2022; 柴若愚等, 2023)。研究已證實, 氧化磷酸化信號通路不僅與ATP 的生成途徑有關, 還是能量代謝的重要樞紐(Linet al, 2020; Suet al,2020)。本研究中, HM 組大黃魚肌肉組織中參與氧化磷酸化信號通路的基因表達顯著上調, 這些基因編碼的產物已證實在呼吸鏈電子傳遞過程中發揮重要作用(Neelet al, 2020), 推測抗流組個體中氧化磷酸化產生的ATP 是其具有較強抗流能力的能量基礎。黏著斑是細胞外基質和細胞內骨架之間聯系的橋梁,對細胞運動有重要作用(Moulderet al, 2010; Liet al,2019; Mackayet al, 2021), 而ECM-受體相互作用信號通路可參與調控細胞生長分化(Kanchanawonget al,2023)。本研究結果顯示, 抗流能力強的大黃魚肌肉組織中黏著斑與ECM-受體相互作用信號通路的相關基因同樣顯著上調。

3.4 抗流分組后累計死亡率分析

應激是影響魚類生命活動的重要影響因素之一(Nitzet al, 2020)。于麗娟等(2014)發現不同水流速度下的中華倒刺鲃 肌肉(Spinibarbus sinensis)自由基代謝會發生改變, 而自由基代謝紊亂會引起氧化/還原系統的失衡, 進而導致魚體的應激損傷。本研究中SM 組大黃魚肌肉組織中泛素介導的蛋白水解、蛋白酶體、內質網中的蛋白加工等與氧化應激相關的信號通路基因明顯上調, 表明水流刺激顯著誘發了SM 組大黃魚的應激反應。婁宇棟等(2019)在對美國紅魚(Sciaenops ocellatus)肝臟轉錄組分析中同樣發現, 高流速下美國紅魚會觸發細胞內的脅迫相關基因PKA、Wnt的上調,這表明過快的水流會引起魚類的氧化應激反應。內質網作為外界環境的壓力感受器(Yooet al, 2017), 可以啟動機體通過自噬或細胞凋亡來清除損傷的細胞。本研究內質網中蛋白加工信號通路被顯著富集, 而其在介導錯誤折疊的蛋白進入胞質中發揮重要作用(Iurlaroet al, 2016), 這也進一步表明了水流的刺激導致了SM組大黃魚機體內產生了應激損傷, 此是導致SM 組累計死亡率明顯高于HM 組的主要原因之一。

4 結論

本研究基于不同抗流能力大黃魚群體間肌肉組織DEGs 的GO 功能與KEGG 通路富集分析, 推測抗流能力強的大黃魚可能是通過改變肌肉收縮的方式及能量代謝來提高機體應對水流的能力。結合死亡率等表型數據的分析, 進一步印證了較強的抗流能力可顯著提高大黃魚在強水流刺激下的成活率, 研究結果為深入挖掘抗流相關基因提供了基礎信息, 并為培育適宜深遠海養殖的大黃魚抗流新品種奠定了理論基礎。