碱度长期胁迫对花鲈肾组织转录组的影响❋

王龙斌， 泮淼军， 王明阳， 王润泽， 李 丽，2， 董双林，2， 李卫东， 田相利 ，2❋❋

(1. 海水养殖教育部重点实验室(中国海洋大学)， 山东 青岛 266003； 2. 青岛海洋科学与技术试点国家实验室， 山东 青岛 266237； 3. 唐山海都水产食品有限公司， 河北 唐山 063000)

新疆是中国盐碱化土壤分布面积较广、土壤积盐较重的地区，其中南疆盐碱水具有碱度高、离子比例不平衡、水质变动大等特点，合理有效地开发利用盐碱地已成为该地区紧迫而艰巨的任务[1]。碱性水通常具有较高的碱度和pH，严重影响水生动物的正常生理过程，例如直接抑制鱼类的氨排泄并增加CO2排泄[2]。尽管如此，仍有一些鱼类可以在高碱的环境中生存[3]，如肯尼亚马加迪湖的马加迪罗非鱼(Alcolapiagrahami)[4]、美国金字塔湖的美洲鲑(Oncorhynchusclarkihenshawi)[5]、我国达里湖的瓦氏雅罗鱼(Leuciscuswaleckii)[6]以及青海湖的裸鲤(Gymnocyprisprzewa-lskii)[7]等。这些鱼类已经进化出特殊的生理机制以适应碱性环境[7-9]。

转录组分析是解释功能基因组学元素和揭示细胞与组织分子机制的有力工具之一[10]。近年来，RNA测序技术(RNA-seq)已广泛应用于基因表达谱的研究，包括水生动物响应盐碱胁迫的通路研究。如马加迪罗非鱼[11]、瓦氏雅罗鱼[6,12]、尼罗罗非鱼(Oreochromisniloticus)[13]、奥尼罗非鱼(Oreochromismossambicusfemale×O.urolepishornorummale)[14]等，这些研究提供了大量与碱度胁迫相关的候选通路和基因。

花鲈(Lateolabraxmaculatus)属鲈形目(Perciformes)鮨科(Serranidae)花鲈属(Latelabrax)，具有生长快、病害少、对盐度适应广等特点。近年来，关于花鲈对盐度的适应能力[15-16]以及耐受机制[17-18]已得到了很好的研究。花鲈对碱度有较高的耐受性，在盐度为10 时，碱度96 h半致死浓度为98.51 mmol/L[19]。但目前国内外对于花鲈耐碱机制的研究很少，仅见于碱度胁迫下Slc4基因家族[20]和热休克蛋白 (Hsp)超家族[21]的相关研究。

在本研究中，我们利用RNA-seq技术探究了不同碱度下花鲈肾组织的基因表达情况，对其响应碱度胁迫相关的通路及基因进行了综合分析，为了解花鲈适应碱性环境的分子生物学机制提供有价值的信息。

1 材料与方法

1.1 样本来源及实验设计

实验所用花鲈是由唐山由海都水产食品有限公司(河北，唐山)提供的人工繁育苗种。暂养一周后随机选取135尾健康、个体大小均匀的花鲈进行养殖实验。

本研究参照南疆阿拉尔地区盐碱水的盐度水平(10左右)，通过添加NaHCO3，设计了不同的碱度梯度。将天然海水与淡水混合成盐度为10的海水，再添加NaHCO3(分析纯)，全部溶解并稳定24 h后使用。根据前期急性胁迫的结果，设置了0 mmol/L(con组)、15 mmol/L(AW15组)和30 mmol/L(AW30组)的碱度梯度，具体实测碱度如表1所示。每组设置3个重复，每个重复15尾鱼，实验鱼体长(15.21±1.56) cm，体质量(60.23±5.54) g，实验时间为2020年9月15日至10月15日，实验周期为 30 d。

表1 碳酸盐碱度设置

