张文馨 潘 霞 沈锡权 徐永健

(宁波大学海洋学院, 应用海洋生物技术教育部重点实验室, 宁波 315211)

海马是一种深受东南亚国家喜爱的传统中药[1]。当前, 因全球野生海马种群的栖息地遭到破坏而导致其资源数量骤减, 但海马的市场需求却与日俱增。因此, 海马养殖业极具潜力, 不仅可以满足市场需求，还能够保护其野生资源[2]。尽管多年来针对海马的繁殖和培育开展了许多研究工作, 但仍然存在诸多养殖技术难题[3], 这些问题主要由养殖环境的胁迫所致。

盐度是鱼类生长发育的重要环境因素之一, 盐度变化对海马的生长发育及生理生化等方面存在很大影响[4]。已有一些研究涉及其他海洋鱼类, 证明盐度变化会影响鱼类的存活和生长[5], 甚至影响鱼类的性别决定和早期发育[6]。而盐度胁迫会引起渗透胁迫, 导致鱼类血浆应激相关激素增多, 能量代谢受到影响, 电解质平衡被破坏, 进而影响生物体内的渗透压稳态和免疫稳态[7]。Lu等[8]的研究发现, 高盐胁迫使黄鳝(Nibea albiflora)肝胰腺抗氧化酶和非特异性免疫酶活性发生显著变化。有报道指出, 长期低盐胁迫显著降低了小龙虾(Cherax quadricarinatus)免疫相关基因(Toll基因和ProPO基因)的表达水平[9]。这些工作对于海马养殖条件的研究都有很好的借鉴。海马属于广盐性鱼类, 可以存活于10‰—32‰盐度的环境中[10], 但海马幼体对盐度变化的适应能力差, 盐度急剧变化会导致渗透胁迫,并增加渗透调节和其他几种生物过程所需的能量需求[11,12], 因此盐度波动对海马幼体的影响较大。为减轻和避免盐度变化带来的不利影响, 本研究以养殖大海马(Hippocampus kuda)幼体为研究对象,采用RNA-Seq技术对盐度胁迫下大海马幼体的肝脏组织进行了转录组测序, 分析其肝脏组织在高/低盐度胁迫下的转录表达情况, 揭示了大海马盐度应激的调控机制, 并从中挖掘响应盐度应激的候选基因, 作为海马养殖过程中的生物标记指示物。

1 材料与方法

1.1 试验材料

本实验用生物大海马幼体由宁波大学海马养殖研究实验室培育。挑选健康无病且活力良好的个体 [个体质量: (0.46±0.07) g, 体长: (5.78±0.42) cm(P＞0.05)], 于实验前置于循环水养殖系统中暂养1周。暂养期间使用盐度25‰、温度(25±0.5)℃、含氧量DO＞5 mg/L过经砂滤后的天然海水, 每日定时投喂 2 次, 摄食 2h 后吸污并换水50%。实验开始后, 各组盐度设置如下: 对照组盐度25‰(Control,CK)、高盐胁迫组31‰(HS-test)和低盐胁迫组17‰(LS-test), 各组水温均保持25℃。急性盐度胁迫 12h 后, 每组随机取 3 尾大海马的肝脏置于同一冻存管中制成一个混样以消除个体误差, 以此每组各设置3个生物学重复[13]。样品组织以液氮速冻,置于-80℃ 冰箱中保存。

1.2 RNA 提取、文库构建和测序

采用 RNAiso plus 试剂盒(大连宝生生物有限公司)提取肝脏总 RNA。采用Multiskan SkyHigh微孔板分光光度计(Thermo Fisher Scientific, 美国)检验RNA样品的浓度(ng/μL)和纯度(A260/A280值), 利用 1%的琼脂糖凝胶电泳检查 RNA 样品的完整性,采用 Bioanalyzer 2100 系统(Agilent Technologies,美国)评估 RNA 的浓度和完整性。取每组2 μg RNA样品, 采用 NEBNext UltraTMRNA 文库制备试剂盒(New England Biolabs, 美国)构建测序文库并在 Illumina 平台上测序, 测序策略为 PE150 [安诺优达基因科技(北京)有限公司]。

