山羊产羔数全基因组关联分析

兰 蓉，朱 兰，姚新荣，王 鹏，邵庆勇，洪琼花*

(1.云南省畜牧兽医科学院，昆明 650224； 2.云南省种羊繁育推广中心，寻甸 650000)

山羊产羔数全基因组关联分析

兰 蓉1，朱 兰1，姚新荣2，王 鹏2，邵庆勇1，洪琼花1*

(1.云南省畜牧兽医科学院，昆明 650224； 2.云南省种羊繁育推广中心，寻甸 650000)

本研究旨在通过对山羊产羔性状的全基因组关联分析(GWAS)，寻找和定位与山羊产羔性状密切关联的新基因。以云南黑山羊6个公羊家系的302只山羊为试验材料，用Illumina 公司Iselect Goat60k芯片技术进行SNP分型，分型结果利用plink V1.07的线性回归模型对山羊产羔数性状进行全基因组关联分析。研究结果表明：在2号染色体上有2个SNPs位点与山羊产羔数达到5%基因组水平显著相关(P<1.48E-6)，分别位于SLC4A10基因的下游和TBR1基因的上游；5个SNPs位点与山羊产羔数达潜在关联(P<2.97E-5)，分别位于1号染色体SENP7基因上游，21号染色体Hypothetical Protein基因的上游，以及28号染色体WDFY4基因和TMEM26基因的上游、BICC1基因的下游。这些基因可作为山羊产羔数性状的相关候选基因，也可为山羊产羔性状的分子机制研究和今后标记辅助选择的开展提供理论基础及新的研究线索。

云南黑山羊；产羔数；全基因组关联分析

山羊产羔数的不断提高一直是山羊生产的重要目标，它关系到生产效率的提高和生产成本的降低。增加山羊的胎产羔数，提高生产效益，促进养羊业向持续环保发展道路的转变，作为一个畜牧大国和人口大国的中国显得尤为重要。但是山羊的产羔性状遗传力较低，是一个受微效多基因控制的复杂数量性状，同时也受一个或多个主基因控制，常规育种手段的运用不能在较短时期内取得明显进展，这在很大程度上阻碍了山羊业的快速发展。准确筛选、鉴定产羔性状主效基因，进而进行分子育种，成为快速改良和提高山羊产羔数的前提条件。

目前已有很多学者对山羊产羔数相关基因做了研究，包括FSHR[1]、gdf-9[2]和NPY-Y1R等基因[3]，而利用目前较为先进的全基因组关联分析(GWAS)技术，开展山羊产羔性状的基因挖掘和定位的研究还未见报导。本研究拟采用GWAS技术对山羊产羔数的分子基础进行研究，以期挖掘和定位与山羊产羔数密切相关的新基因。

1 材料与方法

1.1 试验材料

本研究试验羊只来自云南省种羊繁育推广中心云南黑山羊核心群，由6个公羊家系的302只母羊组成，采集颈静脉EDTA抗凝血2 mL，-20 ℃保存备用，记录试验羊只第二胎产羔数。

1.2 方法

1.2.1 基因组DNA提取及检测 按blood genome DNA Extraction kit试剂盒操作流程提取血样DNA，1×TE缓冲液溶解，用NANO核酸蛋白分析仪检测DAN的浓度和质量。使DNA 的吸光度OD260 nm/OD280 nm为1.7～2.0，浓度大于50 ng·μL-1、体积大于20 μL的样品作为基因分型样品，保存于-20 ℃备用。

1.2.2 SNP基因分型及质量控制 将质量达标的DNA样品送至北京斯克尔基因生物技术有限公司，用美国Illumina 公司的山羊Iselect Goat60k芯片进行SNP分型。

1.2.3 芯片分型流程 包括DNA定量、扩增、片段化、片段化DNA沉淀及回收、芯片杂交、芯片清洗、单碱基延伸及染色、芯片扫描、数据分析、得出SNP分型结果。

