郭軍，曲亮，竇套存，沈曼曼，胡玉萍，王克華

（江蘇省家禽科學(xué)研究所，江蘇省家禽遺傳育種重點實驗室，江蘇 揚(yáng)州225125）

動物遺傳育種研究是一門開放的科學(xué)，對新涌現(xiàn)的方法和理論不斷地兼收并蓄、圓融貫通，尤其是計算技術(shù)[1-2]。遺傳選育是將家畜或家禽按照育種值大小排序，留選育種值高的個體的技術(shù)。育種值的準(zhǔn)確性與計算策略緊密相關(guān)，常用算法主要包括限制最大似然法（restricted maximum likelihood method, REML）和貝葉斯法。PATTERSON 和THOMPSON 將REML 與最優(yōu)線性無偏估計（best linear unbiased predictor, BLUP）方法整合，用來分析非均衡育種數(shù)據(jù)的方差組分，并成為遺傳參數(shù)估計和育種值預(yù)測的標(biāo)準(zhǔn)程序[3-4]。隨著馬爾科夫鏈蒙特卡洛（Markov chain Monte Carlo,MCMC）方法和吉布斯抽樣技術(shù)的發(fā)展，通過條件分布采樣模擬聯(lián)合分布，再通過模擬的聯(lián)合分布直接推導(dǎo)出條件分布，貝葉斯法逐漸成為REML 的備選方法[5]。近年來，積分嵌套拉普拉斯逼近方法（integrated nested Laplace approximations,INLA）被引入貝葉斯分析，與MCMC方法相比較，INLA方法不僅縮減了計算量，還提高了計算準(zhǔn)確度[6]。隨后，HOLAND等[7]和RUE 等[8]基于INLA 方法開發(fā)出遺傳評估R軟件包，促進(jìn)了INLA 方法在畜禽遺傳育種及野生動物生態(tài)研究中的應(yīng)用。

蛋雞體質(zhì)量與飼料消耗有關(guān)，從節(jié)約飼料成本角度考慮，早期選育時可淘汰體質(zhì)量過高的個體。據(jù)ANDERSON等統(tǒng)計，與隨機(jī)交配群體相比，商品蛋雞早期體質(zhì)量持續(xù)下降[9]。管理水平、疾病狀況、季節(jié)變化及自身生理狀況等均是影響蛋雞體質(zhì)量的重要因素。NORRIS等在分析文達(dá)雞4周齡體質(zhì)量遺傳變異時，將批次和性別列入固定效應(yīng)[10]。NAVARRO 等用REML 法估計了4 種肉用型雞6 周齡體質(zhì)量遺傳參數(shù)，將批次、性別及母本周齡納入固定效應(yīng)[11]。MANIATIS等以REML方法對安偉捷育種公司肉雞5 周齡體質(zhì)量進(jìn)行了遺傳評估，遺傳模型固定效應(yīng)包括性別、遺傳組、批次及父本、母本周齡[12]。總之，國外學(xué)者多以REML 算法評估雞早期體質(zhì)量遺傳數(shù)據(jù)，尚少見以貝葉斯算法進(jìn)行評估，罕見以2 種或2 種以上算法評估同一組數(shù)據(jù)?；诖耍驹囼炓缘半uF2資源群體數(shù)據(jù)為基礎(chǔ)，采用REML、MCMC及INLA分析遺傳方差組分，旨在為地方雞種遺傳評估及選育提供支持。

1 材料與方法

1.1 試驗動物

蛋雞資源群體以東鄉(xiāng)綠殼蛋雞（DXBS）、白來航雞（WL）為親本，依據(jù)F2設(shè)計，經(jīng)正反交組建。其中，F(xiàn)1代雛雞1 581 只，包括正交群體（東鄉(xiāng)綠殼蛋雞×白來航雞）552 只、反交群體（白來航雞×東鄉(xiāng)綠殼蛋雞）1 029 只。F2代雛雞3 749 只，包括母雞1 744只、公雞2 005只。第1周雛雞進(jìn)行24 h光照，熱風(fēng)爐加溫；從第2周起，依據(jù)自然光照控制人工補(bǔ)光時間。試驗雞只單獨戴翅號，出雛時經(jīng)頸部皮下注射馬立克氏病疫苗。機(jī)械化喂料、清糞。

