川中黑山羊基因組近交系數(shù)和選擇信號(hào)分析

郭家中，孫學(xué)良，向秋楠，肖鳳敏，張紅平

（四川農(nóng)業(yè)大學(xué)動(dòng)物科技學(xué)院，成都 611130）

我國(guó)地方畜禽品種資源豐富，但群體規(guī)模通常較小，因此需要控制近交從而避免近交衰退[1]。傳統(tǒng)上，個(gè)體的近交系數(shù)（inbreeding coefficients，F(xiàn)）是利用系譜信息進(jìn)行估計(jì)[2]；但這類(lèi)方法主要有兩方面不足：一是需假定祖先個(gè)體是非近交個(gè)體，這往往背離事實(shí)從而低估個(gè)體的近交水平；另一方面，系譜近交系數(shù)反映特定親緣關(guān)系下等位基因同源相同的期望值，無(wú)法準(zhǔn)確反映個(gè)體真實(shí)的近交水平[3]。此外，很多群體系譜記錄不完善甚至缺乏系譜。而利用基于芯片或重測(cè)序的SNP數(shù)據(jù)估計(jì)的近交系數(shù)在全基因組水平上反映了個(gè)體真實(shí)的等位基因同源性的概率或相關(guān)性，克服了基于系譜方法的局限性，這類(lèi)近交系數(shù)統(tǒng)稱(chēng)為基因組近交系數(shù)。基因組近交系數(shù)估計(jì)方法較多，但從原理上主要包括基于純合子過(guò)度的估計(jì)值（excess of homozygosity inbreeding coefficient，F(xiàn)HOM）[4]、長(zhǎng) 純 合 片 段（runs of homozygosity，ROH）占基因組比例（FROH）[5]和基于頻率校正的多位點(diǎn)平均純合度[6]。在實(shí)踐中，第三類(lèi)方法理論上是基于基因組親緣關(guān)系矩陣（genomic relationship matrix，GRM）而獲得近交系數(shù)（FGRM）。目前，關(guān)于牛[7]、豬[8]和綿羊[9]等物種的基因組近交系數(shù)及不同方法的報(bào)道較多，但是在山羊[10]上研究較少。

全基因組選擇信號(hào)研究已成為鑒定影響畜禽重要性狀遺傳變異位點(diǎn)的一種重要方法，并在山羊上已有諸多報(bào)道[10-12]。川中黑山羊（Chuanzhong Black goat，CZ）是我國(guó)優(yōu)良地方山羊品種，又進(jìn)一步分為金堂型和樂(lè)至型，該品種具有體型大、生長(zhǎng)速度快和產(chǎn)羔數(shù)多等優(yōu)點(diǎn)[13]。該品種在西南地區(qū)很受歡迎，被多個(gè)地區(qū)引入作為父本，改良地方黑山羊品種的生長(zhǎng)性狀；但關(guān)于川中黑山羊的遺傳研究相對(duì)較少。我們前期利用9只公羊的重測(cè)序數(shù)據(jù)，對(duì)該品種進(jìn)行了群體遺傳分析[14]；但樣本量太少。因此有必要擴(kuò)大樣本量對(duì)川中黑山羊近交水平、有效群體大小和選擇信號(hào)等群體遺傳特征進(jìn)行更深入的探究。

本研究利用前期收集的41只川中黑山羊（金堂型）的重測(cè)序數(shù)據(jù)獲得了高密度SNP基因分型數(shù)據(jù)，主要開(kāi)展基因組近交系數(shù)、有效群體大小估計(jì)和選擇信號(hào)分析，旨在為川中黑山羊遺傳資源的發(fā)掘利用提供理論參考。

1 材料和方法

1.1 試驗(yàn)群體和高通量重測(cè)序數(shù)據(jù)

