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

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

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

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

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

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

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

1 材料與方法

1.1 試驗(yàn)材料與樣品采集

本試驗(yàn)隨機(jī)選取國(guó)家級(jí)北京油雞保種場(chǎng)2019年隨機(jī)交配保種群體40只個(gè)體，其中公雞17只，母雞23只，采用完全隨機(jī)交配。同時(shí)，對(duì)北京油雞鳳冠、五趾、脛?dòng)鸬韧饷蔡卣鬟M(jìn)行整理，對(duì)近年外貌特征的變化趨勢(shì)進(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，使用兩種方法對(duì)血液DNA提取效果進(jìn)行檢測(cè)：1)瓊脂糖凝膠電泳分析DNA降解的程度以及是否存在RNA污染；2)Nanodrop檢測(cè)DNA的純度，OD260 nm/OD280 nm合格范圍為1.8～2.0。

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

1.3 基因組變異檢測(cè)

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

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

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

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

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

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

(1)

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

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

(2)

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

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

(3)

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

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

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

(4)

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

(5)

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

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

(6)

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

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

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

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

根據(jù)國(guó)家級(jí)北京油雞保種場(chǎng)1979—2019年繁育記錄，統(tǒng)計(jì)本保種群各世代實(shí)際留種公母雞數(shù)，估計(jì)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 北京油雞保種群體外貌特征的變化情況

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

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

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

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

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

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

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

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

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

圖3 ROH長(zhǎng)度及分布統(tǒng)計(jì)圖

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位點(diǎn)計(jì)算的近交系數(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位點(diǎn)計(jì)算的近交系數(shù)離散程度較大。

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

2.3.4 不同算法近交系數(shù)相關(guān)性分析 計(jì)算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)。

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

3 討 論

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

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

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

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

本研究中，F(xiàn)ROH的平均值大于3種基于SNP位點(diǎn)計(jì)算的近交系數(shù)FHOM、FGRM、FUNI，可能是由于FHOM、FGRM、FUNI中某些個(gè)體值為負(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)，且都是基于單點(diǎn)計(jì)算的，與曹愉夏等[38]的研究一致。由于FROH是基于血緣同源的估計(jì)，而FHOM、FGRM、FUNI是基于狀態(tài)同源的估計(jì)，F(xiàn)ROH與FHOM、FROH與FGRM、FROH與FUNI呈中度相關(guān)，在近交系數(shù)估計(jì)時(shí)盡量避免使用FHOM、FGRM、FUNI3種計(jì)算方法。

4 結(jié) 論

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