流速胁迫对美国红鱼的转录特性研究

娄宇栋，冯 建 ，何娇娇，晁 帅 ，王 萍 ，2

(1.浙江海洋大学海洋科学与技术学院，浙江省海洋养殖装备与工程技术重点实验室，浙江舟山 316022；2.宁波大学海洋学院，浙江宁波 315211)

水流是鱼类生长发育的基础环境，鱼类的生长、繁殖、迁徙和种群分布都受到水流的影响[1-3]。现阶段有关水流对鱼类的影响主要集中在游泳机制和游泳速率等方面[4]。随着研究技术和方法的深入，许多学者开始从分子水平研究流速对鱼类的生理影响。MAGNONI，et al[5]在研究虹鳟鱼Oncorhynchus mykiss的肌肉转录组时发现，虹鳟鱼在最大续航游泳速度下游泳时能量和物质表达的相关基因大量参与表达。PALSTRA，et al[6]研究发现，在持续的游泳刺激下，虹鳟幼鱼骨骼肌中白肌肉生长发育的基因出现明显上调。同时PALSTRA，et al[7]还发现，在流速影响下，斑马鱼Danio rerio的胰岛素样生长因子受体、生长激素受体等出现明显差异表达。在其它的鱼类研究中，免疫[8]、低温胁迫[9]、营养代谢机制[10]等大量相关基因被发现并得到证实。这些研究使人们从分子水平上对鱼类的行为机制有了新的认识，同时给鱼类的网箱养殖、增殖放流等提供了一些合理的建议。

美国红鱼Sciaenops ocellatus是我国网箱养殖鱼类中的主要经济鱼种[11]。关于美国红鱼的流速研究已有相关报道[12]。但是，关于流速胁迫对美国红鱼肝脏转录组影响的研究鲜有报道。该试验主要探究了美国红鱼的肝脏转录组在最大续航游泳速度下的变化。本实验结果为美国红鱼的养殖海区选址提供了一定的参考依据，并进一步完善了流速对美国红鱼肝脏转录组的研究。

1 材料与方法

1.1 试验材料与设备

美国红鱼购于舟山岱山县养殖基地，并于实验室暂养2周，以适应实验室环境。挑选规格相近，活力较好的美国红鱼作为实验对象。实验室海水为天然海水，并经过过滤沉淀，用次氯酸钠消毒并曝气。实验过程中，实验水温20～21℃；海水盐度27～28；水中溶解氧大于7.0 m·L-1，实验光源为自然光。实验美国红鱼总计 25 尾，体长 31.5～33.8 cm，体重 501.3～589.0 g。避免测量基础数据对实验鱼的影响，基础数据测量均是实验结束后进行。实验全程用摄像机摄像实时监控，以此确定实验鱼的游泳时间，并避免人为因素干扰。

1.2 美国红鱼肝脏转录组的测序

实验前对实验水槽的参数与性能进行校正并测定该规格美国红鱼的最大续航游泳速度，标定实验流速为0.9 m/s。实验前将美国红鱼放入水槽中进行低流速的适应，适应流速控制在0.2 m·s-1左右，适应时间为1 h[13]。将实验组(T01)的适应流速均匀加速到0.9 m·s-1，通过视频监控软件远程监控实验鱼，避免人为因素影响实验鱼的游泳状态，实验鱼触网20 s以上即视为达到疲劳，结束实验[14]。对照组(T02)美国红鱼规格与实验组基本一致，控制实验流速为0 m·s-1。

1.2.1 采样与处理

实验结束后立即捞出实验鱼，并在无菌操作台迅速解剖获得肝脏，放入液氮中冷冻，保存至-80℃冰箱中。根据Tiangen的要求，提取肝组织的RNA(Trizol法)。RNA完整性检测(1%琼脂糖凝胶电泳法)；RNA纯度检测(NanoPhotometer法)；RNA浓度检测(Qubit?)；RNA完整性的评估(NA Nano 6000 Assay Kit of the Agilent Bioanalyzer 2100)[15]。

1.2.2 构建 RNA 文库

检测样品完毕并合格后，启动RNA文库构建，在Illumina平台HiSeq 2500开始测序[16]。cDNA文库的构建和测序由测序公司完成。

1.2.3 生物信息学分析

为了得到较高质量的数据，需要对所得数据进行过滤从而降低数据噪音。同时去除碱基质量低于10，碱基数低于4和碱基质量低于13，碱基数不高于6的序列；删除无插入片段序列、palyA序列、长度小于15 nt的小片段序列、5′或3′接头污染序、列插入片段过长的序列。确保测序的可靠性。经过上述处理，获取干净序列(clean reads)。

