基于简化基因组测序的红罗非鱼低温体色变异全基因组关联分析

徐鸿飞，陈 诏*，黄彩林，袁宗伟，赵何勇，李 华，杨宾兰，周大颜，苏换换，朱华平

(1.广西壮族自治区水产引育种中心，广西 南宁 530031；2.中国水产科学研究院珠江水产研究所/农业农村部热带亚热带水产资源利用与养殖重点实验室，广东 广州 510380)

【研究意义】红罗非鱼(Red tilapia，Oreochromisspp.)隶属于鲈形目丽鱼科，由罗非鱼属间杂交选育获得，其体色是由色素细胞分泌的色素及遗传、环境、营养和生理等多种因素共同作用的结果。红罗非鱼肉质细嫩、营养丰富[1-2]，体色艳丽，无肌间刺且腹腔无黑膜，近年来深受消费者青睐，在我国和东南亚等国家的养殖面积不断扩大[3]。但在遗传育种和越冬过程中，红罗非鱼的肤色会发生变异，进而制约了其养殖产业的良性发展[4]。简化基因组测序(RAD-Seq)因具有酶切DNA片段大小一致、文库构建简单快速、标签密度易于调节、成本低廉和准确性高等优点已得到广泛应用[5]；而全基因组关联分析(GWAS)在揭示复杂性状的数量遗传变异规律、解析关键基因的遗传调控机制及推动分子辅助育种等方面具有广阔前景[6-9]。因此，基于RAD-Seq和GWAS技术分析红罗非鱼低温体色变异的分子机制，筛选出与体色变异相关的分子标记，对扩大我国红罗非鱼养殖规模及提升其养殖效益具有重要意义。【前人研究进展】传统的遗传育种主要根据不同表型性状直接进行选择[10]，是一个耗时费力的繁杂过程。随着分子生物学和基因组学的发展，鱼类遗传育种研究正经历由传统选择育种、杂交育种到基于基因组信息分子育种的转变[11-12]。Vignal等[13]、Liu和Cordes[14]、全迎春等[15]研究认为，基于基因组分子育种的核心内容是进行基因型和表型多态性研究，而单核苷酸多态性(SNP)因具有分布广、可实现高通量检测等优点已逐渐成为主流分子标记。Tsai等[16]基于高密度SNP序列开展了大西洋鲑(Salmosalar)幼鱼生长性状的GWAS和基因组预测分析。RAD-Seq主要包括复杂度降低的多态序列(CRoPS)测序[17]、限制性酶切位点相关的DNA(RAD)测序[18]和基因分型测序(GBS)[19]，其中RAD分子标记测序技术已在多个物种中得到广泛应用，通过RAD测序能在大多数生物中获得成千上万的SNP分子标记[20-21]。Houston等[22]利用RAD-Seq技术对传染性胰脏坏死病病毒(IPN)感染抗性和敏感的大西洋鲑亲本及14个子代进行基因分型，并构建了遗传连锁图谱。Baird等[23]利用RAD-Seq技术对3种链胞菌属(Neurosporacrassa)进行基因分型，发现13000多个SNPs位点。Carmichael等[24]研究显示，利用RAD-Seq技术检测到鲑鱼虱(Lepeophtheirussalmonis)不同性别群体中存在Lsa101901标记位点，该位点在雌性中为杂合子，在雄性中为纯合子。GWAS则在动植物复杂性状研究中发挥着重要作用，可从全基因组范围内筛选出与目标性状相关联的位点[25]，目前已在大西洋鲑(SalmosalarL.)[6]、虹鳟(Oncorhynchusmykiss)[7]、鲤(Cyprinuscarpio)[8]和鲶(Ictalurusspp.)[9]等多种水生动物中得到应用。【本研究切入点】目前，从分子水平上揭示红罗非鱼体色变异机理的研究主要集中于转录组[3]和红斑数量性状基因座位点绘制[26]，而有关红罗非鱼在低温后体色性状发生变异的分子标记筛选鲜见研究报道。【拟解决的关键问题】在RAD-Seq的基础上，结合GWAS策略在全基因组范围内寻找与红罗非鱼低温处理后体色变异性状相关的SNP位点，旨在筛选出与红罗非鱼低温体色变异相关的遗传标记，为红罗非鱼体色变异的改良与优良品种选育提供基础研究资料。

1 材料与方法

1.1 试验材料