川中黑山羊（金堂型）樣本均來(lái)自于成都市某公司育種場(chǎng)，共包括41只種羊，其中9只公羊、32只母羊。我們前期研究已經(jīng)發(fā)表并公布了上述41個(gè)個(gè)體的短序列基因組重測(cè)序數(shù)據(jù)[14-15]。為進(jìn)行選擇信號(hào)分析，還使用了前期獲得的30只建昌黑山羊的重測(cè)序數(shù)據(jù)[10]并下載了21只野山羊的重測(cè)序數(shù)據(jù)（NCBI accession number:PRJEB3136）[16]。

1.2 高通量測(cè)序數(shù)據(jù)處理和短變異檢測(cè)

參考我們之前使用的分析流程[10]，進(jìn)行重測(cè)序數(shù)據(jù)比對(duì)和短變異檢測(cè)。使用BWA軟件（v 0.7.17）[17]將短序列映射到山羊參考基因組（ARS1[18]，GCA_001704415.1），獲得比對(duì)結(jié)果。使用GATK軟件（v 4.0.5.2）[19]初步檢測(cè)SNPs和Indels并進(jìn)行硬過(guò)濾，再使用VCFtools[20]進(jìn)行群體遺傳學(xué)質(zhì)控，僅保留最小等位基因頻率大于0.05、缺失率小于10%的位點(diǎn)；在進(jìn)行各項(xiàng)遺傳分析前，進(jìn)一步過(guò)濾掉哈迪溫伯格平衡P值小于10-10的雙等位基因SNPs。使用SnpEff軟件（v 4.3）[21]對(duì)短變異進(jìn)行功能注釋。

1.3 有效群體大小和全基因組ROH分析

使用 PopLDdecay（v 3.4.1）軟件[22]分析群體中SNPs之間的連鎖不平衡（linkage disequilibrium，LD）。使用基于連鎖不平衡理論開(kāi)發(fā)的SNeP（v 1.1）[23]和默認(rèn)參數(shù)估算山羊歷史有效群體大?。╡ffective population size，Ne）。使用PLINK（v 1.9）軟件[24]中的“--homozyg”命令檢測(cè)川中黑山羊基因組中的長(zhǎng)純合片段（runs of homozygosity，ROH）。主要參數(shù)值設(shè)置如下：每個(gè)ROH片段最少包含10個(gè)SNPs，最小長(zhǎng)度為100 kb，最低SNP密度為10 kb/SNP，每個(gè)ROH內(nèi)SNP的最大間隔為100 kb；滑動(dòng)窗口大小為50個(gè)SNPs，每個(gè)滑動(dòng)窗口中最多允許1個(gè)雜合位點(diǎn)和5個(gè)缺失位點(diǎn)，滑動(dòng)窗口閾值為0.05。

1.4 基因組近交系數(shù)的估計(jì)

基于ROH的近交系數(shù)（FROH）定義為每個(gè)個(gè)體基因組的ROH總長(zhǎng)度占常染色體基因組長(zhǎng)度的比例（參考基因組ARS1的常染色體總長(zhǎng)度為2 466 191 353 bp）。根據(jù)ROH長(zhǎng)度與世代數(shù)之間的關(guān)系，依據(jù)不同長(zhǎng)度將ROH分成4類(lèi)：0.1～0.2、0.2～0.5、0.5～1.0和＞1 Mb，計(jì)算對(duì)應(yīng)世代的FROH。使用PLINK 軟件[24]中的“--het”命令計(jì)算FHOM。Van-Raden在2008年提出了3種GRM計(jì)算方法[25]，并被廣泛使用。GMAT軟件[26]實(shí)現(xiàn)了VanRaden提出的第一種方法，使用該軟件中的“--grm agrm--fmt 0”命令獲得GRM中對(duì)角線元素，估計(jì)FVR1。使用GCTA[27]（v 1.92）“--ibc”命令計(jì)算FVR2和FUNI。

1.5 選擇信號(hào)鑒定與功能富集分析

