王 珊，李晓林，牛志刚，史洪才

(1.石河子大学 动物科技学院，新疆 石河子 832003；2.农业部草食家畜遗传育种与繁殖重点开放实验室/新疆维吾尔自治区畜牧科学院 生物技术研究中心，新疆 乌鲁木齐 830000)

多浪羊多胎基因全基因组关联分析研究

王 珊1,2，李晓林1，牛志刚2，史洪才2

(1.石河子大学 动物科技学院，新疆 石河子 832003；2.农业部草食家畜遗传育种与繁殖重点开放实验室/新疆维吾尔自治区畜牧科学院 生物技术研究中心，新疆 乌鲁木齐 830000)

利用全基因组重测序技术，对多浪羊产羔性状进行了全基因组关联分析(GWAS)，挖掘了绵羊繁殖性状相关候选基因及分子标记，结果表明：Fst检验结果共筛选出135个显著窗口，可以发现大部分的Fst值都集中在0.35～0.50之间。对筛选到的显著窗口进行基因注释，共得到51个基因。再进行富集分析后发现Fst的最大值为0.8710，且位于NCOA1基因内部。在雌配子发生条目中，多胎组和单胎组共同基因为INHBA。同样，利用ENSEMBLE数据库进行基因注释，共有28个基因与显著窗口有交集，如ACP1、ING5、ARNT等。这些候选基因和分子标记的揭示对绵羊繁殖性状功能基因的挖掘具有参考价值和理论意义。

全基因组关联分析；多浪羊；多胎

高繁殖力是绵羊实现高产的基础，也是影响绵羊养殖经济效益最重要的因素之一。多胎羊在肥羔和羔皮羊的生产中显示出了很大的优越性，在生产效率和经济效益方面多胎羊是单胎羊的2～3倍。因此，研究绵羊的多胎性状，揭示多胎性状的遗传基础，在促进绵羊育种，培育多胎羊品种并利用多胎羊来提高绵羊繁殖力等方面有重要意义。目前，影响绵羊高繁殖力的主要基因有生长分化因子9(GDF9)基因、骨形态发生蛋白15(BMP15)基因和骨形态发生蛋白受体IB(BMPR-IB)基因[1]。褚明星等[2]采用PCR-RFLP技术，检测了生长分化因子9(GDF9)基因和骨形态发生蛋白15(BMP15)基因， 发现BMP15的B2突变(C→T)对小尾寒羊的高繁殖力有显著影响。Demars等[3]利用全基因组关联分析(GWAS)，在BMP15基因上找到了2个影响Gribette羊和Olkuska羊的突变(BMP15T317I和BMP15N337H)，同时发现以上2个突变均不同程度地影响了细胞通路。史洪才等[4]在多浪羊群体中检测BMPR-IB基因时发现了BB、B+和++3种基因型，同时发现B+型多浪羊产羔数均值比++型多0.51只(P<0.01)，从而发现了影响多浪羊的主效基因。然而野生型个体的产羔数仍然达157%，由此判定还有其他控制高繁殖力的基因存在。与传统基因分析方法相比，全基因组关联分析(GWAS)具有周期短、检测效率高等特点。

全基因组关联分析(GWAS)是利用高通量基因分型技术，分析数以万计的覆盖整个基因组的SNP与观测性状的相关性，从而全面揭示出复杂性状的遗传机制和基础，由于操作简便、实施简单，随着高通量测序成本的降低，GWAS在人类疾病以及畜禽经济性状的研究上都表现出巨大的优势。Rubin等[5]利用全基因组关联的方法分析了鸡体质量性状，发现9个SNP位点在5%全基因组显著水平，位点附近包括2个基因。Zhu等[6]利用全基因组分析方法，检测了3个绵羊品种X染色体上的选择信号，利用绵羊50K芯片和整体单体型评分，在3个绵羊品种上找到49、34、55个候选区，长度分别为27.49、16.47、25.42 Mb，这些基因与人类的具有很高的同源性，同时部分与繁殖性能相关。迄今为止，关于新疆绵羊品种的全基因组关联分析还鲜有报道。本研究利用第二代测序技术获得绵羊基因组序列，并从中检测出SNP用以进行GWAS分析，挖掘影响多浪羊繁殖性状的重要功能基因及分子标记。

