全基因组关联分析(genome-wide association study, GWAS)是在全基因组范围内，以数以万计的单核苷酸多态性(single nucleotide polymorphism, SNP)作为分子标记，探究基因变异与性状之间关联的研究。GWAS已成为当今基因定位的主流方法，广泛应用于人类疾病、动植物疾病、动植物性状改良、动植物遗传育种等方面。由于GWAS的样本来源于不同的种群，每个种群独特的遗传历史、交配习惯、繁殖扩张以及随机变异都会使个体之间等位基因频率产生差异^[1-2]，存在群体分层，会导致关联分析过程中检验统计量膨胀，从而降低对数量性状核苷酸(quantitative trait nucleotide, QTN)的检测效力。主成分分析法(principal component analysis, PCA)是当下应对该问题的主要策略^[3-4]。

间断性状，即类似于性别、是否患病等表型，可用二进制或者是不连续数字表示的性状。在进行基因定位时，间断性状和连续性状一样，也需考虑除检验标记外的群体分层、家系结构和隐藏亲缘等混杂因素^[5-6]。然而与呈现正态分布的连续性状用剩余误差作为校正后的表型不同，间断性状的表型形式较为特殊，无法像连续性状一样估计表型的剩余误差，这就需要引入广义线性模型(generalized linear model, GLM)中的logit回归，对间断性状进行GWAS。即使校正了一些固定效应协变量，logit回归仍会受到群体分层的影响导致检验统计量膨胀。与线性模型(linear model, LM)相比较，广义线性模型求解要消耗更多时间和计算机内存。当协变量比较多的时候，GLM还会产生由近似引起的严重偏差问题^[7-8]。

半滑舌鳎(Cynoglossus semilaevis)在动物分类上属于鲽形目(pleuronectiformes)、舌鳎科(cynoglossidae)、舌鳎属，属于近海大型底栖暖温性动物，分布于中国、朝鲜、日本，在中国主要分布在黄海、渤海近海区域^[9]，是目前中国重要的海水养殖鱼类品种之一^[10]。不同性别的半滑舌鳎成体体重差异较大，雌性舌鳎生长速度较快，体重较大，相比之下，雄性舌鳎的体重较轻，生长缓慢。在人工养殖过程中，生理性别是雌性的舌鳎仅占20%~35%^[11]; 此外，半滑舌鳎是性逆转生物，即在发育过程中，生理性别可能发生改变，即由遗传雌性逆转为生理雄性，该现象在昆虫^[12]、爬行动物^[13]、两栖动物^[13]和鱼类^[14]中时有发生。已有研究认为半滑舌鳎在Z染色体上可能存在导致性逆转的基因^[15]。出于经济利益的考虑，在半滑舌鳎养殖过程中控制性逆转的产生，提高雌鱼的比例尤为重要。

半滑舌鳎性逆转性状在数量遗传学中被认为是间断性状，对该类性状QTN的筛选通常采用广义线性回归分析，对于存在群体分层的群体，则可采用考虑主成分(principal component, PC)的广义线性回归模型进行基因定位。为了更好地校正存在于半滑舌鳎性逆转性状中的复杂群体分层，以及提高QTN检测效力，本研究将提出一种间断性状GWAS方法，该方法可以将简单回归模型结果中显著位点的效应值和遗传力的尺度，转化为可解释的广义线性回归模型结果，降低了广义线性回归模型的求解难度和异常解出现的概率。本研究通过使用上述策略对半滑舌鳎性逆转性状进行GWAS，并与直接考虑PC的广义线性回归模型GWAS结果作比较，检验本研究提出的方法对群体分层的校正效果以及QTN检测效力，同时对检测出的影响半滑舌鳎性逆转性状的多态性核苷酸位点进行分析阐述，为以后半滑舌鳎性逆转性状的研究提供理论依据。

