基于SNP芯片數(shù)據(jù)分析不同奶牛場基因組近交系數(shù)及篩選功能性基因

王振宇,張賽博,劉文慧,梁 棟,任小麗,閆 磊,閆躍飛,高騰云,張 震,3*,黃河天*

(1.河南農(nóng)業(yè)大學(xué)動物科技學(xué)院,鄭州 450046;2.河南省奶牛生產(chǎn)性能測定中心,鄭州 450045;3.河南省種業(yè)發(fā)展中心,鄭州 450046)

基因組長純合片段(runs of homozygosity, ROH)一般存在于二倍體生物中,它是親代將單倍型基因中同源相同(identity by descent, IBD)的片段遺傳給子代,并且在子代的基因組中形成連續(xù)性的純合片段[1],即子代從親代繼承了同源的染色體片段,從而導(dǎo)致后代基因組中的純合片段產(chǎn)生并上升到ROH[2]。連鎖不平衡、種群瓶頸、遺傳漂變、近親交配和選擇都可能是引起ROH產(chǎn)生的因素[1,3-4]。不同的群體歷史會產(chǎn)生不同長短的ROH,長片段ROH通常由群體近幾個世代近交產(chǎn)生,短片段ROH通常來自更遠(yuǎn)的祖先[5-7]。因此,通過全基因組ROH特征的檢測,可以了解種群歷史、結(jié)構(gòu)、近交情況。

ROH最早在人類染色體基因組發(fā)現(xiàn),并被認(rèn)為可能對人類健康有重要影響。隨著ROH在人類群體遺傳學(xué)中研究的深入[8-10],不同畜禽的ROH分析研究也逐漸開展[11-13]?；赗OH估計基因組近交系數(shù)已成為利用全基因組信息評估近交的常用方法,即利用ROH計算基因組近交系數(shù)FROH(inbreeding calculated from ROH),它可以準(zhǔn)確計算個體近交系數(shù)?，F(xiàn)已有多項研究證明了基于系譜信息計算的近交系數(shù)要低于真實(shí)的近交系數(shù)。楊湛澄等[14]利用牛54K SNP芯片數(shù)據(jù)對北京地區(qū)2 107頭荷斯坦?；蚪MROH分布進(jìn)行了統(tǒng)計,并計算了基因組近交系數(shù)和系譜近交系數(shù),發(fā)現(xiàn)基于ROH計算的基因組近交系數(shù)能更準(zhǔn)確地反映個體的真實(shí)近交情況。Peripolli等[15]利用770K SNP芯片數(shù)據(jù)比較了2 908頭吉爾牛(Gyr)基于ROH(FROH)、基因組關(guān)系矩陣(genomic relationship matrix,FGRM)、基因組純合子百分比(homozygosity,FHOM)、系譜信息(pedigree,FPED)4種方法計算的近交系數(shù),結(jié)果表明在沒有系譜記錄的情況下,FROH可用作近交估計的替代方法。此外,通過識別群體的高頻ROH片段,鑒定到了與產(chǎn)奶量、乳成分、熱適應(yīng)相關(guān)的基因。Nani和Peagaricano[16]研究發(fā)現(xiàn),基因組ROH與荷斯坦公牛繁殖性狀顯著相關(guān),公牛群體中高度純合的基因組區(qū)域與公牛繁殖性狀呈現(xiàn)負(fù)相關(guān),并在低繁殖力公牛ROH富集區(qū)域鑒定到與精子生物學(xué)和雄性生育能力密切相關(guān)的基因。Liu等[17]利用簡化基因組測序的方法,通過ROH與綜合單倍型評分(integrated haplotype score, iHS)分析,檢測到與上海荷斯坦奶牛群體健康、繁殖、環(huán)境適應(yīng)等有關(guān)的候選基因。通過對全基因組ROH進(jìn)行檢測,可以更準(zhǔn)確地掌握群體的近交程度,幫助研究者在育種實(shí)踐中制定科學(xué)合理的選種選配方案。鑒定全基因組的ROH也可以更好的了解ROH在染色體上的分布規(guī)律,進(jìn)而挖掘可能影響畜禽重要性狀的候選基因[18-20]。

