基于表型和基因組信息評價北京油雞保種群保種情況

王海龍，王 巧，邢思遠，王 杰，李慶賀，鄭麥青，崔煥先，劉冉冉，趙桂蘋，文 杰

(中國農業科學院北京畜牧獸醫研究所 動物營養學國家重點實驗室，北京 100193)

畜禽遺傳資源是生物多樣性的重要組成部分，在農業、經濟和科學研究中發揮著重要作用。我國幅員遼闊，不同地理和氣候條件孕育了豐富的畜禽遺傳資源，約占世界畜禽遺傳資源總量的六分之一。據2012年版《中國畜禽遺傳資源志》統計，我國有777種動物遺傳資源正式命名，包括556個地方品種，其中雞品種116個，北京油雞是其中之一。

北京油雞起源于大約300年前的清代早期，以洼里和大屯兩個地區最為集中。典型的北京油雞具有鳳冠、脛羽和五趾等外觀特征，兼具地方雞種肉、蛋品質優良，耐粗飼、抗病、抗逆性強等優良特性[1]，但生產性能相對于國外品種較差。70年代末，為滿足市場需求，企業主推高產品種，北京油雞數量驟降[2]。為保護北京油雞這一地方雞種，中國農業科學院北京畜牧獸醫研究所等單位開始承擔北京油雞的品種保護工作。2001年北京油雞被農業農村部列為國家級畜禽品種資源重點保護品種。

親緣關系較近個體交配的過程稱為近交。在畜禽育種中，一般使用近交系數評估個體之間的近交情況。近交系數是指當個體由于近交而雜合基因減少時，純合基因或純合子的百分比[3]。近交可以增加種群的純合性，也會增加純合致病基因的概率，導致后代繁殖力下降和表型衰退等情況的發生。群體數量是影響保種情況的重要因素，群體數量較小，更容易導致近交發生。

隨著二代全基因組重測序技術的發展，測序成本不斷降低，重測序技術越來越多的應用到畜禽保種工作中[4-6]?；赟NP標記計算得到基因組近交系數(genomic inbreeding coefficient based on SNPs，FSNP)，成為保種評價的新方法。使用SNP位點信息計算近交系數，可以有效避免系譜信息不完整、誤差較大等特點[7]。連續性純合片段(runs of homozygosity, ROH)[8]是基因組中的一段長純合片段，長度為數百kb到數Mb不等。自Ferencakovic等[9]首次利用ROH進行近交評估以來，ROH已廣泛應用于畜禽近交系數的計算。另有研究發現，FROH的準確性相較于FHOM、FGRM、FUNI計算的結果準確性更高[10]。

本研究以國家級北京油雞保種場隨機交配保種群體為研究對象，基于外貌特征和基因組信息對保種群體的外貌特征、近交系數和有效群體大小進行分析，研究北京油雞的保種情況。

1 材料與方法

1.1 試驗材料與樣品采集

本試驗隨機選取國家級北京油雞保種場2019年隨機交配保種群體40只個體，其中公雞17只，母雞23只，采用完全隨機交配。同時，對北京油雞鳳冠、五趾、脛羽等外貌特征進行整理，對近年外貌特征的變化趨勢進行分析。其中，1979年和1991—1993年的外貌特征記錄來自優質黃羽肉雞品系選育和配套研究論文集(內部資料)，2001年、2005年外貌特征記錄來源于現有群體的早期數據[11]。

翅下靜脈采血，于EDTA抗凝管中，-20 ℃保存，用于提取DNA。

1.2 雞基因組DNA提取

采用酚-仿傳統法提取DNA，使用兩種方法對血液DNA提取效果進行檢測：1)瓊脂糖凝膠電泳分析DNA降解的程度以及是否存在RNA污染；2)Nanodrop檢測DNA的純度，OD260 nm/OD280 nm合格范圍為1.8～2.0。

合格DNA樣品送往北京康普森生物技術有限公司進行二代全基因組重測序，基于Illumina NovaSeq技術測序平臺，利用雙末端測序(Paired-End)的方法，測序策略為Illumina PE150，測序深度為10×。

1.3 基因組變異檢測

