大菱鲆(Scophthalmus maximus)，俗称“欧洲比目鱼”或“多宝鱼”，是中国北部沿海地区重要的海水养殖经济鱼类。大菱鲆具有生长速度快、肉质鲜美以及抗逆性强等诸多优点，其养殖区域迅速扩展、产量产值不断攀升(雷霁霖等, 2002、2005)。然而，随着工厂化养殖规模和集约化程度的扩大，特别是种质资源的退化，成活率低、病害严重和生产缓慢等问题逐渐显现(李杰等, 2019; 王岚等, 2017)。因此，开展大菱鲆相关经济性状的选育研究，培育具有强抗逆性和高生长性能的新品种对大菱鲆产业发展具有重要意义。传统的育种方法是基于目标性状的表型值来开展选育工作(Hulata, 2001)，对于高生长性能品种选育，通常挑选符合标准的个体进行群体混养或是建立多个家系并选择优良的家系作为亲本用于进一步选育。然而，分子生物学和基因组学的快速发展则促使鱼类育种研究由传统选育和杂交育种向基于基因组信息的分子育种转变(Xu et al, 2015)，如以基因定位和数量性状位点(Quantitative trait locus, QTL)作图为基础的分子标记辅助选择(Marker-assisted selection, MAS)育种(Wang et al, 2015; 刘晓菲等, 2019)。利用目标性状与分子标记间的显著关联性，依据标记对后代个体进行全基因组选择或是综合标记与表型信息构建选择指数筛选个体，加速育种过程，提高选育效率。

随着测序成本的下降与高通量基因分型技术的迅猛发展，基于海量单核苷酸多态性(Single-nucleotide polymorphisms, SNP)标记的全基因组关联分析(Genome-wide association studies, GWAS)已取代利用稀疏遗传标记的QTL定位，成为鉴定与目标性状显著关联的遗传变异的主要方法。早期GWAS研究主要针对单一性状表型(Yu et al, 2006; Burton et al, 2007)，而在进行复杂性状的遗传解析且记录有多个相关性状表型值时，考虑多个性状联合信息的多性状GWAS策略通常要优于逐个性状与遗传位点间的关联分析(Kim et al, 2009; Bolormaa et al, 2010)。相较单一性状GWAS分析，表型具有潜在相关性的多性状GWAS分析更具优势，主要是因为：当不同性状间存在遗传相关时，因考虑了单一性状分析中忽略的性状间协方差，多性状分析提高了检验功效(Stephens, 2013)和参数估计的精度(Zhu et al, 2009)；在多性状GWAS分析过程中，数量性状核苷酸(Quantitative trait nucleotides, QTNs)的统计推断只进行一次统计检验，与逐个性状单独分析相比，降低了多重检验造成的误差(Klei et al, 2008)；存在基因多效性时，单个遗传变异与多个性状相关联，使得多性状GWAS分析更具生物学意义(Chavali et al, 2010)；Porter等(2017)综合比较了当下流行的多性状GWAS分析方法，发现大多数现存的多性状GWAS分析方法具有明显相似的统计检验功效，与单一性状GWAS相比，能够大幅增加显著遗传变异的检测效率。

本研究采用实验成本较低的简化基因组(2b- RAD)测序技术进行SNP分型，基于多性状线性混合模型对体重–体长和体长–体宽–尾柄宽共2种多性状相关表型组合进行全基因组关联分析，关联定位控制这些性状的一因多效QTNs，为大菱鲆的体重和体尺性状改良提供理论基础。此外，对逐个性状进行了单一性状线性混合模型GWAS分析，并将其结果与多性状GWAS分析结果进行比较。

