锦鲤基因组数据分析及体色相关基因的筛选

史东杰 胡金有 朱华 张欣 李荣妮 孙砚胜

摘要：为了获得红白锦鲤的基因组信息，筛选与其肤色相关的基因，采用Illumina高通量测序技术对红白锦鲤皮肤组织的基因组进行测序，获得127.23 Gb clean data，Q20碱基比例在95.59%及以上，Q30碱基比例在90.81%及以上，GC含量为37.32%～42.38%，测序错误率为0.07。与鲤鱼基因组序列进行比对的结果显示，比对效率为 96.35%。研究共鉴定了1 048 576个SNPs（单核苷酸多态性），其中3.12百万～5.40百万个SNPs位于短reads比对不到的区域，其中变异位点位于外显子区域的有579 778个SNPs。SNP位点分布于锦鲤的50条染色体上，不包含scaffold（染色体骨架）。经ANNOVAR软件进行功能注释，纯合类型的SNPs数量是574 310个，杂合类型的SNPs数量是474 265个。SNPs位于基因间的数量最多，SNPs位于基因内的外显子区域的多态性最高。通过对8个重要候选基因注释的理解，发现微管蛋白LOC109046532、LOC109049213这2个基因与色素颗粒运输有关。其中基因LOC109046532含有突变，而另1个基因LOC109049213则不含有任何突变。8个候选基因都含有外显子SNP位点，但是没有发现终止密码子突变。

关键词：基因组重测序;锦鲤;体色基因;候选基因

中图分类号： S917  文献标志码： A

文章编号：1002-1302（2019）16-0052-04

收稿日期：2018-05-03

基金项目：北京市财政局、北京市农业农村局观赏鱼产业技术体系北京市创新团队建设专项（编号：BAIC03）;北京市农林科学院项目（编号：KJCX20170101）。

作者简介：史东杰（1985—），女，北京人，硕士，副研究员，主要从事观赏鱼繁育及养殖技术的研究工作。

通信作者：朱 华，博士，研究员，主要从事水产繁殖、养殖以及水产养殖环境水质调控方面的研究与推广工作。

全基因组重测序是对已知参考基因组序列的物种进行不同个体间的基因组测序，并在此基础上对个体或群体进胁迫行差异性分析[1]。近年来，随着测序技术的发展，人们已经在众多水产动物中开展了全基因组测序，目前，鲤鱼（Cyprinus carpio）[2]、大黄鱼（Larimichthys crocea）[3]、半滑舌鳎（Cynoglossus semilaevis）[4]、大西洋鲑（Salmo salar）[5]、鲇鱼（Ictalurus punctatus）[6]、凡纳滨对虾（Litopenaeus vannamei）[7]和牡蛎（Ostrea gigas Thunberg）[8]等的基因组计划已经完成。2011年，由中国水产科学研究院和中国科学院北京基因组研究所共同实施的“鲤鱼基因组计划”成功完成了鲤鱼的全基因组测序，并绘制了鲤鱼基因组框架图谱、基因组物理图谱和高密度连锁图谱，进而利用各方面的资源和数据实现了鲤鱼基因组的基因识别定位和精确的功能注释等。全基因组序列海量数据的获得，为水产基因组辅助育种研究、优良品种的快速培育提供了重要基础。

锦鲤（Cyprinus carpio L.）是经济合作与发展组织（OECD）规定的5种试验生物之一，也是我国主养的观赏鱼类。该鱼隶属于鲤形目（Cypriniformes）鲤科（Cyprinidae）鲤属（Cyprinus）。锦鲤以其雄健的身躯、绚丽的色彩、华丽的斑纹、潇洒的泳姿、温顺的习性而享誉世界，被人们称为“水中活宝石”。该鱼经过几百年的自然分化、基因突变、人工选育，形成了体色艳丽、斑纹丰富、鳞片迥异等十三大品系100余个品种，是目前鲤科鱼类种质资源和基因组资源最丰富的鱼类。本研究通过对锦鲤进行基因组重测序，与鲤鱼进行参考基因组比对，以期找到大量单核苷酸多态性位点（SNP）、拷贝数变异（copy number variation，简称CNV）、插入缺失（insertion/deletion，简称InDel）、结构变异（structure variation，简称SV）等变异信息，分析锦鲤与鲤鱼的遗传多样性，同时研究锦鲤是否有与驯化选择相关的差异位点，并在测序的基础上，筛选出与肤色相关的候选基因。本研究不仅对锦鲤基因组辅助育种研究、体色斑纹定向培育提供了重要基础，而且对鲤科鱼类的基础研究具有重大意义。