1.2 數(shù)據(jù)測定

用MP6001 型電子天平稱量蛋雞每世代第5 周時體質(zhì)量。蛋雞資源群體5周齡體質(zhì)量原始數(shù)據(jù)共有5 405 條，其中親代數(shù)據(jù)640 條，F(xiàn)1代數(shù)據(jù)1 239條，F(xiàn)2代數(shù)據(jù)3 526 條。去除性別不明個體記錄63條（主要是F2代）、翅號重復(fù)的記錄18條和離群值記錄（3倍標(biāo)準(zhǔn)差之外）18條，剩余5 306條。系譜數(shù)據(jù)包含5 355 只雞。以SPSS 21.0 軟件分析世代、性別對體質(zhì)量的影響，確定進(jìn)入固定效應(yīng)的因素。應(yīng)用R軟件包ggplot2，以5周齡體質(zhì)量為縱坐標(biāo)、以世代為橫坐標(biāo)、以分位數(shù)為邊界繪制箱線圖[13]。

1.3 統(tǒng)計分析

資源群體遺傳組系數(shù)由R 軟件包nadiv 計算得出[14]。以加性遺傳效應(yīng)為隨機(jī)效應(yīng)，遺傳組系數(shù)、批次效應(yīng)和性別效應(yīng)為固定效應(yīng)，蛋雞5 周齡體質(zhì)量可寫作為：

式中：yi為第i只雞體質(zhì)量，y=(y1,y2,…,ym)；b0為斜率；b=(b1,b2,…,bn)，為固定效應(yīng)；ziT為指示矩陣；Zu為分子矩陣，其逆矩陣可以通過系譜直接計算得出；ui為加性遺傳效應(yīng)；εi為殘差效應(yīng)。

基于貝葉斯算法，蛋雞5周齡體質(zhì)量可表達(dá)為：

1.4 遺傳評估軟件

WOMBAT：由澳大利亞新英格蘭大學(xué)學(xué)者Karin Meyer 負(fù)責(zé)開發(fā)維護(hù)，下載網(wǎng)址http://didgeridoo.une.edu.au/km/wmbdownload1.php，是著名軟件DFREML的繼任者。該軟件應(yīng)用FORTRAN 語言編寫，專門用于估計高斯性狀遺傳參數(shù)及預(yù)測育種值[15]。WOMBAT 軟件采用REML 方法剖分方差組分，有多種算法供選擇，例如EM算法、AI算法、PX-EM算法。WOMBAT 軟件目前有windows 版本、linux 版本和MAC版本，其優(yōu)勢在于免費、運(yùn)行速度快。本研究用WOMBAT軟件進(jìn)行REML分析。

MCMCglmm：由英國愛丁堡大學(xué)Jarrod Hadfield基于R平臺開發(fā)的軟件包，于2009年供公眾免費使用[16]。從軟件名稱可以看出，該軟件應(yīng)用馬爾科夫鏈-蒙特卡洛方法擬合廣義線性模型，多用于分析高斯性狀、泊松分布性狀，也可用于分析多項分布性狀和零膨脹泊松分布性狀。該軟件支持多性狀模型、隨機(jī)回歸模型以及Hurdle 模型，應(yīng)用領(lǐng)域包括動物遺傳育種、系統(tǒng)進(jìn)化以及薈萃分析。本研究用MCMCglmm軟件包進(jìn)行貝葉斯MCMC分析。

AnimalINLA：由挪威科技大學(xué)Anna Marie Holand等基于R-INLA平臺開發(fā)的適用于動物模型的R 軟件包[7]。R-INLA 平臺是由挪威科技大學(xué)Rue、英國愛丁堡大學(xué)Lindgrey 等開發(fā)的應(yīng)用積分嵌套拉普拉斯逼近方法進(jìn)行貝葉斯計算的軟件包，用于流行病調(diào)查、動態(tài)性狀分析、基因表達(dá)雜合優(yōu)勢分析等[8]。AnimalINLA可用于分析高斯性狀，也可用于分析泊松分布性狀、二項分布性狀和零膨脹泊松分布性狀。本研究用AnimalINLA軟件包進(jìn)行貝葉斯INLA分析。

2 結(jié)果與分析

2.1 蛋雞5 周齡體質(zhì)量表型描述

