滩羊毛色的全基因组关联分析

王晓薇， 马 青

(宁夏农林科学院 动物科学研究所，宁夏 银川 750002)

毛皮颜色不仅是毛皮动物的重要经济性状，还可用于辅助鉴定品种纯度、确定杂交组合等育种工作。滩羊是我国农业农村部重点保护的地方优良品种之一，其所产的二毛裘皮轻薄保暖，在国际上也享有盛誉。滩羊体躯毛色纯白，多数个体头部有褐、黑、黄等斑块。但也有少量滩羊呈黑头或体躯含有黑色斑块，约占总数的0.2%左右[1]。目前生产滩羊二裘皮主要使用白色滩羊。

研究表明，与哺乳动物毛色相关的基因和位点多达378个[2]，这些基因功能多是调控黑色素的生成。黑色素是由黑色素细胞内的黑色素小体产生的一种不溶于水的蛋白质衍生物，广泛存在于高等动物的皮肤中。

在过去的十多年里，随着高通量测序技术和芯片技术的不断发展，全基因组关联分析(genome-wide association study，GWAS)的想法最早在1996年由人类遗传学家为展开人类复杂疾病的遗传研究而提出的[3]，已逐渐成为寻找兴趣性状因果突变位点的基本策略。GWAS是一种利用群体水平连锁不平衡，通过识别大量单核苷酸多态性(SNP)标记，在全基因组范围内定位与兴趣性状变异相关联的分子标记，进而分析这些标记对性状的遗传效应。与传统的标记相比，SNP标记具有更高的基因组覆盖率，因此无需基于候选基因，可直接在全基因组范围内寻找因果关联位点。鉴于GWAS在人类疾病研究上的成功经验，无数动物育种学家也将该技术引入动物经济性状的候选基因的挖掘工作中。在绵羊经济性状中，已有大量关于产毛性状、生长发育性状、肉质性状、繁殖性状、产奶性状的GWAS研究报道，挖掘了大量的候选基因[4]。其中，关于绵羊毛色的研究并不多见。Li等[5]使用50K绵羊芯片对99只芬兰羊的毛色表型进行GWAS分析。他们使用了两种不同的策略：首先把表型分为5类(白毛、灰毛、棕毛、黑毛和黑毛带白点)进行GWAS分析，鉴定到TYRP1、ASIP、MITF共3个候选基因。他们又将表型分为两类(白色、非白毛)仅ASIP基因被鉴定到。但到目前为止，并未见关于滩羊毛色候选基因研究的报道。因此，鉴定滩羊毛色的候选基因不仅有助于解析滩羊毛色性状的遗传机理，还可以辅助滩羊育种。

1 材料与方法

1.1 试验动物及性状

本试验所用的96只滩羊均来自宁夏盐池三个不同滩羊场的核心育种群。根据滩羊毛色不同将所研究的96只滩羊分为3类：全身白毛(简称：全白)、身体全白，头部黑色或带部分黑色斑点(简称：白毛黑头)或身体全白，头部带褐色或黄色斑点(简称：白毛褐或黄头)。

1.2 基因分型及质量控制

从羊颈静脉采取10 mL血样装入抗凝管中，随后采用常规苯酚/氯仿方法提取样本DNA。用美国Affymetrix绵羊600K芯片对96只滩羊进行基因分型。通过iScan扫描芯片结果，并利用Genome Studio软件对原始数据进行分析，从而得到每个样品的SNP分型数据。

得到的SNP数据需进行常规质量控制后才可以进行GWAS分析。我们使用Plink 1.90进行标记质量控制，标准如下：去除检出率低于90%的SNP，去除哈代温伯格平衡检验P值小于0.000 01的SNP，去除最小等位基因频率小于0.01的SNP，去除性染色体上的SNP，最后去除检出率低于90%的个体。

1.3 群体结构分析

对存在亚群的群体进行关联分析会造成大量假阳性结果的出现，因此我们需要了解所研究滩羊群体的亲缘关系及群体结构分布情况。我们使用Plink 1.90[6]对SNP分型数据计算IBS遗传距离矩阵和主成分分析。对于IBS距离矩阵使用R语言包qplots绘制热图进行可视化。在GWAS分析中利用主成分分析校正群体结构的方法最早由Price等[7]提出。由于该方法简单有效，迅速被广大研究人员认可。这里我们提取解释方差最大的前两个主成分绘制散点图查看群体分层情况。