1 材料与方法

1.1 试验材料

试验用红白锦鲤来自观赏鱼产业技术体系北京市创新团队通州综合试验站。从生长状态良好的健康红白锦鲤成鱼上取适量皮肤组织样品（设3个生物学重复），采用TIANamp Genomic DNA Kit（血液/细胞/组织基因组DNA提取试剂盒）进行DNA提取，并通过琼脂糖凝胶电泳、NanoDrop检测和Qubit定量进行DNA样本的检测。取样前，采用MS-222（Sigma，USA）使试验鱼麻醉后安乐死，并根据我国在科学技术方面应用的法律法规人性化对待试验动物。

1.2 试验方法

基因组DNA利用Covaris破碎仪随机打断成长度为 350 bp 的片段，经末端修复和加A尾后，片段两端分别连接接头制备DNA文库。文库构建完成后，先使用Qubit 3.0进行初步定量，随后使用Qseq 100对文库的insert size（插入片段大小）进行检测，insert size符合预期后，使用Q-PCR方法对文库的有效浓度（2 nmol/L）进行准确定量，以保证文库的质量。库检合格后，根据文库的有效浓度及数据产出需求，进行Illumina HiSeq X Ten PE150測序。PE150（Pairend 150 bp）指高通量双端测序，每端各测150 bp。在构建的小片段文库中，insert DNA，即插入片段是高通量测序直接测序的单位。双端测序是对插入片段的两端进行测序的方法，由于插入片段的长度分布已知，双端测序时不仅可以知道片段两端的序列，也能知道这两段序列之间的长度，从而便于后续组装和比对。

对测序获得的reads数据进行质量过滤得到clean reads，用于后续生物信息学的分析。将clean reads与参考基因组进行比对，基于比对结果，使用samtools[9]进行去重复（mark duplicates），使用GATK[10]进行局部重比对（local realignment）、碱基质量值校正（base recalibration）等处理，再使用GATK进行单核苷酸多态性（single nucleotide polymorphism，简称SNP）的小片段插入缺失（small InDel）的检测、过滤，并得到最终的SNP和small InDel的位点集。通过BreakDancer[11]可以得到结构变异（structure variation，简称SV）数据集，其中一般以插入（insertion，简称INS）和缺失（deletion，简称DEL）为主。对SNP、InDel、SV、CNV的检测结果进行注释，从而实现DNA水平差异基因挖掘和差异基因功能注释等。

1.3 数据处理与分析

将下机数据进行过滤，得到clean data，将其与指定的参考基因组进行序列比对，得到mapped data，进行插入片段长度检验、随机性检验等文库质量评估;进行可变剪接分析、新基因发掘和基因结构优化等结构水平分析;根据基因在样品中的表达量进行差异表达分析、差异表达基因功能注释和功能富集等表达水平分析，从而筛选出与体色相关的功能基因。

2 结果与分析

2.1 红白锦鲤基因组重测序数据质量评估

共完成3个样品的基因组重测序分析，通过高通量测序法获得127.23 Gb clean data，Q20碱基的百分比在95.59%及以上，Q30碱基的百分比在90.81%及以上，GC含量为 37.32%～42.38%，测序错误率为0.07%。

2.2 红白锦鲤基因组与参考基因组的比对

将红白锦鲤皮肤样品的clean reads与指定的参考基因组进行序列比对（网址为ftp：//ftp.ncbi.nlm.nih.gov/genomes），比对软件选择BWA，结果显示，比对效率为9635%（表1）。

2.3 红白锦鲤基因组的SNP检测及注释