2 实验结果

2.1 测序数据统计

测序共计获得9.32Gb Clean Data，两组样品Clean Data 均达到 4.54 Gb，Q30 碱基百分比在 92.12% 及以上。De novo拼装后共获得70 148条Unigene，其中长度在1 kb以上的Unigene有12 465条。对Unigene进行功能注释，包括与 NR、Swiss-Prot、KEGG、COG、KOG、GO和Pfam数据库的比对，共获得31 704条Unigene的注释结果。具体的统计信息如图1与表1。

图1 Unigene长度分布图Fig.1 The length distribution map of Unigene

表1 组装结果统计表Tab.1 The assembly results

2.2 SSR分析

分析检测组装得的70 148条基因共得到SSR标记10 214个，同时单碱基(Mono-nucleotide)重复个数至多为6 052个，六碱基(Hexa-nucleotide)重复至少为2个(表2)。同时发现双碱基对与三碱基对也出现重复，其中AC/GT，共有重复1 685个；AGG/CCT，重复404个。

表2 重复碱基的长度和数量的分布Tab.2 Length and number distribution of SSRs based on the number of repeat units

2.3 SNP分析

通过对实验样品的SNP分析检测发现：实验组T01和对照组T02的纯合型数目分别为38 020和29 627，杂合型SNP数目分别为57 835和66 088个(表3)。

表3 SNP数量统计表Tab.3 Statistics of SNP number

2.4 基因表达分析

通过Bowtie5[17]将实验样本所测定的reads与基因库比较，同时和RSEM5[18]对比，并进行表达量水平的评价。基因转录本的丰度表示为FPKM。总体看，基因片段(编码蛋白质)经过测序得到的基因转录本的丰度值跨度很大[19]，详细如图2所示。

图2 两样品FPKM密度分布对比图Fig.2 The comparison of FPKM density in two samples

对样品T01与样品T02进行了差异表达分析，默认将FDR<0.01且差异倍数(两样品间表达量的比值)≥2作为筛选标准[20]。筛选后得到1 773个显著性差异表达基因，包括204个上调基因和969个下调基因。

绘制两组样品基因表达差异的火山图，图上坐标点各表示1个表达基因，图中的横坐标表示两组样品的差异性表达倍数的对数值，随着差异表达倍数的对数值的线性增加，其绝对值也增大；图中纵坐标表示错误发生率的负对数值，差异表达性决定其数值，差异表达越大，数值越大。将差异表达基因进行统计，无显著性表达差异的基因用红色表示(FALSE)，显著性表达的差异基因用绿色表示(TURE)，如图3所示。

图3 样品T01和样品T02间差异表达基因火山图Fig.3 The Differentially expressed gene volcanic map between T01 and T02

2.5 差异表达基因聚类分析

通常，具有相似表达模式的基因其功能上也存在较高的相似性。因而，通过MySQL Cluster，对两组数据的表达模式进行聚类分析。如图4所示，左侧代表试验组T01的基因表达，右侧代表对照组T02的基因表达，行表示基因表达量(用6种不同颜色表示不同表达量FPKM)。此图说明两个样本存在相似表达模式的基因同时其功能上也存在较高的相似性。

图4 样品T01和样品T02间差异表达基因聚类图Fig.4 Cluster analysis of differentially expressedn genes in T01 and T02

2.6 差异表达基因GO功能富集

GO富集的分析结果说明了在流速胁迫下，美国红鱼应对胁迫所表达的生物过程。富集分析发现，生物过程、细胞组分和代谢过程3组，总计高达424个富集到生物过程中，其中有90个富集到细胞组分中，310个富集到分子功能。在实验结果中发现，总计有80个差异基因富集在代谢过程上；氧化还原过程(oxidation-reduction process)和负调节肽链内切酶的活动(negative regulation of endopeptidase activity)显著性富集；在细胞组分上，细胞类别上共有90个差异表达基因富集，其中发现显著性富集的有胞外区 (extracellular region)和细胞外间隙(extracellular space)；此外有310个差异表达基因富集到分子功能类别上，实验发现只有一个内肽酶抑制活性(endopeptidase inhibitor activity)存在显著性富集。

实验样本T01和对照样本T02间差别表达基因以及全部检测的基因在GO二级Term的详细解释结果如图5所示。从图5可以看出，GO各功能位于不同差异表达基因和所有基因两个背景下的处所，从图中发现红色柱体与蓝色柱体具有明显差别的Term也许和差异有关。

图5 差异表达基因的GO分类Fig.5 Gene Ontology classifications of differentially expressed genes

