全州禾花鲤皮肤转录组特征分析

许艺兰，文露婷，黄 姻，陈 忠，覃俊奇，潘贤辉，周康奇，林 勇，邓 潜，武 霞，杜雪松

（1.广西壮族自治区水产科学研究院 广西水产遗传育种与健康养殖重点实验室，广西 南宁 530021；2.广西壮族自治区兽医研究所 广西兽医生物技术重点实验室，广西 南宁 530001）

全州禾花鲤（Cyprinus carpiovar.Quanzhounensis）又名禾花鱼，分类上隶属鲤形目鲤科鲤属，为鲤鱼的一个变种，是我国传统稻田养鱼产业的主要品种之一，养殖区域主要分布在桂林市全州县及周边的兴安县、灌阳县等地，是广西北部地区稻鱼共生系统中形成和保存的特色养殖品种，已成为当地稻渔综合种养产业发展的支柱品种，在乡村产业发展中发挥着重要作用[1]。全州禾花鲤在体色上表现为乌褐色，背部颜色稍深，鳃盖和腹部肌肉半透明，是视觉鉴别的直接依据，具有重要经济价值，但随着全州禾花鲤稻田综合养殖规模的扩大，以及外来鲤鱼群体的进入，全州禾花鲤出现了体色混杂、退化的现象。

RNA-seq 技术是能够快速获得某一组织或器官在特定条件下细胞内所有转录组水平基因信息（mRNA、rRNA、tRNA 及其他ncRNA）的测序技术。该技术能在有参考或无参考基因组的情况下，对转录本测序分析，筛选潜在的差异基因，为解决各种生物学问题提供了丰富的基因信息。随着测序技术成本的降低和生物信息学的发展，转录组技术被广泛应用在水产动物的体色研究中，如湖栖鳍虾虎鱼（Gobiopterus lacustris）[2]、锦鲤（Cyprinus carpio）[3]、刺 参（Stichopus japonicus）[4]、细 锯 脂 鲤（Pristella maxillaris）[5]、红鲫（Carassius auratus）[6]、红鳍笛鲷（Lutjanus erythropterus）[7]等。目前，对全州禾花鲤体色性状的相关研究极少，多集中在生长性状遗传改良和肌肉营养品质等方面[8-12]，体色性状资料匮乏，为保护全州禾花鲤体色性状的特色资源，促进全州禾花鲤养殖可持续发展，亟须开展相关研究。为此，以全州禾花鲤为试验对象，建鲤（Cyprinus carpiovar.Jian）作为对照品种，取皮肤组织进行转录组测序分析，揭示全州禾花鲤皮肤转录组信息的特点，以期为全州禾花鲤的体色遗传机制研究奠定基础，为全州禾花鲤的选择育种和特色性状资源保护提供参考。

1 材料和方法

1.1 试验材料

试验鱼均采自广西农业良种海南南繁育种基地养殖的全州禾花鲤和建鲤群体。其中，全州禾花鲤群体为2017 年从广西壮族自治区全州县才湾镇引进的禾花鲤群体，建鲤群体为2019年从广西安桂水产养殖有限公司引进的商品系。

1.2 试验方法

1.2.1 皮肤样品总RNA 的提取 选取生长状况良好、体色典型的全州禾花鲤（HQ）、建鲤（JH，对照）各3尾，分别标识为HQ-1、HQ-2、HQ-3、JH-1、JH-2、JH-3，其体质量及体长信息见表1。用MS222 快速麻醉后取其背鳍下方到侧线上方区域的皮肤组织。使用Trizol Reagent 法提取皮肤样品的总RNA，通过琼脂糖凝胶电泳、Nanodrop 微量分光光度计、Agilent 2100 等检测样品RNA 的降解与被污染程度、纯度、完整性。

表1 全州禾花鲤和建鲤的体质量与体长Tab.1 Weight and body length of C.carpio var.Quanzhounensis and C.carpio var.Jian

