不同抗流能力大黄鱼(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 通路富集分析, 推测抗流能力强的大黄鱼可能是通过改变肌肉收缩的方式及能量代谢来提高机体应对水流的能力。结合死亡率等表型数据的分析, 进一步印证了较强的抗流能力可显著提高大黄鱼在强水流刺激下的成活率, 研究结果为深入挖掘抗流相关基因提供了基础信息, 并为培育适宜深远海养殖的大黄鱼抗流新品种奠定了理论基础。