1.2.4 质量控制 使用Plink(1.07) 进行质量控制[4]，将MAF小于0.05，检出频率小于0.9，哈迪-温伯格平衡卡方检验将P<1 E-06的SNP位点排除，同时还排除检出率<0.9的样本。

1.3 数据统计分析

1.3.1 表型性状统计分析 用SAS(8.01)中的CORR过程对试验羊只第二胎产羔数进行描述性统计及羊只出生年度与产羔数的相关性分析。

1.3.2 群体结构分析 为了消除连锁不平衡对个体间遗传相关估计的偏差，选用常染色体上独立的SNPs，应用Plink(1.07)全基因组IBS(Identity-by-state)多维度分析[4]，进行MDS分析。SNPs筛选时以25个位点为一个窗口，以5 个位点为步长，r2设为0.2，基于得到的独立SNP点计算个体之间的IBS 距离，应用R(V2.13.1)作图[5]，进而进行群体结构分析。

1.3.3 全基因组关联分析 运用PLINK(V1.07)线性回归模型进行产羔性状的全基因组关联分析。为减少重复检验的假阳性，使用连锁不平衡修正的Bonferroni 修正，以25个位点为一个窗口，5个位点为步长，r2为0.4，得到独立的LD 块和SNP，并据此进行全基因组关联分析值校正。

1.3.4 生物信息分析 应用软件BLAST(2.28+)将潜在关联的SNP位点上下60 bp序列比对到参考基因组中[6-7]，得到SNP所在的染色体位置坐标， 结合山羊参考基因组Annotation-GFF文件，查找SNP位置前后500 K的基因，获得Gene-ID，根据Gene-ID获取基因的编码区(CDS)序列和pep序列，再用基因的pep序列进行Nr库注释，获取基因的注释信息。

2 结 果

2.1 产羔数统计结果

云南黑山羊二胎产羔数平均为2.04只，变异系数为0.75，通过相关统计分析，产羔数与母羊的出生年龄相关系数为0.15，相关不显著(P>0.05)，说明母羊年龄对其产羔数无影响，因而在后面GWAS分析中采用了线性回归模型。

2.2 基因分型数据质控结果

经质控后MAF小于0.05的位点有6 162个，检出频率<0.9的位点有337 个，哈迪-温伯格平衡卡方P<0.000 001的位点有327 个，检出率<0.9的样本、排除染色体位置信息不确定点474个，最终获得295 个样本在45 608 个位点上的分型结果(表1)。

2.3 群体结构分析结果

通过MDS分析，在常染色体可筛选出15 309个独立的SNPs位点，以计算的MDS组分1 和组分2 为坐标，应用R(V2.13.1)绘制群体结构散点图(图1)，6个公羊家系分别用一种颜色表示，由图1可见同一个家系个体几乎聚集在一起，并未均匀分布，表明存在群体分层，同一家系个体间的遗传相关程度较高。

表1 每条染色体上SNP数及对应的Bonferroni校正的5%染色体水平显著阈值

Table 1 Bonferroni-corrected 5% chromosome-wise significance threshold after quality control

染色体Chromosome标记数No.ofSNPs独立的标记数No.ofindependentSNPs5%染色体水平显著阈值5%chromosome-widesignificantthresholdch1293222332.24E-05ch2255019462.57E-05ch3206215023.33E-05ch4216715533.22E-05ch5197013523.70E-05ch6201014333.49E-05ch7194814413.47E-05ch8207015293.27E-05ch9169213313.76E-05ch10189314623.42E-05ch11180912663.95E-05ch12157311684.28E-05ch13147811344.41E-05ch14166512593.97E-05ch1514539885.06E-05ch16143810644.70E-05ch1712709915.05E-05ch1810968535.86E-05ch1911228166.13E-05ch2012959985.01E-05ch2113199565.23E-05ch2210437996.26E-05ch239306487.72E-05ch2411918555.85E-05ch257695658.85E-05ch269487226.93E-05ch278346527.67E-05ch288016467.74E-05ch298676887.27E-05ChX14138266.05E-05Total4560833674

