全基因组关联分析(Genome-wide association study, GWAS)是应用全基因组范围内的大量分子标记(一般为SNP)，将标记基因型结合性状表型进行联合分析，统计每个标记与目标性状之间的关联性大小(一般用P值表示)，鉴定出与目标性状密切相关且具有特定功能和育种潜力的基因位点或分子标记，主要用于物种经济性状相关SNP分子标记以及功能基因的鉴定，从而达到缩短育种周期和提高育种效率的目的，目前已在畜禽等脊椎动物育种中广泛应用(Tavares et al, 2020; Müller et al, 2019; Cui et al, 2016; Zhang et al, 2019)。近年来，随着基因组高通量测序技术的发展及测序成本的降低，GWAS开始应用于水产养殖动物的育种研究，如在大黄鱼(Larimichthys crocea)、鲶鱼(Silurus asotus)、凡纳滨对虾(Litopeneaus vannamei)、龙胆石斑鱼(Epinephelus lanceolatus)、虾夷扇贝(Patinopecten yessoensis)等物种的生长性状关联SNP位点、候选基因的挖掘和鉴定(Zhou et al, 2019; Li et al, 2017; Yu et al, 2019; Wu et al, 2019; Ning et al, 2019)方面应用并取得了一定进展。但是，与陆生脊椎动物相比，GWAS在水产动物育种中的应用尚处于起步阶段。

黄条(Seriola lalandi)属鲈形目(Perciformes)、鲹科(Carangidae)、属(Seriola)，俗称“黄犍牛”、“黄金鲅”，是一种全球广泛分布且具有长距离洄游特性的大洋性经济鱼类。黄条体型大，生长迅速，肉质鲜嫩，具有较优的口感鲜度和较高的营养价值，品质可与三文鱼、金枪鱼媲美，国内外市场消费需求旺盛，是属鱼类中经济价值最高的一种。近年来，全球很多国家纷纷开展黄条的人工繁育和养殖(Sicuro et al, 2016)，目前，已在中国、日本、澳大利亚、新西兰和智利等国形成一定规模的养殖产业。为支撑养殖产业可持续发展，初步开展了黄条种质资源的相关研究，主要包括染色体核型和带型(史宝等, 2017; 刘永山等, 2018)和群体遗传学(Premachandra et al, 2017; Sepulveda et al, 2017; Nguyen et al, 2018a, 2018b)等研究，相关结果为黄条种质资源的科学利用和育种研究提供了遗传信息参考。我国于2017年在黄条人工繁育技术方面取得突破，培育出批量优质苗种，开启了我国黄条人工养殖的序幕。目前，本实验室正在开展黄条的群体遗传特性研究，利用2b-RAD简化基因组测序技术获得了大量SNP分子标记，对亲鱼和F₁代苗种的遗传多样性变动情况进行监测。生长性状在水产养殖业中具有重要的经济价值，直接和产量相关(Li et al, 2017)。黄条体型较大，生长速度快，前期研究发现，体高和体长是影响黄条体质量的2个重要性状，共同决定程度达42.88%(李荣等, 2017)。本研究选择全长和体质量这2个重要生长相关性状为目标，利用获得的SNP分子标记进行全长和体质量性状的GWAS分析，以期找到影响体质量和全长性状的关键SNP分子标记和候选功能基因，为黄条种质资源评价和持续利用提供遗传信息积累和技术依据。

1 材料与方法 1.1 实验材料

选择来自我国黄海种群养殖黄条 F₁代1~2龄成鱼共119尾，用卷尺测量其全长(cm)，精确到1 mm；电子天平测量其体质量(g)，精确到0.1 g。利用高压高温(180℃)灭菌后的手术剪，分别剪取样品的胸鳍，以95%无水乙醇浸润后，在4℃冰箱中保存待用。

1.2 基因组DNA提取与质量检测

使用QIAGEN公司生产的动物基因组DNA提取试剂盒(DP121221)，参照试剂盒使用说明，提取鳍条基因组DNA。用NanoDrop 2000分光光度计(Thermo, 美国)测定基因组DNA浓度，通过1%琼脂糖凝胶电泳检测DNA的完整性，通过A_{260 nm}/A_{280 nm}的比值来判断DNA的质量。将质检合格的DNA浓度稀释至100 ng/µl，于–20℃条件保存备用。

