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

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

(中國農(nóng)業(yè)科學(xué)院北京畜牧獸醫(yī)研究所 動物營養(yǎng)學(xué)國家重點實驗室，北京 100193)

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

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

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

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

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

1 材料與方法

1.1 試驗材料與樣品采集

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

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

1.2 雞基因組DNA提取

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

合格DNA樣品送往北京康普森生物技術(shù)有限公司進(jìn)行二代全基因組重測序，基于Illumina NovaSeq技術(shù)測序平臺，利用雙末端測序(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/)，并去掉未比對上、低質(zhì)量(MQ<4，MQ為mapping的質(zhì)量值)、重復(fù)的Reads。剩余的高質(zhì)量Reads用于后續(xù)分析。

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

1.4 測序數(shù)據(jù)質(zhì)控

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

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

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

(1)

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

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

(2)

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

根據(jù)PLINK軟件計算的結(jié)果，利用公式計算FROH：

(3)

其中∑iLROHi為常染色體上ROH的總長度，Lauto為常染色體的長度，F(xiàn)ROH為基于ROH計算的近交系數(shù)。

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

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

(4)

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

(5)

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

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

(6)

其中，pi同上，xi同上，N為SNP的數(shù)目。

1.7 FROH、FHOM、FGRM和FUNI不同算法近交系數(shù)間的相關(guān)性分析

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

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

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

(7)

(8)

Ft=1-(1-ΔF)t

(9)

其中，Ne為有效群體大小，NS為公雞，ND為母雞，ΔF為近交增量，t為世代數(shù)，F(xiàn)t為t世代近交系數(shù)。

2 結(jié) 果

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

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

圖1 北京油雞外貌特征統(tǒng)計圖

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

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

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

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

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

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

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

圖3 ROH長度及分布統(tǒng)計圖

2.3.3 基因組近交系數(shù)比較 比較基于ROH、純合基因型、基因組關(guān)系G矩陣、聯(lián)合配子之間的相關(guān)性4種不同算法下的基因組近交系數(shù)(FROH、FHOM、FGRM、FUNI)，結(jié)果如表2所示。2019年保種群FROH平均值為(0.079 8±0.017 1)，范圍為0.009 1～0.115 5。3種基于SNP位點計算的近交系數(shù)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結(jié)果，且3種基于SNP位點計算的近交系數(shù)離散程度較大。

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

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

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

3 討 論

3.1 外貌特征變化對保種工作的指導(dǎo)

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

3.2 基因組信息對保種工作的指導(dǎo)

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

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

4 結(jié) 論

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