1.3 转录组的组装, 注释和DEGs分析

使用 Perl 脚本处理原始数据, 以确保后续用于信息分析的数据质量, 包括去除接头污染的 Reads、低质量的 Reads 及含N 比例大于5%的Reads。对过滤后所得到的 Clean Data的质量和数据量等性状进行统计, 包括过滤后的总序列中质量值大于30(错误率＜0.1%) 的碱基数的比例(Clean Q30 bases rate,Q30)、数据量和碱基含量等。采用 Trinity 软件对进行组装, 对得到的全长转录本进行评估, 然后进行开放阅读框(Open reading frame, ORF)预测和功能分析[14]。采用 RPKM法计算基因表达量[15], 采用DEGseq v1.18.0 软件[16]对差异表达基因进行分析。对DEGs(Q≤0.05 且| log2 ratio |≥1)进行基因本体(Gene Ontology, GO)富集分析及京都基因与基因组百科全书(Kyoto Encyclopedia of Genes and Genomes, KEGG)数据库通路富集分析。

1.4 实时荧光定量 PCR (RT-qPCR) 分析

本实验以 β2-微球蛋白(beta-2 microglobulin,B2M)作为内参基因[14], 对10 个差异表达基因的RNA-Seq 结果进行了验证。根据转录组文库中的Unigene 序列, 用 Primer Premier 5 软件设计特异性引物(表 1)。使用Eppendorf Mastercycler ep realplex4(Eppendorf, 德国)和TB Green Premix ExTaqTMⅡ(Tli RNaseH Plus; TaKaRa, 大连)试剂盒对差异表达基因的表达量进行定量测定, 采用 2-ΔΔCt法[17]计算目的基因的相对表达水平。

表1 RT-qPCR所用基因及其引物序列Tab. 1 Primers used for RT-qPCR

2 结果

2.1 大海马肝脏 mRNA 测序结果及从头组装

高通量测序结果显示(表 2), CK、HS-test、CK和 LS-test 组分别获得55510236、58040146和55840746个原始读段(Raw reads)。过滤后分别获得54453888(CK)、56489512(HS-test)和54407308(LS-test)个Clean reads。Q30分别为97.80%(CK)、95.88%(HS-test)和95.84%(LS-test), 均大于95%, 证明数据可靠。Trinity 软件组装高质量的测序数据后(表 3), 共生成154430个转录本, 最大长度为15866 bp, 最小长度为201 bp, 平均长度为1234.69 bp,N50 长度为2318 bp。所得转录本序列进一步组装处理后获得71794个单基因簇(Unigene), 平均长度为820.71 bp, N50为1780 bp。其中68.07%的Unigenes在200—600 bp, 9.96%的Unigenes在600—1000 bp, 1000—2000 bp占10.87%, 11.10%的Unigenes在2000 bp以上。

表2 大海马测序数据的统计汇总Tab. 2 Statistic summary of the sequencing data of H. kuda

表3 组装结果统计表Tab. 3 Statistics of assembly results

2.2 基因功能注释

将获得的 Unigene 数据信息分别在蛋白质家族域数据库(Pfam)、蛋白质序列数据库(Swiss-Prot)、非冗余蛋白数据(Nr)、真核生物蛋白相邻类的聚簇(KOG)、基因本体论(GO)数据库和东京基因与基因组百科全书(KEGG)等数据库比对注释。结果显示(表 4), 37.61%的 Unigenes在Nr 数据库中注释, 占比最多;核酸序列数据库(Nt)注释的Unigenes占 31.19%; 4781 个Unigenes在KEGG数据库中注释, 仅占比6.66%; 24450个Unigenes 在GO数据库中得到了注释, 占比34.06%。GO数据库涵盖3个方面(图 1): 生物过程(Biological process)、细胞组成(Cellular component)和分子功能(Molecular function)。利用Trinotate的注释结果, 统计每个GO条目中注释到的基因, 并根据二级GO条目得到统计结果。生物学过程中, 占比前5的主要类别依次为细胞过程(72.75%)、生物调节(50.84%)、代谢过程(51.46%)、生物过程调节(47.50%)和刺激反应(36.05%); 在分子功能类别中代表性的类别分别是结合功能(63.10%)、催化活性(34.67%)和分子功能调节(7.62%); 细胞组成部分中最主要的类别为分别为细胞(79.42%)、细胞组分(79.26%)、细胞器(63.61%)和膜(46.63%)。