实验所用水体的体积为200 L，盐度为10，水温控制在(20±2) ℃，溶氧(8.4±0.1) mg/L，pH=8.5±0.3。实验期间，每天8:00和18:00各投喂花鲈专用颗粒饲料1次，投喂量约为鱼体质量3%左右，每天换水50%，并检测碳酸盐碱度和pH。pH以DELTA320 型精密pH计测定，碳酸盐碱度以酸碱滴定法(GB/T 9736—2008)测定[22]。

1.2 生长指标的测定

实验开始和结束时测量花鲈的体质量，按以下公式计算体质量特定生长率(SGR)：

SGR=(lnW2-lnW1)/t×100%。

式中：W1为初始体质量(g)；W2为终末体质量；t为实验时间(30 d)。

实验数据以平均值±标准差表示,实验结果采用SPSS 19.0软件对实验数据进行统计分析。先对数据做单因素方差分析，组间若有显著性差异，再用Duncan’s多重比较分析，P<0.05时认为有显著差异。

1.3 RNA提取、建库和测序

实验结束时，每个处理随机取3尾鱼(共9尾)使用MS-222麻醉，解剖以获得肾组织，在液氮中速冻并储存在-80 ℃用于RNA分离。使用Invitrogen公司的TRIzol®Reagent试剂盒提取组织RNA，然后使用Roche公司的DNaseI去除基因组DNA污染，采用Nanodrop 2000、Qubit 2.0和Aglient 2100方法检测RNA样品纯度、浓度和完整性，检验合格的样品送上海美吉生物医药科技有限公司，采用Illumina TruseqTM RNA sample prep Kit方法进行文库构建，基于Illumina Novaseq 6000 测序平台以Paired-end 150 bp双末端测序模式分别对9个样本进行转录组测序。

1.4 序列比对