由图1、表2可知，利用重测序变异检测方法得到的结果显示，以鲤鱼基因组为参考，过滤掉测序深度在10X以下的位点，共鉴定了1 048 576个SNPs，其中312万～540万个SNPs位于短reads比对不到的区域，变异位点位于外显子区域的有579 778个SNPs。SNP位点分布于锦鲤的50条染色体上，不包含scaffold（染色体骨架）。用ANNOVAR软件进行功能注释，结果显示，纯合类型的SNPs数量为574 310个，杂合类型的SNPs数量为474 265个。SNPs位于基因间的数量最多，SNPs位于基因内外显子区域的多态性最高，由此可以看出，与鲤鱼相比，红白锦鲤的变异位点很多，且分布在染色体的各个位置。

2.4 红白锦鲤肤色相关基因注释及SNP分析

通过对8个重要候选基因注释的理解，发现微管蛋白的2个基因LOC109046532、LOC109049213与色素颗粒运输有关。其中基因LOC109046532含有突变，而另1个基因LOC109049213则不含有任何突变。8个候选基因都含有外显子SNP位点，但是没有发现终止密码子突变，详见图2、表3、表4。

3 讨论

在全基因组测序过程中，基因组DNA的提取和检测是关键。通常情况下，DNA的检测主要是通过NanoDrop检测DNA纯度（D260 nm/D280 nm值），用Qubit对DNA浓度进行精确定量[12]。其中D260 nm/D280 nm值在1.8～2.0之间，DNA浓度≥20 ng/μL，总量为1 μg以上的DNA样品被用来建库。在本试验中，红白锦鲤皮肤样品DNA的Q20碱基百分比在95.59%及以上，Q30碱基百分比在90.81%及以上，GC含量为37.32%～42.38%，测序错误率为0.07%，可见样品质量满足建库测序要求，且总量满足2次或者2次以上的建库需要。对测序获得的reads数据进行质量过滤得到clean reads，用于后续生物信息学的分析。将clean reads与参考基因组进行比对，基于比对结果，使用samtools[1]进行去重复（mark duplicates），用GATK[2]进行局部重比对、碱基质量值校正等处理，再使用GATK进行单核苷酸多态性的小片段插入缺失（small INDEL）的检测、过滤，并得到最终的SNP和Small INDEL的位点集。通过BreakDancer[3]可得到结构变异（structure variation，简称SV）数据集，其中一般以插入和缺失为主。并对SNP的检测结果进行注释，实现DNA水平差异基因的挖掘和筛选等。利用基因组比对软件BWA[1]，将过滤后的clean reads比对到参考基因组上，统计比对结果。对于重测序分析而言，比对率以及覆盖度指标能反映样本、建库及测序以及参考序列等的质量。在本试验中，将clean reads与鲤鱼参考基因组序列进行比对，结果显示，mapping率达到96.3%，说明测序样本与鲤鱼参考基因组的相似度很高。

SNP检测主要使用GATK软件工具包[2]。根据clean reads在参考基因组的定位结果，使用SAMtools[3]进行去重复（mark duplicates），使用GATK进行局部重比对、碱基质量值校正等预处理，以保证检测得到的SNP的准确性，再使用GATK进行单核苷酸多态性的检测、过滤，并得到最终的SNP位点集。SNP是通过ANNOVAR软件进行注释的。SNP分布图通过R语言ggplot2包进行绘制展示。在本试验中，将锦鲤测序数据比对到参考基因组上，以分析SNP位点的分布情况，为了使SNP连续显示，过滤去除了测序深度在10X以下的位点，共鉴定了1 048 576个SNPs，其中3.12百万～5.40百万个SNPs位于短reads比对不到的区域，其中变异位点位于外显子区域的有579 778个SNPs。此外，统计结果显示，SNPs位于基因間的数量最多，SNPs位于基因内的外显子区域的多态性最高，由此可以看出，与鲤鱼相比，红白锦鲤的变异位点很多，且分布在染色体的各个位置。此外，没有发现外显子SNP位点含有终止密码子突变，因此SNP位点并没有影响基因的正常编码和表达。可是就目前发现的SNP位点而言，由于鲤鱼基因组缺乏相应的SNP功能注释信息，无法看出SNP位点会对相应基因功能带来何种变化，可能需要进行进一步的功能验证试验。

