大黄鱼是我国养殖规模最大的海水鱼类, 被列为我国八大优势出口养殖水产品之一(农业农村部渔业渔政管理局等, 2021; 韩承义等, 2022)。传统的大黄鱼养殖主要集中在近海海域, 近海网箱养殖过程中, 除陆源输入外, 排泄物及残饵以及养殖区自身污染物的大量排放可引起局部海域氮磷营养盐、有机质含量过高, 造成水质恶化。随着大黄鱼养殖规模的扩大, 病害频发、品质退化、成活率低等问题日益突出, 严重影响了大黄鱼的产业效益(Cao et al, 2015; Ma et al, 2020; 徐鹏等, 2022)。因此, 推动大黄鱼深远海养殖产业发展是实现大黄鱼养殖提质增效、产业转型升级的有效途径, 也是践行大食物观, 向江河湖海要食物, 推进我国深蓝渔业发展的重要战略举措(徐皓等, 2021)。深远海养殖海域营养丰富、水体交换率高、污染物含量少, 配合现代化养殖装备及先进的养殖技术, 可高效获得高品质的养殖大黄鱼(邹国华等, 2021)。然而, 深远海养殖水域浪高流急、海况恶劣, 且易受台风大浪的侵袭, 这对于大黄鱼的养殖是一个巨大的挑战, 要求养殖大黄鱼需要具备较强的抗流能力。目前, 通过国审的4个大黄鱼新品种其选育的性状主要聚焦在生长、条形及耐低温等方面, 这些品种无法应对复杂的海域环境, 特别是较快的水流条件。深远海养殖鱼类新品种现已被列入《2022农业农村产业发展重大技术需求》。因此, 为满足深远海养殖大黄鱼种业发展需求, 需培育抗流能力强适宜深远海养殖的大黄鱼新品种, 以推动我国大黄鱼深远海养殖产业的高质量发展。

鱼类性状的遗传解析可为新品种创制及其育种技术开发提供重要的理论基础。目前已有的针对大黄鱼性状遗传基础解析的研究主要聚焦在疾病(Zhao et al, 2021; Bai et al, 2022; Zhang et al, 2022)、低氧(Mu et al, 2020)、低温(曾霖等, 2022)、生长(Zhou et al, 2019)、性别控制(Luo et 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) ℃, 盐度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 ℃冰箱保存。

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 Hiseq^TM 2500平台上测序。

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标准化, 以差异倍数|log₂ fold 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; ddH₂O 10.5 μL, 总体积25 μL。采用两步法PCR扩增标准程序, 循环参数为反应程序: 95 ℃预变性10 min, 95 ℃变性15 s, 60 ℃退火延伸1 min, 共40个循环。对所得数据进行溶解曲线分析, 用2^{–∆ ∆ Ct}法计算基因的相对表达量。

1.3 数据分析

累计死亡率数据($\bar x$±SD)用SPSS 20.0进行独立样本T检验来检验差异性, P < 0.05表明差异显著。

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%, Q₃₀ > 92.9% (表 2), 表明测序数据质量可用于对后续数据分析。

将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)。

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.3 DEGs的GO富集分析

为进一步探究HM组与SM组DEGs的功能, 对DEGs进行GO功能注释。DEGs在细胞组分(cellular component, CC)、分子功能(molecular function, MF)和生物学过程(biological process, BP)中均有富集。在CC中显著富集在肌钙蛋白复合物、蛋白酶体辅助复合物、肌球蛋白、线粒体内膜等条目; 在MF中显著富集在MAP激酶酪氨酸/丝氨酸/苏氨酸磷酸酶活性条目; 在BP中显著富集在Rho蛋白信号转导的调节条目(图 4)。

2.4 DEGs的KEGG富集分析

本研究中共有24 177个基因在KEGG数据库中匹配, 并将富集程度最显著的20条通路进行散点图绘制。HM组与SM组相比, 上调的DEGs富集到的信号通路主要与能量代谢和信号转导相关, 包括心肌收缩、氧化磷酸化信号通路、黏着斑、ECM-受体相互作用、AGE-RAGE、MAPK、心肌细胞中的肾上腺素能等信号通路, 而下调的DEGs富集的信号通路主要为真核生物核糖体的发生、蛋白酶体、内质网中的蛋白加工、RNA转运以及泛素介导的蛋白水解等(图 5)。

2.5 RT-qPCR验证

随机选取转录组数据中8个DEGs做RT-qPCR验证。RT-qPCR所得到的基因表达情况与转录组获得的基因表达趋势一致, 说明转录组结果可靠(图 6)。

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)。

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聚类分析的结果显示其一致性较好。众所周知, 鱼类表型性状受遗传与环境的双重控制(Jonsson et al, 2019; Raffini et al, 2020; Downie et al, 2021), 本研究中大黄鱼的来源相同, 因此可以断定不同大黄鱼个体间抗流能力的差异主要受遗传因素的影响, SM组个体间的抗流能力存在较大差异进而导致了DEGs的表达模式未能呈现出较好的一致性。

3.2 DEGs的GO富集分析

肌球蛋白、肌钙蛋白复合物是构成肌肉结构的关键组成部分, 肌节是肌肉收缩的基本单位, 参与各种类型的细胞运动和细胞骨架的维持(Drazic et al, 2018; Ochiai et al, 2020)。本研究发现, 不同抗流能力大黄