图1 GO 统计柱状图Fig. 1 GO statistics histogram

表4 各个数据库的注释汇总表Tab. 4 Summary of comments in each database

在 KOG 数据库中, 14428个基因注释到了 25个直系同源组(图 2), 其中信号转导机制(T,23.61%)与一般功能预测(R, 13.54%)代表了最大的类群, 其次是转录(K, 11.73%)、翻译后修饰/蛋白质周转/分子伴侣(O, 8.18%)和细胞内运输/分泌和囊泡运输(U, 6.10%)。

图2 KOG统计图Fig. 2 KOG chart

由图 3 可知,HS-testvs. CK 共获得2740 个差异表达基因(P＜0.05), 其中显著下调495 个, 显著上调 2245 个; LS-testvs. CK 共获得 3715 个差异表达基因, 1854 个显著上调, 1861 个显著下调。高温胁迫组的上调基因数量低于下调基因; 而在低温胁迫组中, 上下调基因数量占比相当。

图3 大海马肝脏转录组中差异表达基因个数统计图Fig. 3 Statistical diagram of differentially expressed genes in liver transcriptome of H. kuda

肝脏转录组中差异表达基因个数统计结果显示(图 4), 两两比对后三组均显著差异表达的基因共计161个, HS-testvs. CK 与 LS-testvs. CK有1260个相同的差异表达基因, 而 1480 个差异基因仅在 HS-testvs. CK显著表达, 共有2455 个差异基因仅在 LS-testvs. CK显著表达。

图4 大海马肝脏转录组中差异基因韦恩图Fig. 4 Venn diagram of differential gene in liver transcriptome of H. kuda

2.3 差异表达基因的 KEGG 富集分析

KEGG 是一个实用程序数据库, 整合了系统功能信息、基因组信息和化学信息, 有助于了解生物系统的高级功能和实用性, 如在分子水平探究细胞、生物和生态系统的信息, 而差异表达基因的Pathway 注释分析有助于在被整合的分子互动网络中进一步解读基因的功能[18], 是国际最常用的生物信息数据库之一。将获得的转录组差异表达基因在KEGG中进行富集比对时, 对 KEGG数据库中的Pathway 信息应用超几何检验进行富集分析, 归纳差异表达基因中显著富集的 Pathway(Q＜0. 05)[14]。其中Q值为多重假设检验校正之后的P值,Q值越小, 表示差异表达基因在该通路中的富集显著性越高[19]。在本研究中, HS-testvs. CK 的差异表达基因共富集到 227 条代谢调节通路, 其中显著富集的通路是类固醇生物合成(Steroid biosynthesis)。LS-testvs. CK 的差异表达基因富集到252条代谢调节通路,其中显著富集的通路有蛋白酶体(Proteasome)、氧化磷酸化(Oxidative phosphorylation)、类固醇生物合成(Steroid biosynthesis)及柠檬酸循环(Citrate cycle)等。

表 5和表 6 分别为HS-testvs.CK和LS-testvs.CK的差异表达基因富集得到的前 20 个 KEGG 代谢通路, 其中高盐胁迫组中代谢途径 (如类固醇生物合成、维生素消化吸收、甘氨酸、丝氨酸和苏氨酸代谢等氨基酸代谢、甘油磷脂代谢、蛋白质运输、磷酸盐和次磷酸盐代谢、不饱和脂肪酸的生物合成、类固醇激素生物合成等)和免疫途径(如PPAR信号通路、抗原加工提呈、NF-κB信号通路、HIF-1信号通路、IgA生成的肠道免疫网络、细胞色素P450对异生素的代谢等)受到影响。低盐胁迫组富集得到的代谢通路主要涉及氨基酸相关代谢有关通路 (如甘氨酸、丝氨酸、苏氨酸、色氨酸、丙氨酸、天冬氨酸和谷氨酸代谢、β-丙氨酸代谢、赖氨酸降解、赖氨酸生物合成)、能量和脂肪酸代谢有关通路 [如类固醇生物合成、氧化磷酸化、柠檬酸循环、碳代谢、淀粉和蔗糖代谢、核糖体生物发生 (真核生物)]和免疫代谢相关通路(如蛋白酶体、抗原加工提呈、抗坏血酸和藻酸盐代谢)等。