1.2.2 cDNA 文库的构建及上机测序 质检合格的皮肤RNA 样品，经过mRNA 纯化、反转录、末端修复、加A 尾、连接adapter 测序接头，筛选200 bp 左右的cDNA，进行PCR 扩增，再使用AMPure XP beads纯化PCR 产物，获得cDNA 文库。cDNA 文库质检后，使用Illumina HiSeq 2500测序仪测序。

1.2.3 转录组测序数据质控与组装 上机测序直接获得的序列是原始序列数据，称为raw data 或raw reads，不能直接用于拼接、组装比对分析，需进行数据质控，即过滤低质量raw reads，获得高质量的clean reads。再使用HISAT 2 软件[13]对质控得到的clean reads与参考基因组比对，即全州禾花鲤、建鲤转录组测序的clean reads，分别与鲤鱼基因组（NCBI登录号：GCF_000951615.1）比对。

1.2.4 转录组数据分析 使用软件GATK 4[14]检测序列中的单核苷酸多态性位点（SNP）和插入缺失（InDel）位点，采用Stringtie 软件[15]计算能够比对上鲤鱼基因组的基因的表达数目。通过DESeq 2 软件[16]筛选2 个组之间差异表达的转录基因。筛选显著差异基因的标准为|log2FC|＞1、FDR＜0.05，其中log2FC＞0，为上调基因；log2FC＜0，为下调基因。以log2FC数值判断差异的大小情况，数值越大，差异倍数越大。利用R 软件包和KOBAS 软件进行差异基因的GO功能分析和KEGG通路分析。

1.2.5 差异表达基因的实时荧光定量PCR（qPCR）验证 采用qPCR 验证转录组数据的准确性，随机选 取8 个 显 著 差 异 基 因mapk14b、krt13、Myh7、alas2、tcf7l1a、mknk2、PSPH、ddit4，以eef1a基因作为内参基因[17-18]，使用Primer 6.0 软件设计基因引物（表2），每个基因设置3 个技术重复和3个试验重复。按照R223 逆转录试剂盒（诺维赞公司）说明书操 作，构 建20 μL 的 反 应 体 系：2×ChamQ SYBR qPCR Master Mix 10 μL、PCR Forward Primer（10 μmol/L）0.4 μL、PCR Reverse Primer（10 μmol/L）0.4 μL、cDNA 模板4.0 μL、ddH2O 5.2 μL。PCR 反应程序：预变性95 ℃90 s；变性95 ℃5 s，退火60 ℃15 s，延伸72 ℃20 s，共40 个循环。最后根据eef1a基因表达量标准化，按照2-△△Ct方法处理数据。

表2 实时荧光定量PCR引物序列Tab.2 Primer sequence of real-time quantitative PCR

2 结果与分析

2.1 全州禾花鲤皮肤转录组测序数据质量

使用Illumina HiSeq 2500 测序平台进行转录组测序，获得raw reads 37 131 866～58 798 394条，其中clean reads 数目占99.82%以上，测序碱基质量值达到Q30 以上水平的占94.96%以上，数据过滤后GC含量为51.00%～51.70%，全部可定位到基因组上的reads 数量占reads 比例75.03%以上（表3）。可见全州禾花鲤、建鲤的转录组测序质量较高，碱基组成相对平衡，可用于后续分析。

表3 全州禾花鲤和建鲤的皮肤转录组测序数据过滤统计结果Tab.3 Statistics of sequencing data filtering of C.carpio var.Quanzhounensis and C.carpio var.Jian

2.2 全州禾花鲤皮肤转录组基因类型统计与新基因注释

利用转录组数据鉴定新基因，将新基因合并比对上的已知基因进行类型统计，结果显示，全州禾花鲤和建鲤皮肤样本中检测到的基因数占参考基因组基因总数的62.98%～72.44%。利用Stringtie 软件重构转录本，在测序结果中挖掘到新基因10 537个，从各样本中挖掘到的新基因占总数的76.79%～84.32%。与公共数据库比对注释，其中Nr数据库比对注释8 975 个Unigene，注释比例最高，达85.17%；GO 数据库注释3 714 个Unigene，占35.24%；KEGG数据库注释3 466个Unigene，占32.89%。