2.4 全基因组关联分析结果

使用连锁不平衡修正的Bonferroni 进行修正，得到独立的LD 块和33 674个独立的SNPs用于GWAS分析，据此得到Bonferroni 修正的达5%基因组显著水平的P值为1.48E-6(0.05/33674)，基因组潜在关联阈值为2.97E-5(1/33674)。同时根据每条染色体上的独立SNPs 数目，得到每条染色体水平的显著性阈值(表1)，通过对比，发现达到潜在基因组水平显著差异的SNPs位点共7个(图2，表2)。图2和表2显示：2个SNPs位点与山羊产仔数达到5%基因组水平显著相关(P<1.48E-6)，分别位于SLC4A10基因的下游和TBR1基因的上游；5个位点达潜在关联(P<2.97E-5)，分别位于1号染色体SENP7基因上游，21号染色体Hypothetical Protein基因的上游，以及28号染色体WDFY4和TMEM26基因的上游、BICC1基因的下游。

图1 群体结构多维度分析图Fig.1 Population structure identification with multidimensional scaling analysis

横坐标为染色体号，30 代表X 染色体。纵坐标为关联分析P 值的-log10 结果。黑色实线为达基因组水平5%显著阈值线，黑色虚线为达基因组潜在关联阈值线The scale on the X-axis represents ID of chromosomes.X chromosome is represented by 30.The scale on the Y-axis is the -log10(P-value) score of association analysis.The black solid line is drawn at -log10(1.48E-6) to show those significant at the genome-wide 5% threshold，and the black dashed line indicates genome-wide significance of suggestive association图2 产羔数全基因组关联分析曼哈顿图Fig.2 Manhattan plot of genome-wide association analysis for lambing number

3 讨 论

全基因组关联分析(Genome-wide association study，GWAS)是一种利用遍布于生物体基因组中数以百万的分子标记(主要是单核苷酸突变，即SNP)，并借助统计学工具对影响某些复杂性状如疾病易感性或重要经济性状的遗传变异进行鉴定和分析的一种新策略，其基本思想是直接检测基因本身或基因附近微小区域的SNP 标记与性状表型信息的关联来实现基因的精细定位。目前，GWAS技术已在奶牛、绵羊、猪、鸡等畜种中得以应用[8-12]，挖掘出许多与重要经济性状相关的新基因，包括繁殖、遗传抗性、产乳性状、生长性状等重要经济性状相关的新基因，而有关山羊产羔性状的GWAS研究则未见报导。

表2 产羔数Bonferroni 校正的5%染色体水平显著的SNPs位点信息

Table 2 Bonferroni-corrected 5% chromosome-wide significant SNPs for lambing number

标记SNPs染色体Chromosome物理位置BP基因型Allele最近基因NearestgeneP值P-valuesnp20467-scaffold202-4458644234235387T/CSLC4A10(D)8.57E-07snp20468-scaffold202-4490513234267256T/GTBR1(U)8.57E-07snp11435-scaffold1416-3795782840521920T/CWDFY4(U)4.07E-06snp45731-scaffold627-4874627144576529A/GSENP7(U)4.49E-06snp54840-scaffold838-31946002814948614T/CTMEM26(U)1.50E-05snp25958-scaffold2688-2609162165919086T/CHypotheticalProtein(U)1.86E-05snp54781-scaffold838-7287692812482783A/CBICC1(D)1.92E-05

D.SNP位于基因的下游； U.SNP位于基因的上游

D.SNP located in the downstream of the gene；U.SNP located in the upstream of the gene