使用基因組比對軟件BWA(v0.7.17)[12]中BWA-MEM算法，將過濾后的Clean Reads比對參考基因組(ftp://ftp.ensembl.org/pub/release-101/fasta/gallus_gallus/dna/)，并去掉未比對上、低質量(MQ<4，MQ為mapping的質量值)、重復的Reads。剩余的高質量Reads用于后續分析。

在比對到參考基因組序列的基礎上，通過突變分析軟件GATK(v4.1.7.0)[13]檢測全基因組中所有SNPs位點，過濾得到高可信度的SNP數據。

1.4 測序數據質控

本研究使用PLINK(v1.90)[14]軟件對全基因組測序數據生成的文件進行質量控制。質控標準如下：SNP檢出率大于0.9，最小等位基因頻率大于0.05，哈迪-溫伯格平衡P值大于10-6，同時只保留常染色體，避免性別影響。經過質控后，共剩余40只個體的6 252 214個SNPs用于后續分析。使用BEAGLE.18May20.d20[15]軟件對缺失基因型進行填充。

1.5 基于基因組信息估計有效群體大小(Ne)

本研究使用PLINK軟件對質控數據進行連鎖不平衡修剪。把r2<0.1默認為兩個SNPs不連鎖。經過連鎖不平衡修剪后共剩余43 167個SNPs位點，再使用SNeP(v1.1)[16]軟件估計北京油雞歷史世代的有效群體大小，公式[17]：

(1)

使用NeEstimator(v2.1)[18]軟件的基于連鎖不平衡方法估計當前世代的有效群體大小[19]。

1.6 基于基因組信息估計近交系數(F)

(2)

其中，l為組成ROH的最少SNP數目，a為設定的檢測到假陽性ROH的百分率，本研究中設定為0.05，ns為個體的SNP數目，ni為樣本數，het為SNP的平均雜合度。本研究中由此公式計算得出組成ROH的最少SNP數目為153個。

根據PLINK軟件計算的結果，利用公式計算FROH：

(3)

其中∑iLROHi為常染色體上ROH的總長度，Lauto為常染色體的長度，FROH為基于ROH計算的近交系數。

使用PLINK軟件對質控后數據計算基于純合基因型的近交系數(FHOM)和基于聯合配子之間相關性計算的近交系數(FUNI)，公式如下：

FHOM=(O-E)/(L-E)

(4)

其中，O為個體觀測純合基因型的數量，E為期望純合基因型的數量，L為基因型的SNP數目。

(5)

其中,pi為第i個SNP位點等位基因A的初始頻率，xi是基因型中a的個數，當SNP標記的基因型為AA、Aa或aa時，xi分別為0、1或2。

使用GCTA(v1.93.2beta)軟件構建基因組關系G矩陣，再使用R(v3.6.3)軟件對構建的G矩陣計算基于基因組關系G矩陣的近交系數(FGRM)[21]，公式如下：

(6)

其中，pi同上，xi同上，N為SNP的數目。

1.7 FROH、FHOM、FGRM和FUNI不同算法近交系數間的相關性分析

使用R的PerformanceAnalytics包對2019年北京油雞保種群體FROH、FHOM、FGRM和FUNI4個近交系數進行相關性分析。

1.8 基于各世代實際留種數量估計近交系數

根據國家級北京油雞保種場1979—2019年繁育記錄，統計本保種群各世代實際留種公母雞數，估計2019年北京油雞保種群有效群體大小(Ne)、近交增量(ΔF)、近交系數(Ft)，公式[22]如下：

(7)

(8)

Ft=1-(1-ΔF)t

(9)

其中，Ne為有效群體大小，NS為公雞，ND為母雞，ΔF為近交增量，t為世代數，Ft為t世代近交系數。

2 結 果

2.1 北京油雞保種群體外貌特征的變化情況

對1979年、1991—1993年、2001年、2005年、2017年、2018年和2019年的外貌特征進行統計，結果見圖1。脛羽百分比始終在97%以上，表現穩定。鳳冠、五趾的百分比在2005—2017年都有不同程度的下降，鳳冠百分比從100%下降到93%，五趾百分比從27%下降到15%。2018年鳳冠百分比上升到98%，五趾百分比上升到27%，并保持穩定。

圖1 北京油雞外貌特征統計圖

2.2 基于各世代實際留種數量估計近交系數

根據國家級北京油雞保種場1979—2019年各世代實際留種公、母雞數量(共38個世代)，估計有效群體大小和近交增量，結果如表1所示。北京油雞保種群38個世代平均近交增量為0.001 9，38個世代繁育后，2019年保種群近交系數為0.070 3。

表1 北京油雞保種群有效群體大小、近交系數估計

2.3 基于基因組信息的北京油雞保種效果分析

2.3.1 有效群體大小估計 基于連鎖不平衡方法估計2019年北京油雞保種群體的有效群體大小為193。使用SNeP估計北京油雞保種群的有效群體大小，結果如圖2所示。發現有效群體大小隨著世代的減小逐漸呈現平緩下降的趨勢，98世代之前的有效群體大小為595，13世代之前的有效群體大小176。北京油雞保種群有效群體大小從45世代前開始下降加快，北京油雞世代間隔大約為1年，70年代 剛好是北京油雞群體數量開始減少的時期，與現實情況相吻合。

圖2 北京油雞歷史世代有效群體大小變化圖

2.3.2 ROH統計 對40只北京油雞進行分析，共檢測出7 101個ROH。分別統計ROH不同長度所占比例(圖3A)、不同染色體ROH數量所占比例(圖3B)、不同染色體上ROH的平均長度(圖3C)和ROH數量與長度相關統計(圖3D)。較短的ROH(0～0.5 Mb)所占比例最多，約占總數的75%，且ROH數量隨ROH長度的增加逐漸減少。常染色體ROH數量各不相同，ROH分布不均勻。每條染色體ROH的數量大體隨染色體長度的增加而增加。其中1號染色體ROH數量最多，約占ROH總數的19%，16、25、30、31、32號染色體ROH數量最少，而不同染色體上ROH的平均長度沒有顯著區別。在不同個體中，ROH的數量與ROH覆蓋的長度也有很大的不同。隨著ROH數量的增加，ROH的總長度也在增加。在這一群體中，ROH數量最多的個體有232個ROH，總長度約為110 Mb，ROH數量最少的個體有34個ROH，總長度約為8.7 Mb。

圖3 ROH長度及分布統計圖

2.3.3 基因組近交系數比較 比較基于ROH、純合基因型、基因組關系G矩陣、聯合配子之間的相關性4種不同算法下的基因組近交系數(FROH、FHOM、FGRM、FUNI)，結果如表2所示。2019年保種群FROH平均值為(0.079 8±0.017 1)，范圍為0.009 1～0.115 5。3種基于SNP位點計算的近交系數FHOM、FGRM、FUNI平均值分別為(0.060 5±0.039 8)、(0.066 2±0.034 7)、(0.063 7±0.035 4)，范圍分別為-0.021 8～0.176 1、-0.003 2～0.162 3、0.000 6～0.164 3，分布范圍明顯大于FROH結果，且3種基于SNP位點計算的近交系數離散程度較大。

表2 基因組信息不同算法計算的近交系數

2.3.4 不同算法近交系數相關性分析 計算2019年保種群體FROH、FHOM、FGRM和FUNI4種不同算法所得近交系數值間的相關性，結果如圖4所示。從圖4中上三角可以看出，FROH與FGRM顯著相關(P<0.01)，相關系數為0.45。其他各組兩兩之間相關極顯著(P<0.001)，相關系數都大于0.5。從圖4中下三角可以看出，FHOM與FGRM，FHOM與FUNI以及FGRM與FUNI之間存在較高線性相關。

對角線：近交系數的頻率分布直方圖，圖中曲線是百分位曲線。對角線上方：不同近交系數之間的Spearman相關系數。＊.P<0.05，＊＊.P<0.01，＊＊＊.P<0.001。對角線以下：不同近交系數之間的散點圖，散點位置由上側和右側F值共同確定，圖中曲線是趨勢線

3 討 論

3.1 外貌特征變化對保種工作的指導

北京油雞保種群體脛羽、鳳冠、五趾等品種特有外貌特征得到了良好保留。2005—2017年北京油雞五趾、鳳冠百分比有所下降，自2017年起，繁育過程中適當增加具有鳳冠、五趾特征的留種比例，至2018年鳳冠、五趾特征百分比上升，恢復到2017年前的水平，并保持穩定。對1979—2019年各世代的外貌特征百分比進行統計，這一數據為國家級北京油雞保種場利用品種特有外貌特征指導保種工作提供了依據。

3.2 基因組信息對保種工作的指導

傳統計算近交系數的方式是基于系譜記錄估計的，但經常遇到系譜信息不完整和錯誤等問題，且忽略初代個體間近交系數也會影響計算準確性[23]。使用系譜信息是基于親緣關系的血緣同源概率估計期望值，未考慮減數分裂過程中基因重組是隨機事件，使用基因組信息估計的是個體間實際近交情況[4]，更能反映保種群體的真實狀況。FHOM、FGRM、FUNI是基于狀態同源估計的[4, 24]，無法區分血緣同源與狀態同源，且某些個體近交系數為負值，因此，分別使用上述幾種方法估計近交系數都是不夠準確的。FHOM、FGRM、FUNI是通過單點計算得到，結果會受等位基因頻率估算的影響，而FROH直接反映了自身長純合片段占參考基因組長度的比例，受DNA測序樣品的質量影響更小[25]，利用ROH估計近交系數可以解決其他近交系數估計存在的缺陷[26-28]。同時，在人類[29-30]、牛[31-32]、豬[33-35]研究中也發現，基于ROH估計近交系數準確性最高。Peripolli等[36]的研究中，系譜計算的近交系數與FROH相關性最強，提出FROH可以作為代替系譜評估保種群體近交系數的有效方式。

本研究中，FROH的平均值大于3種基于SNP位點計算的近交系數FHOM、FGRM、FUNI，可能是由于FHOM、FGRM、FUNI中某些個體值為負數導致的近交系數偏低，造成的結果不準確，這一結果與Alemu等[37]和Shi等[33]的研究一致。相關性分析發現，FHOM與FGRM，FHOM與FUNI以及FGRM與FUNI具有較高的線性相關，且都是基于單點計算的，與曹愉夏等[38]的研究一致。由于FROH是基于血緣同源的估計，而FHOM、FGRM、FUNI是基于狀態同源的估計，FROH與FHOM、FROH與FGRM、FROH與FUNI呈中度相關，在近交系數估計時盡量避免使用FHOM、FGRM、FUNI3種計算方法。

4 結 論

北京油雞1979—2019年40年來外貌特征明顯且百分比保持穩定，經過40余年保種工作后FROH為0.079 8，近交系數增長緩慢，說明國家級北京油雞保種場隨機交配保種群體的保種工作是切實可行的。將來對北京油雞隨機交配保種群體每年隨機選取一定數量的個體進行二代全基因組重測序檢測，有利于以后對保種狀況進行動態監控，隨時調整保種工作方案。