1 材料与方法

1.1 试验动物

试验动物选自新疆喀什市疏附县其木布拉克合作社羊场的300只多浪羊，这个群体中选择繁殖力性状表型差异大的个体，即产羔数差异较大的母羊个体30只用于全基因组重测序，其中，连续2年产2胎及以上的母羊15只，记为多羔组；连续2年产单胎的母羊15只，记为单羔组。

1.2 基因组DNA的提取及质控

按常规酚-氯仿抽提法提取多浪羊血液组中DNA，利用分光光度计和凝胶电泳检测DNA的纯度和浓度。质量控制包括DNA的完整性(>10 kb)，DNA纯度OD260 nm/OD280 nm值应在1.8～2.0之间，最后DNA浓度应大于50 ng/μL。产品合格后送北京诺禾致源生物信息科技股份有限公司测序。

1.3 基因型数据的质控

使用Qualitystrim-1.5.3软件进行初步质控, 标准为最小的片段平均质量分数(Minimum read average quality)大于20，最短片段长度(Minumum read length)大于50，片段中未测出碱基数(Minimum N bases included)小于3，样本检出率(call rate)大于0.90。

1.4 选择信号检测

本研究首先使用Fst法来检测基因组中受到选择的区域，计算公式为：

其中，Nk=p1k(q2k-q1k)+p2k(q1k-q2k)，

Dk=p1kq2k+q1kp2k。

1.5 Pool-GWAS

本研究利用基于池测序数据的全基因组关联分析(Pool-GWAS)，探究了SNP与多浪羊繁殖力的相关性。首先将SNP信息整理成文本文件，然后利用Popoolation 2软件中的Fisher精确性检验对多浪羊的繁殖力性状进行全基因组关联分析。得到初步结果后，再利用FDR校正所得P值，根据校正后的P值来判断显著性，显著性水平为0.05。

1.6 基因注释与富集分析

利用Ensemble数据库进行基因注释，标准为：当某个基因与显著SNP的上下游50 kb范围内有1个以上碱基重叠时认为该基因与性状显著关联。

为了了解受到选择的基因的生物学功能，对受选择区域内的基因进行富集分析，利用David在线平台对基因的分子功能、生物学过程、细胞组分进行富集分析。然后从富集分析的结果中挑选与绵羊繁殖性状有关的基因。

2 结果与分析

2.1 群体SNP

本研究利用Illumina公司高通量测序平台对多浪羊基因组序列进行测序，并比对出SNP的基因型，多胎池共得到13812250个SNP，单胎池也得到了13633661个SNP，其中有10447450个SNP为2个池共同拥有。不同池的各染色体数目和平均间距见表1。

2.2Fst值统计

将SNP利用Fst检验，本研究得到了多浪羊样本群体全基因组范围内的一系列Fst值，各染色体上Fst的分布情况见图1。

将所有Fst值按降序排列之后，取前0.1%的位点为显著位点，最大值为0.8832，阈值为0.3437，共筛选出135个显著窗口。通过观察基因组Fst密度曲线图，可以发现大部分的Fst值都集中在0.35～0.50之间(图2)。

2.3 基因注释与富集分析

将筛选到的显著窗口利用ENSEMBLE数据库平台进行基因注释，共得到51个基因，说明这些基因在多胎池和单胎池个体之间出现了分化。将基因在DAVID在线平台上进行富集分析。采取滑动窗口的计算方式对基因组Fst值又进行了一次计算，以500 bp为窗口大小，筛选极端值的方法也是取前0.1%的值，然后对筛选到的区域在ENSEMBLE上做了基因注释，接着利用在线工具DAVID对通过基因注释得到的基因进行富集分析，在富集分析得到的结果中筛选了与繁殖有关的显著条目，包括有性生殖、生殖行为、卵泡发育、配子发生、排卵周期、卵母细胞成熟等等(表3)。

表1 多胎池和单胎池各染色体上的SNP信息统计

