豹纹鳃棘鲈差异流速下肝脏转录组分析

高 进，王永波，刘金叶，郭一兰，符书源,3

1.海南热带海洋学院/热带海洋生物资源利用与保护教育部重点实验室，海南 三亚 572022

2.海南省海洋与渔业科学院，海南 海口 571126

3.海南省热带海水养殖工程技术研究中心，海南 海口 571126

豹纹鳃棘鲈 (Plectropomus leopardus) 俗称东星斑，是一种广泛分布于热带和亚热带海域的暖水性礁栖鱼类[1]。因色彩艳丽、肉质细嫩、经济价值较高，豹纹鳃棘鲈的市场前景广阔，养殖发展潜力巨大[2]。生长性状为水产养殖鱼类的主要经济性状，主要受基因、环境以及两者间相互作用的影响[3]。阐明生长发育过程中遗传因子与环境因子协同作用的调控机理，将为豹纹鳃棘鲈生长相关性状的进一步遗传改良提供理论基础。迄今，国内外针对豹纹鳃棘鲈的研究多集中于人工繁育[4]、养殖技术[5]、营养[6]、免疫[7]和生态[8]等方面。目前豹纹鳃棘鲈的主要养殖模式有工厂化流水养殖和网箱养殖两种，工厂化流水养殖由于具有可控性强、水体稳定和管理方便等优势而成为主要养殖模式[9]，流速是该养殖模式下影响豹纹鳃棘鲈生长的重要环境因子[10-12]。然而，流速差异对豹纹鳃棘鲈生长发育的分子作用机制尚未见报道。

受多变的水流速度甚至湍流的影响，自然界中的鱼类特别是洄游性和礁栖性鱼类，通常需要选择合适的栖息环境以平衡适应环境的能量消耗与自身生长发育之所需[10]。随着陆基水产养殖技术的快速发展，流速作为鱼类养殖环境中一个极为重要的环境因子被广泛研究。Ogata和Oku[11]开展了流速对日本牙鲆 (Paralichthys olivaceus) 幼鱼生长性能的影响实验，结果表明流速能够改变日本牙鲆的体脂肪含量。Merino等[12]在研究差异流速对北美牙鲆 (P.californicus) 幼鱼生长的影响时发现，差异流速会显著影响其的饲料转化率和生长速度。对循环流水水产养殖系统中的大菱鲆 (Scophthalmus maximus) 幼鱼进行差异流速下生长性能和水质的研究发现，在4种流速养殖实验中，大菱鲆的比生长速率会随着流速的增加显著提高，而氨氮总量和弧菌总数对鱼体的影响及其在水体中的累积却随之下降，进而改善了大菱鲆的养殖环境[13]。此外，流速同样影响其他水生生物的生长发育过程，在差异流速下对全缘马尾藻 (Sargassum integerrimum) 幼孢子的体质量、叶绿素a含量、超氧化物歧化酶(Superoxide dismutase, SOD)、过氧化氢酶 (Catalase,CAT) 和蛋白质浓度等生理指标进行比对研究，发现流速对其体质量增长和各种生化指标均产生了显著影响[14]。因此，研究流速对水生生物的影响有助于提高养殖经济效益和促进流水健康养殖的可持续发展。

基于二代高通量测序的RNA-Seq是分析不同组织在特定生理状态下基因表达情况的技术，由于转录组会随着生长发育阶段和外界环境变化而改变，其被广泛应用于生长发育、应激生理和抗病免疫等生物学过程中与某些特定功能相关的差异表达或未知基因的发掘研究[15-17]。随着测序技术的日臻成熟和成本下降，RNA-Seq被应用于多种水产动物的关键基因网络调控研究[18-19]，如在光周期影响尼罗罗非鱼 (Oreochromis niloticus) 脑组织[20]、卵形鲳鲹 (Trachinotus ovatus) 盐度胁迫[21]和幼体克氏海马 (Hippocampus kelloggi) 温度胁迫[22]研究中的应用均表明，该技术已成为阐明水产动物环境适应性分子作用机制的重要手段，为开展目标性状调控机制和育种实践研究提供了参考。Wang等[23]利用RNA-Seq对豹纹鳃棘鲈的体色差异进行转录分析，共筛选到38个参与体色差异调控的候选基因。Mekuchi等[24]采用RNA-Seq和代谢组方法对豹纹鳃棘鲈的代谢调控机制进行研究，结果表明支链氨基酸在豹纹鳃棘鲈禁食期间的能量代谢中发挥重要作用，为豹纹鳃棘鲈的代谢机制解析提供了基础数据支持。本文基于RNA-Seq技术，比较了高、低流速下豹纹鳃棘鲈肝脏的基因表达差异，探究了豹纹鳃棘鲈对于不同流速水文环境的适应性分子调节机制，为开展豹纹鳃棘鲈“陆海接力”、网箱养殖和增殖放流等提供基础理论依据。

