李 青， 何 斌， 陈群利， 游 萍

(贵州工程应用技术学院生态工程学院，贵州毕节 551700)

鲫鱼()隶属于鲤形目鲤科鲫属，由于其营养价值高、味鲜肉嫩、生长速度快、杂食性和适应性强等特点，是我国重要的淡水养殖经济鱼类。鲫鱼多生活于水域的中下层，冬天低温条件下可以忍受缺氧数月，具有极强的低氧耐受性，被公认为是低氧耐受力极强的鱼类之一。此外，鲫鱼具有二倍体、三倍体和四倍体3种不同的染色体倍性类型，不同类型具有不同的繁殖方式。二倍体通过两性繁殖方式产生雌雄同体后代；四倍体通过雌核发育产生全雌后代；而三倍体鲫鱼同时存在这2种繁殖方式，即通过精子激活雌核发育成全雌后代和雌雄交配产生雌雄同体三倍体后代，鲫鱼是研究鱼类进化发育基因组学和低氧胁迫适应分子机制的适当模型。

近年来，随着测序技术的不断发展，在研究目标物种全基因组数据未知的情况下，转录组测序已成为一种短期、低成本得到海量基因数据的高效方法。基因的表达具有组织特异性，研究生物过程相关组织中基因的表达模式，可以为进一步阐明其分子机制提供基础依据。目前，关于鱼类多倍体起源、遗传多样性和耐低氧机制尚未完全了解。血液可以反映低氧应激时机体代谢水平和器官机能状态的变化，肝脏作为动物体最大的代谢器官，低氧环境会导致鱼类肝脏中氧化氢酶、谷胱甘肽过氧化物酶和超氧化物歧化酶等活性上升。此外，研究表明，低氧可能导致鱼类产生强烈的应激反应而影响其生殖、生长和发育等。

因此，本研究通过RNA-Seq测序技术对鲫鱼肝脏、血液和卵巢组织进行转录组测序分析，旨在发掘鲫鱼低氧胁迫相关调控功能基因，以期为水产养殖过程中鱼类适应低氧胁迫策略提供参考；此外，开发一批分子标记，以期为鲫鱼种质资源评价、分子育种和群体遗传多样性分析等方面提供基础信息。

1 材料与方法

1.1 样品采集

试验用鲫鱼为个体大且体表无损伤的成鱼(体质量为120～150 g，体长为 20～25 cm)，于2019年购自贵州省毕节市水产品交易市场，实验室充气暂养1周后，尾部取血后，分别取肝脏和卵巢组织，经液氮速冻后于-80 ℃保存、备用。

1.2 总RNA提取

分别取肝脏、血液和卵巢组织各150 mg，采用TRIzol法提取其总RNA。琼脂糖凝胶电泳分析RNA降解程度及是否有污染，Nanodrop检测RNA纯度，Qubit2.0对RNA浓度进行精确定量，Agilent 2100精确检测RNA的完整性。将5尾个体相同组织样品mRNA等量混合，用于cDNA文库的构建。

1.3 建库测序和拼接组装

使用NEBNext® UltraRNA Library Prep Kit方法构建鲫鱼肝脏、血液和卵巢转录组文库，库检合格后，进行Illumina HiSeq高通量测序。得到的原始数据经CASAVA碱基识别分析转化为原始测序序列(raw data)，对原始测序序列过滤，得到clean reads。采用Trinity对clean reads进行拼接，转录本序列信息以FASTA格式储存。

1.4 基因功能注释

测序数据采用无参考基因组分析方法，为获得全面的基因功能信息，进行七大数据库的基因功能注释，包括非冗余蛋白序列数据库(Non-Redundant Protein Sequence Database，NR)，核酸序列数据库(Nucleotide Sequence Database，NT)，蛋白质家族的集合数据库(Protein Family，Pfam)，蛋白相邻类的聚簇数据库(Cluster of Orthologous Groups of Proteins，COG)，真核生物蛋白质同源簇数据库(euKaryotic Orthologous Groups，KOG)，Swiss-Prot蛋白质序列数据库(Manually Annotated and Reviewed Protein Sequence Database，Swiss-Prot)，京都基因和基因组百科数据库(Kyoto Encyclopedia of Genes and Genomes，KEGG)和基因本体数据库(Gene Ontology，GO)。

1.5 SNP、InDel和SSR分析

使用Samtools和PicardTools工具对比，将结果进行染色体坐标排序，去掉重复的reads，通过变异检测软件GATK2分别进行单个核苷酸变异多态性(single nucleotide polymorphisms，SNP)和核苷酸插入或缺失(insertion-deletion，InDel)，并对原始结果进行过滤，过滤掉质量值小于30，距离小于5的SNP。采用MISA 1.0对Unigene进行重复序列标记(simple sequence repeats，SSR)检测，并对不同SSR类型在基因转录本的密度分布进行统计。