在我國,北京[14]、上海[17]、寧夏[21]基于荷斯坦牛群體基因組ROH估算群體近交系數(shù)、檢測與經(jīng)濟(jì)性狀相關(guān)候選基因及選育過程中的選擇信號等的研究,為中國荷斯坦奶牛育種提供了重要數(shù)據(jù)參考。然而,通過基因組ROH信息估計不同牧場荷斯坦奶牛群體近交水平和檢測群體選擇特征的研究仍然較少。本研究旨在利用奶牛150K SNP芯片數(shù)據(jù)對河南省7個奶牛場荷斯坦牛進(jìn)行全基因組ROH檢測,計算ROH的長度、頻率、數(shù)目和分布以及基因組近交系數(shù)FROH,比較不同牧場荷斯坦牛基因組近交程度,并在高頻ROH區(qū)域注釋與荷斯坦牛經(jīng)濟(jì)性狀相關(guān)的候選基因。以期為詳細(xì)了解河南省荷斯坦牛群體基因組ROH分布特征及基因組近交程度,為牧場今后選種選配提供參考。也可通過ROH富集區(qū)域鑒定一些與奶牛經(jīng)濟(jì)性狀相關(guān)的基因,為奶牛標(biāo)記輔助選擇提供候選基因信息,為奶牛場科學(xué)選種選配提供指導(dǎo)。

1 材料與方法

1.1 試驗動物

根據(jù)系譜、生產(chǎn)數(shù)據(jù)記錄的完整性,篩選出7個存欄量在150～5 000頭的規(guī)?；翀?按存欄量10%的比例抽取牧場核心群個體進(jìn)行血液樣本采集,最終共采集了900頭荷斯坦牛。具體樣本分布情況詳見表1。

表1 不同奶牛場荷斯坦牛ROH長度和數(shù)量Table 1 The mean length and number of runs of homozygosity(ROH) in Holstein of different dairy farms

1.2 SNP芯片分型及數(shù)據(jù)質(zhì)量控制

采集尾椎靜脈血,提取DNA,利用GGP Bovine 150K芯片進(jìn)行基因分型。用PLINK(v1.90)[22]對原始數(shù)據(jù)進(jìn)行質(zhì)控,設(shè)定條件:1)SNP檢出率大于95%;2)個體檢出率大于99%;3)最小等位基因頻率大于0.01;4)哈迪-溫伯格平衡P值大于10-6;5)保留常染色體數(shù)據(jù)。

1.3 群體結(jié)構(gòu)及連鎖不平衡分析

基于SNP信息,使用GCTA(v1.93)軟件[23]對900頭荷斯坦牛群體進(jìn)行主成分分析(principal component analysis, PCA)。采用PopLDdecay(v3.42)軟件[24]計算每個牧場的連鎖不平衡(linkage disequilibrium, LD)程度,并使用軟件自帶的Plot_MultiPop.pl腳本繪制LD衰減曲線圖。

1.4 ROH檢測及基因組近交系數(shù)的計算

ROH檢測使用PLINK軟件[22],使用滑動窗口的方法對常染色體進(jìn)行檢測,具體檢測參數(shù)如下:1)滑動窗口閾值使用0.05;2)滑動窗口設(shè)置50個SNPs位點(diǎn);3)每一個滑動窗口中允許丟失的基因型為5個;4)每一個滑動窗口中允許的雜合子數(shù)目為1個;5)組成ROH的SNP的最大間隔為1 Mb;6)組成ROH的SNP的最低密度為每50 kb 1個SNP;7)ROH片段的最小長度設(shè)為500 kb;8)每個ROH至少由50個SNPs組成。

利用ROH計算近交系數(shù)(FROH),公式如下:

其中,∑LROH為常染色體上ROH片段長度之和,Lgenome為常染色體基因組物理長度之和(2.49 Gb)。

