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

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

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

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

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

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

1 材料與方法

1.1 試驗動物

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

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

1.2 SNP芯片分型及數據質量控制

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

1.3 群體結構及連鎖不平衡分析

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

1.4 ROH檢測及基因組近交系數的計算

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

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

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

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

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

2 結 果

2.1 SNP質控結果及群體遺傳結構和連鎖不平衡分析

在質控后每個個體保留了96 789個SNPs位點,相鄰SNPs之間的平均距離為25.72 kb,以供后續(xù)分析。圖1A顯示了7個牧場荷斯坦牛群體的PCA分析結果。從圖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 群體遺傳結構及連鎖不平衡Fig.1 Population genetic structure and linkage disequilibrium

2.2 ROH數目、長度及分布的統計

由表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的數目。

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

2.3 基因組近交系數評估

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

表2 基于ROH的不同奶牛場的近交系數(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,以確定統計閾值,本研究選取頻率大于29.78%作為高頻率的ROH區(qū)域閾值。共檢測到8個高頻區(qū)域,并通過Ensembl數據庫對ROH中的高頻區(qū)域進行基因注釋,共注釋到79個基因,見表3。其中,14號染色體上22.78～23.38 Mb位置的區(qū)域,80%的個體都在該區(qū)域內發(fā)生ROH片段,并注釋到3個基因。利用KOBAS對注釋到的基因進行KEGG通路富集分析,結果見表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基本統計分析

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

3.2 基于ROH的基因組近交系數

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

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

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

4 結 論

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