1.6 数据处理

采用RESM软件，以Trinity拼接得到的转录组作为参考序列(ref)，将每个样品clean reads往ref上做mapping，得到每个样品比对到每个基因上的readcount数目，使用每百万fragments中来自某一基因每千碱基长度的fragments数目(expected number of fragments per kilobase of transcript sequence per millions base pairs sequenced，FPKM)表示基因的表达水平。使用DEGSeq软件分析不同组织样品基因表达差异，从而筛选出差异基因，然后使用K-means和SOM进行差异基因聚类分析，获得3种组织中差异基因表达量的聚类模式。为了更好地对差异基因的功能进行研究，分别使用GOSeq和KOBAS软件对3个组织间差异基因进行GO富集和KEGG富集分析，同时将差异基因根据上调或下调分别进行富集分析。

2 结果与分析

2.1 数据质量评估与序列组装分析

经过数据质量控制，肝脏、卵巢和血液组织分别得到20 332 138、20 209 531、20 451 910个clean reads，平均GC含量为46.21%。将3个组织的reads片段拼接组装，共得到219 192条转录本(transcript)，平均长度为1 088 bp，最长的片段长度为17 286 bp，最短的片段长度为181 bp，N50和N90(将拼接转录本按照长度从长到短排序，累加转录本的长度，到不小于总长50%～90%的拼接转录本的长度)分别为2 096 bp和397 bp。共得到127 801条unigene，平均长度为735 bp，最长的片段长度为17 286 bp，最短的片段长度为201 bp，N50和N90分别为1 288 bp和283 bp。转录组拼接组装的数据已提交至BioProject数据库(BioProject ID：PRJNA735422)。

2.2 功能基因注释与分类

2.2.1 功能注释与相似性 通过BLAST与相应数据库比对，最终获得有注释信息的unigene数量为117 414个，约占总unigene数量的91.87%。比对到NT数据库的unigene数量最多，为116 885条，占总unigene数量的91.45%，其次，为NR(39.44%)、KO(17.23%)、KOG(16.21%)、Swiss-Prot(6.52%)、GO(0.08%)，而没有unigene比对到Pfam数据库。在以上7个数据库中至少1个数据库注释成功的unigene有117 414个，占总unigene数量的91.87%。

以NR数据库为例，对鲫鱼unigene序列的相似性进行分析，在斑马鱼()中相似序列匹配的比例(73.1%)最高，其次为墨西哥丽脂鲤()(52)、虹鳟()(3.4%)、鲤鱼()(1.2%)和罗非鱼()(1.2%)。匹配序列的相似度都大于40%，其中相似度在80%～95%占的比例(45.1%)最高，其次为相似度介于95%～100%(29.6%)、60%～80%(21.1%)、40%～60%(4.2%)。值在0区间内的unigene数量最多，占总体的19.2%，值介于10～10之间的unigene数量占总体的15.0%，值介于10～10之间的unigene数量最少，占总体的10.5%。

2.2.2 unigene的GO注释分类 由图1可知，鲫鱼3个组织的转录组共有105条unigene在GO数据库生物学过程(biological process)、细胞组分(cellular component)和分子功能(molecular function)三大类42个功能中找到对应。其中，生物学过程获得的注释信息最多，分子功能获得的注释信息最少。在生物学过程类中，获得的注释信息多集中于细胞学过程、单一有机体过程、生物调节、代谢过程、生物过程调节、应激反应和信号等功能。细胞组分类中，获得的注释信息多集中于细胞、细胞组件、大分子复合物和细胞器等方面。分子功能类中，在绑定和催化活性中分布较多，在酶调节器活性中分布最少。GO数据库注释分类有助于从整体上大致了解目的组织全部基因产物的功能，为进一步探究基因的功能提供参考。

2.2.3 unigene的KOG注释分类 由图2可知，KOG数据库可以对基因产物进行直系同源分类，结果显示共有20 725条unigene得到注释，依据功能分为26个直系同源功能类型。其中，参与一般功能预测和信号传导机制功能的最多，其次为翻译后修饰、蛋白质折叠与分子伴侣，转录，细胞内转运、分泌与囊泡运输。

2.2.4 unigene的KEGG注释 由图3可知，通过与KEGG数据库比对，鲫鱼unigene归属于A(细胞过程)、B(环境信息处理)、C(遗传信息处理通路)、D(代谢)和E(有机系统通路)五大类；有机系统通路注释比例最高，而遗传信息处理通路所占比例最低。其中注释数量最多的代谢途径是信号传导途径(3 856个)，注释数量大于1 000的为内分泌系统(1 776个)，免疫系统(1 769个)，运输和分解代谢途径(1 426 个)，细胞群落(1 356个)，折叠、排序和降解途径(1 191个)，神经系统(1 190个)，信号分子与相互作用途径(1 086个)。