本研究所用半滑舌鳎实验群体育成于河北黄骅。2013年3月，随机选取6尾亲本雌鱼和11尾表型雄鱼作为亲本，混合养殖，亲本间随机交配繁殖，产生后代。半滑舌鳎性逆转通常发生在孵化后的前90 d，为避免低温性逆转现象发生，对子代幼鱼进行持续90 d、恒温22 ℃的养殖。90 d之后，随机选取268尾半滑舌鳎进行DNA提取、亲代分析、遗传和表型性别检测，去除缺失数据，115尾雌性遗传基因被用于GWAS。采集这115只雌性半滑舌鳎的鱼鳍和性腺，并用鱼鳍提取DNA，提取鱼鳍组织DNA的具体操作步骤在Jiang等^[16]的文章中有详细的描述。采用Chen等^[17]提出的方法确定遗传性别，利用遗传雌性的生殖腺组织切片进行表型性别鉴定，得到表型数据。基因型数据由上海欧易生物技术有限公司测序分析得到，使用Illumina HiSeq2500平台对115尾遗传型雌性的半滑舌鳎的鱼鳍DNA进行单端测序，采用2b-RAD方法对测序结果进行基因分型，质量控制后，最终获得共17618个SNP标记。

考虑群体分层时，通常以亲缘关系矩阵的主成分去校正由群体分层导致的分层效应，根据广义线性模型理论，连接函数建立起了间断性状表型向量y (0或1)的期望和被检验标记效应的关系：

式中，μ是y的期望，$Xb$是固定效应项，β_l为当前检验标记的加性遗传效应，z为所对应的指示变量向量。

通常，为了简化计算，事先估计出$X\hat b$并将它作为已知的回归项，此时只需要用${y^*} - X\hat b$代替y求解标记的效应：

式中，w=μ(1–μ)是权重，${y^{\bf{*}}}{\bf{ = }}z{\beta _{\rm{l}}}{\bf{ + }}\frac{{y - \mu }}{w}$是由表型值和回归系数组成的新的因变量。

采用卡方统计量，推断SNP与间断性状关联：

这个统计量服从自由度为1的卡方分布。

考虑群体分层时，通常以亲缘关系矩阵的特征向量作为基因组简单线性回归模型的协变量去校正由群体分层导致的分层效应，简单回归模型变为：

式中，β_c是简单回归模型求出的效应值，β_c可以直接用最小二乘法求解，e~$N(0,\sigma _e^2)$是回归的误差项。采用卡方统计量，推断SNP与间断性状关联：

这个统计量服从自由度为1的卡方分布。

当广义线性模型的回归系数多，群体小或标记等位基因频率较低时，广义线性模型很难收敛，理论和实践证明，GLM与LM具有相同的统计效力^[18]，所以本研究先采用(4)中的简单回归模型对每个标记进行检验，并以Bonferroni矫正阈值为标准，筛选候选QTN。假设检测到了q个候选QTN，将这q个标记指示变量z逐个代入考虑群体分层的广义线性模型中重新求解β_l，将新求解出的β_l代替简单回归模型求解出的β_c，这样标记的效应值就从简单回归尺度转化到了广义线性回归尺度。此时，总检测到QTN的遗传力是$\frac{{\sum\limits_1^q {{\mathop{\rm Var}\nolimits} ({z_i}{\beta _{{\rm{l}}i}})} }}{{\sum\limits_1^q {{\mathop{\rm Var}\nolimits} ({z_i}{\beta _{{\rm{l}}i}})} + 1}}$，其中第i个QTN的遗传力是$\frac{{{\mathop{\rm Var}\nolimits} ({x_i}{\beta _{{\rm{l}}i}})}}{{\sum\limits_1^q {{\mathop{\rm Var}\nolimits} ({x_i}{\beta _{{\rm{l}}i}})} + 1}}$。

在GWAS中，Q-Q图是算法统计性质的直接体现；曼哈顿图常用于显著QTN的定位和筛选，分布于水平线即5% Bonferroni校正阈值线(2.838× 10^–6)之上的点即为显著QTN; 基因组控制(genomic control, GC)值可以直接显示是否存在群体分层，GC值接近1说明基因组控制效果较好。