本研究使用ROH島、iHH12和XP-EHH 3種方法鑒定川中黑山羊基因組選擇信號(hào)；其中，前兩種方法分別基于群體內(nèi)的長(zhǎng)純合片段和單倍型分布特征，第3種方法則是比較群體間擴(kuò)展單倍型純合度的分化程度?；赑LINK獲得的ROH結(jié)果，使用R包detectRUNS[28]進(jìn)行ROH島檢測(cè)，鑒定標(biāo)準(zhǔn)包含3個(gè)指標(biāo)：①個(gè)體間共享的ROH出現(xiàn)的頻率大于50%（即在41個(gè)個(gè)體中至少21個(gè)個(gè)體攜帶相同的ROH片段）；②ROH島內(nèi)至少包含10個(gè)SNPs；③ROH島長(zhǎng)度最短為1 000 bp。使用selscan軟件[29]計(jì)算iHH12，主要流程如下：使用Beagle軟件[30]對(duì)SNPs進(jìn)行填充和定相，然后計(jì)算原始的iHH12值，再使用norm模塊對(duì)原始值進(jìn)行歸一化處理。最后利用自編R腳本，以10 kb窗口和10 kb步長(zhǎng)沿染色體滑動(dòng)計(jì)算每個(gè)窗口內(nèi)的iHH12平均值。將所有窗口平均值按從高到低進(jìn)行排序，前0.5%的窗口作為選擇信號(hào)的候選窗口。使用selscan軟件計(jì)算XP-EHH，主要步驟類(lèi)似于iHH12計(jì)算。由于在分析中，將川中黑山羊作為試驗(yàn)群體、野山羊作為參考群體，XP-EHH值為正值則代表某區(qū)域在川中黑山羊中受到正選擇。因此，在XP-EHH值從高到低的分布中，選擇前0.5%的窗口作為選擇信號(hào)的候選窗口。使用BEDTools[31]對(duì)3種方法鑒定的候選區(qū)域分別進(jìn)行基因注釋?zhuān)?種方法共享基因被定義為正選擇基因。

使用R包c(diǎn)lusterProfiler（v 4.4.1）[32]對(duì)正選擇基因進(jìn)行GO功能富集分析，顯著富集的GO條目篩選標(biāo)準(zhǔn)為P＜0.05。在分析中，山羊全基因組范圍的基因功能注釋信息來(lái)自AnnotationHub（編號(hào)為“AH101444”）。

2 結(jié)果與分析

2.1 全基因組SNPs和Indels檢測(cè)

在41只川中黑山羊常染色體基因組中，共檢測(cè)到14 043 333個(gè)雙等位基因SNPs、70 478個(gè)復(fù)等位基因SNPs和1 197 402個(gè)Indels。變異注釋結(jié)果表明，位于基因間區(qū)和內(nèi)含子區(qū)的SNPs的比例最高，分別為45.38%和44.15%，而外顯子區(qū)的SNPs僅占0.92%。類(lèi)似地，內(nèi)含子區(qū)和基因間區(qū)的Indels的比例也是最高，分別為45.89%和42.95%。

2.2 連鎖不平衡和有效群體大小分析

由圖1a可知，當(dāng)SNPs之間的物理距離為10 bp時(shí)，r2平均值（0.51）最高；隨后LD迅速衰減，當(dāng)SNPs之間的距離增加到1 000 bp時(shí)，r2等于0.2。圖1b可知，川中黑山羊有效群體大小持續(xù)縮減，在999世代前Ne值為5 696只，而13世代前Ne值為190只。

圖1 川中黑山羊連鎖不平衡和有效群體大小分析Figure 1 Neand linkage disequilibrium decay in Chuanzhong Black goats

2.3 基于5種方法的川中黑山羊基因組近交系數(shù)估計(jì)值

在41只川中黑山羊基因組上共檢測(cè)到47 831個(gè)ROH，ROH在1～29號(hào)染色體上均有分布；其中，1號(hào)染色體的ROH數(shù)量最多（3 287），27號(hào)染色體上ROH數(shù)量最少（742）。最長(zhǎng)的ROH位于18號(hào)染色體（29 071 104～35 631 696 bp）長(zhǎng)度為6.56 Mb，該區(qū)域包含CDH8、TK2和CMTM3等26個(gè)基因。