1.3 文库构建与测序

将≥200 ng的各样品基因组DNA采用IIB型限制性内切酶BsaXI进行酶切，酶切产物分别加入5组不同的接头，使用T4 DNA Ligase连接，然后PCR扩增连接产物，最后根据5组接头信息，将5个标签按顺序串联，连接产物添加barcode序列，混库，使用Illumina Hiseq测序平台对混合好的文库进行Paired-end测序。

1.4 数据分析 1.4.1 表型数据分析

使用R语言中fivenum()函数，对黄条体质量和全长表型数据进行描述性统计，计算其平均值、中位数、最小值、最大值、标准差和变异系数。

1.4.2 测序数据分析与SNP分型

Illumina HiSeq测序平台得到的原始图像数据文件经碱基识别转化为Raw Reads，过滤删除含有接头序列的Reads，得到Clean Reads，过滤删除含有N碱基比例大于8%的Reads，过滤删除低质量Reads(质量值低于Q30的碱基超过15%)；利用Pear (Zhang et al, 2014)软件(V0.9.6)将成对的Clean Reads拼接，提取出各样品对应的Reads，过滤删除不含酶切识别位点的Reads后，得到各样品的Enzyme Reads；利用电子酶切从参考基因组中提取含有酶切识别位点的标签，作为参考序列，利用SOAP软件将各样品的Enzyme Reads比对到参考序列上，主要参数为-r0–M4–v2 (-r0指唯一比对；–M4指最优比对；–v2指比对允许2个错配)，对比对到相同标签的reads聚类，得到unique标签深度，选择样品深度 > 3×且深度 < 500、标签长度为27 bp的标签，利用SOAP软件(V 2.21) (Li et al, 2008)将测序数据比对到参考序列，利用最大似然法(ML)进行位点的分型(Fu et al, 2013)，过程中使用的RAD分型软件包(RAD typing)，包含10余个软件组分，覆盖了从数据预处理至最终分型结果输出的全过程。

1.4.3 全基因组关联分析

使用EMMA eXpedited (EMMAX)高效混合模型(Kang et al, 2010)，通过方差分量方法进行SNP分子标记和表型性状的全基因组关联分析，所用模型：

$y = Xb + Ga + e$

式中，y为表型值；X为固定效应关联矩阵，b为固定效应向量，G为通过SNP标记计算得到的关系矩阵，a为随机加性遗传方差的参数，e为剩余效应的向量。

每个SNP位点能得到1个关联值P。对GWAS给出的P值划定2条显著性水平线，其中1条经Bonferroni校正P=0.05/N来确定全基因组显著性阈值(Bonferroni, 1936)，N为SNP标记的个数，2个性状经Bonferroni校正后显著关联阈值–lgP=5.726；另一条使用R软件包中的p.adjust()函数计算得到经FDR校正后的阈值，体质量性状潜在显著关联阈值–lgP= 4.091，全长性状潜在显著关联阈值–lgP=4.413，挑选Scaffold长度的前30使用R软件包的qqman绘制曼哈顿图，绘制QQ图对关联分析进行评价，判断关联分析结果是否可靠。

1.5 候选基因鉴定及功能分析