白來航雞5 周齡公雞體質(zhì)量平均值為268.08 g（變異系數(shù)為10.49%），母雞體質(zhì)量平均值為236.76 g（變異系數(shù)為9.86%）；東鄉(xiāng)綠殼蛋雞公雞體質(zhì)量平均值為166.99 g（變異系數(shù)為16.38%），母雞體質(zhì)量平均值為147.26 g（變異系數(shù)為16.34%）；F1代公雞體質(zhì)量平均值為272.52 g（變異系數(shù)為11.72%），母雞體質(zhì)量平均值為241.89 g（變異系數(shù)為12.10%）；F2代公雞體質(zhì)量平均值為252.48 g（變異系數(shù)為13.56%），母雞體質(zhì)量平均值為225.08 g（變異系數(shù)為12.62%）。白來航雞經(jīng)過系統(tǒng)選育，5周齡體質(zhì)量均勻度好于同周齡的東鄉(xiāng)綠殼蛋雞。F2代蛋雞基因出現(xiàn)分離，而F1代蛋雞雖處于基因雜合狀態(tài)，但基因不表現(xiàn)分離，因而F1代蛋雞體質(zhì)量均勻度好于F2代蛋雞。各世代5 周齡體質(zhì)量分位數(shù)分布見圖1。鑒于蛋雞資源群體親本、F1代與F2代5 周齡體質(zhì)量平均值的差異有統(tǒng)計學(xué)意義（P＜0.01），遺傳模型固定效應(yīng)宜包含品種及批次效應(yīng)。t檢驗結(jié)果表明，公雞體質(zhì)量顯著高于同世代、同品種母雞（P＜0.01），因此，動物模型固定效應(yīng)宜包含性別因素。

圖1 東鄉(xiāng)綠殼蛋雞（DXBS）與白來航雞（WL）資源群體5周齡體質(zhì)量箱線圖Fig.1 Boxplot on the body mass for Dongxiang blue-shelled layers (DXBS) crossbred with White Leghorn layers(WL)at the five-week-old age

2.2 蛋雞5 周齡體質(zhì)量遺傳分析

蛋雞5 周齡體質(zhì)量加性遺傳方差、殘差及遺傳力如表1 所示。為減少計算量，對數(shù)據(jù)標(biāo)準(zhǔn)化之后再進(jìn)行方差組分剖分。3 種算法所得遺傳力相近，REML方法獲得的估計值置信度與貝葉斯方法的置信區(qū)間相符。蛋雞5周齡體質(zhì)量遺傳力高于0.2，屬于中等偏高遺傳力。為進(jìn)一步評價這3 種算法，用各軟件獲得的育種值進(jìn)行Spearman 秩相關(guān)分析。由表2可知，INLA方法預(yù)測的育種值排序與REML方法預(yù)測的結(jié)果基本一致，而MCMC方法預(yù)測的育種值排序與其余2種方法預(yù)測的結(jié)果存在差異。

2.3 計算用時分析

由于育種數(shù)據(jù)樣本量大，因而計算用時也是考量計算策略的重要指標(biāo)。3種算法分別運(yùn)行同一組數(shù)據(jù)，計算機(jī)配置是Windows 7系統(tǒng)、英特爾i7處理器（3.50 GHz）、8G 內(nèi)存。結(jié)果表明，REML 方法用時1 min，INLA方法用時8 min，MCMCglmm方法用時81 min。MCMCglmm 方法運(yùn)行時間與迭代次數(shù)有關(guān)，為了完成蛋雞5 周齡體質(zhì)量遺傳分析數(shù)據(jù)收斂，迭代次數(shù)設(shè)置為100萬次重復(fù)，每100次重復(fù)抽樣1 次，舍棄前5 萬次迭代重復(fù)。如圖2 所示：馬爾科夫鏈平均值為0.59，表明后驗分布達(dá)到穩(wěn)定狀態(tài)，也表明數(shù)據(jù)分析已完成收斂；密度圖呈現(xiàn)單一光滑單峰，表明遺傳評估結(jié)果可靠。

表1 5周齡蛋雞體質(zhì)量加性遺傳方差、殘差及遺傳力Table 1 Additive genetic variance, residual variance and heritability on body mass of the five-week-old layer