图1 各染色体上Fst值的统计

图2 基因组Fst值密度曲线

表3 500 bp选择信号基因富集条目汇总

2.4 基因的筛选和主要功能

选择信号中Fst的最大值为0.8710，位于3号染色体32024500～32025000位置，经过基因注释后发现该窗口位于一个基因内部，即NCOA1基因。在雌配子发生条目中(表3)，多胎组和单胎组分别富集到8和11个基因，两组间只有1个相同基因：INHBA基因。同样利用ENSEMBLE数据库进行基因注释，共有28个基因与显著窗口有交集，如ACP1、ING5、ARNT等基因，以上部分基因的位置和部分功能见表4。

3 讨论

本研究基于全基因组重测序数据信息，分别对多浪羊多胎和单胎群体进行了分析。在样本群体中分别挑选出多胎和单胎母羊，进行混池测序，利用Fst方法检测选择信号。结果共得到51个基因。通过文献资料发现这些基因可能和多浪羊的一些重要经济性状的生物学过程相关，例如IL1RN、IL36A等，它们均属于干扰素基因家族，它们的主要功能是产生干扰素来调控细胞因子，进而调控细胞和组织的生长和分化[7]。

表4 基因的位置及其主要功能

通过观察基因组Fst密度曲线图，可以发现大部分的Fst值都集中在0.35～0.50之间，说明多胎组和单胎组的遗传结构在很大程度上分化并不是非常严重，考虑到本研究的样本来自于同一羊场的同一群体，因此这一结果符合实际情况。以500 bp为窗口大小，筛选极端值的方法也是取前0.1%的值，经过基因注释和富集分析筛选了NCOA1、INHBA等一些候选基因。其中NCOA1在动物多种组织中表达[8]，通过为其他辅酶因子提供结合位点而形成稳定的起始前复合物，激活或抑制下游基因的表达。Lango等[9]认为NCOA1基因是导致哺乳动物繁殖性能增加的重要原因。Lindberg等[10]研究表明：NCOA1 基因上突变与梅山猪的产仔性状密切相关，同时与母猪的排卵率显著相关。张玉龙等[11]基于INHBA基因设计了两对引物，结果在1042 bp处发现了一处突变，最小二乘法分析显示与徐淮山羊没有相关性。而Leyhe等[12]报道：当绵羊品种平均产羔数增加时，绵羊INHBA基因座TaqI A等位基因频率也增加。故结合本研究的结果，NCOA1、INHBA基因可以作为多浪羊高繁殖力候选基因用于后续研究。

利用500 bp窗口找到的选择信号，经过基因富集分析后，重点关注了多浪羊的卵泡发育与成熟、卵子发生、胚胎发育以及生殖激素等相关条目。在多胎池中，以上条目里富集到了一些共同的基因，比如REC8、RHOXF1、BMP等基因家族，因此可以认为这些基因受到了选择。同样，在单胎池中，则存在BCL2L1、BAX、CCNB1、EREG、CDC25B这些共有的基因，但是这并不能说明以上这些基因的变异与多浪羊繁殖力的差异有直接的关系。不过，多胎池中的这些基因几乎都是与生长发育有关的促发育基因，比如BMP4基因，该基因属于BMP蛋白家族的一员，是一种生长分化因子，能够促进卵泡的发育。而在单胎池中，虽然有像EREG、CDC25B这样能促进发育和细胞分化的常规促发育基因，但是也发现了编码细胞凋亡因子的基因，如BCL2L1、BAX，因此也许可以推断出多浪羊的高繁殖力有可能是由于该类编码细胞凋亡因子的基因在与生殖有关的生物学过程中的表达受到了抑制而导致的。后续应通过测序等方法检测以上基因是不是控制多浪羊产羔数的主效基因，或与主效基因是否存在连锁。这些基因所对应的SNPs位点并未达到全基因组显著水平，这可能是因为统计方法、样本含量等多种因素造成的，但它们同样是与绵羊体重性状相关且有意义的候选基因。由于本研究样本含量不是很大，涉及的品种也较单一，所获得的结论只是初步的，因此还需要扩大样本数和增加绵羊品种数对遗传标记与产羔性能之间的相关性做进一步研究。

