多浪羊多胎基因全基因組關聯分析研究

王 珊，李曉林，牛志剛，史洪才

(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