1.4 全基因组关联分析模型

本研究表型为分级性状，我们采用Logistic回归对全部标记进行单点扫描。分析模型如下：

link=Xβ+Zkγk+e。

其中P为显著性阈值，m为标记个数。

1.5 基因注释

我们将检测到的显著标记定位到绵羊参考基因组(Ovi_arise_V3.1)上，以寻找候选基因。考虑到基因上下游的调控作用，以及潜在的连锁不平衡效应，我们使用显著标记前后100 kb的物理区间来筛选候选基因。若显著标记前后100 kb的物理区间与多个基因重叠，则挑选距离该标记最近的基因，若没有与任何基因重叠，则认为没有注释到基因。最后，我们将注释到的基因与已发表的文献进行对比，看是否已经报道过该基因有相关功能，以进一步佐证我们的研究结果。

2 结果与分析

2.1 表型分布及基因分型质控

96只滩羊毛色分布为全身白毛33只、白毛黑头28只、白毛褐或黄头33只，三种毛色滩羊分布基本均匀。基因分型结果共包含604 727个SNP，去除检出率小于90%的SNP 151 158个，MAF小于0.01的SNP 27 450个，哈代温伯格平衡小于0.000 01的SNP 12 650个，性染色体的SNP 4 815个，总计剩余408 654个SNP用于后续分析。

2.2 群体结构检验

图1为96只滩羊的遗传距离热图。遗传关系热图利用颜色深浅不同展示了每一只羊与其余羊的遗传距离，颜色越浅代表遗传距离越大，两者的亲缘关系也就越小。从图1中，我们可以看出所研究滩羊样本大部分个体之间存在一定的遗传关联，少部分个体亲缘关系很近。这表明试验样本在不同家系间选取比较均匀，并没有形成明显的几个亚群。这有助于提高后续分析的准确性。图2展示了不同毛色滩羊的群体结构分布，从总体上来看，大部分滩羊都是聚在一起的，少部分滩羊有偏离群体的趋势。因为这里没有其他品种的羊作参照，无法看出这些偏离群体的羊是否真的具有非常明显的群体分层。我们还将不同毛色的羊用不同颜色的点表示出来，以便更好地观测不同毛色羊群的分布情况。从不同毛色羊的分布情况来看，不同毛色的滩羊是混杂在一起的，因此滩羊并没有因为毛色不同构成不同的亚群。

图1 IBS遗传距离热图Fig.1 The heatmap plot of IBS genetic distances

2.3 全基因组关联分析及候选基因注释

为了更好地将SNP的显著性情况可视化，我们利用显著性检验P值的负对数作y轴绘制曼哈顿图(图3)。本研究采用Bonferroni法校正多重比较，因此SNP的显著性阈值为1.22×10-7。在所有的候选标记中，共找到5个P值小于1.22×10-7的SNP。

这5个SNP都位于14号染色体上14.2 Mb处，因此该区域可能存在影响滩羊毛色的基因。表1展示了这5个SNP注释基因的结果，其中SNP AX-173779077位于MC1R基因内，SNP AX-173572102位于TCF25基因下游550 bp处，其余3个SNP都位于TCF25基因内。

3 讨论

随着基因芯片技术的发展，给动物育种工作者提供了更好的方式去探索基因与家畜经济性状之间的关系。全基因组关联分析更是成为了解分析家畜经济性状遗传机理的强力工具。尽管已有绵羊毛色性状全基因组关联分析的报道[5， 8]，并鉴定到了一些候选基因(例如MC1R等)，我们仍然关注这些候选基因是否对滩羊群体毛色有同样的作用。本研究通过对滩羊毛色性状进行全基因组关联分析，共鉴定到5个显著SNP，这5个SNP均位于14号染色体上14.21～14.24 Mb的区间内。我们根据这5个SNP成功挖掘到了两个毛色候选基因MC1R和TCF25。其中MC1R基因是哺乳动物毛色性状的明星基因，已在猪[8]、马[9]、羊[10]、犬[11]、鼠[12]等多个物种中均发现与毛色相关。对于TCF25基因与毛色相关的报道就很少见，TCF25基因可能仅对绵羊或是滩羊的毛色调控起到作用。考虑到TCF25基因与MC1R基因靠得很近(仅3.57 kb)，也可能是由于SNP位点之间的连锁不平衡，才鉴定到TCF25基因。