2.3 全州禾花鲤皮肤转录组的突变检测

对全州禾花鲤、建鲤的皮肤RNA 测序序列进行SNP、InDel 变异位点检测，共检测到1 359 144 个变异位点，包括1 278 906 个SNP 位点和80 238 个InDel位点。将SNP位点按功能类型分类表明，同义突变（Synonymous SNV）类型占的比例最大（60.16%），其次是非同义突变（Nonsynonymous SNV，31.67%）（图1A）。对SNP 突变位点的位置进行统计，结果显示，突变位置在外显子区域的数量最多（35.04%），其次是基因间区（27.02%）、突变内含子区（14.07%）、3′UTR区（12.84%）（图1B）。

图1 全州禾花鲤与建鲤的转录本SNP功能分类（A）和位置（B）Fig.1 SNP functional classification（A）and location（B）of transcripts of C.carpio var.Quanzhounensis and C.carpio var.Jian

2.4 全州禾花鲤皮肤转录组的差异表达基因分析

2.4.1 差异表达基因统计 基于表达量信息，利用R 语言开展主成分分析（PCA）显示，全州禾花鲤样本与建鲤样本表现出各自聚集的分布，第一主成分对样品差异的贡献值为70.3%（图2A）。

图2 全州禾花鲤与建鲤皮肤转录组PCA分析（A）和显著差异基因的火山图（B）Fig.2 PCA analysis（A）and volcano map of significantly different gens（B）screened from C.carpio var.Quanzhounensis and C.carpio var.Jian

以|log2FC|＞1、FDR＜0.05 为标准筛选显著差异表达基因，共筛选到11 717 个显著差异基因。与建鲤皮肤中的转录表达基因相比，在全州禾花鲤皮肤中有5 133 个上调基因，6 584 个下调基因，各占总差异基因数目的43.81%、56.19%（图2B）。

2.4.2 差异表达基因的GO 分析 对所有差异表达转录本进行GO 分类注释（图3），应用超几何检验，查找与整个基因组背景相比显著富集的GO 条目（图4）。GO 分类注释结果显示，生物学过程（Biological process）富集在25 个二级分类，细胞组分（Cellular compoment）富集在19 个二级分类，分子功能（Molecular function）富集在12 个二级分类。在细 胞 过 程（Cellular process）、单 体 过 程（Singleorganism process）、代谢过程（Metabolic process）、细胞（Cell）、细胞部分（Cell part）、细胞器（Organelle）、结合（Binding）、催化活性（Catalytic activity）等主要二级分类中，下调基因数量显著高于上调基因数量。应用超几何检验，查找与整个基因组背景相比显著富集的GO 条目127个（P＜0.05），其中显著性排在前20 的GO 条目包括细胞内（GO：0005622）、细胞内部分（GO：0044424）、细胞器（GO：0043226）、细胞内细胞器（GO：0043229）、细胞（GO：0005623）、细胞部分（GO：0044464）、膜结合细胞器（GO：0043227）、细胞内膜结合细胞器（GO：0043231）、细胞器包膜（GO：0031967）、包膜（GO：0031975））等GO 条目（图4）。显著性排在前20 的GO 条目中，生物学过程（Biological process）分类GO 条目占55%（11 条），细胞组分（Cellular component）分类GO 条目占40%（8条），分子功能（Molecular function）分类GO 条目仅占5%（1 条）。以上结果表明，与建鲤皮肤相比，禾花鲤皮肤组织转录组差异在生物学过程和细胞组分方面更明显。

图3 全州禾花鲤与建鲤皮肤转录组中差异表达基因的GO 富集分类结果Fig.3 Histogram of GO enrichment classification of differentially expressed genes of C.carpio var.Quanzhounensis and C.carpio var.Jian

图4 全州禾花鲤与建鲤皮肤转录组中差异表达基因的GO富集结果Fig.4 GO enrichment circle map of differentially expressed genes of C.carpio var.Quanzhounensis and C.carpio var.Jian