根据表 5和表 6富集通路中的有关基因, 结合所查阅的相关文献, 最终筛选得到Gst、Hsp70、Hsp90、Sod、Bcl-2、Gadd45α、Tcrβ、Tap2和Traf3等免疫相关基因,Fadsd6、Fas、Sqle、Cyp51、Elovl6和Slc27a6等脂肪酸代谢相关基因,Vlcad、Pdha1、Mdh1、Idh3b、G6pd和Sdhd等能量代谢相关基因, 及Gldc、Atp6v1e1、Sms、Fadh、Asl、Ass1和Glud1等与氨基酸代谢以及渗透转运有关的基因在盐度应激中发挥重要作用。

本研究共筛选出了Gst、Hsp70、Hsp90、Sod、Bcl-2、Gadd45α、Fadsd6、Fas、Sqle、Cyp51、Elovl6、Pdha1、Mdh1、Idh3b、G6pd、Sms、Asl和Glud1 等对盐度胁迫敏感的候选基因。

2.4 实时荧光定量 PCR 验证

为验证RNA-Seq检测结果的准确性, 通过查阅文献共选取 10个差异表达基因 (表 1)进行验证。以B2M作为内参基因, 采用 RT-qPCR 相对定量方法, 检测这些差异表达基因在试验组与处理组中的表达量差异。将 RT-qPCR 检测结果与转录组分析结果进行比较, 结果表明, 10个差异表达基因的RNA-seq 相对表达水平与 RT-qPCR 相对表达趋势基本一致 (图 5和图 6), 证明了转录组分析结果的可靠性。

图5 高盐胁迫组 (HS-test) RT-qPCR 及 RNA-Seq 的比较分析Fig. 5 Comparative analysis of RT-qPCR and RNA-Seq in HStest

图6 低盐胁迫组 (LS-test) RT-qPCR 及 RNA-Seq 的比较分析Fig. 6 Comparative analysis of RT-qPCR and RNA-Seq in LS-test

3 讨论

水环境中盐度的变化能够对水产动物的行为和生理生化等方面产生十分复杂的影响, 从行为至分子水平的影响均十分严重, 且影响程度存在物种差异性[20]。盐度变化一般会直接破坏水产动物的电解质平衡, 导致渗透压改变, 此外还会影响氨基酸代谢、能量代谢、免疫代谢和应激激素水平等方面[21]。有研究表明, 鱼类和甲壳类水产动物的鳃[22,23]、肌肉[24]和肝脏[25,26]等组织器官都会因盐度变化而受到影响。本研究结果也表明, 高盐和低盐对大海马幼体造成的影响是不同的 (表 5和表 6): 在暴露于低盐时, 主要涉及与氨基酸代谢 (大部分基因上调)、能量代谢 (大部分基因上调)和免疫功能 (大部分基因上调)相关的途径发生改变; 而暴露于高盐时, 免疫功能 (大部分基因上调)和脂质代谢 (大部分基因下调)相关的途径发生了明显变化。

3.1 盐度对大海马免疫功能的影响

盐度胁迫会造成水产动物体内氧化损伤, 其中氧化应激下的酶活性变化及其相关机制已有许多报道[22,27]。在本研究中的低盐胁迫下, 蛋白酶体途径显著富集, 且富集于此途径的基因均显著上调(如26s蛋白酶体非ATP酶调节亚单位11A、B等)。蛋白酶体途径可以响应感染、热休克或氧化损伤等多种细胞应激, 此途径具有重要且复杂的功能。热休克蛋白 (HSP)可以识别错误折叠或未折叠的蛋白质, 通过蛋白酶体的水解作用降解蛋白[28]。本研究发现, 低盐胁迫会使大海马幼体体内产生氧化损伤, DEGs中抗氧化基因 (如Sod和Gst)和热休克蛋白基因 (如Hsp70和Hsp90)的显著上调证实了这一点。此外, 表 5和表 6的结果发现了多个免疫相关通路, 如抗原加工提呈、NF-κB信号通路和IgA生成的肠道免疫网络等, 经分析发现高/低盐胁迫组中多数免疫相关基因显著上调, 如Tcrβ、Tap2、Traf3和Gadd45α等。在水体盐度变化时, 大海马免疫防御系统中的免疫分子转录增加, 推测是为环境应激可能导致的疾病做准备[29]。