如圖2a所示，41只川中黑山羊FROH值范圍為0.06～0.21，平均值為0.12。在群體水平上，當(dāng)前世代川中黑山羊所累積的近交水平主要來(lái)自于250～500世代（FROH0.1-0.2Mb）和100～250世代（FROH0.2～0.5Mb）。相關(guān)分析表明，當(dāng)前世代個(gè)體間FROH的變異與50～100世代（FROH0.5-1.0Mb，r=0.94，P＜2.2×10-16）相關(guān)性最強(qiáng)。如圖2b所示，在41只川中黑山羊中，F(xiàn)VR1近交系數(shù)最小值為-0.03，最大值為0.30，平均值為0.21。FVR2近交系數(shù)最小值為-0.02，最大值為0.27，平均值為0.18。類(lèi)似地，F(xiàn)UNI近交系數(shù)最小值為-0.03，最大值為0.26，平均值為0.18。而FHOM近交系數(shù)最小值為0.07，最大值為0.27，平均值為0.19。如表1所示，除FROH與FVR1和FVR2之間無(wú)顯著相關(guān)，其余不同方法獲得的川中黑山羊近交系數(shù)之間均正呈現(xiàn)顯著性正相關(guān)（P＜0.05）。其中，F(xiàn)UNI和FVR1兩種近交估計(jì)值的線性相關(guān)最高（r=0.983，P=2.2×10-16），而FUNI和FHOM兩種估計(jì)值的相關(guān)性也較高（r=0.893，P=4.0×10-15）。

圖2 基于5種方法的川中黑山羊基因組近交系數(shù)估計(jì)值Figure2 Summary of estimated genomic inbreeding coefficients in Chuanzhong Black goats using five methods

表1 基于5種方法的川中黑山羊基因組近交系數(shù)的皮爾遜相關(guān)Table 1 The Pearson’s correlations between the genomic inbreeding coefficients in Chuanzhong Black goats based on five methods

2.4 川中黑山羊選擇信號(hào)的檢測(cè)

如圖3a所示，基于iHH12、ROH島和XP-EHH統(tǒng)計(jì)量在川中黑山羊中分別鑒定到1 218（iHH12＞7.71）、88（ROH共享率大于50%）和1 227（XP-EHH＞3.10）個(gè)離群值窗口并分別注釋到361、164和324個(gè)基因（以Ensemble ID為標(biāo)準(zhǔn)）。其中，NCAPG（chr6:37 858 170～37 903 004 bp）、LCORL（chr6:37 905 295～38 068 616 bp）、ESR1（chr9:76 096 964～76 376 135 bp）、KIT（chr6:70 711 312～70 794 841 bp）等67個(gè)基因被3種方法均檢測(cè)到，被定義為川中黑山羊的正選擇基因。

圖3 川中黑山羊選擇信號(hào)和正選擇基因Figure 3 Summary of selection signals and positively selected genes in Chuanzhong Black goats

上述67個(gè)正選擇基因顯著富集在269項(xiàng)GO條目中（P＜0.05），其中在生物學(xué)過(guò)程、細(xì)胞組分和分子功能條目上各富集到175、50和44項(xiàng)條目。表2展示了顯著富集的前10個(gè)GO條目，其中最顯著富集的條目是上皮運(yùn)輸（GO:0070633，transepithelial transport）和平滑肌細(xì)胞分化調(diào)節(jié)（GO:0051150，regulation of smooth muscle cell differentiation）生物學(xué)過(guò)程；這些過(guò)程包含ABCG2、AHCYL1、KIT和MED28共4個(gè)基因。另外，NCAPG基因顯著富集在有絲分裂染色體凝聚信號(hào)通路（GO:0007076，mi-totic chromosome condensation，P=0.024）。

表2 川中黑山羊正選擇基因顯著富集的前10個(gè)GO條目Table 2 Top ten enriched GO terms for the positively selected genes in Chuanzhong Black goats

2.5 NCAPG-LCORL選擇信號(hào)區(qū)域重要變異的鑒定