2.4.3 差异表达基因的KEGG分析 将筛选的差异表达基因进行KEGG 通路富集分析，共富集到336个通路，其中显著富集的通路60 个（P＜0.05），包括氧化磷酸化途径（Oxidative phosphorylation）、细胞周期途径（Cell cycle）、代谢途径（Metabolic pathways）、DNA 复制途径（DNA replication）、氨基酸的生物合成途径（Biosynthesis of amino acids）等（图5）。对显著富集的通路进行网络分析并绘制网络图，发现网络图以代谢途径（ko01100：Metabolic pathways）、p53信号通路（ko04115：p53 signaling pathway）和细胞周期（ko04110：Cell cycle）通路为核心分成两部分（图6），表明试验组间转录组水平的调控差异主要集中在代谢活动和细胞分裂有关的生理过程。此外，精氨酸和脯氨酸代谢（ko00330：Arginine and proline metabolism）、氧 化 磷 酸 化（ko00190：Oxidative phosphorylation）、FoxO 信 号 通 路（ko04068：FoxO signaling pathway）、产热（ko04714：Thermogenesis）等通路相联的通路相对较多，在后续研究中应该关注。

图5 全州禾花鲤与建鲤皮肤转录组中差异表达基因的KEGG富集差异结果Fig.5 KEGG enrichment difference bubble chart of differentially expressed genes of C.carpio var.Quanzhounensis and C.carpio var.Jian

图6 全州禾花鲤与建鲤皮肤转录组中差异表达基因的KEGG富集分析网络Fig.6 KEGG enrichment analysis network of differentially expressed genes of C.carpio var.Quanzhounensis and C.carpio var.Jian

2.5 全州禾花鲤皮肤转录组差异表达基因验证

选取8 个差异表达基因进行定量分析，验证已获得的转录组数据的可靠性。结果显示，8 个差异表达基因的qPCR 定量结果与RNA-seq 结果在基因相对表达水平的上下调趋势上是一致的（图7），说明转录组测序数据准确，可信度高，能用于基因信息分析。

图7 全州禾花鲤与建鲤皮肤转录组中的差异表达基因qPCR结果分析Fig.7 Analysis of qPCR results of differentially expressed genes of C.carpio var.Quanzhounensis and C.carpio var.Jian

3 结论与讨论

3.1 全州禾花鲤皮肤转录组特征分析

全州禾花鲤作为一种地方性稻田养殖鲤鱼种群，其特殊的体色被市场作为确定经济价值的判断依据。深入解析全州禾花鲤体色形成机制，有利于保持其特色性状及后续的遗传育种研究。前期研究发现，该种群特殊体色的形成基础是体表鸟嘌呤结晶缺失和黑色素含量升高[19]。鱼类体色性状及其形成调控机制是受到广泛关注的研究领域，利用二代转录组测序技术，对斑马鱼（Danio rerio）[20]、红鳍笛鲷（Lutjanus erythropterus）[12]和锦鲤（Cyprinus carpio）[21]等品种皮肤转录组特征分析为体色调控机制研究提供了大量线索。本研究针对全州禾花鲤体色性状进行皮肤转录组特征分析，获得了raw reads数目37 131 866 ～58 798 394条，并发现新基因10 537 个，丰富了鲤鱼皮肤转录组数据，对基因表达种类和丰度信息的揭示为探索与禾花鲤体色相关的生物学问题提供了直接信息。同时，本研究挖掘1 359 144 个变异位点，包括1 278 906 个SNP 位点（94.1%）和80 238 个InDel 位点（5.9%），与三色锦鲤皮肤中变异位点相似[21]，但突变位点的类型、分布情况显著不同于与禾花鲤亲缘关系较远的黄尾鲴（Xenocypris davidiBleeker）、银 鲴（Xenocypris argentea）等物种[22-23]。SNP 是基因组上单个核苷酸突变造成的遗传标记，被广泛应用在水产动物的性状关联研究，关键基因中的SNP 位点可能与重要经济性状产生关联，如暗纹东方鲀（Fugu obscurusAbe）Leptin基 因[24]、大 口 黑 鲈（Micropterus salmoides）HSC70-1基因SNP 多态性位点[25]与其生长性状相关联。本研究获得的SNP 位点，可作备选标记位点用于全州禾花鲤体色、生长速度等经济性状研究，以及后续种群遗传结构、遗传多样性、品种鉴定、分子辅助育种等研究。

