全基因组关联分析(genome-wide association study, GWAS)为全基因组范围内单个核苷酸变异所引起的 DNA 序列多态性即单核苷酸多态性(single nucleotide polymorphism, SNP)的检测提供了强有力的工具，它能鉴别与复杂性状相关联的数量性状核苷酸(quantitative trait nucleotides, QTNs), 实现目标性状的基因定位。目前，大多数的GWAS模型都是针对单个性状进行遗传位点和表型间的关联分析。当对个体进行多尺度表型测量时，共同的环境效应和连锁不平衡使得表型间具有潜在相关性，存在一个基因同时影响多个性状的情况，即一因多效性。因此当不同性状间存在遗传相关时，可以对多个潜在相关性状进行联合分析，由于考虑了性状表型间的协方差，与单性状分析相比会显著提升统计检验的效力，且具有一因多效性的QTN能够解释多个性状间的遗传亲缘相关，结果更具有生物学意义^[1,2,3]。

线性混合模型(linear mixed model, LMM)由于能够综合考虑群体分层，家系结构和复杂的亲缘关系等混杂因子已成为GWAS中广泛使用的分析模型。在多性状GWAS中，多变量线性混合模型(multivariate linear mixed model, mvLMM)可高效定位一因多效QTN, 提高多个相关性状中单个性状的检测效率^[4]。此外，mvLMM在估计组织间基因表达的遗传力^[5], 理解进化模型和辅助动物遗传育种中也发挥着重要作用^[1,6,7,8,9]。对于多个相关表型性状，在使用约束最大似然法估计多性状表型间的方差-协方差矩阵时，mvLMM的计算复杂度将比LMM增加性状个数倍。基于mvLMM的精确关联分析需要逐个SNP重复求解多性状混合模型，这对具有大量SNP和中等规模群体而言是不切实际的。GEMMA使用无效模型中相同方差-协方差矩阵替换每个SNP间不同的协方差矩阵，并在计算复杂度相当低的情况下估计mvLMM中的协方差矩阵，这大大降低了计算时间的消耗和计算机内存运行的负担。然而在更大规模群体的多性状GWAS分析中，对更多相关联表型性状进行更为复杂的矩阵运算时，GEMMA将表现出相对较低的关联统计检验效率。表型正交化方法是对这些相关表型进行降维的一种策略，通过对所分析的相关性状的表型协方差矩阵执行谱变换，对原始表型性状加权可以获得相互独立的“超表型”, 然后逐个性状进行关联分析^[10], 加和单个性状检验的卡方统计量，最后检验多个关联性状的一因多效QTN。

大菱鲆(Scophthalmus maximus, Linnaeus, 1758)是一种广泛分布在包括波罗的海、黑海和地中海的欧洲大西洋沿岸的经济鱼类^[11]。因其生长速度快，耐冷水能力强，味道鲜美，是世界上养殖最广泛的商业鱼类之一。我国自1992年引入以来，大菱鲆已经发展成为我国北方海水养殖产业的主要养殖种类^[12]。因此，探索大菱鲆的生长和耐低温等重要经济数量性状的关联基因或QTNs, 对深入解析其数量性状遗传机制以及提高标记辅助育种效率有着重要的意义。本研究通过比较GEMMA和表型正交化方法对大菱鲆体质量和体尺性状多性状全基因组关联分析结果，来评价表型正交化方法的QTN检测效率，为水产动物育种在多性状GWAS的研究提供更为高效便捷的方法思路。