实验用鱼共585条，取自29个全同胞家系组成的大菱鲆育种群体，4月龄时，对群体中所有大菱鲆个体注射电子标记，并提取鳍条组织DNA用于后续的2b-RAD测序。随后，将标记后的大菱鲆随机混养在2个6 m×6 m×1.5 m具有循环水养殖系统的水泥池中，控制水温在5℃~24℃，每天定时投喂商用饵料2次至饱食。实验用鱼按孵化日期至表型测量时的生长天数为275~1001 d不等。对实验用鱼进行11次不定期表型测量，每次测量前利用50 mg/ml MS-222鱼用安定剂麻醉待测大菱鲆个体，避免应激反应造成的鱼体损伤。利用电子秤称量每个个体的体重(Body weight, g, BM)性状，同时在统一的参考标尺下，用数码相机由固定高度向下垂直拍摄相应个体的体尺性状。根据拍摄图形，利用ImageJ软件标定每条大菱鲆的体长(Body length, cm, BL)、体宽(Body width, cm, BW)和尾柄宽(Caudal peduncle width, cm, CPW)共3个体尺性状的表型值。采用简化基因组2b-RAD高通量标记分型技术对具有表型测量记录的个体进行SNP分型，共获得30049个多态SNP分子标记，参考大菱鲆基因组(ASM318616v1)建立多态标记物理图谱。挑选生长周期为473 d左右(前后皆不超过5 d)大菱鲆个体的表型观测值作为多性状GWAS分析的表型值，使用PLINK v1.9 (http://www.cog-genomics.org/plink2)对相应样本个体的基因型数据进行质量控制。剔除低于90%最小检出率的个体以及最小检出频率低于95%、最小哈代温伯格平衡为1.0×10^–6、最小等位基因频率小于3%和方差变异大于0.05的SNPs，最终得到441个样本的23988个SNP标记用于全基因组关联分析。

本研究使用GEMMA软件(Zhou et al, 2014)中的多性状线性混合模型(Multivariate linear mixed model, mvLMM)进行多个性状表型值和SNPs标记的全基因组关联分析，所用模型为：

式中，Y为n×t维表型矩阵，n是样本个数，t是分析性状个数；W为n×c维协变量(非遗传固定效应，如性别、年龄等)矩阵，A为c维相应系数行向量，c是包含截距项在内的协变量个数；x为n维当前检验SNP的基因型指示变量列向量，β为t维当前检验SNP的加性遗传效应行向量；G为n×t维不包括当前检验SNP的剩余多基因效应矩阵，G~MN(0, V_g$ \otimes $K)，K是由全基因组SNP标记构建的基因组亲缘关系矩阵，V_g是t×t维剩余多基因方差–协方差矩阵；E为n×t维误差项矩阵，E~MN(0, V_e$ \otimes $I)，V_e是t×t维误差方差–协方差矩阵，I是单位矩阵。

最终，全基因组上每个SNP都能得到1个一因多效Wald统计量和相对应的统计概率P值。此外，还利用因子谱分解线性混合模型(Factored spectrally transformed linear mixed model, FaST-LMM) (Lippert et al, 2011)逐个性状进行单一性状的GWAS分析，检测与各个性状相关联的QTNs。使用R语言绘制曼哈顿和Quantile-Quantile (QQ)图，同时，统计用于判别群体分层影响大小的膨胀系数或称基因组控制(Genomic control, GC)值，在实际研究中GC值被定义为所有SNP标记的卡方统计量均值(Price et al, 2010)。

将质量控制后的基因型数据用PLINK软件(Chang et al, 2015)处理为分析所需的格式，再利用单性状约束最大法估计法(Restricted maximum likelihood, REML)逐个估计各性状的遗传参数。查找关联分析得到的显著SNP位点在全基因组上所处的物理位置，于大菱鲆全基因组(GCA_003186165.1)上选择距离其最近的候选基因进行注释分析。

选取生长天数为473 d左右的441条大菱鲆的体重、体长、体宽和尾柄宽4个性状进行联合GWAS分析。4个性状的表型频数分布见图 1。从图 1可以看出，各个性状的表型值基本都符合正态分布，具有一定的可靠性，适合后续的GWAS分析。各性状的原始表型的统计分析和遗传力估计值见表 1。此外，还对原始以及校正了性别、池子和测定日期等非遗传固定效应后的表型值进行了性状间的相关分析，发现各性状表型值间具有较强的相关性(表 2)。

图 1(Fig. 1)

图 1 4个性状的频数分布 Fig.1 Frequency distribution histograms of four traits

表 1(Tab.1)

表 1 大菱鲆体重和体尺性状的描述性统计 Tab.1 Descriptive statistics of body weight and morphological traits in turbot 参数Parameter 性状Trait 体重 Body weight (g) 体长 Body length (cm) 体宽 Body width (cm) 尾柄宽 Caudal peduncle width (cm) 最大值 Maximum 541.60 26.71 9.16 3.62 最小值 Minimum 74.40 13.69 4.45 1.58 平均值 Mean 245.72 18.69 6.53 2.53 标准误 SE 1.77×10–1 1.04×10–3 1.51×10–3 7.78×10–4 变异系数 VC 0.32 0.10 0.10 0.14 遗传力 Heritability 0.72 0.58 0.59 0.48

表 1 大菱鲆体重和体尺性状的描述性统计 Tab.1 Descriptive statistics of body weight and morphological traits in turbot

表 2(Tab.2)

表 2 大菱鲆体重和体尺性状相关程度 Tab.2 Degree of correlation for body mass and morphological traits in turbot 相关系数 Correlation coefficient 体重 Body weight 体长 Body length 体宽 Body width 尾柄宽 Caudal peduncle width 体重 Body weight 1 0.759 0.801 0.662 体长 Body length 0.834 1 0.825 0.772 体宽 Body width 0.807 0.832 1 0.729 尾柄宽 Caudal peduncle width 0.784 0.842 0.787 1 注：下三角中的值为性状原始表型间相关程度，上三角中的值为校正了非遗传固定效应后表型间的相关程度 Note: Values in lower-triangular matrix represent degrees of correlation among original phenotypes, and values in up-triangular matrix mean degrees of correlation among phenotypes which were adjusted by non-genetic fixed effects

表 2 大菱鲆体重和体尺性状相关程度 Tab.2 Degree of correlation for body mass and morphological traits in turbot

根据质量控制结果，441条大菱鲆的23988个SNPs被用于体重和体尺性状的联合GWAS分析，分别进行了BM-BL和BL-BW-CPW两个性状组合的多性状线性混合模型以及逐个性状的单一性状线性混合模型GWAS分析。对BM-BL两个性状的多性状GWAS分析后(图 2和表 3)发现，超过基因组显著水平(5% Bonferroni校正阈值，即2.084×10^–6)的QTN位点共9个，其中，各有1个显著QTN位于1、10和20号染色体上，即SNP_1_31495825、SNP_10_ 6893888和SNP_20_18773114，分别有2个显著QTNs位于3 (SNP_3_21601589与SNP_3_24007357)、5(SNP_5_25635891与SNP_5_26888833)和22(SNP_22_ 6190492与SNP_22_6380029)号染色体上。而在相同显著水准下，BL-BW-CPW 3个体尺性状的GWAS分析则在10号染色体上检测到2个显著的QTN位点SNP_10_2724296和SNP_10_6893888。2组多性状GWAS分析的GC值分别为BM-BL(1.035)和BL- BW-CPW(1.044)。

图 1(Fig. 1)

图 2 大菱鲆体重与体尺多性状GWAS曼哈顿和QQ图 Fig.2 Manhattan and Q-Q plots of multi-trait GWAS on body mass and morphological traits in turbot

表 3(Tab.3)

表 3 大菱鲆体重和体尺性状显著关联的QTNs信息 Tab.3 Information of QTNs significantly associated with body mass and morphological traits in turbot 性状 Trait SNP位点 SNP ID 等位基因 Allele 染色体 Chr 位置 Position P值 P value 附近基因 Nearby gene 区域 Region (bp) 多性状线性混合模型Multivariate Linear Mixed Models BM-BL SNP_1_31495825 G/A 1 31, 495, 825 1.78×10–6 GRIK2 31268951~31431331 SNP_3_21601589 C/G 3 21, 601, 589 1.09×10–6 cyp21a2 21601235~21604716 SNP_3_24007357 G/T 3 24, 007, 357 1.96×10–6 novel gene* 23968869~23970083 SNP_5_25635891 A/G 5 25, 635, 891 2.64×10–7 tshz3b 25624917~25661383 SNP_5_26888833 C/T 5 26, 888, 833 8.85×10–7 adamts17 26823608~26965578 SNP_10_6893888 A/G 10 6, 893, 888 3.91×10–7 SYT8 6878968~6899153 SNP_20_18773114 G/A 20 18, 773, 114 1.88×10–8 myt1la 18751768~18792541 SNP_22_6190492 T/C 22 6, 190, 492 7.73×10–7 psmd4a 6191068~6196105 SNP_22_6380029 T/C 22 6, 380, 029 4.98×10–7 novel gene* 6365881~6382399 BL-BW-CPW SNP_10_2724296 A/C 10 2, 724, 296 2.83×10–7 novel gene* 2717495~2721226 SNP_10_6893888 A/G 10 6, 893, 888 1.48×10–6 SYT8 6878968~6899153 单性状线性混合模型Single-trait Linear Mixed Model BM SNP_3_21601589 C/G 3 21, 601, 589 1.61×10–7 cyp21a2 21540626~21547626 SNP_5_26888833 C/T 5 26, 888, 833 5.94×10–7 adamts17 26823608~26965578 SNP_10_2724296 A/C 10 2, 724, 296 5.75×10–7 novel gene* 2717495~2721226 SNP_10_6893888 A/G 10 6, 893, 888 7.32×10–8 SYT8 6878968~6899153 BL SNP_6_6047375 G/A 6 6, 047, 375 1.42×10–6 grm4 5928125~6045591 BW SNP_8_3870447 G/A 8 3, 870, 447 1.98×10–6 BFAR 3863424~3871475 CPW – – – – – – – 注: *表示未在大菱鲆全基因组中进行功能注释的新基因 Note: * represents novel gene which has not found any gene function annotation in genome of S. maximus

表 3 大菱鲆体重和体尺性状显著关联的QTNs信息 Tab.3 Information of QTNs significantly associated with body mass and morphological traits in turbot

对体重和体尺的逐个性状单一性状GWAS分析结果(图 3和表 3)显示，在BM性状中，于3和5号染色体上各检测到1个显著的QTN位点，即SNP_3_21601589和SNP_5_26888833，10号染色体上检测到2个显著的QTN位点SNP_10_2724296和SNP_10_6893888。BL和BW性状分别在62和8号染色体上检测到1个显著的QTN位点SNP_6_ 6047375和SNP_8_3870447。而在CPW性状中则没有检测到显著的遗传变异位点。各性状单一性状GWAS分析的GC值分别为1.026 (BM)、1.021 (BL)、1.034 (BW)和1.031 (CPW)。

图 1(Fig. 1)

图 3 大菱鲆体重与体尺单一性状GWAS曼哈顿和QQ图 Fig.3 Manhattan and Q-Q plots of single-trait GWAS on body mass and morphological traits in turbot

以大菱鲆多性状与单一性状GWAS获得的显著SNP位点为探针，根据其在大菱鲆全基因组的位置，向上下游寻找距离最近的候选基因。在本研究检测到的17个QTNs附近共找到12个候选基因(表 3)，分别为GRIK2 (Glutamate receptor, ionotropic, kainate 2)、cyp21a2 (Cytochrome P450, family 21, subfamily A, polypeptide 2)、tshz3b (Teashirt zinc finger homeobox 3b)、adamts17 (ADAM metallopeptidase with thrombo- spondin type 1 motif, 17)、SYT8 (Synaptotagmin 8)、myt1la (Myelin transcription factor 1-like, a)、psmd4a (Proteasome 26S subunit, non-ATPase 4a)、grm4 (Glutamate receptor, metabotropic 4)、BFAR (Bifunctional apoptosis regulator)和3个新基因。对这些候选基因进行GO注释分析发现，细胞组分主要涉及核膜、细胞质膜及细胞核等；分子功能主要涉及离子通道活性、激素活性、氧化还原酶活性、肽酶活性、核苷酸结合、蛋白结合及G蛋白偶联受体活性等；而生物学进程则主要涉及离子转运、糖皮质激素生物合成过程、氧化还原过程、细胞信号转导、基因表达调控、蛋白酶解、胞外分泌、转录调控和G蛋白偶联受体信号通路等。其中，3个新基因在大菱鲆基因组中未有相关注释信息。

GWAS作为一种定位影响重要经济性状分子标记的有效方法，是以全基因组SNPs直接关联复杂性状表型，广泛应用于动植物分子标记开发和辅助选择育种等方面(Santana et al, 2015; He et al, 2017)。随着诸多水产动物基因组测序工作相继完成以及测序成本的下降，水产动物生长和抗病等复杂性状的GWAS分析多有研究报道。Gutierrez等(2015)对大西洋鲑(Salmo salar)在性成熟时的生长速率性状进行GWAS分析，检测到多个在大西洋鲑代谢和生长发育过程中发挥重要作用的候选基因。Zhou等(2019)利用GWAS揭示了半滑舌鳎(Cynoglossus semilaevis)抗病的遗传机制，共定位了与抗病性状显著关联的33个SNP位点，结合基因表达和甲基化分析，检测到数个影响半滑舌鳎抗病力的候选基因。Jiang等(2019)则对罗非鱼(Tilapia)的耐盐性状进行GWAS和QTL定位研究，为罗非鱼耐盐性遗传机理揭示及进一步的功能分析奠定了基础。目前，大菱鲆生长和体尺相关性状的GWAS分析鲜有应用，但QTL定位研究已在大菱鲆生长(Enrique et al, 2011)和气单胞菌(Aeromonas salmonicida)耐药性(Silvia et al, 2011)等相关性状中开展。

现有的绝大多数GWAS研究都是关联全基因组上SNP位点与单个性状表型值，即使当多个相关表型性状存在时，往往也会选择逐个性状分别进行单一性状GWAS分析(Stephens, 2013)。在多个相关表型性状条件下，单一性状关联定位方法忽略了QTNs对多个表型性状的共同影响。相较于单一性状关联分析方法，多性状GWAS方法具有更高的统计功效和参数估计及定位精确度。通过比较大菱鲆体重和体尺多性状GWAS与各个性状逐个进行单一性状GWAS的定位结果发现，在相同的显著水准下，前者检测到的QTN数多于后者(表 3)。在BM-BL性状组合的多性状GWAS分析中，共检测到9个QTNs，而在相对应的BM和BL单一性状GWAS分析中，BM中检测到4个QTNs，BL中仅定位到1个QTN，且BM中检测到的4个QTNs中有3个与多性状GWAS分析得到的QTNs重复。同样，BL-BW-CPW性状组合多性状GWAS分析检测到的QTNs也多于3个性状各自进行单一性状GWAS分析得到的QTNs。除了能够检测到更多的QTNs，由于多性状GWAS分析检测到的QTN是一因多效QTN，相较于单一性状GWAS分析检测到的QTN在生物学和实际应用层面更有意义。

在大菱鲆全基因组上寻找距离每个QTN最近的基因，共找到12个候选基因，其中，有9个已知的功能基因及3个新基因。这些影响大菱鲆体重和体尺性状的候选基因在全基因组分布较为分散，没有明显集中区域。本研究为首次对大菱鲆体重和体尺性状进行全基因组关联分析，检测到在大菱鲆相关研究中尚未见报道的显著的已知功能基因。相关研究表明，GRIK2为离子型受体，它偶联离子通道并形成受体通道复合物，且与细胞缺血缺氧损伤有关(张冬梅等, 2010)。Eachus等(2017)和Weger等(2018)的研究均表明，cyp21a2在斑马鱼(Brachydanio rerio var)幼鱼的糖皮质激素生物合成中发挥着重要作用。Erickson等(2011)以斑马鱼为实验动物模型，发现tshz3b可能调控后脑中的Hox功能，而Hox功能则被认为与形态发生和器官形成有关联(Gair et al, 2003)。由此可见，该基因可能参与大菱鲆体形态发育过程。Myt1la可能在有丝分裂结束阶段发挥重要作用(Nakajima et al, 2008)，BFAR则能与p75NTR的蛋白相互作用并抑制p75NTR信号转导，而p75NTR可与高亲和力受体TrkA协同作用促进细胞增殖或与细胞内配体结合诱导细胞凋亡(李红梅等, 2011)。本研究找到的候选基因可能都是大菱鲆体重和体尺性状的重要候选功能基因，它们在大菱鲆生长发育过程中的影响有待进一步的功能实验验证。