1 材料与方法

1.1 实验材料

实验所用豹纹鳃棘鲈为同一繁育批次的平均体质量100 g、平均体长约20 cm的6月龄幼鱼，养殖于海南海研热带海水鱼类良种场 (110°68'00''E,19°37'12''N) 的工厂化循环水养殖车间。豹纹鳃棘鲈在循环水养殖车间中的流速一般控制在约0.1 m·s−1，海区养殖的流速介于0.2～0.8 m·s−1，在主要海区养殖的平均流速约为0.4 m·s−1。因此，本研究将规格一致的豹纹鳃棘鲈幼鱼分别在正常流速 (0.1 m·s−1,Low flow velocity, LFV) 和高流速 (0.4 m·s−1, High flow velocity, HFV) 下养殖5个月，在LFV和HFV组中分别选取相同规格的实验鱼各3尾，取肝脏组织，液氮速冻后−80 ℃保存备用。

1.2 组织学观察

将实验鱼麻醉后沿排泄孔前部剪开，打开腹腔，切取肝脏组织块，用4%多聚甲醛液固定，常规脱水，石蜡包埋，Leica RM 2135型切片机切片，厚度为5 μm，常规H.E染色和油红O染色，中性树胶封片，显微镜下观察并摄影。

1.3 文库构建和测序

取冷冻的豹纹鳃棘鲈实验用鱼肝脏组织，按照TRIzol Reagent操作流程提取总RNA，利用DNaseI (RNaseFree) 去除RNA酶污染和残留的基因组DNA。分别用1%的琼脂糖凝胶电泳、Nanodrop和Qubit RNA检测提取的总RNA的完整性、纯度和浓度，判定其是否符合下一步文库构建的质量要求。分别取高、低流速组中实验鱼各1.0 μg的RNA构建两个RNA样品池，在两个样品池中各取1.5 μg总RNA构建cDNA文库。利用AMPure XP beads纯化cDNA，再进行末端修复、加A尾并连接测序接头，经片段大小选择后通过PCR富集得到cDNA文库，使用Q-PCR方法对已构建文库的有效浓度 (>2 nmol) 进行准确定量。基于北京百迈客生物科技有限公司的 Illumina NovaSeq 高通量测序平台，以双末端测序模式完成库检合格cDNA文库的转录组测序。

1.4 测序数据质控和组装