图1 所示是半滑舌鳎性逆转性状考虑复杂群体分层的广义线性回归模型的GWAS结果。观察Q-Q图可以得知，1个PC和3个PC的Q-Q图前半段均未贴合理论线，5个PC的Q-Q图的前半段能够很好地贴合理论线，后半部分上扬；而且1个PC时的GC值为2.37, 3个PC时的GC值为1.64, 5个PC时的GC值为1.01, GC值和Q-Q图结果都显示，1个PC和3个PC时的群体结构未得到很好的校正，存在假阳性。考虑5个PC的广义线性回归模型的GC值接近1，基因组控制效果较好，观察曼哈顿图可见，这时在Z染色体有2个QTN超过阈值线。

半滑舌鳎性逆转性状使用基因组简单线性回归模型尺度转换后的GWAS结果如图2所示，此方法把简单回归结果中显著位点的效应值和遗传力的尺度转化为可解释的广义线性回归结果。可以看出，1、5、30个PC的Q-Q图前面部分都能较好地贴合理论线，尾部上扬，说明3种情况都能对群体结构较好校正。计算GC值能直接了解基因组控制效果，1个PC的GC值为1.1, 5个PC的GC值为1.1, 30个PC的GC值为1.01，可以看出考虑30个PC的简单回归尺度转换方法对于群体分层的校正效果最好。由曼哈顿图可看出，检测到6个QTN, 5个在Z染色体上，1个在W染色体上；同时还检测到13号染色体的1个接近阈值的位点，以及Z染色体上3个接近阈值的位点。

半滑舌鳎性逆转性状考虑30个PC，使用简单回归模型进行尺度转换，然后进行GWAS，检出的QTN位点的详细信息如表1所示。表中QTN遗传力说明各个位点能解释性逆转性状遗传变异的百分比，6个P值超出阈值的SNP位点即为检测到的与半滑舌鳎性逆转性状相关的QTN，其余4个接近阈值的位点，可作为候选QTN。

图1

图1 　半滑舌鳎性逆转性状广义线性回归模型GWAS曼哈顿和Q-Q图A. 1个PC, B. 3个PC, C. 5个PC Fig. 1 　Manhattan and QQ plots of generalized linear regression model GWAS on sex reversal traits in Half-smooth ongue soleA. 1 PC, B. 3 PCs, C. 5 PCs

图2

图2 　半滑舌鳎性逆转性状简单回归模型尺度转换GWAS曼哈顿和Q-Q图A. 1个PC, B. 5个PC, C. 30个PC. Fig. 2 　Manhattan and QQ plots for scale transformation of simple regression model GWAS on sex reversal traits in alf-smooth Tongue soleA. 1 PC, B. 5 PCs, C. 30 PCs.

表1

表1 　检测到的与半滑舌鳎性逆转性状相关联的QTN信息 Tab. 1 　QTN information associated with sex reversal trait detected in Cynoglossus semilaevis 染色体 chromosome 位置 location SNP编号 SNP ID 基因名称 gene name Lm效应 Lm effect Lm标准误 Lm standard error Glm效应 Glm effect Glm标准误 Glm standard error QTN遗传力/% QTN heritability P Z 9793913 M65579 Si:deky-193c22.1 −0.3032 0.0466 −1.2138 0.1864 0.189 8.24E-09* Z 8643646 M65437 DAPK1 −0.4751 0.0735 −1.9015 0.2942 0.198 1.01E-08* Z 5833328 M65164 ADGRD2 −0.2998 0.0538 −1.1999 0.2155 0.240 4.06E-07* Z 6676874 M65244 FBXL17 −0.2540 0.0477 −1.0168 0.1910 0.236 1.08E-06* Z 7256721 M65301 DMXL1 0.4827 0.0915 1.9322 0.3663 0.275 1.31E-06* W 1218194 M63786 LOC103396896 0.5462 0.1067 2.1863 0.4270 0.275 2.40E-06* Z 1633247 M64758 LOC103397473 −0.4116 0.0879 −1.6475 0.3516 0.270 1.26E-05 Z 9232660 M65511 无 0.3734 0.0811 1.4946 0.3245 0.275 1.70E-05 13 2140529 M40452 CD27 −0.9204 0.2006 −3.6840 0.8029 0.227 1.82E-05 Z 781540 M64676 MALT1 −0.2973 0.0657 −1.1902 0.2631 0.252 2.31E-05 注：*表示该SNP位点与性逆转性状显著相关(P<2.838×10–6). Note: * represent that the SNP location is significantly correlate with sex reversal trait (P<2.838×10–6).