本研究MDS分析发现试验样本存在群体分层。已有的人类疾病的全基因组关联分析发现，由于祖先的遗传差异引起疾病和对照间等位基因频率的差异造成的群体分层现象，导致标记和性状间的伪关联[13]，利用全基因组标记的遗传信息对样本遗传相关进行估计，比系谱信息可更真实的反映个体间的相关度，因此本研究GWAS分析时使用了MDS分析产生的组分1作为协变量，以校正群体结构对关联分析结果的影响[14]，对应的全基因组膨胀系数由1.14～2.15减小为1～1.05，可以认为修正后得到的结果没有受到群体分层等混杂因素的影响。

本研究在生物信息学分析中均是用差异位点与附近基因的编码区进行比对才确认SNP位点是位于该基因的上游或下游，因此本研究发现的与产羔数相关的基因均是位于功能基因上。

本研究还发现，在28号染色体上有3个SNPs位点与产羔数存在潜在关联，分布在12～40 Mb，这个区域内有WDFY4、TMEM26及BICC3个候选基因，但对这3个基因的相关功能研究极少，在山羊中更是未见报导，所以本研究发现，它们与山羊的产羔数存在潜在的关联还需要做大样本、更为深入的研究，以确保得到的结论真实可靠。

本研究结果可作为山羊产羔性状的相关候选区域和候选基因，也可为山羊产羔性状的分子机制研究和今后的标记辅助选择的开展提供新的研究线索和理论基础。但本研究由于经费及时间所限，只采用了302只云南黑山羊作为试验样本，因此所获候选区域及基因仍然需要进行更大样本的重复验证，以确保今后在分子育种中应用的准确性和效率。

4 结 论

本研究鉴定出位于2号染色体上的SLC4A10和TBR1基因为影响山羊产羔性状的重要候选基因；28号染色体12～40 Mb的区域内分布的WDFY4、TMEM26和BICC1基因对山羊产羔数有潜在的影响。

[1] 朱 吉，杨仕柳，欧阳叙向，等.山羊FSHR基因5′端多态性与产仔数性状的关系 [J].草食家畜，2007， 134(1)：12-16. ZHU J，YANG S L，OUYANG X X，et al.Relationship between polymorphism of follicle stimulating hormone receptor gene and litter size traits in goats.[J].Grass-FeedingLivestock，2007， 134(1)：12-16.(in Chinese)

[2] 任丽娜，刘思国，胡大伟.崇明白山羊gdf-9基因多态性与产仔率的关联分析 [J].畜牧与兽医，2012， 44(10)：25-29. REN L N，LIU S G，HU D W，et al.Polymorphism of growth differentiation factor 9 gene and its relationship with litter size of Chong-ming white goats.[J].AnimalHusbandry&VeterinaryMedicine，2012，44(10)：25-29.(in Chinese)

[3] 储明星，冯 涛，赵有璋，等.神经肽Y-Y1受体基因(NPY-Y1R)多态性及其与山羊繁殖性能关系 [J].农业生物技术学报，2009，17(6)：960-966. CHU M X，FENG T，ZHAO Y Z，et al.Polymorphism of neruropeptide Y-Y1 receptor gene(NPY-Y1R) and its relationship with reproductive performance in goats [J].JournalofAgricultureBiotechnology，2009，17(6)：960-966.(in Chinese)

[4] PURCELL S，NEALE B，TODD-BROWN K，et al.PLINK：a tool set for whole-genome association and population-based linkage analyses [J].AmJHumGenet，2007，81(3):559-575.

[5] DEVELOPMENT R，CORE T R.A language and environment for statistical computing R Foundation for Statistical Computing [M].Vienna，Austria，2011.

[6] http://goat.kiz.ac.cn/GGD/download9.htm.

[7] DONG Y，XIE M，JIANG Y，et al.Sequencing and automated whole-genome optical mapping of the genome of a domestic goat (Caprahircus) [J].NatBiotechnol，2013，31(2)：135-141.