2.3 不同组织差异表达基因分析

FPKM是每100万个fragments中来自某一基因每千个碱基长度的fragments数目，其同时考虑了测序深度和基因长度对fragments计数的影响，是目前最为常用的基因表达水平估算方法。因此，笔者所在课题组将readcount数进行了FPKM转换，以不同样品作为横坐标，不同样品表达量(FPKM+1)的对数值作为纵坐标，绘制3个组织基因表达量的FPKM箱线图，由图4可知，3种组织表达量依次为血液>卵巢>肝脏，表明基因在不同组织中转录水平具有明显差异。

进一步比较差异表达基因数目，由表1可知，卵巢vs血液差异表达基因数量最多(9 726个)，肝脏vs卵巢次之(8 390个)，肝脏vs血液最少(6 142个)，且卵巢vs血液上调基因数量最多，肝脏vs卵巢上调基因数量最少，与之相反，肝脏vs卵巢下调基因数量最多，卵巢vs血液下调基因数量最少。

表1 差异表达基因数目(DEGs) 个

由图5可知，卵巢vs血液和卵巢vs肝脏差异表达基因模式相似，多集中于细胞组分中的胞内、细胞内组分、细胞器、细胞内细胞器和细胞质，和生物学过程中的初级代谢过程、细胞大分子代谢过程、细胞蛋白代谢过程、细胞组分组织或生物合成；在分子功能方面，仅在核糖体结构组成部分有少量差异表达基因，而肝脏vs血液差异表达基因数量最少，多集中于细胞组分中的蛋白质复合体、中间丝状体和中间丝状体细胞骨架，在生物学过程中的气体运输和氧气运输，分子功能中的氧气转运子活性、氧气结合、铁离子结合、亚铁血红素结合和四吡咯结合方面有少量差异表达基因(图5-C)。差异基因的富集分析，为进一步筛选组织特异性表达基因及探究组织特异功能分子调控机制奠定了基础。

在生物体内，不同基因相互协调行使其生物学功能，通过Pathway显著性富集能确定差异表达基因参与的最主要生化代谢途径和信号转导途径。由图6可知，KEGG是有关Pathway的主要公共数据库，依据KEGG代谢途径数据库，找出差异基因相对于所有有注释的基因显著富集的Pathway，绘制差异基因KEGG富集散点图，KEGG富集程度通过Rich factor、Qvalue和富集到此通路上的基因数量来衡量。其中，Rich factor指差异表达的基因中位于该Pathway条目的基因数量与所有有注释基因中位于该Pathway条目的基因总数的比值。Q-value是做过多重假设检验校正后的P-value，Q-value的取值范围为[0，1]，越接近于0，表示富集越显著，笔者所在课题组挑选富集最显著的20条Pathway在图6中展示，不足20条的则全部展示。由图6可知，肝脏vs卵巢和卵巢vs血液中，差异表达基因富集数量在细胞周期和核糖体Pathway都较多；卵巢vs血液和肝脏vs血液中，差异表达基因富集数量在内质网蛋白质加工Pathway富集较多，肝脏vs卵巢在错配修复Pathway富集数量最少，卵巢vs血液在脂肪酸延长Pathway富集数量最少，肝脏vs血液在药物代谢-细胞色素P450 Pathway富集数量最少。

2.4 SNP和SSR统计结果分析

统计分析结果显示，在鲫鱼卵巢组织中检测到的SNP位点最多(348 795个)，其次为血液(343 516个)，肝脏组织中最少(307 740个)。其中，肝脏和卵巢组织编码区SNP数量多于非编码区，而血液中非编码区SNP数量多于编码区数量。由表2可知，3个组织中同义突变数量都多于非同义突变数量。

表2 SNP统计

利用MISA从鲫鱼127 801条unigene中共检测到48 808个SSR位点，其中，具有SSR位点的unigene 32 769个，SSR位点大于1的unigene 10 260个。不同类型SSR出现频率不同，由图7可知，单碱基、二碱基和三碱基重复类型所占比例较高，其他3种重复类型所占比例较少。其中，单碱基重复类型中以9～12次重复数目的SSR最多，其次为 13～16次重复数目的SSR。二碱基重复类型中以 5～8次重复数目的SSR最多，9～12次重复数目的次之。三碱基、四碱基、五碱基和六碱基重复类型中均以5～8次重复数目的SSR最多。且发现二碱基重复基元中，AC/GT含量最高，其次为AT/AT，三碱基重复基元中，AAT/ATT含量最高，ATC/ATG次之。

3 讨论与结论