表1 　检测到的与半滑舌鳎性逆转性状相关联的QTN信息 Tab. 1 　QTN information associated with sex reversal trait detected in Cynoglossus semilaevis

本研究结果表明，简单回归尺度转换方法和广义线性回归模型都可以很好地校正关联分析中的群体分层，无论是评价指标中的Q-Q图还是基因组控制值，都显示只要考虑作为协变量的主成分数目恰当，就可以很好地控制全基因组关联分析中的假阳性。比较两种方法检出QTN的个数可知，基于简单回归尺度转换的简单回归模型比直接考虑PC的广义线性回归模型的QTN的检出效力更高，且使用简单回归尺度转换方法检测到的QTN高于已发表的方法，不但检测到了Z染色体上的QTN，还检测到了W染色体的一个QTN。

Z染色体上位置为9793913的SNP位于si:deky-193c22.1基因上，未有研究者对该基因的功能进行探究。Z染色体上位置为8643646的SNP在DAPK1 (death-associated protein kinase 1)基因上，DAPK1与细胞的死亡有关，Raval等^[19]研究发现，DAPK1的表达缺失或降低是人类慢性B淋巴细胞白血病发病的原因；Siqueira等^[20]在不同性别牛胚胎对群落刺激因子2反应的研究中证实，DAPK1的表达有受性别影响的趋势。Z染色体上位置为5833328的SNP在ADGRD2 (adhesion G protein-coupled receptor D2)的基因上。Z染色体上位置为6676874的SNP在FBXL17 (F-box and leucine-rich repeat protein 17)的基因上，Jiang等^[16]有关半滑舌鳎性逆转性状的研究发现，FBXL17基因上位置为6676874的SNP与半滑舌鳎性逆转性状有关，本研究所应用的方法也检测到了这一位点。Z染色体上位置为7256721的SNP，在DMXL1 (Dmx like 1)基因上，人类DMXL1基因编码一个11 kb的RNA，开放阅读框3027个氨基酸，主要具有调节功能^[21]。W染色体上位置为1218194的SNP，是唯一检测到的W染色体上的QTN，在LOC103396896基因上，还没有研究者对该基因功能进行深入的研究，但是这个基因位于W染色体上，有可能参与到性逆转性状的调控中。

除了检出的5个在Z染色体上的QTN和1个在W染色体上的QTN之外，还检测到Z染色体上的另外3个QTN，但未超过阈值，13号染色体上也检测出了一个接近阈值的QTN，这几个位点如下：Z染色体上位置为1633247的SNP，在LOC103397473基因上；Z染色体上位置为9232660的SNP，未在基因编码区内；Z染色体上位置为781540的SNP，在malt1 (MALT paracaspase 1)基因上，MALT1是在抗原受体介导的淋巴细胞活化中起关键作用的信号蛋白^[22]; 13号染色体上检出的位点位置是2140529，在CD27基因上，CD27是T细胞的共刺激分子，可支持naïve T细胞的抗原特异性扩增，若删除或改变编码CD27的基因，免疫应答会发生延迟，特异性T细胞数量减少^[23]。这些位点虽未超过阈值，但可作为候选QTN，为后续有关半滑舌鳎性逆转性状以及其他性状的全基因组关联分析研究奠定基础。

运用基于尺度转换的简单回归模型所检测出的与半滑舌鳎性逆转性状相关的位点，除了Z染色体上位置为6676874的SNP与半滑舌鳎性逆转性状有关之外^[16]，其余位点到目前为止没有报告或有明显的证据表明与性别决定的各种途径直接相关，但它们仍值得被关注，因为这些位点基本都在Z染色体或者W染色体上。众所周知，Z染色体和W染色体是半滑舌鳎的性染色体，雌性为ZW，雄性为ZZ，这两条染色体上的基因以及他们的连锁基因对于性别决定都可能具有潜在的影响，任一位置的突变都有可能对决定性别的基因产生影响。另外4个接近阈值的位点可作为候选QTN，为研究者们研究半滑舌鳎性逆转性状提供了候选研究位点。