White，全白毛滩羊；White black，白毛黑头滩羊；White yellow，白毛黄头滩羊。White， Pure white Tan sheep; White black， Tan sheep with white hair and black head; White yellow， Tan sheep with white hair and yellow head.图2 PCA群体结构图Fig.2 The plot of population structure by PCA

图3 GWAS结果曼哈顿图Fig.3 Manhattan plot of GWAS results

表1 显著SNP注释基因结果Table 1 The annotation of significant SNPs

MC1R(黑色素皮质激素受体1)基因属于G蛋白耦合受体-黑色素皮质激素受体家族成员之一。MC1R基因与α-促黑色素细胞激素(α-MSH)和肾上腺皮质激素结合可以活化酪氨酸，活化酪氨酸不足时，可以促进细胞合成真黑色素，活化酪氨酸过量时，细胞会合成褐黑素[13-14]。Deng等[15]报道了云南黑骨绵羊的毛色与MC1R基因内的碱基突变相关。Fontanesi等[16]利用161只马萨斯绵羊研究发现，MC1R基因上的3个错义突变会导致马萨斯绵羊毛色变异。付冬丽[10]对10个中国地方品种羊的毛色与MC1R基因的关联性进行了研究，在MC1R基因内发现了5个可能与毛色性状相关的单碱基突变。何军敏等[17]对220只哈萨克羊的毛色(黑色、棕色、青灰色和白色)与MC1R基因内的突变位点进行关联分析，发现MC1R基因突变位点在白色和棕色背毛群体中，只存在BB基因型。因此，可以认为MC1R基因的突变会影响哈萨克羊的毛色。景炅婕等[18]对MC1R基因在不同毛色(纯黑、纯白、黄褐和黑斑白)的太行山羊皮肤中的表达量进行了比较研究，发现MC1R基因在纯黑个体中表达量最高，在黄褐色个体中表达量次之。他们的研究亦证明了MC1R基因与太行山羊的毛色性状有关。综上，大量的研究都表明了MC1R基因与绵羊和山羊的毛色性状有关。我们的试验再次证明了这一事实，也证明了MC1R基因与宁夏盐池滩羊的毛色性状相关。

转录因子25(Transcription Factor 25，TCF25)基因是在胚胎发育中很重要的转录因子[19]，属于基本-螺旋-环-螺旋(bHLH)结构。Ishida等[20]对貂的MC1R和TCF25基因测序并进行单倍型分析，指出两种存在连锁不平衡。Li等[21]利用随机森林法对42只Manchega绵羊和Rasa Aragonesa绵羊(18只黑色，24只白色)进行GWAS分析，发现有86%黑色羊在位点s26449发生碱基突变，该位点距离TCF25基因419 bp，这揭示了TCF25基因与毛色存在潜在关联。Kim等[22]对3种韩牛群体(棕色韩牛、斑纹韩牛、济州黑牛)进行基因组选择信号检测，利用群体分化指数(Fst法)检测斑纹牛基因组中的受选择区域，定位到11个候选基因，其中包含MC1R和TCF25。在人类发色研究中，Hysi等[23]研究者对290 891个欧洲人的五类头发颜色(金色、红色、浅棕、深棕、黑色)进行全基因组关联分析，该研究定位到一个显著位点靠近TCF25基因。除此之外，Morgan等[24]用同样来自欧洲的343 234个个体进行全基因组关联分析，该研究发现近200个影响发色的候选基因，同时揭露人类发色之所以会有黑色到深棕色，再到浅棕色直至金色的梯度变化，正是由这200个基因遗传差异的叠加而造成的，其中就包括TCF25和MC1R基因，并且发现MC1R基因中突变位点导致红发所能解释的遗传方差仅约为70%，MC1R内突变位点与其他显著位点产生的共同效应可解释90%的遗传方差。