2.7 差异表达基因KEGG注释

生物体的统一性由大量差别基因配合调节，彼此作用。不同基因行驶不同功能，各司其职。我们进行了基因差异性表达的Pathway注释分析去认识各个基因的差异功能。将得到的Pathway和KEGG(Kyoto Encyclopedia of Genes and Genomes)数据库进行对比。实验数据表明，共有1 096个基因数差异表达；发现423条基因被注释到KEGG通路上。其中富集到新陈代谢通路、环境信息处理和细胞过的差异表达基因比较多，分别为136、69和67条。分析这423条代谢通路发现有4条代谢通路显著富集(Qvalue≤0.05)，分别参与甘氨酸、丝氨酸和苏氨酸代谢(glycine，serine and threonine metabolism)(p=0.000 9)、色氨酸代谢(tryptophan metabolism)(p=0.0036)、视黄醇的代谢(retinol metabolism)(p=0.003 6)和类固醇类的生物合成(steroid biosynthesis)(p=0.030 3)，具体如图 5 所示。

2.8 与运动生理相关的差异表达基因

对差异表达基因的KEGG代谢通路富集分析，得到了美国红鱼应对流速胁迫的一些相对基因。如Hedgehog 信号通路(Hedgehog signaling pathway)中 Hh、PKA、CK1、Wnt等基因表现为上调、RNA 降解(RNA degradation)中的Dcp1、EDC3基因表现为下调，而TOB基因则表现为上调；甾类激素生物合成(steroid hormone biosynthesis)过程中的 17beta-estradiol 17-dehydrogenase酶、testosterone 17beta-dehydrogenase(NADP+)酶会表现下调，葡萄糖醛酸转移酶(glucuronosyltransferase)等酶表达量会表现上升、组氨酸的代谢(histidine metabolism)过程中的 aldehyde dehydrogenase(NAD+)、组氨酸解氨酶(histidine ammonia-lyase)、天冬酰转移酶(aspartoacylase)表达量会表现上升。通过GO功能富集分析与KEGG信号通路富集分析，得到了美国红鱼在流速胁迫下与运动和疲劳有关的基因，同时对比了在胁迫下与正常情况下这些基因在肝脏中表达的不同。具体实验结果如表4。

表4 根据GO和KEGG信号通路富集分析游泳相关的差异表达基因Tab.4 According to GO and KEGG signaling pathways to analyzed the swimming related differentially expressed genes

图6 前20个差异表达基因信号通路富集Fig.6 Scatter chart of top 20 pathway enrichment of differentially expressed genes

2.9 转录组测序数据的验证

为了验证RNA-seq分析的表达谱，对11个基因进行了相关mRNA水平上qPCR测定。这些基因的表达数据和通过RNA-seq和qPCR检测到的转录亚型见表5和图7。这些基因和转录本的表达，均通过qPCR得到证实。经过Spearman的rho测试发现qPCR数据之间显著相关(R2=0.979 7)。根据qPCR及RNA测序的数据，上调和下调基因均非常吻合。

表5 RT-PCR和RNA-Seq的对比验证结果Tab.5 Comparisons between RNA-Seq data and RT-PCR results

3 讨论与分析

本实验发现，美国红鱼的游泳能力很强，与王萍等[21]的研究结果基本一致。在分子方面，美国红鱼无参考基因[22-23]，从基因上研究流速胁迫对其的影响机制仍不明确。本文采用Nanordrop、Qubit2.0、Aglient2100 方法对 RNA 样品进行分析。测序总共获得9.32Gb Clean Data，各样品 Clean Data 均达到 4.54 Gb，Q30碱基百分比在92.12%及以上。De novo组装后共获得70 148条Unigene。其中长度在1 kb以上的Unigene有12 465条。本次样品测序效果较好，文库构建成功可用，同时进行基因结构注释、基因表达量分析和基因功能注释等一系列生物信息分析，为今后有关美国红鱼功能基因研究提供了宝贵的资源与数据库。

图7 用RT-PCR验证的结果Fig.7 Validation of RNA-seq data using qPCR

相关文献发现，在其它胁迫下，如温度、盐度、pH、重金属等会使鱼类的基因表达产生巨大的差异[24-25]。本试验中，通过差异基因的分析发现，在流速胁迫的下，美国红鱼试验组有大量显著性差异表达的基因。经过筛选后获得显著差异表达的基因有1 773个，其中上调的基因数量为204个，下调的基因数目为969个。这些差异表达的基因与许多重要的过程和途径相关，如糖的代谢、RNA剪接、蛋白质的分解代谢、核糖体合成、剪接体等。该研究结果对鱼类应对流速胁迫的特异信号网络的进一步探究具有十分重要的意义。