本研究以来自29个全同胞家系共585尾大菱鲆组成的育种群体作为验证表型正交化方法性能的实际资料。于4月龄时，对养殖的育种群体中所有大菱鲆个体注射电子标记，并提取鳍条组织DNA用于标记测序。随后，在自然光周期下，采用循环水系统将已标记的个体混养在6 m×6 m× 1 m的水泥池中，每天投喂两次商用饵料至饱食，水温控制在5~24 ℃。自注射电子标记起，用电子秤定期称量每个个体体重，并在固定参考标尺下定期用数码相机从一定高度向下垂直拍摄个体的形态性状。到36月龄时，绝大多数个体均测量9次左右。每次测量前，利用50×10^–6 MS-222鱼用安定剂将待测个体麻醉，避免处理压力。根据拍摄图像，利用Image J软件获得不同日龄体长和尾柄宽的大菱鲆体尺表型测量值。采用简化基因组2b-RAD高通量标记分型技术，共获得 30049个多态性SNP标记，参考大菱鲆基因组(ASM318616v1)建立多态标记物理图谱。在其中挑选生长周期为473 d左右的大菱鲆个体的表型值，使用PLINK v1.9 (http://www.cog-genomics. org/plink2)对与表型个体所对应的基因型数据进行严格的质量控制，去除低于90%最小检出率的个体，以及最小检出频率低于95%、最小哈代温伯格平衡为1.0×10^–6、最小等位基因频率小于3%和方差变异大于0.05的SNPs, 最终得到 585个样本的21993个SNP标记。

为了去除性别、池子、测量日期的固定效应影响，利用R语言中的GLM模型对表型数据进行校正：

式中，y为个体表型测量值，Sex为性别，Pool为养殖池，Date为测量日期，e为剩余效应，本研究把剩余效应作为校正后的表型值Y用于后续的全基因组关联分析。

按下列公式计算：

其中，Y为n×d维的多性状表型值矩阵，n为分析的个体数，d为分析的性状数；b为d维当前检验SNP的加性遗传效应行向量；z为n维当前检验SNP的基因型指示变量列向量；G和E分别是n×d维的随机多基因效应和剩余误差矩阵。通常，G和E服从多变量正态分布：$G \sim N({\bf{0}},\;K \otimes {V_g})$, $E \sim $$N({\bf{0}},\;{I_n} \otimes {V_e})$。${V_g}$和${V_e}$分别是d×d维的多基因和剩余协方差矩阵，I_n是单位矩阵，K和$ \otimes $分别表示基因组亲缘关系矩阵和克罗内克积，满足$E(Y) = Z \otimes b$, $V(Y) = K \otimes {V_g} + {I_n} \otimes {V_e}$。

用R表示d个性状的表型协方差矩阵。将R谱分解为$R = {U_R}{S_R}U_R^T$, 令${Y^*} = Y{U_R}$,${b^*} = b{U_R}$,${G^*} = G{U_R}$和${E^*} = E{U_R}$则模型(1)可被转化为下式：

由$R{ = }{V_g} + {V_e}$可知，$V({G^*}{\bf{ + }}{E^*}) = U_R^T({V_g} + {V_e}){U_R} = $$U_R^{\rm{T}}R{U_R} = {S_R}$。令$V({G^*}{\bf{ + }}{E^*}) = ({S_R} - {I_d}) + {I_d}$, 则模型(2)变成相互独立的：

式中，$l = 1,\;2,\; \cdots ,\;d$, $g_l^*\~N(0,{s_R} - 1)$与$e_l^*\~N(0,1)$。

为了消除个体间的相关，把基因组亲缘关系矩阵谱分解$K = {U_k}{S_k}U_k^T$, ${S_k}$和${U_k}$是相对应的特征值和特征根矩阵。令 ${y'_l} = {({S_K} + {I_n})^{{\textstyle{1 \over 2}}}}U_k^Ty_l^*$, $z' = {({S_K} + {I_n})^{{\textstyle{1 \over 2}}}}U_k^Tz$, ${g'_l} = $${({S_K} + {I_n})^{{\textstyle{1 \over 2}}}}U_k^Tg_l^*$以及${e'_l} = {({S_K} + {I_n})^{{\textstyle{1 \over 2}}}}U_k^Te_l^*$, 然后把模型(4)转换为一般线性回归模型：

经最小二乘法求解，我们使用检验统计量推断标记的效应：

式中，${\hat b_l}$是估计的遗传效应，$V({\hat b_l})$是当前检验SNP的方差。

对d个独立的模型(3), 我们加和服从d自由度卡方分布的d个统计量(5)去检验一因多效QTN:

本研究选取生长天数为473 d左右的585尾大菱鲆的体重、体长、尾柄宽3个性状，分别采用GEMMA和表型正交化两种方法进行多性状GWAS分析。图1为大菱鲆3个生长性状的表型值正态分布图，从中可以看出，各个性状的表型值都符合正态分布，适合进行后续的多性状GWAS分析。表1为各性状的原始表型值最大值、最小值、平均值、标准误差、变异系数的统计分析和遗传力估计参数。表2为去除性别、池子、测量日期等固定效应后的表型值进行了性状间的协方差分析，发现所研究的3个大菱鲆生长性状具有较强的相关性。

图1

图1 　大菱鲆体重和体尺性状的的正态分布图 Fig. 1 Normal distribution of body mass and morphological traits of turbot

表1

表1 　大菱鲆体重和体尺性状的描述性统计 Tab. 1 Descriptive statistics of body mass and morphological traits in turbot 参数  parameter 体重/gbody mass 体长/cmbody length 尾柄宽/cmcaudal ped un cle width 最大值 maximum 676.40 27.01 3.90 最小值 minimum 74.40 13.69 1.58 平均值 mean 274.94 19.51 2.67 标准误 SE 88.69 1.85 0.37 变异系数 VC 0.32 0.10 0.14 遗传力 heritability 0.48 0.81 0.18

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

表2

表2 　大菱鲆体重和体尺性状间的相关 Tab. 2 Degree of correlation for body mass and morphological traits in turbot 相关系数correlation  coefficient 体重body weight 体长body length 尾柄宽caudal pedu n cle width 体重body weight 1 0.844 0.834 体长body length 0.844 1 0.805 尾柄宽 caudal peduncle width 0.834 0.805 1

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

在全基因组关联分析中，通常使用曼哈顿图和QQ图来反映所分析的所有标记的显著与否以及假阳性和假阴性的控制，图2显示，585个个体，21993个SNP大菱鲆体重和体尺性状数据进行GEMMA和表型正交化分析得到的基因定位结果。

由图2可知GEMMA方法超过基因组显著水平(5% Bonferroni校正阈值，2.273451×10^–6)的QTNs共11个位于1、5、6、8、12、16、20、21、22 号染色体上，而表型正交化方法共检测到14个QTNs位于1、3、5、6、8、12、16、20、21、22号染色体上。其中与GEMMA方法分别在1(SNP_1_31688963), 5(SNP_1_19132118), 8(SNP_ 8_14134711), 12(SNP_12_14910921), 16(SNP_16_ 11311257), 20(SNP_20_18773114), 21(SNP_21_ 808333), 22(SNP_22_6190492)号染色体上共检测到相同的QTNs有9个（表3）。表4中记录了GEMMA和表型正交化方法在分析不同性状个数时的运算时间。结果表明，在相同的内存条件下，表型正交化方法的运算时间远远低于GEMMA, 具有更高的关联检测效率。

动植物生长性状的研究分析能够直接有助于生产实践，因此这些生长性状是重要的经济性状，探索与大菱鲆重要经济性状关联的QTNs、候选基因，对解析水产动物的遗传机制，提高基因辅助育种效率有着重要意义。随着二代测序技术的迅猛发展，全基因组测序成本大大降低，有关水产动物经济性状的GWAS研究逐渐增多。Gutierrez等^[13]使用混合模型方法发现了与大西洋鲑性成熟时生长速率、年龄相关的主效QTN位点。Gonzalez-Pena 等^[14]使用单变量模型对虹鳟10、13 月龄的体重进行GWAS分析，检测到虹鳟在omy5染色体上检测到一个SNP位点可分别解释表型1.4%和1%变异。Jiang等^[15]使用QTL-seq的方法在距LG18染色体23.0Mb的位置上定位到与罗非鱼耐盐性相关的QTNs, 揭示了罗非鱼的耐盐性渗透调节遗传机制。目前，有关大菱鲆经济性状的定位分析大都采用QTL定位的方法，Sánchez-Molano等^[16]通过QTL筛选和微卫星搜索发现了与生长性状相关的共11个QTL。

图2

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

表3

表3 　大菱鲆体重和体尺性状显著关联的QTNs信息 Tab. 3 Information of QTNs significantly associated with body mass and morphological traits in turbot 性状trait SNP位点SNP ID 等位基因allele 染色体Chromosome 位置position -log10P 附近基因nearby gene 区域/bpregion GEMMA BM-BL-CPW SNP_1_31688963 C/A 1 31688963 8.788 Zgc:77056 29724588-30027306   SNP_5_19132118 G/A 5 19132118 10.787 novel gene* 19137312-19237777   SNP_5_27008012 A/G 5 27008012 6.143 novel gene* 27022392-27082668   SNP_6_21034563 T/G 6 21034563 6.141 novel gene* 20953184-21013608   SNP_8_14134711 G/C 8 14134711 7.834 novel gene* 14121245-14133456   SNP_12_14910921 G/T 12 14910921 5.930 Slc1a6 14903868-14913454   SNP_12_20732502 G/A 12 20732502 5.710 lpp 20675041-20804711   SNP_16_11311257 C/T 16 11311257 6.679 novel gene* 11336381-11348325   SNP_20_18773114 G/A 20 18773114 10.610 Myt1la 18751768-18792541   SNP_21_808333 T/C 21 808333 6.623 Zgc112356 804720-812315   SNP_22_6190492 T/C 21 619492 7.620 novel gene* 617712-621986 Fast-LMM   SNP_1_30018909 G/A 1 31688963 5.737 VPS13B 29724588-30027306   SNP_1_31036272 G/A 1 31688963 5.723 novel gene* 31063130-31066902   SNP_1_31688963 C/A 1 31688963 9.300 Zgc:77056 31688749- 31690775   SNP_1_31719031 C/A 1 31719031 5.781 novel gene* 31717594-31718666   SNP_3_20633231 C/A 3 20633231 5.701 novel gene* 20633231   SNP_5_19132118 G/A 5 19132118 11.502 novel gene* 19137312-19237777   SNP_6_4721035 T/G 6 4721035 6.223 novel gene* 4808453-4971126   SNP_6_21034563 A/G 6 21034563 6.220 novel gene* 20953184-21013608   SNP_8_14134711 G/C 8 14134711 8.380 novel gene* 14121245-14133456   SNP_12_14910921 G/T 12 14910921 5.679 Slc1a6 14903868-14913454   SNP_16_11311257 C/T 16 11311257 7.653 novel gene* 11336381-11348325   SNP_20_18773114 G/A 20 18773114 11.773 Myt1la 18751768-18792541   SNP_21_808333 T/C 21 808333 7.850 Zgc112356 804720-812315   SNP_22_6190492 T/C 22 619492 8.402 novel gene* 617712-621986 注：*表示未在大菱鲆全基因组中进行功能注释的新基因. Note: * represents novel gene which has not found any gene function annotation in the genome of Scophthalmus maximus.

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

表4

表4 　GEMMA和表型正交化方法多性状关联检验运算时间的比较 Tab. 4 Comparison in runtimes (min) between GEMMA and phenotypic principal component methods for association tests 方法 method 性状个数 number of traits 2 3 4 5 6 GEMMA 0.94 2.11 4.44 9.70 32.74 表型正交化phenotypic 0.19 0.28 0.36 0.49 0.57

表4 　GEMMA和表型正交化方法多性状关联检验运算时间的比较 Tab. 4 Comparison in runtimes (min) between GEMMA and phenotypic principal component methods for association tests

大多数动植物GWAS研究主要针对单个性状进行QTN定位分析，有关水产动物生长性状的多性状GWAS研究尚处在初步阶段。Bolormaa等^[17]使用主成分分析和一系列双变量分析对奶牛状进行了多性状GWAS, 发现在统计错误率没有增加的情况下，多性状GWAS比单性状GWAS具有更高的统计效力。到目前为止，基于表型正交化的大菱鲆多性状GWAS还未见报道，我们的研究通过对相关表型正交化，使得多变量混合模型关联检验变成独立的单变量检验，与 GEMMA 相比，在相同的实验群体下表型正交化方法通过对性状的单独关联检验不但增加了所分析性状个数，且保证了QTN的检测效率，甚至比GEMMA的检测效率更高。这为分析其他水产动物性状的多性状GWAS研究提供一个有价值的工具。

在大菱鲆全基因组上寻找距离每个QTN最近的基因，结果找到VPS13B、Zgc:77056、Slc1a6、Myt1la、Plxdc2、Zgc112356共6个已知候选基因。相关研究表明，Vps13B是顶体生物发生所必需的^[18]。Slc1a6是转运体基，因为代谢提供氨基酸，保持细胞外氨基酸水平，具有高的非耦合阴离子电导作用^[19]。Myt1la髓磷脂转录因子是一种神经特异性的含锌指的DNA结合蛋白，在神经系统的发育中发挥作用^[20]。Plxdc2是I型跨膜蛋白，与巢蛋白和丛状蛋白有一定的同源性，它在发育中的神经系统中以高度离散和动态的模式表达，在各种模式中心有显著的表达^[21]。这些候选基因可能在大菱鲆生长发育过程中发挥着重要的作用，它们的筛选定位极大地提高了对大菱鲆遗传结构的认识。