将筛选到的具有关联性的SNP位点上下游1 Mb范围内的碱基序列与GenBank数据库中已有的黄条参考基因组(https://www.ncbi.nlm.nih.gov/genome/? term=Seriola+lalandi+dorsalis)进行序列比对，使用SnpEff软件(Version 4.3T) (Cingolani et al, 2012)对得到的SNP位点进行注释，以确定SNP位点在基因元件的位置、对氨基酸的变化影响等找到距离SNP位点最近的基因。将所有关联基因与KEGG数据库比对，进行Pathway分析，并用超几何分布检验的方法计算每个Pathway条目中基因富集的显著性，公式如下：

$p = 1 - \mathop \sum \limits_{i = 0}^{m - 1} \frac{{\left({\begin{array}{*{20}{c}} M \\ i \end{array}} \right)\left({\begin{array}{*{20}{c}} {N - M} \\ {n - i} \end{array}} \right)}}{{\left({\begin{array}{*{20}{c}} N \\ n \end{array}} \right)}}$

式中，N为所有基因中具有KEGG注释的基因数目，n为N中差异表达基因中具有的KEGG注释的基因数目，M为所有基因中注释为某特定KEGG的基因数目，m为注释某特定KEGG的差异表达基因的数目。计算的结果会返回一个富集显著性的P值，小的P值表示基因在该Pathway中出现富集，当P≤0.05表示显著富集。

2 结果 2.1 表型性状描述性统计

根据本研究所用黄条样品的体质量和全长性状的描述性统计结果(表 1)可知，体质量范围为228~ 3000 g，变异系数为58.12，全长范围为27~67 cm，变异系数为19.49。统计分析显示，黄条体质量和全长的表型数据呈近正态分布。

2.2 SNP分型

对2b-RAD简化基因组测序数据按照以下指标进一步过滤。剔除所有样品中低于80%个体可以分型的位点；剔除MAF低于0.05的位点，剔除等位基因大于2的位点。最终，测序获得26665个SNP位点进行GWAS分析。

2.3 全基因组关联分析

利用R软件包分别绘制黄条体质量和全长性状GWAS分析的QQ图(图 1和图 2)，比较观测值(纵轴)和期望值(横轴)的一致性，散点分布与阈值实线在前端基本贴合一致，代表SNP假阳性较低，关联分析结果可靠。黄条体质量与全长性状GWAS分析的曼哈顿图分别见图 3和图 4。根据全基因组关联分析软件EMMAX方差分量法的计算和Bonferroni校正P值，共筛选到17个与体质量性状显著相关的SNP位点(表 2)。通过NCBI和Ensembl数据库比对，当体质量性状的显著性值P=3.49×10^–10时，有1个SNP位点，在基因组上找到1个基因LOC111655074；当P=1.20× 10^–8时，鉴定出15个显著关联SNP位点，找到15个关联基因，分别是ptpn14、LOC111651007、LOC111655681、LOC111658683、LOC111659677、gpr139、LOC111668306、LOC111670411、LOC111673333、LOC111645011、dag1、LOC111645458、poglut1、LOC111646889和LOC111648303；当P=6.29×10^–8时，定位了1个SNP位点，找到1个关联基因LOC111657860。全长性状未找到经Bonferroni校正后的阈值线以上的显著关联位点，但经FDR校正后找到全长性状潜在显著关联SNP位点12个(表 3)，可能的候选基因12个。

2.4 候选基因生物信息学分析

本研究共鉴定黄条17个体质量性状显著关联SNP位点，12个全长性状潜在显著关联SNP位点。通过NCBI和Ensembl数据库比对，分别找到17个体质量性状关联基因和12个全长性状关联基因。通过与KEGG数据库比对，进行Pathway分析，发现体质量和全长性状共5个关联基因参与的6条相关代谢通路(表 4)，主要与营养代谢、生长代谢、免疫调控、生殖调控等生物学过程密切相关；这些与体质量和全长性状关联的SNP位点和相关基因参与多个生物学进程，具有重要的调控作用，在今后遗传选育中可能具有较大的应用价值。

3 讨论

黄条作为一种全球分布的大洋性经济鱼类(Swart et al, 2016)，因其高营养和经济价值在国际上诸多国家开始了苗种繁育和养殖(Whatmore et al, 2013)。有学者分析了全球鱼养殖发展现状，认为种苗早期成活率低、遗传育种研究迟滞以及专用高效饲料缺乏是目前鱼养殖产业发展的主要瓶颈(Sicuro et al, 2016)，加之过度捕捞和环境退化加剧，近年来，全球黄条种质资源下降迅速，有必要加强黄条种质资源与遗传育种相关研究。目前，有关黄条遗传育种的研究较少，新西兰已开展了黄条人工选育研究并取得了一定进展(Symonds et al, 2014)。另外，日本学者Ohara等(2005)通过微卫星标记建立了黄条和五条(Seriola quinqueradiata)的遗传连锁图谱，为其遗传育种提供了重要工具，但未进行经济性状的关联分析。澳大利亚学者Nguyen等(2018b)利用SNP分子标记构建了黄条的高密度遗传连锁图谱，并通过全基因组关联分析，鉴定出6个与体质量性状显著关联的位点，为体质量性状控制的遗传基础解析和下一步育种应用提供了依据。本研究首次利用2b-RAD简化基因组测序技术获得了黄条SNP位点26665个，并在此基础上进行全长、体质量性状的GWAS分析，筛选到17个与体质量性状显著关联的SNP位点、12个与全长性状潜在显著关联的SNP位点，为黄条全长与体质量性状遗传基础解析和遗传育种提供了潜在可用的分子标记与候选基因。

GWAS依赖于连锁不平衡(Linkage Disequilibrium, LD)检测目标物种群体的遗传变异与性状之间的关联，然后通过统计基因型和表型的关联性大小筛选出影响显著的遗传变异，定位影响表型性状的重要数量性状位点(QTL)和候选基因，确定其遗传机制(陶林等, 2019; Naha et al, 2016)。本研究根据鉴定出与体质量和全长性状关联的SNP位点，在每个SNP上、下游1 Mb序列范围内扫描，共挖掘到17个体质量性状显著关联基因、12个全长性状潜在显著关联基因，但这些基因尚未在黄条研究中报道。已有研究表明，脊椎动物包括鱼类的生长过程是一个由多基因共同参与的级联式网络调控过程，除生长轴相关因子外，其他相关的基因、激素、代谢产物等通过内分泌、旁分泌的通路对生长起到调控作用，成为潜在的育种靶位点(Zhou et al, 2019)。Wang等(2019)通过克隆获得了黄条生长相关关键基因并验证了其在早期生长发育过程中的表达调控作用。但本研究未能关联体质量、全长性状至GH、IGF等经典生长调控功能基因。本研究定位到的所有黄条体质量、全长性状关联基因中，通过KEGG或文献检索，发现部分基因在其他动物生长性状调控过程中起重要作用。如adhesion G protein-coupled receptor B2-like、G protein-coupled receptor 139等均为G蛋白偶联受体家族成员，主要参与有丝分裂、肌肉收缩、离子通道调控、基因转录等生理过程(Raise et al, 2015)。Myosin-10-like和four and a half LIM domains protein 2-like都是重要的细胞骨架调控蛋白，参与信号转导、细胞黏附和细胞增殖(Matthias et al, 2006; Johannessen et al, 2006; Pi et al, 2008; Woolner et al, 2008)。Dag1通过ECM-receptor interaction通路发挥作用，细胞与细胞外基质相互作用导致直接或间接控制细胞粘附、迁移、分化、增殖和凋亡等。LIM结合结构域基因在鸡(Gallus gallus)脚重性状中显著关联并被证明是影响脚重的重要基因(陈则东等, 2016)。在已报道的其他鱼类生长性状全基因组关联分析研究中，生长性状大多关联到生长轴以外的重要功能基因，这些生长轴以外的基因都通过相关通路对生长起到间接调控作用(Zhou et al, 2019; Li et al, 2018; Yu et al, 2019; Wu et al, 2019; Ning et al, 2019)，提示除经典生长轴功能基因外，其他生长因子也会在鱼类生长性状选育中成为重要的靶向分子标记。

本研究发现了29个与黄条体质量、全长性状显著或潜在显著关联的SNP位点，并定位到29个显著关联基因。这些基因在其他物种中主要参与细胞分化增殖、免疫、组织器官发育等相关的代谢通路调控过程，因此，可能会在今后黄条选育中发挥重要的作用，但其具体的生理功能和遗传育种应用价值尚需研究和验证。全基因组关联分析(GWAS)已广泛应用于脊椎动物经济性状相关SNP位点筛选与鉴定，是进行QTL精细定位的重要工具，但在水产动物特别是鱼类上的研究报道较少。今后，随着鱼类基因组研究的深入发展和技术进步，GWAS技术有望在养殖鱼类经济性状解析和遗传育种中得到更加广泛的应用(陈松林等, 2019)。