鱼群体间的DEGs主要富集在肌球蛋白、肌钙蛋白复合物、肌节和Rho蛋白信号转导的调控等与肌肉收缩功能有关的条目, 这表明了大黄鱼肌肉组织中肌球蛋白、肌钙蛋白复合物与肌节等重要组分的生理功能与其抗流能力的强弱之间存在密切的关联。这一点在斑马鱼(Danio rerio)和金头鲷(Sparus aurata)的研究中也得到了证实, 其游泳能力与骨骼肌肌球蛋白和肌钙蛋白基因的表达密切相关(Van Der Meulen et al, 2006; Palstra et al, 2020)。当然, 鱼类肌肉组织中重要组分的生理功能与其相关基因的表达不仅与遗传因素相关, 与环境因素同样密切相关。Moya等(2019)指出持续性的游泳会使金头鲷(Sparus aurata)幼鱼肌纤维特征发生明显变化, 使虹鳟(Oncorhynchus mykiss)肌肉收缩发育相关的多个基因表达发生改变(Palstra et al, 2013), 此外, Totland等(1987)研究表明较强的水流条件会使大西洋鲑鱼(Salmo salar)白肌中单位面积的纤维数量减少, 单个肌肉纤维的平均面积增加。本研究中, 实验所用大黄鱼的养殖环境因素是一致的, 因此可以确定, 遗传因素的差异影响了大黄鱼个体间肌肉组织中重要组分的生理功能进而决定了其抗流能力的强弱。

3.3 DEGs的KEGG通路分析

心肌收缩对鱼类维持正常的生理活动具有重要作用(Singh et al, 2016; Kubly et al, 2019; Nabeebaccus, 2022), 周盈颖(2001)发现心肌收缩受心肌细胞中肾上腺素能信号通路的调控。本研究发现, 在抗流能力强的大黄鱼肌肉组织中心肌收缩相关基因的表达显著上调, 表明大黄鱼心肌收缩功能与其抗流能力间存在正相关关系。心输出量和血氧含量是鱼类游泳能力的重要指标, 持续游泳中鱼类通过心输出量和血氧含量的改变影响能量代谢和氧气利用, 来维持细胞稳态。Castro等(2013)发现大西洋鲑鱼的游泳训练可以通过激活特定基因来诱导心脏的肥大和增生参与对心脏生长的调节。反过来, 心脏生长可以潜在地增加心输出量, 促进氧气和营养物质的输送。Li等(2023)利用临界游泳速度对筛选得到的快速和慢速游泳的大黄鱼幼鱼进行转录组分析发现差异基因主要富集在氧化磷酸化、心肌收缩等途径, 说明游泳过程中大黄鱼通过增加心输出量影响能量代谢和氧气利用, 来维持细胞稳态。同样地, 具有较强游泳能力的蓝鳍金枪鱼(Thunnus orientalis), 强大的心肌收缩能力是其在快速游泳状态下保持心室充盈、维持心输出量的重要生理基础(Ciezarek et al, 2020)。除此之外, 较强的水流条件也会引起鱼类新陈代谢的显著变化, 进而引起机体能量需求的增加(袁喜等, 2012; Schreck et al, 2016; Yu et al, 2022; 柴若愚等, 2023)。研究已证实, 氧化磷酸化信号通路不仅与ATP的生成途径有关, 还是能量代谢的重要枢纽(Lin et al, 2020; Su et al, 2020)。本研究中, HM组大黄鱼肌肉组织中参与氧化磷酸化信号通路的基因表达显著上调, 这些基因编码的产物已证实在呼吸链电子传递过程中发挥重要作用(Neel et al, 2020), 推测抗流组个体中氧化磷酸化产生的ATP是其具有较强抗流能力的能量基础。黏着斑是细胞外基质和细胞内骨架之间联系的桥梁, 对细胞运动有重要作用(Moulder et al, 2010; Li et al, 2019; Mackay et al, 2021), 而ECM-受体相互作用信号通路可参与调控细胞生长分化(Kanchanawong et al, 2023)。本研究结果显示, 抗流能力强的大黄鱼肌肉组织中黏着斑与ECM-受体相互作用信号通路的相关基因同样显著上调。

3.4 抗流分组后累计死亡率分析

应激是影响鱼类生命活动的重要影响因素之一(Nitz et al, 2020)。于丽娟等(2014)发现不同水流速度下的中华倒刺鲃肌肉(Spinibarbus sinensis)自由基代谢会发生改变, 而自由基代谢紊乱会引起氧化/还原系统的失衡, 进而导致鱼体的应激损伤。本研究中SM组大黄鱼肌肉组织中泛素介导的蛋白水解、蛋白酶体、内质网中的蛋白加工等与氧化应激相关的信号通路基因明显上调, 表明水流刺激显著诱发了SM组大黄鱼的应激反应。娄宇栋等(2019)在对美国红鱼(Sciaenops ocellatus)肝脏转录组分析中同样发现, 高流速下美国红鱼会触发细胞内的胁迫相关基因PKA、Wnt的上调, 这表明过快的水流会引起鱼类的氧化应激反应。内质网作为外界环境的压力感受器(Yoo et al, 2017), 可以启动机体通过自噬或细胞凋亡来清除损伤的细胞。本研究内质网中蛋白加工信号通路被显著富集, 而其在介导错误折叠的蛋白进入胞质中发挥重要作用(Iurlaro et al, 2016), 这也进一步表明了水流的刺激导致了SM组大黄鱼机体内产生了应激损伤, 此是导致SM组累计死亡率明显高于HM组的主要原因之一。

4 结论

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