表2 REML、MCMC及INLA預(yù)測的育種值相關(guān)性分析Table 2 Correlation analysis on estimated breeding values with REML,MCMC and INLA approaches

圖2 5周齡蛋雞體質(zhì)量遺傳力軌跡圖和密度圖Fig.2 Trace and density plots of heritability on body mass of the five-week-old layer

3 討論與結(jié)論

貝葉斯法由于數(shù)學(xué)性質(zhì)優(yōu)越，在統(tǒng)計領(lǐng)域被廣泛應(yīng)用。對于多層次模型，如畜禽遺傳評估中的動物模型，可選擇貝葉斯法或REML方法。先前已有學(xué)者比較了用這2 種方法處理動物模型的優(yōu)缺點[17]：REML方法優(yōu)勢在于計算用時短，貝葉斯法優(yōu)勢在于可以計算參數(shù)的誤差并且有校準(zhǔn)點。本研究以真實的育種數(shù)據(jù)檢驗新近出現(xiàn)的INLA 方法、MCMC 方法及REML 方法，聚焦3 種算法在家禽遺傳評估中的應(yīng)用性能。結(jié)果表明，蛋雞資源群體5周齡體質(zhì)量遺傳力性狀屬于中等偏高遺傳力。INLA 方法與REML 方法獲得的遺傳參數(shù)相近，預(yù)測的育種值排序秩相關(guān)系數(shù)較高，表明這2 種方法等效。MANIATIS等在分析安偉捷公司肉雞體質(zhì)量時得到近似結(jié)果，即INLA 和REML 方法準(zhǔn)確度高于MCMC 方法[18]。然而，MATHEW 等應(yīng)用REML、MCMC、INLA方法對春大麥、水稻農(nóng)藝性狀進(jìn)行遺傳評估時，未發(fā)現(xiàn)這3 種算法在準(zhǔn)確度上存在明顯差異，只是MCMC方法計算用時較長[19]?？紤]到數(shù)據(jù)結(jié)果及樣本量可能影響分析效率，幾組試驗結(jié)果不足以證明一種算法優(yōu)于另外一種算法。

雞早期體質(zhì)量是家禽遺傳評估研究熱點。本研究表明，蛋雞5 周齡體質(zhì)量遺傳力為0.49～0.59，表明該群體有足夠多的加性遺傳變異供早期體質(zhì)量選擇。KAPELL等以動物模型對安偉捷公司4個肉雞品系5 周齡體質(zhì)量的遺傳評估結(jié)果表明，其遺傳力分別為0.36、0.39、0.32、0.40，略低于本研究結(jié)果，但同屬于中等遺傳力[20]。GAYA 等以REML 算法評估了羅斯肉雞38日齡體質(zhì)量，發(fā)現(xiàn)其遺傳力為0.40[21]，本研究結(jié)果與此相近。前文提到的MANIATIS 等以ASReml 軟件評估5 周齡雞體質(zhì)量時發(fā)現(xiàn)，其遺傳力為0.21[12]。由于本研究所用試驗素材為F2分離群體，其基因庫蘊(yùn)含豐富的遺傳變異資源，因而獲得的雞早期體質(zhì)量遺傳力高于其他團(tuán)隊研究結(jié)果。

本研究針對蛋雞資源群體5 周齡體質(zhì)量性狀，應(yīng)用REML、MCMC、INLA 3 種算法開展遺傳評估的結(jié)果表明：蛋雞早期體質(zhì)量受世代、親代品種及遺傳組系數(shù)影響，動物模型宜考慮將這些因素納入固定效應(yīng)，否則將影響遺傳力估值的準(zhǔn)確性及育種值排序；蛋雞早期體質(zhì)量屬于中等偏高遺傳力性狀，宜用個體選育獲取遺傳進(jìn)展。3種算法相比較，MCMC方法雖然能給出遺傳參數(shù)估計值置信區(qū)間，但計算用時過長，結(jié)果準(zhǔn)確性不如另外2 種算法；INLA方法與REML方法所得遺傳參數(shù)及育種值排序結(jié)果相近，但REML 方法計算用時較短。總之，隨著計算機(jī)硬件發(fā)展及算法效率的提高，阻礙貝葉斯法在畜禽育種中應(yīng)用的計算問題終將被克服，越來越多的畜禽遺傳評估工作將會使用貝葉斯法。