3.2 全州禾花鲤皮肤转录组差异表达基因富集分析

全州禾花鲤的特殊体色表型与兴国红鲤、瓯江彩鲤、锦鲤等品种的体色类型具有显著差异，从遗传角度分析，该群体应负载着有别于其他鲤鱼品种的与体表色素细胞的形成、分化和迁移过程具有紧密联系的体色调控多样性或基因突变资源。分析不同类型的皮肤组织转录组，可发现对体色具有显著影响的信号通路或基因。以体表色素细胞种类不同为变异特征的皮肤转录组研究表明，黑色素生物合成通路、酪氨酸代谢途径、Wnt和MAPK 信号通路影响皮肤的色素沉着[26]，与黑色素合成、酪氨酸代谢、蝶呤代谢等过程有关的一系列基因在不同体色类型皮肤组织中差异表达[27]。对色素细胞种类间的转录组分析显示，不同种类色素细胞的生成可能与嘌呤合成途径、糖酵解途径及磷酸戊糖途径有关[28]。以鸟嘌呤结晶缺失为体色变异性状的皮肤转录组研究表明，MPKA 通路、肌动蛋白细胞骨架调节通路、嘌呤代谢通路等可能与鸟嘌呤结晶缺失有关[29]。

本研究分析2 种鲤的皮肤转录组信息显示，11 717 条转录本表达量显著不同，并在60 条KEGG通路和127 个GO 条目显著富集，显示2 种皮肤在转录水平的重要差异。在氧化磷酸化、糖酵解/糖异生、能量代谢、戊糖磷酸途径、嘌呤代谢等KEGG 通路中，全州禾花鲤差异表达基因显著下调，可能是建鲤皮肤中分布有较多含有大量鸟嘌呤结晶的虹彩细胞，而糖酵解、氧化磷酸化等通路在鸟嘌呤合成中发挥重要作用[30-31]。在对显著富集KEGG 通路的网络分析中，代谢活动和细胞分裂有关的生理过程作为两大网络核心与众多KEGG通路关联且相对独立，表明全州禾花鲤和建鲤皮肤的物质和能量代谢存在差异，细胞周期调控过程也存在差异。在GO 富集分析中，最显著富集的条目多数与细胞有关，集中在细胞、细胞过程、细胞结构、细胞骨架等条目，提示全州禾花鲤和建鲤皮肤细胞分化、构成等方面的差异，也可能提示全州禾花鲤鸟嘌呤结晶的缺失与虹彩细胞的增殖分化有关[21，32]。以上结果与大颚细锯脂鲤（Pristella maxillaris）[33]、草金鱼（Carassius auratus）[29]及斑马鱼[20]透明突变体（虹彩细胞缺失/鸟嘌呤结晶缺失）皮肤中的结果差异较大，不能形成较为一致的富集结果，可能提示该性状的形成和调控机制多样。在此基础上开展深入研究，有可能发现和阐释新的体色调控机制。此外，皮肤是鱼类机体的物理屏障和非特异性免疫屏障，其生理功能对个体生存和适应环境极为重要[34]。全州禾花鲤是用于稻田养殖的鲤鱼群体，其养殖环境和运输条件与其他养殖形式不同，接受外源刺激更为复杂。本研究中，全州禾花鲤与建鲤皮肤的转录产物集合在表达量上差异明显，这是否会造成全州禾花鲤皮肤屏障、免疫功能等经济性状相关的生理过程差异值得深入研究，其结果将有助于全州禾花鲤性状鉴定、养殖性能评估和遗传改良工作的开展。