利用fastp对测序获得的原始数据进行质量评估，然后用SeqPrep(https://github.com/jstjohn/SeqPrep)和Sickle(https://github.com/najoshi/sickle)去除测序接头序列、低质量读段、N率较高序列及长度过短序列。通过TopHat2[23]将质控后的数据与花鲈参考基因组(http://gigadb.org/dataset/100458)进行比对，获得用于后续分析的mapped data。

1.5 表达量差异分析

使用RSEM V1.2.15[24]对mapped data进行计数，然后通过FPKM进行数据标准化。利用DESeq2 R包[25]筛选差异表达基因，将P-adjust<0.05且|log2FC|≥1作为筛选标准。将花鲈肾组织在三个处理下两两对比得到的基因集(con_vs_AW15，con_vs_AW30，AW15_vs_AW30)中进行比较分析。

1.6 GO富集和KEGG通路富集分析

使用软件Goatools(https: //github.com/tanghaibao/GOatools)对差异基因进行GO富集分析，采用Bonferroni方法对P值进行校正，当经过校正的P值(FDR)小于0.05 时，认为此GO功能存在显著富集。利用KOBAS软件(v2.0.12)对差异基因进行KEGG通路富集分析，以判定差异基因主要影响的生物学功能和通路。

1.7 qRT-PCR验证

选择9个差异表达基因，使用Primer Premier 6.0软件[26]设计引物(见表2)，并由生工生物工程(上海)股份有限公司合成，进行实时荧光定量PCR验证。使用诺唯赞公司HiScript®III All-in-one RT SuperMix Perfect for qPCR试剂盒将总RNA逆转录合成cDNA，再以逆转录产物为模板、以18S ribosomal RNA为内参基因，使用诺唯赞公司的ChamQTMSYBR®Color qPCR Master Mix试剂盒，在ABI PRISM 7500型荧光定量PCR仪上进行表达量分析。总体系20 μL的反应体系如下：稀释10倍的cDNA模板2.0 μL，2×ChamQ SYBR Color qPCR Master Mix 10 μL，正反向引物各0.4 μL，RNase-Free ddH2O 7.2 μL。反应程序：95 ℃ 30 s，95 ℃ 10 s，60 ℃ 30 s，40个循环。最后在60～95 ℃制作熔解曲线。每组样品的每个基因均重复3次，采用2-ΔΔCt法[27]分析基因的相对表达量。使用SPSS 19.0软件计算RNA-seq结果和qRT-PCR结果之间的相关性系数。

表2 RT-qPCR引物信息

2 结果

2.1 生长指标

养殖实验结果见表3。由表3可以看出，碱度严重影响花鲈生长，随着碱度升高，花鲈的终末体质量、特定生长率显著降低。在15 mmol/L的中碱度下，花鲈的生长速率显著降低。而在30 mmol/L的高碱度下，花鲈的终末体质量显著低于初始体质量。

表3 不同碱度对花鲈生长性能的影响

2.2 测序数据统计与参考序列比对分析

所有样品的RNA浓度均高于200 ng/μL，RNA总量8 μg，均符合建库要求。对9个转录组测序结果进行分析可知，所有样品获得的Clean reads数量在46 266 138～55 723 922之间，利用SeqPrep软件过滤掉不合格的序列后，所有样品得到的有效读数均占原始读数总数量的97%以上(见表4)。表4中，Q20、Q30分别指测序质量在99%和99.9%以上的碱基占总碱基的百分比，在得到的9个转录组中，Q20均在96%以上，Q30均在91%以上。利用HISAT2软件将有效读数比对到参考基因组，结果发现有79.71%～82.12%读数成功映射到参考基因组。这些结果表明，花鲈肾组织各处理组测序量大，获得的有效读数数量多，转录组质量高，符合进一步生物信息学分析的要求。

表4 花鲈肾组织转录组高通量测序基本数据统计

2.3 表达量差异分析

使用基于负二项分布的DESeq2软件对raw counts进行分析，基于一定的标准化处理和筛选条件获得比较组间表达差异的基因，参数如下：P-adjust<0.05 & |log2FC| ≥1。分析得到的差异表达基因结果表明，在不同的碱度环境下，花鲈发生显著差异表达基因的数量也不同(见表5)。与con组相比，AW15组有142个基因上调，104个基因下调；AW30组有2 037个基因上调，1 882个基因下调；与AW15相比，AW30组有1 318个基因上调，964个基因下调。对不同碱度下差异表达基因进行综合分析，结果如图1所示，与con组相比，AW15组与AW30组分别有68和2 189个特有的差异表达基因；与AW15组相比，AW30组有577个特有的差异表达基因；3个基因集共有的差异表达基因为49个。

表5 花鲈肾组织差异表达基因统计结果

图1 花鲈肾组织con_vs_AW15组、con_vs_AW30组和AW15_vs_AW30组显著差异表达基因的维恩图

2.4 GO和KEGG富集分析

为确定可能参与碱度胁迫的生物过程或通路，分别对3个基因集(con_vs_AW15，con_vs_AW30，AW15_vs_AW30)中的显著差异基因进行GO和KEGG富集分析。

GO富集结果显示，差异基因被富集到生物过程、细胞组分和分子功能3个类别体系中。其中最丰富的类别是生物过程，其次是分子功能和细胞组分。图2显示了肾中3个基因集前20个显著GO术语。可以看出，con_vs_AW15中气体运输、氧结合、血红蛋白复合物富集程度最高，con_vs_AW30中水解酶活性的正调控、GTPase活性的正调控、GTPase活性的调节富集程度最高，而AW15_vs_AW30中DNA复制启动、蛋白质磷酸化、丝氨酸家族氨基酸代谢过程富集程度最高。其中，GTPase活性的正调控、水解酶活性的正调控、水解酶活性的调节、GTPase活性的调节、丝氨酸家族氨基酸代谢过程、α-氨基酸代谢过程、细胞氨基酸代谢过程等在con_vs_AW30和AW15_vs_AW30中均显著富集。

图2 花鲈肾组织con_cs_AW15(A)、con_vs_AW30(B)和AW15_vs_AW30(C)组差异表达基因的GO富集分析

KEGG分析表明，在3个基因集中，分别有254、344、343条通路与碱度胁迫有关。图3显示了肾con_vs_AW30和AW15_vs_AW30两个基因集前20个显著KEGG通路。可以看出，con_vs_AW30中细胞周期、叶酸一碳库、过氧化物酶体显著富集，AW15_vs_AW30中DNA复制、癌症中的蛋白多糖、轴突导向显著富集。其中，细胞周期、DNA复制、叶酸一碳库、上皮细胞的细菌侵袭、趋化因子信号通路、轴突导向、肌动蛋白细胞骨架的调节、错配修复、PPAR信号通路和癌症中的蛋白多糖在con_vs_AW30和AW15_vs_AW30均显著富集。

图3 花鲈肾组织con_vs_AW30(A)组和AW15_vs_AW30(B)组差异基因KEGG散点图

2.5 花鲈肾组织碱度胁迫响应的关键基因分析

图1表明con_vs_AW15、con_vs_AW30、AW15_vs_AW30这3个基因集中有49个共有差异基因，这些基因是花鲈肾组织碱度胁迫响应的关键基因，其中大多随着碱度升高而降低(见表6)。

表6 花鲈肾组织碱度胁迫响应的关键基因

2.6 qRT-PCR验证

为了确认Illumina RNA-seq结果的准确性和可重复性，选择9个代表性基因，通过实时定量PCR(qRT-PCR)验证不同处理组基因的表达水平。进一步将所选基因的相对表达量与RNA-seq测序结果进行了相关性分析(见图4)。可以看出，虽然所选基因在表达变化幅度上存在一定差异，但qRT-PCR检测结果与RNA-seq测序结果中基因的表达趋势基本一致(R2=0.959 4)，这表明了Illumina 测序获得的转录组数据以及生物信息学分析的可靠性。

图4 花鲈肾组织RNA-seq产生的9个差异表达基因进行qRT-PCR验证

3 讨论

作为水生动物的重要器官，鱼类肾脏主要具有维持渗透压和离子平衡的功能，在鱼类的渗透调节过程中起着重要作用[28]。研究发现，高碱度会使鱼类肾组织发生不同程度的退化，而碱度越高肾萎缩程度越大[29]，肾功能被明显抑制[30]。在本研究中，我们首次研究了不同碱度胁迫下花鲈肾组织的有参转录组数据，探讨了碱度对花鲈肾组织的影响，获得了9 602个差异基因，富集分析后得到大量碱度胁迫相关通路以及对碱度敏感的基因，为花鲈耐碱机制研究提供了新视角。

3.1 GO富集分析

GO数据库是一个国际化的基因功能标准分类体系，通过GO功能显著性富集分析能确定差异表达基因行使的主要生物学功能[31]。本研究对花鲈肾脏转录组数据GO富集分析表明，相较于对照组，碱度15 mmol/L处理组差异基因GO术语显著富集于氧气和血红蛋白相关术语，说明中碱度可能主要影响肾组织氧气的运输。相类似，何强等[32]发现，随着碱度上升，瓦氏雅罗鱼的耗氧率呈先下降后上升的趋势，Lykkeboe等[33]发现，罗非鱼(Tilapiagrahami)血红蛋白的氧亲和力随着pH升高而降低，这些研究结果与本研究基本相一致。但高碱度对鱼类的影响机制与中碱度可能有所差异。例如，Xu等[6]发现高碱度下瓦氏罗雅鱼肾组织GO术语显著富集于转运蛋白活性和酶调节活性，高碱度导致肾组织跨膜转运体活性明显增加。在本研究中，高碱度(30 mmol/L)对花鲈肾的影响主要集中在蛋白质水解(水解酶活性的正调控、催化活性的正调控)和氨基酸代谢过程(丝氨酸家族氨基酸代谢过程、α-氨基酸代谢过程、细胞氨基酸代谢过程等)中，这可能表明肾组织在高碱度下能量消耗巨大，能量需求的增加导致肾脏蛋白质水解和氨基酸代谢过程活跃。

3.2 KEGG富集分析

KEGG是一个集成了基因组、化学和系统功能信息的生物系统数据库，KEGG提供了一个参考知识库，通过通路映射过程将基因组与生命联系起来[34]。在本研究中，不同碱度下花鲈肾脏KEGG富集分析表明，相较于对照组，中碱度的15 mmol/L处理组对花鲈肾组织的影响主要集中在氨基酸合成与代谢相关通路(苯丙氨酸、酪氨酸和色氨酸生物合成、氨基酸生物合成、半胱氨酸和蛋氨酸代谢等)，而高碱度(30 mmol/L)对花鲈的影响则更加广泛，主要包括细胞有丝分裂(细胞周期、DNA复制)、免疫与疾病相关通路(例如上皮细胞的细菌侵袭、癌症中的通路)以及代谢(叶酸一碳库)。从现有的研究看，Cheng等[35]发现，高碱度会使尼罗罗非鱼肌肉中不饱和脂肪酸含量升高，Zhao等[36]发现不饱和脂肪酸的生物合成通路在尼罗罗非鱼肝组织AW40_vs_AW60中显著富集。除了脂肪酸，氨基酸代谢在渗透调节中也起着关键作用[37]，例如丙氨酸、天冬氨酸和谷氨酸代谢[38-40]以及缬氨酸、亮氨酸和异亮氨酸代谢等[41]。这些研究结果表明，脂肪酸与氨基酸的合成与代谢可能是碱度胁迫影响的关键通路之一。而在碱度胁迫下，奥尼罗非鱼鳃组织富集通路与免疫、错配修复、核苷酸切除修复、DNA复制以及细胞周期显著相关[14]，这与本研究结果基本一致。因此，推测高碱度可能损伤花鲈肾细胞，抑制细胞有丝分裂。值得注意的是，本研究发现花鲈肾中轴突导向和肌动蛋白细胞骨架调节相关基因在高碱度下显著富集，这在其他鱼类相关研究中从未发现。已有研究表明，在肾病中肾小球轴突导向通路显著富集，肌动蛋白细胞骨架失调随着足突消失和粘连丧失而从基底膜脱离而动态重组[42]，说明这两条通路的富集可能与肾组织疾病相关。相类似，Chen等[41]发现癌症中的通路在盐度和碱度胁迫下叶尔羌高原鳅(Triplophysayarkandensis)鳃中均富集，这也与本研究结果相一致。本研究发现花鲈肾脏癌症中的蛋白多糖在con_vs_AW30和AW15_vs_AW30两个基因集中均显著富集，而该通路的富集被认为是肾细胞癌的标志[43]。上述结果表明，中、高碱度均会不同程度地影响花鲈免疫能力，增加了花鲈患病风险。

3.3 花鲈肾组织碱度胁迫响应关键基因分析

环境胁迫会导致组织细胞氧化应激、蛋白质损伤、DNA双链断裂和氧化碱修饰[44-45]。本研究发现，随碱度升高花鲈肾组织生长相关基因表达显著下调，例如fyn、ccna2、cks1、arhgap19、melk、cenpo、pbk、map3k15、rrm1等。fyn是Src家族的成员，fyn的抑制与细胞生长减少有关[46]。ccna2基因编码细胞周期蛋白，有细胞周期调节剂的作用[47]。cks1在调节细胞周期中发挥作用，Cks1蛋白的缺失可导致细胞周期停滞在G2期，有丝分裂受阻[48]。arhgap19的低表达水平则会影响早期有丝分裂期间的细胞形状，导致严重的染色体分离缺陷[49]。而melk基因在细胞周期控制方面发挥着重要作用，melk的表达水平与影响细胞周期进程[50]。cenpo基因编码间期着丝粒复合体，其编码的蛋白质在整个细胞周期中定位于着丝粒，是有丝分裂期间双极纺锤体组装、染色体分离和检查点信号传导所必需的[51]。pbk基因编码的丝氨酸-苏氨酸激酶在有丝分裂中活跃[52]。map3k15基因编码MAP3K蛋白，在细胞凋亡、应激反应和各种疾病中起重要作用，从鱼类相关研究看，三刺鱼(Gasterosteusaculeatus)肾脏中MAP3K15的表达与盐度呈显著负相关[53]。rrm1基因编码核糖核苷酸还原酶的大型催化亚基，对细胞周期S期的DNA复制以及多个DNA修复过程非常重要[54]。本研究发现上述相关基因在中高碱度下显著下调，表明碱度严重抑制了肾细胞的有丝分裂。因此，碱度作为一种应激源，可能抑制了花鲈肾细胞的有丝分裂，影响了细胞的生长与分裂，甚至可导致细胞凋亡。

一般认为，在盐度胁迫下水生生物为了维持体内平衡和渗透平衡，需要合成转运蛋白和酶，这一过程会消耗大量能量[55]。本研究发现随着碱度升高，花鲈肾组织能量代谢通路相关基因(acsl、gaphd)表达显著上调。酰基辅酶A合成酶(ACSL)是负责脂质代谢初始步骤的关键酶[56]，gapdh基因编码的3-磷酸甘油醛脱氢酶是能量代谢中的一种重要酶，在脊椎动物中具有多种细胞调节作用，增强糖酵解对于提供升高水平的 ROS 清除剂丙酮酸可能很重要[57]。我们推测随着碱度升高，花鲈肾组织离子转运与免疫调节会消耗大量能量，因此导致能量代谢相关基因表达上调。

不同类型的应激会导致鱼类免疫和抗病能力发生变化，但具体应激反应取决于应激源的强度及其持续时间。如果应激源是慢性和长期的，免疫反应显示出抑制作用，鱼类患病的几率就可能增加[58]。本研究发现，随着碱度升高，花鲈肾组织免疫相关通路基因(ceacam1、clec4m和b4galt5等)表达显著下调。ceacam1基因编码免疫球蛋白相关糖蛋白，在先天性和适应性免疫反应的调节中发挥作用[59]。clec4m基因编码的C型凝集素是固有免疫系统中的重要模式识别受体[60]。b4galt5基因编码β-1,4半乳糖基转移酶，参与多种免疫细胞表型修饰、募集和迁移[61]。白细胞介素 4 诱导1(il4i1)编码的分泌型 L-苯丙氨酸氧化酶是一种免疫调节剂，主要与机体免疫调节功能有关[62]。crp基因编码C-反应蛋白，参与胆固醇结合过程，功能与宿主防御相关[63]。已有研究报道显示，碱度会导致尼罗罗非鱼免疫基因表达下调[36]和大鳞鲃(Luciobarbuscapito)免疫抑制[64]。本研究结果表明，中高碱度导致花鲈肾组织免疫能力明显下降，可能导致患病风险增加。当然，碱度下调花鲈免疫反应的详细机制尚不清晰，需要进一步研究确认。

4 结语

本研究通过高通量测序探讨了花鲈肾组织转录组对碱度胁迫的应答，获得了大量碱度胁迫相关通路，并在花鲈肾组织中鉴定出49个持续变化基因。研究发现，中、高碱度对花鲈肾组织影响有所不同，其中，中碱度主要影响花鲈脂肪酸与氨基酸的合成与代谢，高碱度主要影响其有丝分裂过程及免疫能力。随碱度升高，花鲈生长受到抑制，能量代谢显著增强，机体免疫能力下降。本研究获得了大量花鲈肾组织与碱度胁迫相关的通路和基因，为理解水生动物碱度适应的分子机制提供了新见解，同时也为花鲈在南疆盐碱水域渔业中的潜在应用提供了科学依据。