由圖3b可知，山羊6號(hào)染色體NCAPG-LCORL座位所在區(qū)域在全基因組范圍內(nèi)顯示出最高的iHH12值（iHH12值=48.52），表明該區(qū)域在川中黑山羊中是一個(gè)強(qiáng)烈的正選擇區(qū)域。如圖4a所示，該區(qū)域在川中黑山羊和野山羊之間的平均加權(quán)Fst值為0.39，遠(yuǎn)高于全基因組Fst的平均值（0.16）。另外，在川中黑山羊中該區(qū)域Tajima’sD（平均值=-0.53）也遠(yuǎn)低于全基因組平均值（1.31）。在NCAPGLCORL座位內(nèi)共檢測(cè)到435個(gè)SNPs和63個(gè)indels；其中包括位于NCAPG基因第6外顯子的1個(gè)錯(cuò)義突變（c.858A＞G，p.Ile286Met）和位于LCORL基因內(nèi)的3個(gè)錯(cuò)義突變（c.4397C＞T，p.Ala1466Val；c.1433A＞G，p.Asn478Ser；c.1298A＞G，p.Tyr433Cys）。費(fèi)希爾精確檢驗(yàn)表明，c.858A＞G（P=1.36×10-5）和c.1298A＞G（P=6.08×10-5）位點(diǎn)在川中黑山羊和野山羊群體中的基因型頻率分布存在顯著性差異。但只有c.858A＞G基因型分布在川中黑山羊和建昌黑山羊群體中（參考型等位基因頻率=93.33%）存在顯著性差異（P=4.68×10-12）。而63個(gè)indels中包括兩個(gè)移碼突變位點(diǎn)（c.1615_1619delTTAAA，p.Leu539fs；c.828dupA，p.Ser277fs），它們均位于LCORL基因的第7外顯子內(nèi)。費(fèi)希爾精確檢驗(yàn)顯示，移碼突變位點(diǎn)c.1615_1619delTTAAA的基因型頻率分布在川中黑山羊（參考型等位基因頻率=13.16%）和野山羊群體（參考型等位基因頻率=100%）存在顯著性差異（P=5.89×10-11）。另外，在川中黑山羊和建昌黑山羊群體中（參考型等位基因頻率=93.33%）該位點(diǎn)的基因型分布也存在顯著性差異（P=4.3×10-12）。

圖4 川中黑山羊基因組的NCAPG-LCORL選擇信號(hào)深入分析Figure 4 In-depth analysis of the selection signal NCAPG-LCORL locus in Chuanzhong Black goats

3 討論

基因組水平的畜禽遺傳研究依賴于全基因組范圍的分子標(biāo)記，本研究利用重測(cè)序技術(shù)獲得了川中黑山羊群體的高密度SNP圖譜。與低密度的SNP芯片（例如，山羊52 k SNP芯片[33]）技術(shù)相比，利用短序列高通量測(cè)序在山羊[10]、綿羊[34]、牛[35]、豬[36]和雞[37]等物種上均可獲得千萬(wàn)級(jí)的SNP基因型數(shù)據(jù)，從而顯著提高了全基因組關(guān)聯(lián)研究等多種遺傳分析的成功率和精確性。與其他物種的注釋結(jié)果相似，在川中黑山羊群體中絕大部分SNPs位于基因間區(qū)和非編碼區(qū)，而外顯子區(qū)的變異占比極低。

有效群體大小本質(zhì)上反映了群體的遺傳多樣性豐富程度。本研究發(fā)現(xiàn)，在最近1 000世代內(nèi)川中黑山羊Ne隨著時(shí)間持續(xù)減少，這和我國(guó)其他山羊品種的變化特征一致[10,38-39]。導(dǎo)致上述變化的主要原因是，我國(guó)大部分地方品種生長(zhǎng)性能不突出、養(yǎng)殖效益低，導(dǎo)致養(yǎng)殖戶的不斷退出，群體規(guī)模逐漸減少。另一方面，由于研究樣本往往僅來(lái)自保種群或育種群，這些群體初始世代包含的種羊數(shù)量偏少。另外，川中黑山羊最近世代的Ne高于我國(guó)其他地方品種。在未來(lái)的川中黑山羊遺傳改良工作中，應(yīng)保持甚至增加家系數(shù)量，防止遺傳多樣性的喪失。