试验所用200尾红罗非鱼由广西壮族自治区水产引育种中心桂虹罗非鱼良种场提供，平均体重50.11 g/尾，其亲本为广西壮族自治区水产引育种中心桂虹罗非鱼良种场的普通生产群体。试验红罗非鱼在室内长方形水泥池(7.6 m×2.3 m×0.8 m)暂养7 d后开始降温试验，暂养期间水温保持25 ℃，正常喂食，每2 d换水1次。暂养7结束后，以每天降低1 ℃的速率进行缓慢降温，降至15 ℃时维持7 d。7 d后选取体色变化明显和体色不变的红罗非鱼各30尾，利用100 mg/L MS-222(美国Sigma MO公司)进行麻醉处理，刮去鳞片后取其皮肤组织投入液氮中速冻，-80 ℃保存备用。

1.2 基因组DNA提取与检测

根据传统的苯酚∶氯仿∶异戊醇法对红罗非鱼皮肤组织样品进行基因组DNA提取，并利用酶标仪和Qubit Fluorometer检测基因组DNA浓度，采用1.0 % 琼脂糖凝胶电泳检测基因组DNA的完整性，用微量分光光度计测定其质量，确定基因组DNA纯度，并置于-20 ℃冰箱保存备用。

1.3 RAD文库构建及测序

RAD文库的构建及测序委托深圳华大基因股份有限公司完成。操作步骤如下：取1 μg基因组DNA，使用EcoR I(NEB)进行酶切。对得到的酶切片段进行5′端修复和磷酸化处理，3′端加poly(A)，连接P1 Index接头，用琼脂糖凝胶电泳筛选300～500 bp目的DNA，参照QIAquick PCR Purification Kit试剂盒说明进行纯化回收，对回收产物进行末端修复、末端加poly(A)及P2接头连接，然后对连接产物进行PCR扩增。PCR产物经琼脂糖凝胶电泳后，切取目的片段，参照QIAquick Gel Extraction Kit(QIAGEN)试剂盒说明进行纯化回收，适量Elution Buffer溶解，完成文库构建。构建好的RAD文库使用Agilent 2100 Bioanalyzer和ABI StepOnePlus Real-Time PCR System进行质量和产量检测，合格文库在Illumina HiseqTM4000平台上进行测序。

1.4 数据过滤与SNP分型