在生物过程中，富集差异基因最多的为代谢过程，显著性富集的类别是氧化还原过程和负调节肽链内切酶的活动，内肽酶抑制活性。研究发现，在环境胁迫下，植物的氨基酸和糖类的代谢受到显著影响[26]。在相关实验中发现美国红鱼的游泳疲劳主要与肝糖原相关[12]。这些差异表达基因的发现，能从分子层次说明美国红鱼在持续游泳200 min内主要进行有氧呼吸，但是高强度的持续游泳，鱼体内的的糖类得不到有效补充，从而发生了糖的异生作用。氧化还原过程直接决定细胞衰老和影响蛋白质的功能[27-29]。细胞组分方面，胞外区和细胞外间隙的显著性富集说明美国红鱼在运动过程中，伴随着剧烈的物质交换，物质跨膜运输频繁。内肽酶[30-31]在植物中决定着植物叶片的老化，在动物中对动物的造血系统[32]、呼吸功能[33]等有着一定的调节作用。本实验中发现在分子功能类别中，内肽酶抑制活性基因有显著性富集，因此推测在流速胁迫下，美国红鱼的造血系统和呼吸系统受到一定的影响，具体的机理尚未明确。应激信号的检测及其细胞内转导是在各种环境胁迫下生物的适应和生存的关键。关于鱼类应对流速胁迫信号并触发细胞内的反映机制的信息是极为稀缺的。在本研究中，Hedgehog信号通路中Hh、PKA、CK1、Wnt等基因表现为上调、RNA降解(RNA degradation)中的Dcp1、EDC3基因表现为下调，而TOB基因则表现为上调。在老鼠的研究中，TOB基因特异性表达于大鼠的海马和小脑，其与大鼠的运动协调性有关[34]。

在早期分的化斑马鱼肌细胞中，Hh信号能控制快慢肌之间的命运选择[35]。除了调节早期二元态(慢与快纤维类型)细胞命运的决定，Hh信号还诱导这些谱系中不同类型的细胞生成(肌肉祖细胞、慢肌纤维的子集、内侧快肌纤维和快肌纤维的子集)[36]。在美国红鱼持续游泳过程中，推测Hh信号将诱导慢纤维的生成。同时肌节分化需要神经管中的Shh和Wnt信号，Wnt信号在胚胎生肌过程中发挥重要作用，且Wnt与Shh相互合作诱导肌节形成[37]。在成熟成骨细胞中，Hedgehog信号通路活性增强将上调细胞内甲状旁腺激素相关蛋白的表达量，通过PKA等信号调控轴促使环磷酸腺苷应答元件结合蛋白与RANKL基因增强子结合，促进RANKL基因表达，促使破骨前体细胞分化为成熟的破骨细胞，进而增强骨吸收活动[38]，从而进一步提高美国红鱼的游泳能力。

甾类激素是一类脂溶性激素，它们在结构上都是环戊烷多氢菲衍生物，普遍存在于动植物细胞之中[39]。甾类激素生物合成过程中的17beta-estradiol 17-dehydrogenase酶、testosterone 17beta-dehydrogenase(NADP+)酶会表现下调，说明美国红鱼在游泳过程中，糖的氧化受到抑制，迫使血液中的糖持续升高，另一方便又促进蛋白质等非糖物质转化为糖(糖的异生作用)。葡萄糖醛酸转移酶(UGT)是重要的II相代谢酶，在大鼠和水牛[40-41]的研究中发现，肝片吸虫感染的肝脏中UGT的活性均有不同程度的下降。在鱼类中，却是鲜有研究。本文中的(UGT)等酶表达量会表现上升组，其中具体原因还需进一步验证。组胺1型受体[42]主要分布于内皮和平滑肌等多种细胞，对于血管扩张有一定的促进作用[43]。本实验中发现在组氨酸的代谢过程中的Aldehyde dehydrogenase(NAD+)、组氨酸解氨酶、天冬酰转移酶表达量会表现上升，说明在组氨酸脱羧酶的作用下，组氨酸脱羧形成组胺，在高强度的游泳过程中，迫使血管急速扩张，体内物质交换加速，血液循环加快。

综上所述，本实验发现流速胁迫对美国红鱼肝脏转录组存在很大的影响。进一步的分析发现很多与运动生理相关的编码如RNA剪切、核糖体真核生物核糖的体合成、蛋白质的分解代谢、剪接、蛋白酶和RNA的运输的重要组成的蛋白质的基因。具体作用机理还需要进一步深入研究讨论。