盡管相比于系譜近交系數(shù)，基因組近交系數(shù)更能反映個(gè)體真實(shí)的近交水平；究竟哪種方法更適應(yīng)于有效群體含量較小的畜禽品種尚無(wú)統(tǒng)一答案[6-7,40]。因此，本研究采用了5種常用方法對(duì)川中黑山羊近交系數(shù)進(jìn)行估計(jì)和比較?？傮w上，基于ROH（＞100 kb）估計(jì)的川中黑山羊的FROH值與瑞士山羊的近交程度接近[41]，而低于建昌黑山羊的近交水平[10]。F.Bertolini等[42]則利用SNP芯片獲得的長(zhǎng)ROH（＞1 Mb）估計(jì)了全世界117個(gè)山羊群體的近交程度，并依據(jù)FROH將近交水平劃分為低（FROH＜0.1）、中（0.1＜FROH＜0.2）和高（FROH＞0.2）3類(lèi)。據(jù)此，總體上川中黑山羊處于中等程度的近交水平；但如果排除掉川中黑山羊基因組中短ROH（＜1 Mb），則川中黑山羊的近交水平較低。與建昌黑山羊的結(jié)果類(lèi)似，川中黑山羊個(gè)體近交系數(shù)之間的差異主要?dú)w因于最近50～100世代的近交。另外，川中黑山羊FROH和FHOM之間存在較高的相關(guān)性，這與其他家畜上的結(jié)果[7,43]一致。而利用GRM矩陣所估計(jì)川中黑山羊的FUNI、FVR1和FVR2值之間相關(guān)較高，主要是因?yàn)檫@些方法的基本原理相同，都是均基于頻率校正的多位點(diǎn)平均純合度。

川中黑山羊具有體型大、生長(zhǎng)速度快和產(chǎn)羔數(shù)多等優(yōu)點(diǎn)，但關(guān)于上述性狀的遺傳基礎(chǔ)鮮有報(bào)道。本研究表明，6號(hào)染色體的NCAPG-LCORL座位是川中黑山羊中最強(qiáng)烈的選擇信號(hào)。NCAPG基因編碼非SMC凝聚素Ⅰ復(fù)合亞基G，該基因除了在有絲分裂和減數(shù)分裂中調(diào)節(jié)染色體的穩(wěn)定和壓縮外，還在腫瘤發(fā)生中扮演重要作用[44]。LCORL基因編碼配體依賴性核受體共抑制因子樣蛋白，該基因最初被認(rèn)為是一個(gè)在精子細(xì)胞中表達(dá)的轉(zhuǎn)錄因子[45]。NCAPG和LCORL在多個(gè)動(dòng)物基因組中均是彼此相鄰，故被合稱(chēng)為NCAPG-LCORL座位?；诙鄠€(gè)群體研究表明，NCAPG-LCORL座位與牛生長(zhǎng)（例如，采食量）、體型（例如，體重）和繁殖性狀均呈現(xiàn)顯著性關(guān)聯(lián)，具有一因多效性[46]。綜合馬[47]、犬[48]和豬[49]等其他物種的研究，NCAPG-LCORL被公認(rèn)為是影響動(dòng)物體型大小（體高、體長(zhǎng)等）的一個(gè)重要遺傳座位?；?個(gè)群體，我們發(fā)現(xiàn)山羊NCAPG和LCORL基因均存在大效應(yīng)突變位點(diǎn)，但其他物種上LCORL基因內(nèi)的變異位點(diǎn)與動(dòng)物體型性狀的關(guān)聯(lián)效應(yīng)更大[46]。具體哪個(gè)基因更可能是影響山羊體型的候選基因，以及上述突變位點(diǎn)的具體效應(yīng)和作用機(jī)制值得擴(kuò)大樣本量進(jìn)一步探究。