表5 HS-test vs. CK差异基因富集的前20个KEGG代谢通路Tab. 5 Top 20 KEGG metabolic pathways enriched for HS-test vs. CK differential genes

3.2 盐度对大海马幼体氨基酸代谢的影响

游离氨基酸、离子通道和水通道蛋白是已知的3种盐度胁迫效应物[30]。据报道游离氨基酸在细胞体积变化和渗透调节过程中具有重要功能, 且大多数水生物种在环境盐度变化时, 将体内的游离氨基酸积累在细胞内[31]。在本研究中, 表 6的前20个KEGG通路中低盐胁迫组富集了多个氨基酸代谢通路, 如甘氨酸、丝氨酸、苏氨酸、色氨酸、丙氨酸、天冬氨酸和谷氨酸代谢、赖氨酸降解和生物合成及β-丙氨酸代谢等代谢途径。除表 6外, 本研究还富集到了半胱氨酸和蛋氨酸代谢、组氨酸代谢、精氨酸生物合成及精氨酸和脯氨酸代谢等通路。以上通路中富集的基因多数显著上调, 表明氨基酸代谢的增加可能是大海马幼体应对低盐胁迫的一个重要调节途径, 类似结果也在牡蛎 (Crassostrea gigas)[32]、尼罗罗非鱼 (Oreochromis niloticus)[33]及凡纳滨对虾 (Litopenaeus vannamei)[26]等水产动物中被发现。但研究发现不同物种中重要的游离氨基酸种类并不一样[32], 因此可以将氨基酸的定性及定量作为今后研究的重要内容之一, 将氨基酸代谢相关的潜在标志基因与氨基酸含量相联系, 进一步了解盐度胁迫下大海马幼体渗透调解过程中氨基酸代谢的机制。

表6 LS-test vs. CK差异基因富集的前20个KEGG代谢通路Tab. 6 Top 20 KEGG metabolic pathways enriched for LS-test vs. CK differential genes

3.3 盐度对大海马幼体渗透调节的影响

肝脏是鱼类代谢的主要器官, 同时可以为渗透调节作用提供能量[34]。关于渗透调节过程中能量代谢的研究主要集中在鱼类的生理生化反应、氧气消耗和氨排泄等方面, 而在分子水平上的研究较少。水产动物的能量代谢途径在响应盐度应激时起关键作用[35], 有研究发现在高盐胁迫下斑点鲈(Lateolabrax maculatus)肝脏中能量代谢相关的许多基因显著上调[25]。而本研究在低盐胁迫下大部分能量代谢的关键基因 (如Vlcad、Pdha1、Mdh1、Idh3b和G6pd等)显著上调, 说明大海马幼体肝脏在低盐胁迫时的能量代谢过程可能承担了重要的角色。推测大量能量用于合成关键蛋白、离子调节和渗透调节等过程, 以维持体内的稳态和渗透平衡[36]。

3.4 盐度对大海马幼体脂肪酸代谢的影响

水产动物的脂肪酸代谢在盐度胁迫下也承担了多种功能, 一种是提供能量维持离子转运并调节渗透压, 另一种是调整细胞膜的结构[26]。在本研究中高盐胁迫下大海马幼体肝脏中富集了多个脂肪酸合成和降解途径 (如不饱和脂肪酸的生物合成、类固醇生物合成、类固醇激素生物合成、脂肪酸代谢、脂肪酸降解和甘油磷脂代谢等), 且通路中的大部分基因显著下调, 如Fadsd6、Cyp51、Sqle、Fabp和Slc27a6等脂肪代谢基因在高盐胁迫下均显著下调, 由此推测大海马幼体体内可能有足够的能量用于抵抗高盐胁迫。