1. 引言
刀鲚(Coilia nasus)属鲱形目(Clupeiformes)鳀科(Engraulidae)鲚属(Coilia)，俗称刀鱼、河刀 [1]。刀鲚主要分布在我国黄渤海和东海一带 [2]，辽宁境内的辽河、鸭绿江及大洋河亦有分布 [3]。刀鲚味道鲜美，长久以来都有“长江三鲜”的美誉，是重要的经济型鱼类，但近年来，由于水利工程兴建以及环境污染等不利因素，导致刀鲚的生存繁殖环境遭到严重的破坏。特别是近10年间，由于不科学的过度捕捞导致刀鲚种群量逐年递减，现已被世界自然保护联盟(IUCN)列为濒危物种 [4]。辽宁境内水域是我国刀鲚分布的最北线，有关其研究相对较少，仅见有繁殖生物学、遗传多样性及活体运输等方面的报道 [5] [6] [7]：利用DNA分子技术探究长江刀鲚、湖鲚的遗传多样性 [8]；利用转录组测序技术对刀鲚嗅觉上皮组织进行比较 [9]。本文以大洋河流域的洄游型刀鲚为研究对象，对来自不同盐度水域刀鲚鳃的转录组进行了比较分析，探究在转录组水平上刀鲚鳃部组织对盐度变化的响应机制，为今后研究大洋河刀鲚的生活习性和环境适应性机制提供研究数据和理论参考。
2. 材料与方法
2.1. 样品采集
2020年8月用流刺网在大洋河下游不同河段采捕刀鲚，本试验取河口咸水区(盐度为12‰，样品标记为DE)和石山桥河段淡水区(盐度为0.2‰，样本标记为SG)大小相近的鲜活刀鲚各5尾，刀鲚全长和体重分别为(27.52 ± 0.81) cm、(68.49 ± 6.96) g。剪去鳃部组织置于RNAfixer (Bioteke北京)保存液中带回实验室备用。
2.2. 总RNA提取与检测
取保存的鳃部组织100 mg左右，采用TRizol法提取total RNA。利用Thermo Nanodrop2000对所提RNA的浓度和纯度进行检测，并用1.5%琼脂糖凝胶电泳检测RNA完整性。采用Agilent2100测定RIN值，单次建库要求RNA量1 ug，浓度 ≥ 50 ng/µL，OD260/280介于1.8~2.2之间。
2.3. 文库构建
利用带有Oligo (dT)的磁珠与ployA进行A-T碱基配对，总RNA中分离出mRNA，加入fragmentation buffer，可以将mRNA随机断裂成300 bp左右的小片段。通过逆转录酶加入六碱基随机引物(random hexamers)，以mRNA为模板反转合成一链cDNA，随后进行二链合成，形成稳定的双链结构。双链的cDNA结构为粘性末端，加入End Repair Mix将其补成平末端，随后在3’末端加上一“A”碱基，用于连接Y字形的接头。
2.4. 转录组测序、质控与组装
本实验采用Illumina HiSeq2500测序平台完成转录组测序，Illumina测序基于循环可逆终止技术工作，流程如下：
1) 技术代表目前新兴的基因组DNA提纯后被随机打断。这一步可以通过物理方法完成，如声波法、剪切法，或者雾化法，通常进一步通过长度分选随机打断的DNA片段。在两端都接上接头。
2) 单链DNA片段共价连接到流动细胞通道的表面。
3) 加入DNA聚合酶和未标记的脱氧核苷酸产生固相“桥扩增”，其中模板DNA使两端连接到通道表面形成U形环。
4) 双链桥生成。双链分子变性，然后继续扩增以形成高度簇集的模板DNA。
5) 加入四个标记的可逆终端(包含引物和DNA聚合酶)。在给定的循环中，一个可逆终端只能被加入一个模板。在特殊的不能延长的碱基处会产生链终止。
6) 在激光的激发下，第一个碱基的身份被记录。
7) 在第二个循环中，可逆终端被去除(保护)。所有四个标记的可逆终端和聚合酶再次被加入流细胞中。这个循环被重复。
之后对获得的测序数据进行质量控制(QC)，之后利用生物信息学手段对转录组数据进行分析。其中使用fastx_toolkit_0.0.14软件对对每一个样本的碱基质量、碱基错误率以及碱基分布进行分析；然后使用软件SeqPrep (https://github.com/jstjohn/SeqPrep)和Sickle (https://github.com/najoshi/sickle)去除低质量reads得到高质量的质控数据(clean data)；最后由Trinity (https://github.com/trinityrnaseq/trinityrnaseq/wiki)软件将clean data进行从头组装。
2.5. 数据处理
将测序过程中的图像信号经CASAVA碱基识别(Base Calling)转换成文字信号，并将其以fastq格式储存起来作为原始数据。根据index序列区分各个样本的数据，以便进行后续分析。单次运行能产生数十亿级的reads。将测序获得的reads与unigenes通过Bowtie 2.3.5进行比对，然后结合比对结果，经由RSEM 1.2.2软件的分析得到表达量水平的结果估计。最终输入不同样本的拷贝数(read counts)信息构成的矩阵，从而得到表达量信息差异分析结果。
2.6. 功能注释
通过BLAST 软件将该次转录组测序获得的所有转录本与六大数据库(NR, Swiss-Prot, Pfam, COG, GO和KEGG数据库)进行比对，获得在各数据库的注释信息，并使用HMMER与Pfam软件对各数据库注释情况进行统计和分析；并对SNP功能区域和SSR进行统计。
3. 结果与分析
3.1. 刀鲚鳃组织RNA检测结果
提取得到的RNA样品具有完整清晰的28S、18S和5S带型，A260/280值为2.02，28S:18S 值为1.80，RIN值为8.0，说明本次提取得到的刀鲚组织样品RNA质量较好，可以用于制备测序文库。
3.2. 测序数据与组装结果
利用Trinity将拼接过滤后得到的reads片段进行聚类及进行拼接组装，共得到70,964条unigenes和100,022条转录本(transcript)，平均长度分别为1061.73 bp和1140.11 bp。transcript与unigenes的N50 (重叠群序列累加后长度超过转录组总长度一半时的重叠群序列长度)分别为2227和2209 (表1)。其中长度在0~500 bp之间的unigenes数量为36,872条，占总数量的52%。长度在501~1000 bp为12,779条，总unigenes数量的18%；1000 bp长度以上的unigene数量均不超过10% (表2)。
Table 1. Statistical table of assembly results
表1. 组装结果统计表
Table 2. Sequence length distribution of unigenes
表2. Unigenes序列长度分布情况
3.3. Unigene功能注释统计
将测序获得的所有转录本在各数据库注释的整体情况进行统计，结果表明，有22,342个unigenes获得了GO注释，占全部unigenes的31.48% (表3)。GO注释的功能分析由3大部分组成，可以对基因和基因产物按照其参与的BP (Biological Process，生物过程)、MF(Molecular Function，分子功能)及CC (Cellular Component，细胞组分)方面进行分类注释(图1)。在这三个大分支下面又分很多小层级，功能上的细分更
Table 3. Annotation result statistics of unigenes
表3. Unigenes注释概况统计表
有助于从整体上了解全部基因产物的功能分类。本次实验重点关注了生物过程这一分支，在有关生物过程下属的层级的注释中(图2)：细胞进程(celluar process)共9214条占29.22%；生物调节(biological regulation)共5646条，占17.90%；代谢过程(metabolic process)共5127条，占16.26%，NCBI_NR(NCBI非冗余蛋白库)为综合数据库，可通过比对查看本物种转录本序列与相近物种的相似情况，以及同源序列的功能信息。如图3所示，刀鲚(Coilia nasus)的unigene序列大西洋鲱(Clupea harengus)相似度最高，为58.68%，另外还有19.60%的序列被注释到其他物种之中。通过与KEGG数据库比对，可获得某基因或转录本可能参与的具体生物学通路情况，这些信息有助于从系统水平解读基因的生物学功能。如图4所示，6大代谢途径中，参与生物体系统(Organismal Systems)的信号通路的unigene所占比例最高，遗传信息处理(Genetic Information Processing)所占比例最低。通过对unigene的注释和蛋白差异表达分析，为后续了解基因功能和解释表型差异提供数据基础。
Figure 2. GO (BP) analysis pie charts
图2. GO (BP)分析饼图
3.4. 表达量统计
RSEM，Kallisto，Salmon是三款常用的转录组定量分析软件，本次实验采用RSEM进行转录组的定量，并完成转录组和组装结果的比对。三款软件的输出结果包含表达量TPM或FPKM (仅RSEM)信息和read counts (比对到基因上的reads个数)信息。不同样本的read counts信息构成的矩阵可以用于差异分析的输入，表达量信息用于后续的样本聚类等分析。如图5所示，样本DE_1的表达量最高；样本SG_3表达量最低，同盐度样本中，样本之间的表达量变化较大，总体来看不同盐度的样本之间，随着盐度的变化的总体表达量无明显的差异。
注：横坐标为样本名称，纵坐标为样本表达量(log10 TPM+1)
Figure 5. Statistical line chart of expression quantity
图5. 表达量统计折线图
3.5. 样本间Venn分析
Venn 分析展示样本间或组别间共有和特有表达的基因/转录本，可简单呈现样本间相关性，同一组别样本表达基因/转录本的数目不应差别很大。从图6可以看出，DE组的表达unigene数目为34,768；SG组的表达unigene数目为23,546；两组共有的组的表达unigene数目为20,626，总体上开看DE组的unigene的总表达量数量要高于SG组。
Figure 6. Expression quantity Venn diagram
图6. 表达量Venn图
4. 讨论
全基因组测序技术的发展和测序数据免费公开极大地改变了生物学的研究方式，测序和其他基因组学数据具有推动农业、环境科学和生态学的巨大潜力。近些年得益于高通量测序技术的快速发展和普及，基因组测序的水平有了显著的提高 [10] [11]，对于刀鲚这类非模式生物而言，NGS (高通量测序)是一个方便且可以快速得到其基因序列的方法。本文通过对大洋河不同水域刀鲚样本的鳃组织进行转录组测序，得到可信的unigenes
1) 其中长度在0~500 bp之间的unigenes数量为36872条，占总数量的52%。长度在501~1000 bp为12779条，占总unigenes数量的18%。测序样本的Mapped read数分别为SG-1:42641542、SG-2:38980922、SG-3:44444358、SG-4:41879374、SG-5:41038376、DE-1:44473856、DE-2:43234696、DE-3:43123052、DE-4:47355458、DE-5:40274730，序列在NR、Swiss-Prot、Pfam、COG、GO和KEGG公共数据库中进行注释。
2) NR结果表明刀鲚的uingene序列与大西洋鲱(Clupea harengus)的最为接近，达到了58.68%，符合二者同属鲱形目的动物学分类。其次有5.59%的unigene在齿鲱(Depticeps clupeoides)中得到注释，说明虽同属鲱形目，但是二者之间差异巨大。另外还有19.60%的unigene被注释到其他物种中。
3) 在GO注释的生物过程unigenes表明，细胞进程、生物调节、代谢过程这三类unigene占据很大比重，这些注释的unigenes均是参与调控和维持细胞正常生理活动的基因，这三类的高占比，在一定程度上体现了刀鲚的鳃细胞代谢旺盛，增殖快，新陈代谢速率快等特点。这些基因既保证了细胞的正常运转，也可能与调控刀鲚鳃细胞的增殖、分化密切相关。另外GO的注释结果中有相当部分的未知功能的unigenes没有得到注释，说明对于刀鲚unigenes的研究还有还大的空间。
4) KEGG通路注释分类结果中，生物体系统(Organismal Systems)的信号通路unigene数量最多，表明在生长发育及生命活动过程中的代谢活动非常旺盛，这些信号通路起着至关重要的生理调控作用；另外参与环境信息处理(Environmental Information Processing)的信号通路的unigene数也很多，说明刀鲚鳃组织细胞对环境因子的敏感程度很高。某些硬骨鱼类鳃组织的酶活性与盐度的变化有一定关系 [12]，鳃丝表皮细胞对维持渗透压等生理功能有重要的调节作用 [13]，可能与刀鲚适应海水盐度或者水温的频繁变化与维持自身生理系统稳定和提高环境自适应能力的unigene信号通路数量占比高，具有一定的联系。
在生物体内，不同基因产物之间通过有序的相互协调来行使其具体的生物学功能，基因表达量的变化一定程度上体现了细胞代谢的强弱，某些酶的表达对鳃的生理功能有着重要的作用，例如在某些硬骨鱼类中鳃丝Na+/K+-ATP酶活力会随着盐度的变化而变化，而某些硬骨鱼类从半海水过渡到淡水或者海水的过程中，其鳃丝Na+/K+-ATP酶活力没有明显的变化；HEIJDEN [14] 等认为，盐度对鱼类鳃丝Na+/K+-ATP酶活力影响并不显著。通过对转录组定量的分析可以看出DE_1的表达量最高，SG_3的表达量最低，从总体来看DE组的表达量大于SG组的表达量(图6)，但不同盐度之间的表达量差异并不与盐度呈现相关性，从结果来看盐度对刀鲚鳃部基因的表达量的影响并不显著。大洋河刀鲚属于洄游性鱼类，不同于其他广盐性鱼类，大洋河刀鲚活动范围仅限于近海，区别于其他鱼类，刀鲚在生存环境上的差异不如其他鱼类那样显著。目前针对环境因子对刀鲚鳃组织基因表达和生理功能的影响的研究较少，对于刀鲚鳃组织基因表达蛋白通路的作用机制尚不明确，今后对于盐度对刀鲚鳃组织转录组影响的研究还有很大的空间。
5. 展望
目前关于转录组的报道和研究主体还是以四大家鱼、鲆鲽类、大西洋鲑等市场经济型鱼类为主 [15] [16] [17] 关于刀鲚及其种属的研究目前尚不多见。刀鲚作为一种名贵鱼类，其经济价值潜力巨大。日后旨在探究在洄游过程中，不同盐度对刀鲚基因表达及生理功能的影响，还应在各个层次水平开展更为深入的研究，为刀鲚的资源养护提供更多的理论基础，提升完善整个刀鲚的理论体系，为日后的发展创造更好的理论条件。
基金项目
辽宁省农业农村厅项目(JH20-210000-39754)、辽宁省农科院项目(2021HQ1918)及辽宁省科技厅项目(2021JH2/10200031)。
NOTES
*通讯作者。