[8] LI C，SUN D，ZHANG S， et al.Genome wide association study identifies 20 novel promising genes associated with milk fatty acid traits in chinese Holstein [J].PLoSONE，2014，9(5)：e96186.

[9] KIJAS J W，SERRANO M，MCCULLOCH R，et al.Genome wide association for a dominant pigmentation gene in sheep [J].JAnimBreedGenet，2013，130(6)：468-475.

[11] JUNG E J，PARK H B，LEE J B，et al.Genome-wide association study identifies quantitative trait loci affecting hematological traits in an F2 intercross between Landrace and Korean native pigs [J].AnimGenet，2014，45(4)：534-541.

[12] SUN Y F，ZHAO G P，LIU R R，et al.The identification of 14 new genes for meat quality traits in chicken using a genome-wide association study [J].BMCGenomics，2013，14：458-468.

[13] PRICE A L，PATTERSON N J，PLENGE R M，et al.Principal components analysis corrects for stratification in genome-wide association studies [J].NatGenet，2006，38(8)：904-909.

[14] 顾晓荣.影响鸡体重性状QTL的全基因组连锁和关联研究[D].北京：中国农业大学，2011. GU X R.Whole genome linkage and association studies on QTLs effecting body weight in chicken[D].Beijing：China Agricultural University，2011.(in Chinese)

[15] DAMKIER H H，AALKJAER C，PRAETORIUS J.Na+-dependent HCO-3import by the slc4a10 gene product involves Cl-export [J].JBiolChem，2010，285(35)：26998-27007.

[16] DEMARCO L A，ESPINOSA F，EDWARDS J，et al.Involvement of a Na+/HCO-3cotransporter in mouse sperm capacitation [J].JBiolChem，2003，278(9)：7001-7009.

[17] TIENTHAI P，JOHANNISSON A，RODRIGUEZ-MARTINEZ H.Sperm capacitation in the porcine oviduct [J].AnimReprodSci，2004，80(1-2)：131-146.[18] BEDOGNI F，HODGE R D，ELSEN G E，et al.Tbr1 regulates regional and laminar identity of postmitotic neurons in developing neocortex [J].ProcNatlAcadSciUSA，2010，107(29)：13129-13134.

(编辑 郭云雁)

Genome-wide Association Study of Lambing Number in Goat

LAN Rong1，ZHU Lan1，YAO Xin-rong2，WANG Peng2，SHAO Qing-yong1，HONG Qiong-hua1*

(1.YunnanInstituteofAnimalScienceandVeterinary，Kunming650224，China；2.SheepBreedingPromotionCenterinYunnanProvince，Xundian650000，China)

The aim of this study was to find and map novel genes related to lambing number in goat by genome-wide association study(GWAS).A total of 302 Yunnan Black goats from 6 families were genotyped by Illumina Goat60K iSelect array.The GWAS on goat lambing number was performed.The Results showed that 2 SNPs located at the downstream ofSLC4A10 and upstream ofTBR1 on chromosome 2 were significantly associated with lambing number (P<1.48E-6) at 5% genome-wide level.Five SNPs were suggestive associated with lambing number (P<2.97E-5).These 5 SNPs were located on chromosome 1，21，28，respectively，includingSENP7，Hypothetical Protein，WDFY4，TMEM26 andBICC1 genes.These genes and SNPs offer essential information for understanding the molecular mechanisms of lambing number and provide foundation for the application of marker-assisted selection in breeding program of goat.

Yunnan Black goat；lambing number；genome-wide association study

10.11843/j.issn.0366-6964.2015.04.006

2014-07-02

国家现代农业产业技术体系建设专项(CARS-39)；云南省科技计划项目(2014BB014)

兰 蓉(1969-)，女，云南人，副研究员，硕士，主要从事家畜分子遗传育种研究，E-mail:rtlankitty@163.com

*通信作者：洪琼花，研究员，E-mail:yxh7168@126.com

S827；S813.3

A

0366-6964(2015)04-0549-06