对测序原始数据(Raw data)中的低质量序列和接头序列进行剪切和过滤得到有效数据(Clean data)。以BWA为比对工具，将测序数据与参考基因组(https://www.ncbi.nlm.nih.gov/assembly/GCA_001858045.3)进行在线比对分析；然后应用GATK进行SNP位点筛选和基因分型，并采用BGI自主研发的软件对其进行注释和统计分析。对多个样品的SNP位点筛选条件为：同一SNP位点上SNP缺失率<50 % ；同一SNP位点上每个样品支持该位点的Reads数≥5；同一SNP位点的碱基类型数≥2。

1.5 GWAS分析

应用一般线性模型(GLM)和混合线性模型(MLM)算法进行GWAS分析，并用R语言绘制Q-Q图(Quantile-Quantile Plot)和曼哈顿图。其中，曼哈顿图中的-log(P)阈值选取α=0.05/SNP的个数，图中的虚线为y=-lg(α)；Q-Q图主要用来估计数量性状观测值与预测值间的差异，获得的数量性状数据均为正态分布数据。Q-Q图的X和Y轴分别代表SNP标记位点的预测-lg(P)值和观察-lg(P)值。

2 结果与分析

2.1 测序数据质量及比对结果

通过Illumina HiSeqTM4000平台测序得到的Raw data经Barcode拆分、过滤获得Clean data。在所有Reads中质量值大于20的碱基占总Reads长度的比例均大于97.44 % ，质量值大于30的碱基占总Reads长度的比例均大于94.07 % ，因此可判断测序数据质量值高，可用于后续的信息分析。以尼罗罗非鱼(Oreochromisniloticus)的基因组序列作为参考基因组，应用BWA软件将所有测序数据与参考基因组进行比对，结果发现所有样本的比对率均介于99.34 %～99.45 % ，说明测序数据质量可靠，可进行后续数据分析和试验操作。

2.2 SNP分型结果

通过BWA与参考基因组进行比对分析后，利用GATK软件进行SNP位点筛选和基因分型，结果共得到1823 228个SNPs位点。根据这些SNPs位点的完整度进一步筛选，得到46 889个高质量SNPs位点可用于GWAS。

2.3 群体结构分析结果

通过Admixture软件分析红罗非鱼的群体结构，根据交叉验证错误率的估值确定最优分群数(K)，结果(图1)表明，K均在0.6以上，说明所分析的样品来自同一个原始祖先，即样品不存在明显的分群，适合进行后续的关联分析。从混合线性模型下群体的Q-Q图(图2)可看出，实际观测值与预测值在前期基本相符，在后期产生分离，实际观测值曲线逐渐翘起，说明所选用的分析模型正确，关联分析结果可靠。

2.4 GWAS分析结果

采用混合线性模型对红罗非鱼低温处理后体色的变异性状进行GWAS分析，并绘制曼哈顿图。从图3可看出，GWAS分析共筛选出24个与红罗非鱼低温体色变异性状显著相关的SNPs位点(P<5.15E-08)，其中有15个SNPs位点位于1号染色体上，有9个SNPs位点位于12号染色体上，且具有成簇分布的特点。取SNPs位点上、下游各1000 bp序列，采用BLAST程序与尼罗罗非鱼基因组进行在线比对，获得每个SNP位点周围的基因，结果显示除5个未注释到基因及蛋白功能重复的基因外，共发现19个已知功能的基因，这些基因主要与信号转导、合成和免疫等功能相关(表1)。

表1 与红罗非鱼低温体色性状显著关联SNPs位点基因的注释情况Table 1 Gene annotation of SNPs significantly associated with body color trait in red tilapia after low temperature treament

3 讨 论

SNP作为一种新兴的分子标记，在草鱼(Ctenopharyngodonidella)[27]、尖吻鲈(Latescalcari)[28]、大口黑鲈(Micropterussalmoides)[29]和牙鲆(Paralichthysolivaceu)[30]等多种鱼类中得到广泛应用，是一种有效的分子辅助育种手段，也是优良品种选育的基础。本研究通过对降温条件下体色变化明显与体色不变各30尾红罗非鱼进行RAD-Seq，共筛选到24个与红罗非鱼低温体色变异性状显著相关的SNPs位点，这些SNP位点主要分布在1号染色体和12号染色体上(在其他染色体上均未发现)，且表现出成簇分布的特点，与Li等[31]、Huang等[32]对其他物种基因组的关联分析结果相似，可能与低温体色变异相关的SNP位点主要集中在部分染色体上有关。

鱼类体色主要通过色素细胞及其所含色素而形成，且受神经系统和内分泌系统调节，能随着环境的改变而产生适应性变化。Li等[26]研究发现，马来西来红罗非鱼遗传性红斑形成相关的微卫星主要位于5号染色体和15号染色体上，虽然在这2条染色体上发现2个与色素合成相关的基因(SLC45A2和ASIP)，但这2个基因未在黑罗非鱼和红罗非鱼的皮肤显著表达，说明自身色素合成机制不一定能直接影响红罗非鱼体色变异，而外界环境刺激也可能对其体色变化产生影响。本研究通过对低温条件下与体色变异性状相关SNP位点附近的基因进行GWAS分析，除5个SNPs位点未注释到的基因外，共发现19个已知功能的基因，这些基因主要与信号转导、生物合成和免疫等功能相关，与Zhu等[4]在不同体色马来西亚红罗非鱼皮肤组织中的转录组测序结果较一致。此外，所注释到的基因均与信号转导、物质合成及免疫等相关，究其原因可能是温度降低对红罗非鱼是一种强烈的外界刺激。红罗非鱼在应对外界温度降低刺激时引起各种生理反应，导致其需调动大量的基因来应答外界刺激，以减少刺激对自身的影响，使与信号转导、合成及免疫应答等功能相关的基因表现较活跃。本研究在所注释到的基因中未找到直接与红罗非鱼体色变异相关的基因，可能与用于RAD-Seq的红罗非鱼群体数量较少有关，因此，下一步还需在更大的红罗非鱼群体中进行筛查。

4 结 论

给予RAD-Seq和GWAS分析共筛选出24个与红罗非鱼低温体色变异性状显著相关的SNPs位点，并注释到19个已知的功能基因，且均与信号转导、生物合成及免疫相关，可作为后续开展红罗非鱼体色相关分子辅助育种研究的参考依据。