为确保后续转录组分析的准确性，需要对原始下机数据 (Raw Data) 进行过滤，包括去除含有接头、N (无法确认的碱基信息) 比例大于10%以及低质量 (Qphred≤10的碱基数占整条read的50%以上) 的reads，经一系列严格质控步骤获得高质量测序数据 (Clean Data)。利用HISAT2软件[25]将测序reads高效比对到豹纹鳃棘鲈参考基因组 (https://db.cngb.org/search/project/CNP0000859)，用StringTie软件[26]将比对到基因组上的reads组装为转录本，实现基因或转录本的表达估计。

1.5 基因功能注释

为获取全面的基因功能注释信息，比对到参考基因组上的基因和转录本以参考基因组注释为准，而未被注释的转录区则被发掘为该物种的新基因和新转录本。使用BLAST软件[27]将新基因与Swiss-Prot、GO、COG、KOG、Pfam和KEGG数据库分别进行序列比对，预测新基因的氨基酸序列，通过HMMER软件[28]与Pfam数据库比对获得注释信息。

1.6 差异基因GO和KEGG富集

基于R语言GOseq程序包的Wallenius noncentral hyper-geometric数学模型，对筛选的差异表达基因进行GO富集分析[29]。以KEGG数据库中Pathway为单元，使用KOBAS 2.0软件[30]的超几何检验方法检测差异表达基因中显著性富集的Pathway。

1.7 差异表达基因筛选

将组织样品中比对到参考基因组的reads数目和转录本长度进行归一化，采用FPKM (Fragments Per Kilobase of transcript per Million fragments mapped) 法[31]计算转录本或基因表达量。利用DESeq2软件[32]对各组样本中的基因表达进行差异显著性分析，获取高低流速组之间的差异表达基因集，具体标准为Fold Change≥2且FDR (False Discovery Rate)<0.05。

1.8 实时荧光定量 PCR (qPCR) 验证

利用qPCR对14个差异表达基因在不同流速下豹纹鳃棘鲈肝脏中的表达水平进行定量验证 (引物见表1)。RT-qPCR的反应体系共15 μL，其中包含 7.5 μL 2×Power Green qPCR Mix (东盛生物，广州)，0.6 μL 上下游引物 (10 μmol·L−1)，1.5 μL cDNA，4.8 μL ddH2O。对每个样品进行3次技术重复，并将反应置于Roche LightCycler 96 (罗氏，瑞士) 定量仪上运行。反应程序为：95 ℃ 3 min；95 ℃15 s，60 ℃ 15 s，72 ℃ 30 s，40 个循环，CT通过溶解曲线分析PCR产物的特异性。以β-actin为内参基因 (F: CACCACAGCCGAGAGGGA；R: TCTGGGCAACGGAACCTCT)，采用 2−ΔΔCT方法[33]分析定量基因相对表达水平。使用R语言对结果进行统计，P<0.01表示差异具有统计学意义。

表1 本研究所用引物及其序列Table 1 Sequences of primers used in this study

1.9 SNP检测

采用GATK软件对HISAT2比对结果进行单核苷酸变异 (Single Nucleotide Polymorphisms, SNP)位点检测。根据SNP位点在参考基因组上的位置以及参考基因组上的基因位置信息，利用SnpEff软件对检测到的SNPs进行注释和变异影响预测。SNP位点按照碱基替换方式不同分为转换和颠换2种类型，而按照等位数目差异则可以划分为纯合型和杂合型。

2 结果

2.1 不同流速下肝脏组织观察

实验鱼肝脏H.E染色显示 (图1) ，LFV和HFV两组鱼肝脏组织的形态差异不大，其最外层均有由一层扁平上皮和结缔组织构成的浆膜层组成，浆膜层上均附有大的脂肪颗粒，肝细胞呈不规则椭圆形排列，细胞核大而圆，一般为一个，胞质丰富，多呈嗜酸性，胞质内含大量脂肪，占据肝细胞的大部分空间，细胞核被空泡状脂滴挤到一侧。处于肝脏结缔组织分支之间的部分为肝实质，其主要成分是肝细胞和窦状隙。彼此相连的肝细胞排列组成放射状的细胞索，而细胞索又互相连接形成网状结构，网状间隙里则是由内皮细胞、星形的枯否氏细胞和脂肪细胞组成的肝血窦。实验鱼肝脏油红O染色显示 (图2) ，LFV和HFV两组实验鱼肝脏组织中脂肪含量差异明显，LFV实验组肝脏中的脂肪含量明显高于HFV实验组。

图1 不同流速下豹纹鳃棘鲈肝脏组织 (H.E染色)Fig.1 Liver tissues (H.E staining) of P.leopardus with different flow velocities

图2 不同流速下豹纹鳃棘鲈肝脏组织 (油红O染色)Fig.2 Liver tissues (Oil red O staining) of P.leopardus with different flow velocities

2.2 测序数据统计及质量分析

本研究分别构建了差异流速组中各3尾豹纹鳃棘鲈实验用鱼的cDNA文库，并进行RNA-Seq分析。Raw Reads经严格质量控制后获得50.20 Gb 的Clean reads，且各样品的Q30碱基占比均≥94.80%，GC含量介于50.44%～51.54%，具体统计数据见表2。测序数据整体质量可靠性高，能够满足后续分析的要求。

表2 豹纹鳃棘鲈肝脏文库测序数据量统计Table 2 Output statistics of P.leopardus liver library sequencing

2.3 转录组数据与参考基因组的比对

利用HISAT2软件将Clean reads与豹纹鳃棘鲈参考基因组进行序列比对，获取reads在参考基因组或是相应基因中的位置信息。从比对效率即比对到参考基因组上的reads占Clean reads百分比的统计结果 (表2) 来看，各样品的reads与参考基因组的比对效率在92.77%～93.78%，且唯一比对位置的reads占比≥85.54%，基于参考基因组的比对和组装能够满足相应分析的需求。

2.4 差异基因GO功能分类和KEGG富集分析

GO (Gene Ontology) 数据库是一个结构化的标准生物学注释系统，能够对各物种进行基因和蛋白功能的限定和描述。将所有与参考基因组比对获得的基因与GO数据库进行注释，共有12146个基因被归属到56个类别，其中高低流速组差异表达基因 (DGEs) 有1124个(图3)。GO分类结果显示，生物学过程共涉及23个功能分类，参与细胞过程(Cellular process, 530)、单组织过程 (Single-organism process, 460) 和代谢过程 (Metabolic process,430) 的差异基因最多；细胞成分涉及16个功能分类，细胞 (Cell, 447)、细胞组分 (Cell part, 440) 和膜(Membrane, 386) 的差异基因最多；分子功能涉及16个功能分类，结合 (Binding, 581)、催化活性(Catalytic activity, 442) 和转运活性 (Transporter activity, 72) 差异基因最多。

图3 豹纹鳃棘鲈基因GO分类图Fig.3 Gene ontology (GO) assignment of genes of P.leopardus

将基因序列与KEGG数据库中的通路进行比对注释，结果显示573个差异表达基因参与了154条KEGG通路 (图4)，其中最显著的PPAR信号通路内包括差异表达基因20个，其次显著的氨基酸生物合成通路内包括差异表达基因21个，涉及其中的基因富集效果都极显著 (P<0.01)。此外，包含精氨酸合成、初级胆汁酸合成、碳代谢、甘油酯代谢和脂肪细胞因子信号通路在内的7个最为显著的通路均直接或间接涉及脂肪合成代谢及沉积、相关酶的活性、基因表达以及信号转导的调节。为了明确豹纹筛棘鲈不同流速下肝脏组织之间的差异性，本研究制备了肝脏组织切片进行组织学观察 (图1和图2)，其观察结果与RNA-Seq分析中筛选到最为显著的差异基因富集通路有相当高的一致性。

图4 豹纹鳃棘鲈差异表达基因KEGG富集图Fig.4 KEGG enrichment of differentially expressed genes of P.leopardus

2.5 差异基因表达分析

对6个转录测序样品比对到参考基因组上的基因采用FPKM法定量，利用DESeq2进行样品组间的差异表达分析，获取高低流速组两个生物学条件之间的差异表达基因 (Different expression gene,DEG) 集。在LFV-HFV的DEG分析过程中，共筛选到DEGs 1 977个，999个基因在低流速组鱼肝脏上调，978个基因在高速组鱼肝脏上调。在不同流速下豹纹鳃棘鲈肝脏差异表达基因中筛选显著富集通路包含的部分基因，如cpt1a、cyp7b1、irs2、socs3、adpgk、kmt5c和acsbg2等 (表3)。本研究选择RNA-Seq分析中的14个差异表达基因，通过RT-qPCR技术对RNA-Seq结果进行验证，结果显示挑选的基因表达模式与RNA-Seq分析一致，表明RNA-Seq分析结果可信 (图5)。

表3 豹纹鳃棘鲈不同流速显著差异通路中部分候选基因表达模式Table 3 Expression pattern of partial candidate genes of P.leopardus in significant pathways with different flow velocity

图5 RT-PCR验证RNA-Seq结果Fig.5 Validation of RNA-Seq data by using RT-qPCR

2.6 SNP注释分析

单核苷酸多态性 (Single nucleotide polymorphisms, SNP) 是指在基因组上由单个核苷酸变异形成的遗传标记，变异位点丰富，是一种重要的辅助育种标记。对豹纹鳃棘鲈肝脏转录组数据进行SNP的筛选和注释将为豹纹鳃棘鲈分子标记的开发提供基础数据，有助于其目标性状的分子遗传机制解析和新品种选育。对各样品中注释的SNP位点数目、转换类型和颠换类型占比进行统计 (表4) 发现，各项统计指标在样品间差距较小，表明基于本研究数据的SNP信息挖掘结果可靠性高，具有进一步分析的意义。SNP突变类型统计结果 (图6) 显示，主要的突变类型为A(G)/G(A)和C(T)/T(C)，占据所有突变类型的65%以上。

表4 SNP位点统计表Table 4 Output statistics of SNP loci

图6 SNP突变类型统计分布图Fig.6 Statistical distribution of type of SNP mutations

3 讨论

3.1 流速变化对鱼类生长的影响

水流速度是鱼类生活环境中最为重要且复杂的生态因子之一，主要通过刺激鱼类感觉器官来影响其新陈代谢和生理生态反应，最终影响鱼类的生长发育、肉纤维组成、体营养成分和免疫功能[34]。娄宇栋等[35]开展了美国红鱼 (Sciaenops ocellatus) 的流速胁迫肝脏转录组实验，筛选出大量参与各类重要生物学过程和代谢途径的显著性差异表达基因，对美国红鱼应对流速胁迫的分子响应机制有了新的认识。而拉氏鱥 (Rhynchocypris lagowskii) 幼鱼生长、非特异性免疫能力及脂肪酸组成受流速差异影响的研究[36]也表明，相较于低流速，高流速会增加鱼体能量的消耗，多不饱和脂肪酸含量更高。目前，豹纹鳃棘鲈在工厂化流水养殖中的流速多控制在0.1 m·s−1左右，而海区网箱养殖的平均流速在0.4 m·s−1左右。对豹纹鳃棘鲈养殖过程的长期观测发现，流速差异并不会导致体长或体质量等生长性状表型产生直接变化。然而，肝脏油红O染色结果 (图2) 显示，光学显微镜下LFV组 (图2-d) 的脂滴细胞定位区域和密度明显高于HFV组 (图2-c)，这与RNA-Seq分析筛选到大量参与脂肪合成代谢相关显著通路的结果相对应。基于上述豹纹鳃棘鲈不同流速环境中的养殖情况，本研究开展了高、低差异流速的养殖实验，通过肝脏组织转录组测序分析差异流速组个体之间的基因表达量变化，探究豹纹鳃棘鲈对流速变化适应性的分子调节机制。

3.2 豹纹鳃棘鲈肝脏RNA-Seq分析

鱼类肝脏是维系其生理机能最重要的器官之一，在鱼脑 (各种神经内分泌因子)-脑垂体 (分泌生长激素)-肝脏 (产生类胰岛素生长因子) 轴的调控过程中发挥重要作用[37]。而胆固醇合成和脂肪酸氧化也主要在肝脏中进行，且经由肠道吸收的营养物质能够被肝脏转换，从而促进机体生长发育和物质代谢[38-39]。由于肝脏在动物机体活动中发挥着重要作用，基于肝脏组织的RNA-Seq被广泛应用于各种目标性状和环境对机体影响的分子生理机制研究[40-42]。本研究利用 0.1 和 0.4 m·s−1两种流速下豹纹鳃棘鲈的肝脏转录组数据进行差异基因表达分析，LFV-HFV上调基因富集最显著的生物学过程为心脏发育和微管相关进程，而LFV-HFV下调基因富集最显著的生物学过程则为基因表达的节律调控和三羧酸循环。鱼类心脏是形成较早且发挥重要功能的器官之一，其发育过程会直接影响血管的发生和生成[43]。微管存在聚合和解聚的动力学特性，在细胞形态维持、信号转导及物质输送等过程中具有重要作用[44]。参与表达节律调控的基因在转录过程中和转录后的调控网络中发挥重要作用，通过节律调控提高生物机体对环境变化的适应性[45]。三羧酸循环是糖类、脂类和氨基酸的最终代谢通路，是它们之间相互联系的枢纽。差异表达基因富集最显著的PPRA信号通路则在调控皮下脂肪与肌内脂肪差异沉积等方面发挥作用[46]，结合组织切片结果 (图2) ，本研究表明流速差异可能导致豹纹鳃棘鲈脂肪的代谢形成产生显著变化。此外，本研究还检测到了446928个潜在的SNP位点。这些差异表达基因和SNP位点的筛选为豹纹鳃棘鲈基因资源挖掘提供了基础数据。

3.3 豹纹鳃棘鲈差异流速适应性相关基因

为探究豹纹筛棘鲈不同流速环境下适应性分子机理，本研究在RNA-Seq分析中通过KEGG分析，从豹纹筛棘鲈肝脏转录组中筛选出多个可能涉及差异流速适应性的显著通路 (图5) 和参与其中的候选基因 (表3)。在上述显著通路和候选基因中，氨基酸生物合成通路在控制脂肪细胞的分化过程中发挥重要作用，氨基酸调控脂肪代谢的机制研究结果[47]发现，亮氨酸、异亮氨酸、缬氨酸和甲硫氨酸这4种必需氨基酸是脂肪合成代谢过程中起重要调控作用的营养调控因子。精氨酸合成通路对脂类代谢具有调节作用[48]，精氨酸经过一氧化氮合酶的催化生成NO，继而促进葡萄糖和脂肪酸的氧化分解和脂肪细胞脂解作用，抑制脂肪组织沉积。初级胆汁酸合成通路主要涉及胆酸、鹅脱氧胆酸及相应结合型胆汁酸，初级胆汁酸在一定条件下可转化为次级胆汁酸，胆汁酸合成是胆固醇和脂代谢的主要途径，其主要功能是促进脂质和脂溶性维生素的消化和吸收[49-50]。甘油酯代谢过程则主要涉及甘油磷脂代谢、磷脂酸合成、脂肪酸的合成与降解。Hao等[51]对一种特殊类型的短链脂肪酸 (丁酸) 进行深入研究发现，cpt1a活性经丁酸盐可增强进而促进脂肪酸氧化，cpt1更被认为是影响脂肪酸氧化的决定因素。山羊cpt1a基因克隆表达分析则表明，cpt1a的基因表达量与各肌肉组织肌内脂肪含量均显著正相关，其可能在脂肪沉积过程中发挥重要调控作用[52]。cyp7b1在胆汁酸合成过程中发挥关键作用[53]，作为胆汁的重要成分，胆汁酸在脂肪代谢中起重要作用。irs2是胰岛素促进肝糖原合成和一致葡萄糖输出的重要中间因子[54]，而miR-33a在绵羊前体脂肪细胞分化研究则表明miR-33a靶向irs2抑制绵羊体脂肪细胞分化和脂滴沉积[54]。kmt5c是一种赖氨酸甲基转移酶，Zhao等[55]在棕色脂肪前体细胞分化过程中敲除kmt5c发现，kmt5c会抑制脂肪细胞产热功能，且kmt5c可以通过调控ucp1基因表达控制白色脂肪棕色化。D'Andre等[56]对鸡脂肪发育基因的鉴定和描述研究发现，acsbg2基因与腹脂重和腹脂率显著相关，在脂肪形成过程中起重要作用。综上，笔者推测豹纹鳃棘鲈肝脏在差异流速养殖环境中可能通过调节脂类合成代谢来调控流速变化适应性的作用机制。

4 结论

本文在转录水平初步探讨了差异流速情况下豹纹鳃棘鲈个体之间相关基因表达水平的变化情况，通过肝脏组织学观察、量化基因表达水平以及差异表达基因分析，对豹纹鳃棘鲈流速变化适应性机制开展了探索研究。本研究不仅为豹纹鳃棘鲈经济性状的分子遗传机制解析和新品种选育提供了理论依据，也为其进一步的分子标记开发、良种选育和养殖海区选址奠定了基础。