鱼类细胞形态变化、定向运动、胞内物质（如色素颗粒）与“器官”的移迁（有丝分裂、减数分裂中的染色体极向移动）都与微管蛋白的聚合与解聚相关[13]。微管是由微管蛋白亚基组装而成的，每个微管蛋白亚基都是由2个非常相似的球状蛋白（α-微管蛋白和β-微管蛋白）结合而成的异二聚体，这种α-β二聚体是微管组装的基本结构单位[14]。鱼类体色的重要调控机制之一是通过微管蛋白对色素颗粒的靶向运输[15]。在本试验中，1、5号基因为微管蛋白基因，与色素颗粒运输有关。鲤鱼基因组gff的基因信息全部是由美国国立生物技术信息中心（NCBI）网站上Gnomon预测软件进行预测的结果，因此该基因组并没有完整、真实的数据来进行支撑。因此，由BLAST得到的这8个候选基因的名称都是以其在染色体上的位置进行命名的，至于其功能也是由预测软件进行功能注释的。

参考文献：

[1]Altshuler D，Pollara V J，Cowles C R，et al. An SNP map of the human genome generated by reduced representation shotgun sequencing[J]. Nature，2000，407（6803）：513-516.

[2]水 科. 鲤鱼全基因组序列图谱绘制完成[N]. 中国渔业报，2014-10-13（A03）.

[3]陈小明，李佳凯，王志勇，等. 基于简化基因组测序的大黄鱼耐高温性状全基因组关联分析[J]. 水生生物学报，2017，41（4）：735-740.

[4]刘 峰. 半滑舌鳎经济性状的遗传评估及基因组选择初步研究[D]. 上海：上海海洋大学，2015：37-40.

[5]Davidson W S，Koop B F. ICSASG international collaboration. Sequencing the Atlantic salmon （Salmo salar） genome the old fashioned way[R]. Plant & Animal Genomes XIX Conference，2011，San Diego，CA，USA：33-41.

[6]Liu J. Strategies for efficient assembly and annotation of the catfish whole genome sequence[R]. Plant & Animal Genomes XIX Conference，2011，San Diego，CA，USA：49-53.

[7]张晓军. 中国甲壳动物学会第十一届年会暨学术研讨会论文摘要集[C]//中国海洋湖沼学会甲壳动物学分会，中国动物学会甲壳动物学分会，2011：18-19.

[8]Zhang G F，Guo X M，Li L，et al. The oyster genome project：an update[C]// Ninth International Marine Biotechnology Conference. Qingdao，China，2010：371-379.

[9]Li H，Durbin R. Fast and accurate short read alignment with Burrows-Wheeler transform[J]. Bioinformatics，2009，25（14）：1754-1760.

[10]McKenna A，Hanna M，Banks E，et al. The genome analysis toolkit：a MapReduce framework for analyzing next-generation DNA sequencing data[J]. Genome Research，2010，20（9）：1297-1303.

[11]Wang K，Li M，Hakonarson H. ANNOVAR：functional annotation of genetic variants from high-throughput sequencing data[J]. Nucleic Acids Research，2010，38（16）：e164.

[12]莫惠棟，顾世梁. 基因组长度的估计方法[J]. 科学通报，2000，45（13）：1414-1419.

[13]尹云厚. 中药复方制剂对缺氧大鼠微管蛋白和驱动蛋白表达影响的研究[D]. 长春：中国人民解放军军需大学，2003：156-158.

[14]Hirokawa N，Takemura R. Kinesin superfamily proteins and their various functions and dynamics[J]. Experimental Cell Research，2004，301（1）：50-59.

[15]薛继鹏. 三聚氰胺、氧化鱼油和脂肪对瓦氏黄颡鱼生长和体色的影响[D]. 青岛：中国海洋大学，2011：125-128.