鲫鱼不仅是一种重要的水产养殖物种，而且还是研究进化基因组学和生理性适应机制的模式生物，近年来，国内外对其开展了大量研究，获取了不同组织(嗅觉上皮、脑、肌肉、肝脏、肾脏和皮肤)的转录组信息。为进一步丰富其基因组数据，本研究对其肝脏、血液和卵巢组织进行了高通量转录组测序。

3.1 转录组质量

本研究利用Trinity软件对所得的reads片段进行组装，得到127 801条unigene，所得unigene的长度为201～17 286 bp，平均长度为735 bp。其N50为1 288 bp，N50值越大说明组织得到的长片段越多，组装效果越好。Q30值在80%以上就认为测序质量非常可靠，本研究中鲫鱼肝脏、卵巢和血液转录组Q30值分别为90.42%、89.59%和90.21%，均大于80%，表明本研究构建的转录组数据库准确可信，可以为后续鲫鱼基因克隆及功能基因验证提供基础数据。

3.2 数据库注释与差异表达基因

数据库注释结果显示，鲫鱼3个组织转录组拼接组装的127 801条unigene在NR、NT、KO、Swiss-Prot、GO、KOG和KEGG公共数据库中均得到注释，有助于进一步深入了解基因的功能。其中GO、KOG和KEGG注释结果表明，肝脏、卵巢和血液组织中除了参与一般功能预测外，参与信号传导机制功能的基因数量最多，其次是内分泌系统和免疫系统，这可能与肝脏肩负免疫防御功能、血液是体液调节的联系媒介，机体分泌的激素、酶和维生素等物质也是依靠血液传递才能发挥对代谢的调节作用，及与类淋巴细胞的免疫防御功能、卵巢主要承担雌性生殖与内分泌两大重要功能等密切相关。

不同组织的转录组表达基因量往往存在差异，如红鳍东方鲀鳔的转录组表达基因量高于鳃的转录组表达基因量，日本七鳃鳗肝脏组织转录组表达基因量则低于血液组织转录组表达基因量，本研究中鲫鱼肝脏、血液和卵巢组织的基因表达量分布显示，血液的表达量最高，这可能与血液是体液调节的联系媒介，参与机体许多生命活动的调节过程有关。差异表达基因数量显示，卵巢和血液差异表达基因数量最多，肝脏和血液差异表达基因数量最少，这可能是由于卵巢主要参与雌性生殖的调控，与血液在功能上重合的较少，二者表达的基因存在较大差异。肝脏是机体的代谢中心，而调节机体代谢过程的一些物质主要依靠血液传递才能发挥对代谢的调节作用，所以肝脏和血液在功能上有较多的重叠，二者表达的基因也相似。一些差异表达基因则与组织的特异功能息息相关，如血液除了传递调节代谢的物质外，还主要承担着气体运输、氧气结合和运输、铁离子结合和亚铁血红素结合的特异功能。

KEGG注释结果显示，肝脏vs卵巢和卵巢vs血液差异主要体现在细胞周期和核糖体Pathway，推测这些活动相关Pathway主要与卵巢组织中雌性生殖细胞的产生密切相关。卵巢vs血液和肝脏vs血液差异主要体现在内质网蛋白质加工Pathway，推测其主要与血液中运输的蛋白质物质加工过程相关。此外，肝脏和卵巢可能共同参与错配修复Pathway，卵巢和血液共同参与脂肪酸延长Pathway，肝脏和血液共同参与药物代谢-细胞色素P450 Pathway。

3.3 分子标记

SNP和SSR是利用转录组数据开发最多的2类标记，本研究依据转录组数据发现肝脏组织中307 740个、卵巢组织中348 795个和血液组织中343 516个SNP位点，此外，利用MISA软件在鲫鱼3个组织中共检测到48 808个SSR位点。与传统方法相比较，转录组数据可以直接反映基因的表达情况，利用转录组测序数据挖掘和开发与性状相关的优势基因型与分子标记更省时、高效，为选育提供遗传依据。因此，本研究将为鲫鱼今后的多态性检测、群体遗传多样性分析及分子鉴定等方面打下基础。

本研究采用高通量测序技术，对鲫鱼的肝脏、血液和卵巢转录组进行测序和分析。经组装后最终获得有注释信息的unigene共117 414条，不同组织间比较结果显示，卵巢与血液差异表达基因数量最多，肝脏与血液的差异表达基因数量最少。此外，在鲫鱼肝脏、血液和卵巢组织中分别检测到 307 740、343 516、348 795个SNP位点，共检测到 48 808 个SSR位点。研究结果为进一步克隆和挖掘鲫鱼功能基因、多态性检测及群体遗传多样性分析以及探究鲫鱼耐低氧分子机制等方面研究奠定了基础。