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

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

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

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

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

1 材料和方法

1.1 試驗群體和高通量重測序數(shù)據(jù)

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

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

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

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

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

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

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

1.5 選擇信號鑒定與功能富集分析

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

使用R包clusterProfiler（v 4.4.1）[32]對正選擇基因進行GO功能富集分析，顯著富集的GO條目篩選標準為P＜0.05。在分析中，山羊全基因組范圍的基因功能注釋信息來自AnnotationHub（編號為“AH101444”）。

2 結(jié)果與分析

2.1 全基因組SNPs和Indels檢測

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

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

由圖1a可知，當(dāng)SNPs之間的物理距離為10 bp時，r2平均值（0.51）最高；隨后LD迅速衰減，當(dāng)SNPs之間的距離增加到1 000 bp時，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ù)估計值

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

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

圖2 基于5種方法的川中黑山羊基因組近交系數(shù)估計值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 川中黑山羊選擇信號的檢測

如圖3a所示，基于iHH12、ROH島和XP-EHH統(tǒng)計量在川中黑山羊中分別鑒定到1 218（iHH12＞7.71）、88（ROH共享率大于50%）和1 227（XP-EHH＞3.10）個離群值窗口并分別注釋到361、164和324個基因（以Ensemble ID為標準）。其中，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個基因被3種方法均檢測到，被定義為川中黑山羊的正選擇基因。

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

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

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

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

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

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

3 討論

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

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

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

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