1.5 高頻ROH區(qū)域候選基因鑒定

使用R語言統(tǒng)計每個SNP在奶牛群體中參與組成ROH的次數(shù)占樣本數(shù)的比例,并將前1%的SNPs區(qū)域作為高頻的ROH區(qū)域。基于高頻ROH區(qū)段的物理位置,并通過生物數(shù)據(jù)庫Ensembl[25]中的BioMart模塊與牛參考基因組(Bos_taurus.ARS-UCD1.2)進(jìn)行比對,檢索基因,然后依據(jù)NCBI(https://www.ncbi.nlm.nih.gov/)、GeneCards(https://www.genecards.org/)網(wǎng)站及文獻(xiàn)查詢基因功能。運(yùn)用KOBAS(http://bioinfo.org/kobas/)[26]在線數(shù)據(jù)庫對注釋到的基因進(jìn)行KEGG通路富集分析,當(dāng)P<0.05時,則表示顯著富集。

2 結(jié) 果

2.1 SNP質(zhì)控結(jié)果及群體遺傳結(jié)構(gòu)和連鎖不平衡分析

在質(zhì)控后每個個體保留了96 789個SNPs位點(diǎn),相鄰SNPs之間的平均距離為25.72 kb,以供后續(xù)分析。圖1A顯示了7個牧場荷斯坦牛群體的PCA分析結(jié)果。從圖1可以看出,7個牛場主要分為了5個亞群。采用PopLDdecay分別計算各牧場群體的成對r2值,用于比較不同荷斯坦牛群體的LD水平(圖1B)。LD分析顯示,7個牧場奶牛群體LD衰減的順序為:H7>H4&H5>H2&H3&H6>H1。

A.主成分分析圖;B.LD衰減圖。H1～H7代表牧場編號A.Principal component analysis of Holstein cattle population;B. LD decay of Holstein cattle population. H1-H7 represents pasture number圖1 群體遺傳結(jié)構(gòu)及連鎖不平衡Fig.1 Population genetic structure and linkage disequilibrium

2.2 ROH數(shù)目、長度及分布的統(tǒng)計

由表1可以看出,在7個牧場荷斯坦牛群體中共檢測出55 908個ROH,ROH的平均長度為4.23 Mb,范圍在1.90～14.07 Mb。其中H6號牛場ROH平均長度最小,為3.27 Mb;H2號牛場ROH平均長度最大為4.49 Mb。在0～5 Mb長度上,ROH總體比例占76.21%,其中H1、H6牧場ROH比例較大(83.70%、84.30%),其余牧場ROH比例范圍為73.33%～76.52%;在5～10 Mb長度上,ROH總體比例占15.14%,其中H1、H6牧場ROH比例較小(10.26%、10.67%),其余牧場ROH比例范圍為14.89%～17.06%;在>10 Mb長度上,ROH總體比例占8.64%,其中H1、H6牧場ROH比例較小(6.03%,5.04%),其余牧場ROH比例范圍為7.61%～9.61%。圖2展示了常染色體上不同長度ROH的數(shù)目。

圖2 染色體上不同長度ROH的數(shù)目Fig.2 Number of ROH with different length on chromosome

2.3 基因組近交系數(shù)評估

不同牧場荷斯坦牛群體基于ROH的近交系數(shù)及變化范圍見表2。全群中基于ROH的基因組FROH范圍為0.021～0.447,近交系數(shù)平均值為0.106,標(biāo)準(zhǔn)差為0.040。其中H2號牧場平均FROH最高(0.123),H7號牧場平均FROH最低(0.082),其他牧場分別為0.112、0.114、0.109、0.108、0.103。在個體層面中,FROH最低的個體出現(xiàn)在H7號牛場中(0.021),FROH最高的個體出現(xiàn)在H4號牛場中(0.447)。

表2 基于ROH的不同奶牛場的近交系數(shù)(FROH)Table 2 Inbreeding coefficient (FROH) of different dairy farms based on ROH

2.4 高頻ROH區(qū)域及候選基因鑒定與注釋、富集

圖3展示了在1～29號染色體上組成ROH的SNPs占群體的百分率。通過選擇組成ROH中前1% SNPs,以確定統(tǒng)計閾值,本研究選取頻率大于29.78%作為高頻率的ROH區(qū)域閾值。共檢測到8個高頻區(qū)域,并通過Ensembl數(shù)據(jù)庫對ROH中的高頻區(qū)域進(jìn)行基因注釋,共注釋到79個基因,見表3。其中,14號染色體上22.78～23.38 Mb位置的區(qū)域,80%的個體都在該區(qū)域內(nèi)發(fā)生ROH片段,并注釋到3個基因。利用KOBAS對注釋到的基因進(jìn)行KEGG通路富集分析,結(jié)果見表4。分析得出79個基因顯著富集于酮體的合成與降解(synthesis and degradation of ketone bodies)、纈氨酸、亮氨酸和異亮氨酸降解(valine, leucine and isoleucine degradation)、丁酸代謝(butanoate metabolism)、Ras信號通路(ras signaling pathway)等11個信號通路。

表3 荷斯坦牛高頻ROH區(qū)域及候選基因Table 3 High-frequency ROH regions and candidate genes in Holstein cattle

表4 高頻ROH區(qū)域基因的KEGG通路富集分析(P<0.05)Table 4 KEGG pathway enrichment analysis of genes in high-frequency ROH regions(P<0.05)

圖3 ROHs中SNPs百分比曼哈頓圖Fig.3 Manhattan plot of SNPs percentages in ROHs

3 討 論

3.1 荷斯坦牛群體基因組ROH基本統(tǒng)計分析

不同育種目標(biāo)及選擇強(qiáng)度會引起不同荷斯坦牛群體中ROH數(shù)目、長度及分布情況的差異[5-6,27]。Kim等[7]通過比較3個北美荷斯坦牛群體在產(chǎn)奶性狀不同選擇強(qiáng)度下基因組ROH的變化,揭示了總體ROH頻率和分布方面的顯著差異,結(jié)果顯示群體內(nèi)ROH平均長度約為6 Mb,小于5 Mb的ROH片段數(shù)目占總片段數(shù)目的53%。而與Kim等[7]的研究結(jié)果相比,本研究中荷斯坦牛群體ROH平均長度為4.23 Mb,小于5 Mb的ROH片段的數(shù)目占總片段數(shù)目的76.21%。另外對比不同牧場群體,小于5 Mb的ROH片段數(shù)目所占比例也有差異。在基因組ROH長度上,Marras等[28]利用50K SNP芯片對5個意大利公牛品種進(jìn)行ROH分析,結(jié)果表明相較于其他品種,乳用品種荷斯坦牛和意大利布朗牛的平均ROH長度更大(3.6、3.9 Mb),其中荷斯坦牛群體的ROH平均長度與本研究的結(jié)果相近。在牧場群體方面,H1和H6號牧場群體在小于5 Mb的ROH片段數(shù)目占總片段數(shù)目最高(83.70%、84.30%),而大于10 Mb的ROH片段數(shù)目占總片段數(shù)目比例最低(6.03%、5.04%)。研究顯示,較近世代的共同祖先會造成長ROH片段的形成,短的ROH來源于關(guān)系較遠(yuǎn)的共同祖先[7,29]。此外,各個牧場奶牛群體ROH平均長度、變化范圍也有差異,這與不同牧場奶牛群體來源以及選配過程中使用不同國別的冷凍精液有關(guān)。因此,本研究基于對不同牧場群體基因組ROH的數(shù)目、長度及分布的研究,評估群體近交情況,為牧場今后的選種選配提供參考。

3.2 基于ROH的基因組近交系數(shù)

目前,ROH常用來計算個體近交系數(shù),且具有較高的準(zhǔn)確性[15,30-33]。本研究中,河南荷斯坦牛群體總平均FROH(0.106)與寧夏[21](0.101)、北京[14](0.007～0.312)荷斯坦牛群體FROH相近,與上海[17]荷斯坦牛群體(0.363)相差較大。上海與北京作為我國的南、北奶牛養(yǎng)殖業(yè)的代表地區(qū),由于選育目標(biāo)、強(qiáng)度、氣候等因素的影響,群體近交程度出現(xiàn)差異,河南地理位置上屬于中原地區(qū),在奶牛育種策略和群體近交情況上與北方更相近。近交水平在一定程度上也可以反映牧場選種選配管理狀況。在牧場選配管理上,由表2可以看到,H1、H2、H3號牧場平均FROH較高(0.112、0.123、0.114),H7號牧場平均FROH較低(0.082),不同牧場之間的差異側(cè)面反映出這些牧場在選配過程中對群體近交問題的管理程度;在牧場規(guī)模上,H1、H2、H3號牧場規(guī)模較小,群體數(shù)量較少,平均FROH較高(0.112、0.123、0.114),H4號牧場規(guī)模較大,群體數(shù)量多,平均FROH較低(0.109)。此外,在H4號牧場中有些個體的FROH明顯較高(>0.285),最大FROH達(dá)到0.458,反映出該牧場在個體選種選配過程中未充分考慮近交問題。因此,通過對近交系數(shù)的計算可以了解不同牧場群體近交狀況,從而在實(shí)際選種選配工作中能更有效的避免近交,減少經(jīng)濟(jì)損失。

3.3 基因組高頻ROH區(qū)域的候選基因分析

本研究在高頻ROH區(qū)域中共鑒定到了79個基因,其中包含與奶牛經(jīng)濟(jì)性狀有關(guān)的基因,如AKAP3、C5H12orf4、CAPN3、ARL15、XKR4、CRBN、IL5RA等。5號染色體上AKAP3、C5H12orf4、FGF6基因與體型、體高有關(guān)[34-36]。10號染色體上CAPN3基因與胴體、繁殖性狀相關(guān)[37-38]。CHST14基因與妊娠維持和胎兒生長直接相關(guān)[39]。22號染色體上IL5RA基因影響牛奶蛋白質(zhì)組成[40]。此外還有一些基因與繁殖、生長等性狀有關(guān),如FGF10基因參與調(diào)節(jié)胎兒卵泡生成[41]。值得注意的是,14號染色體上22.78～23.38 Mb區(qū)域是ROH頻率最高的區(qū)域,80%的個體都在該區(qū)域內(nèi)發(fā)生ROH片段(圖3)。發(fā)現(xiàn)該區(qū)域與寧夏[21]荷斯坦牛群體高頻區(qū)域(21.61～24.99 Mb)高度重合,這可能與不同地區(qū)育種目標(biāo)及選擇強(qiáng)度有關(guān),并隨著選育的推進(jìn),在基因組中出現(xiàn)相近的長純合區(qū)域。這個高頻區(qū)域注釋到TGS1、LYN、CHCHD7基因,這些基因與生長、胴體相關(guān)性狀[42-43]和飼料效率有關(guān)[35,44-45]。因此,本研究在ROH富集區(qū)域鑒定的基因可以為荷斯坦奶牛分子育種提供候選基因信息。

4 結(jié) 論

本研究對河南省荷斯坦牛群體進(jìn)行全基因組ROH檢測與分析,發(fā)現(xiàn)ROH在不同牧場群體中的數(shù)目、長度及頻率存在差異,基于ROH計算的近交系數(shù)范圍在0.082～0.123,反映出不同牧場近交水平存在差異,這有助于了解河南省荷斯坦牛群體近交程度,為牧場選育過程中避免近交提供指導(dǎo)。在全基因組范圍內(nèi)檢測到8個高頻ROH富集區(qū)域,共篩選出79個與奶牛經(jīng)濟(jì)性狀相關(guān)的基因,如AKAP3、C5H12orf4、CAPN3、ARL15、XKR4、CRBN、IL5RA等,可作為奶牛分子育種中進(jìn)行標(biāo)記輔助選擇的候選基因。