山羊產羔數全基因組關聯分析

蘭 蓉，朱 蘭，姚新榮，王 鵬，邵慶勇，洪瓊花*

(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