[1] 柳淑芳,姜运良,杜立新.BMPR-IB和BMP15基因作为小尾寒羊多胎性能候选基因的研究[J].遗传学报,2003,30(8):755-759.

[2] 褚明星,桑林华,王金玉，等.小尾寒羊高繁殖力候选基因BMP15和GDF9的研究[J].遗传学报,2005,32(1):38-45.

[3] Demars J, Fabre S, Sarry J, et al. Genome-wide association studies identify two novelBMP15 mutations responsible for an atypical hyperprolificacy phenotype in sheep[J]. Plos Genetics, 2013, 9(4): e1003482.

[4] 史洪才,高志英,牛志刚，等.新疆多浪羊FecB突变检测及与产羔数的关系[J].农业生物技术学报,2011,19(2):330-334.

[5] Rubin C J, Zody M C, Eriksson J, et al. Whole-genome resequencing reveals loci under selection during chicken domestication[J]. Nature, 2010, 464(25): 587-591.

[6] Zhu C, Fan H, Yuan Z, et al. Detection of selection signatures on the x chromosome in three sheep breeds[J]. International Journal of Molecular Sciences, 2015, 16(9): 20360-2374.

[7] Dinarello C A. Biologic basis for interleukin-1 in disease[J]. Blood, 1996, 87(6): 2095-2147.

[8] Misiti S, Schomburg L, Yen P M, et al. Expression and hormonal regulation of coactivator and corepressor genes[J]. Endocrinology, 1998, 139(5): 2493-2500.

[9] Lango A H, Estrada K, Lettre G, et al. Hundreds of variants clustered in genomic loci and biological pathways affect human height[J]. Nature, 2010, 467(7317): 832-838.

[10] Lindberg M K, Weihua Z, Andersson N, et al. Estrogen receptor specificity for the effects of estrogen in ovariectomized mice[J]. J Endocrinal, 2004, 174(25): 167-178.

[11] 张玉龙,孙伟,苏锐,等.徐淮山羊PRLR和INHBA基因的多态性及其与产羔数的关联分析[J].中国畜牧杂志,2014,50(9):14-18.

[12] Leyhe B, Hiendleder S, Jaeger C, et al. Pronounced differences in the frequency ofTaqⅠ βA-inhibin alleles between sheep breeds with different reproductive performance[J]. Animal Genetics, 1994, 25(1): 41-43.

(责任编辑:曾小军)

Whole Genomic Association Analysis of Polyembryonic Gene in Duolang Sheep

WANG Shan1,2, LI Xiao-lin1, NIU Zhi-gang2, SHI Hong-cai2

(1. College of Animal Science and Technology, Shihezi University of Xinjiang, Shihezi 832003, China; 2. Key Opening Laboratory of Genetic Breeding and Reproduction of Livestock, Ministry of Agriculture / Biotechnology Research Center, Xinjiang Academy of Animal Sciences, Urumqi 830000, China)

The whole genomic association analysis (GWAS) of lambing traits of Duolang sheep was conducted by using the whole genomic resequencing technique, and the candidate genes and molecular markers related to the breeding traits of sheep were excavated. TheFsttest screened out 135 significant windows, and theirFstvalue was mostly concentrated within 0.35～0.50. A total of 51 genes were obtained by the gene annotation to the selected significant windows. The further enrichment analysis found that the maximum value ofFstwas 0.8710, and it was located within the geneNCOA1. In the female gametogenic clauses, the common gene of polyembryonic group and unigerminal group wasINHBA. Through using ENSEMBLE database for gene annotation, a total of 28 genes (such asACP1,ING5,ARNT, etc.) had intersection with the significant window. The disclosure of these candidate genes and molecular markers has reference value and theoretical significance for the excavation of functional genes related to sheep reproductive traits.

Whole genomic association analysis; Duolang sheep; Polyembryony

2016-11-16

王珊，女，新疆乌鲁木齐人，硕士研究生，研究方向：绵羊分子育种。

S827

A

1001-8581(2017)05-0077-05