張夢華，羅漢鵬，魏 趁，胥 磊，王 丹 ，張曉雪，譚世新，李紅燕，黃錫霞*，王雅春*

（1.新疆農(nóng)業(yè)大學(xué)動(dòng)物科學(xué)學(xué)院，新疆烏魯木齊 830052；2.中國農(nóng)業(yè)大學(xué)動(dòng)物科學(xué)技術(shù)學(xué)院，北京 100193；3.新疆天山畜牧生物工程股份有限公司，新疆昌吉 831100）

奶牛育種的核心是在現(xiàn)有優(yōu)良品種中系統(tǒng)和持續(xù)地開展選育工作，實(shí)現(xiàn)群體的持續(xù)遺傳改良[1]。在奶牛群體持續(xù)性遺傳改良中，種公牛是影響牛群遺傳品質(zhì)最主要的因素，據(jù)統(tǒng)計(jì)，種公牛對奶牛的遺傳改良貢獻(xiàn)率可達(dá)總遺傳進(jìn)展的75% 以上，所以一個(gè)奶牛群體的生產(chǎn)性能水平多取決于與配公牛的遺傳水平[2-3]。培育和選擇優(yōu)秀種公牛，采用冷凍精液配種，是提高奶牛群體遺傳改良的根本措施[4]，種公牛選擇方法可以根據(jù)不同來源的評定信息劃分為系譜選擇法[5]、本身信息選擇法[5]、后裔測定法[5-6]和基因組選擇法[7]。后裔測定是目前評定種公牛遺傳價(jià)值最可靠的方法。

相比于其他家畜的育種技術(shù)，奶牛育種技術(shù)的研究和應(yīng)用都走在了世界前列[8]。但是，由于規(guī)?；膛ｐB(yǎng)殖成本的增加，對于奶牛生產(chǎn)性能的遺傳基礎(chǔ)要求變得越來越高，這給奶牛育種提出了更高要求[9]。而我國大部分奶牛養(yǎng)殖場管理人員對于奶牛育種方法及技術(shù)的認(rèn)識和應(yīng)用依然非常落后，尤其是對奶牛后裔測定工作缺乏深刻認(rèn)識，致使積極性不高，嚴(yán)重影響了奶牛育種工作的開展和相關(guān)技術(shù)的發(fā)展。在實(shí)際育種工作中，核心是育種公司的選種和牛場的選種及選配[10]。而種公牛的后裔測定工作則是選擇優(yōu)秀種公牛的主要環(huán)節(jié)。與配種公牛的好壞直接關(guān)系到新疆地區(qū)中國荷斯坦牛品種改良的遺傳進(jìn)展和良種化程度的高低，因此對種公牛進(jìn)行后裔測定是非常必要和緊迫的工作，對加快新疆地區(qū)中國荷斯坦牛群體遺傳改良進(jìn)展具有積極作用。

國內(nèi)外的奶牛群體遺傳改良?xì)v史證明，具有系譜、生產(chǎn)性能等多種信息資料的牛群，其群體遺傳水平提高的速度遠(yuǎn)高于無信息牛群[11-12]?；谛陆貐^(qū)生產(chǎn)性能測定積累的大量生產(chǎn)性能測定（Dairy Herd Improvement，DHI）記錄和較為完整的系譜數(shù)據(jù)（生產(chǎn)性能數(shù)據(jù)屬于有重復(fù)記錄的測定類數(shù)據(jù)，所以這類性狀的遺傳評估多采用測定日模型[13-15]），本研究先通過非遺傳因素分析確定遺傳評估模型中的固定效應(yīng)，再通過泌乳曲線擬合確定測定日模型中的子模型（泌乳曲線模型），最后運(yùn)用測定日模型進(jìn)行種公牛遺傳參數(shù)估計(jì)。

1 材料與方法

1.1 數(shù)據(jù)來源 數(shù)據(jù)來自于新疆地區(qū)2010—2015 年使用過新疆天山畜牧生物工程有限公司種公牛站凍精的千余家牛場。選取出生于2001—2010 年的1 557 頭中國荷斯坦種公牛的54 085 頭女兒，共計(jì)922 753 條DHI測定日記錄數(shù)據(jù)，建立數(shù)據(jù)集。

1.2 數(shù)據(jù)處理

1.2.1 系譜整理 數(shù)據(jù)集的系譜包括牛號、父號、母號、祖父號、祖母號、外祖父號、外祖母號。遺傳評估需要根據(jù)系譜所提供的親緣關(guān)系構(gòu)建A 陣（分子親緣關(guān)系矩陣），基于新疆地區(qū)DHI 記錄中提供的基本系譜信息，有表型記錄的個(gè)體要至少向上追溯三代，剔除系譜記錄不完整個(gè)體的全部數(shù)據(jù)。初步整理后，用于方差組分估計(jì)群體及育種值估計(jì)群體各胎次記錄的系譜情況見表1。

表1 數(shù)據(jù)集系譜統(tǒng)計(jì)

1.2.2 數(shù)據(jù)整理 根據(jù)數(shù)據(jù)結(jié)構(gòu)及數(shù)據(jù)的完整性進(jìn)行以下篩選：1）產(chǎn)奶量1.0～50 kg、乳脂率1.50%～9.00%、乳蛋白率1.50%～6.00%[16]；2）刪除父號及產(chǎn)犢月齡缺失的記錄；3）刪除泌乳天數(shù)小于5 d 和大于305 d 的記錄；4）刪除在泌乳天數(shù)90 d 之前無記錄的個(gè)體；5）刪除連續(xù)2 次測定日間隔超過70 d 的個(gè)體；6）刪除測定次數(shù)小于6 次的個(gè)體；7）刪除后代女兒頭數(shù)低于20頭的公牛后代[17]。

數(shù)據(jù)集效應(yīng)劃分：測定場按照參測牛場劃分為1 032個(gè)水平，將測定日期進(jìn)行重新編號后與測定場合并為場-測定日效應(yīng)；產(chǎn)犢年份效應(yīng)根據(jù)產(chǎn)犢年份自2010—2015 年劃分為6 個(gè)水平，產(chǎn)犢季節(jié)效應(yīng)根據(jù)新疆特殊的氣候條件，以3、4、5 月為春季，6、7、8 月為夏季，9、10、11 月為秋季，12、1、2 為冬季，劃分為4 個(gè)水平[18]，將產(chǎn)犢年份和產(chǎn)犢季節(jié)合并為產(chǎn)犢年季效應(yīng)；泌乳天數(shù)效應(yīng)根據(jù)計(jì)算的泌乳天數(shù)篩選5～305 d 的記錄，每30 d一個(gè)水平，劃分為10 個(gè)水平；并按照1 胎、2 胎、≥3胎次分為3 個(gè)子集分別進(jìn)行計(jì)算。

1.3 統(tǒng)計(jì)分析

1.3.1 混合動(dòng)物模型 運(yùn)用SAS8.1 軟件的GLM 過程，使用混合動(dòng)物模型進(jìn)行最小二乘分析[19]，分析胎次、測定場、產(chǎn)犢年份、產(chǎn)犢季節(jié)和泌乳天數(shù)效應(yīng)對日產(chǎn)奶量和DHI 測定指標(biāo)的影響。模型為：

式中，Yjklmno為性狀觀察值；u為群體均值；Pj為第j個(gè)胎次效應(yīng)；Hk為第k個(gè)測定場效應(yīng)；Yl為第l個(gè)產(chǎn)犢年份效應(yīng)；Sm為第m個(gè)產(chǎn)犢季節(jié)效應(yīng)；Dn為第n個(gè)泌乳天數(shù)效應(yīng)；ejklmno為隨機(jī)殘差。

1.3.2 泌乳曲線擬合 運(yùn)用SAS 8.1 統(tǒng)計(jì)分析軟件對日產(chǎn)奶量進(jìn)行描述性統(tǒng)計(jì)分析，使用NILN（非線性回歸模型）計(jì)算出模型參數(shù)。泌乳曲線的擬合采用Legendre多項(xiàng)式回歸模型[20]：

上述模型中，a、b、c、d、e、f、g均為模型的參數(shù)，t為泌乳天數(shù)，yt表示在第t天的產(chǎn)奶量。

1.3.3 隨機(jī)回歸模型 本研究中隨機(jī)回歸測定日模型的方差組分和育種值估計(jì)均采用DMU（Derivative Free Multivariate）軟件[21]。采用的隨機(jī)回歸模型包括固定效應(yīng)、固定回歸、隨機(jī)回歸效應(yīng)。固定效應(yīng)包括：場-測定日和產(chǎn)犢年季效應(yīng)；以固定回歸擬合產(chǎn)犢月齡效應(yīng)，固定回歸采用4 階Legendre 多項(xiàng)式[22-23]；隨機(jī)效應(yīng)包括：加性遺傳和永久環(huán)境效應(yīng)。

隨機(jī)回歸模型矩陣形式如下：y=Hc+X+Z+W+e

其中，c 和β分別表示場-測定日、產(chǎn)犢年季和產(chǎn)犢月齡固定效應(yīng)（回歸系數(shù)）向量；a 表示加性隨機(jī)回歸系數(shù)向量；p 表示永久環(huán)境效應(yīng)隨機(jī)回歸系數(shù)向量；H、X、Z 和W 是對應(yīng)效應(yīng)的關(guān)聯(lián)矩陣；A、P 分別是加性效應(yīng)和永久環(huán)境效應(yīng)回歸系數(shù)（協(xié)）方差矩陣，假設(shè)殘差同質(zhì)。

1.4 計(jì)算5～305 d 總遺傳力和估計(jì)育種值 使用隨機(jī)回歸測定日模型Legendre 多項(xiàng)式子模型進(jìn)行遺傳參數(shù)和方差組分的估計(jì)則需要經(jīng)過進(jìn)一步的運(yùn)算才能得到遺傳力和估計(jì)育種值，求解過程如下[24]：

1）求相應(yīng)階數(shù)z 矩陣（不同階數(shù)對應(yīng)固定的z 矩陣）；

2）求隨機(jī)回歸系數(shù)間的VCV 矩陣（DMU 結(jié)果文件給出），包括加性回歸系數(shù)的方差協(xié)方差矩陣∑a，永久環(huán)境效應(yīng)回歸系數(shù)的方差協(xié)方差矩陣∑p 和殘差e；

4）個(gè)體育種值的計(jì)算，EBV=Z'C a，a 為DMU 估計(jì)得到的個(gè)體的各階次Legendre 多項(xiàng)式系數(shù)，即每個(gè)個(gè)體的5～305 d 加性效應(yīng)Legendre 正交多項(xiàng)式函數(shù)。將個(gè)體一個(gè)泌乳期內(nèi)5～305 d 的估計(jì)育種值進(jìn)行累加，即可獲得個(gè)體某性狀在整個(gè)泌乳期內(nèi)的估計(jì)育種值。對個(gè)體泌乳期估計(jì)育種值取平均值，即可獲得該個(gè)體相應(yīng)泌乳期內(nèi)的平均育種值。

1.5 育種值排隊(duì) 利用DMU 軟件的DMU4 模塊估計(jì)大群育種值。依據(jù)本研究中涉及的性狀，計(jì)算采用現(xiàn)行的中國奶牛性能指數(shù)2（CPI2，China Performance Index 2）[25]對中國荷斯坦公牛的綜合育種值進(jìn)行排隊(duì)。

其中，Milk 表示產(chǎn)奶量；Fatpct 表示乳脂率；Propct 表示乳蛋白率；SCS 表示體細(xì)胞評分；δmilk表示產(chǎn)奶量的標(biāo)準(zhǔn)差；δfatpct表示乳脂率的標(biāo)準(zhǔn)差；δpropct表示乳蛋白率的標(biāo)準(zhǔn)差；δSCS表示體細(xì)胞評分的標(biāo)準(zhǔn)。

2 結(jié)果與分析

2.1 描述性統(tǒng)計(jì) 通過數(shù)據(jù)篩選，最后用于計(jì)算的樣本量為77 6759 條測定日記錄。如表2 所示，新疆地區(qū)中國荷斯坦牛平均日產(chǎn)奶量為25.74 kg，乳脂量為0.897 2 kg，乳蛋白量為0.823 7 kg，體細(xì)胞評分為3.75。

表2 日產(chǎn)奶量和DHI 測定指標(biāo)的描述性統(tǒng)計(jì)

2.2 最小二乘分析結(jié)果 由表3 可知，胎次、測定場、產(chǎn)犢年份、產(chǎn)犢季節(jié)和泌乳天數(shù)均對產(chǎn)奶量、乳脂量、乳蛋白量和體細(xì)胞評分有極顯著影響，所以，遺傳參數(shù)估計(jì)模型中需要考慮以上效應(yīng)。

表3 日產(chǎn)奶量、乳脂量、乳蛋白量和體細(xì)胞評分的方差分析表

2.3 Legendre 多項(xiàng)式回歸模型擬合泌乳曲線的效果 利用Legendre 多項(xiàng)式回歸模型擬合新疆地區(qū)中國荷斯坦牛泌乳曲線，通過計(jì)算，一胎、二胎、三胎和所有胎次的泌乳曲線擬合度分別為0.975 8、0.980 4、0.969 1、0.940 1（表4）。圖1 為利用Legendre 多項(xiàng)式回歸模型擬合中國荷斯坦牛不同胎次的泌乳曲線，從圖中可以看出，中國荷斯坦牛從開始泌乳，產(chǎn)奶量迅速增加，在泌乳的第40～60 d 出現(xiàn)產(chǎn)奶高峰，一直維持到第100 天左右，之后開始緩慢下降，這一泌乳特征與該品種的實(shí)際泌乳規(guī)律一致。所以，利用Legendre 多項(xiàng)式回歸模型作為遺傳參數(shù)估計(jì)模型的子模型是可行的。

圖1 Legendre 多項(xiàng)式回歸模型擬合中國荷斯坦牛不同胎次泌乳曲線

表4 Legendre 多項(xiàng)式回歸模型各胎次泌乳曲線擬合參數(shù)及擬合度

2.4 最優(yōu)隨機(jī)回歸模型遺傳參數(shù)結(jié)果 由表5 可以看出，新疆地區(qū)中國荷斯坦牛產(chǎn)奶量5～305 d 總遺傳力在0.18～0.22 之間，乳脂量的5～305 d 總遺傳力在0.12～0.15之間，乳蛋白量的5～305 d 總遺傳力在0.14～0.19 之間，體細(xì)胞評分的5～305 d 總遺傳力在0.06～0.10 之間。

表5 不同性狀各胎次5～305 d 方差組分及遺傳參數(shù)

2.5 公牛育種值排隊(duì) 表6 所示為新疆地區(qū)排名前10位的后裔測定種公牛，冠軍公牛是65109051，CPI2 為4 439.86，國別為中國。其父親HOUSAM17058140 的GLPI（Genomic Lifetime Performance Index）為1 807（R2為96%），母親為60635835，為來自美國的優(yōu)秀荷斯坦牛。中國奶牛數(shù)據(jù)中心對該公牛的官方描述為：進(jìn)口美國優(yōu)秀胚胎移植，父親對其女兒的體軀結(jié)構(gòu)、泌乳系統(tǒng)、肢蹄等均有明顯的改良效果，其女兒的平均日產(chǎn)奶量為29 kg，平均乳脂率為3.87%，平均乳蛋白率為3.23%，平均體細(xì)胞數(shù)為18.2 萬個(gè)/mL。

表6 公牛后裔測定排名前10 位的種公牛

3 討 論

本研究所使用的測定日模型包括3 個(gè)部分：固定效應(yīng)、固定回歸效應(yīng)和隨機(jī)回歸效應(yīng)。首先，需要利用混合動(dòng)物模型進(jìn)行最小二乘分析，從而確定模型中需要剖分出哪些顯著或者極顯著的固定效應(yīng)，或者通過劃分?jǐn)?shù)據(jù)集、分別計(jì)算的方式避免或降低固定效應(yīng)的影響，以提高模型的準(zhǔn)確性。本研究中胎次、測定場、產(chǎn)犢年份、產(chǎn)犢季節(jié)和泌乳天數(shù)均對產(chǎn)奶量、乳脂量、乳蛋白量和體細(xì)胞評分有極顯著影響，所以在模型中考慮場-測定日和產(chǎn)犢年季效應(yīng)作為固定效應(yīng)的同時(shí)，按照不同的胎次劃分?jǐn)?shù)據(jù)集。此外，在奶牛群體初次進(jìn)行遺傳分析時(shí)，常存在系譜信息記錄不完整的問題，進(jìn)一步完善系譜深度可以挖掘群體內(nèi)的遺傳聯(lián)系，有望提高遺傳力的估計(jì)值，并增加育種值估計(jì)的可靠性[24]。

其次，測定日模型中的子模型有很多（Legendre多項(xiàng)式回歸模型[14]、Wood 模型[26]、Wilmink 模型[27]、Nelder 逆多項(xiàng)式模型[28]、Ali-Schaeffer 模型[28]等），其作用是利用子模型來計(jì)算群體泌乳曲線，那么就需要根據(jù)當(dāng)前數(shù)據(jù)集進(jìn)行群體泌乳曲線擬合后選定模型擬合度高于0.95 的作為測定日模型的子模型較為適宜。本研究選用了Legendre 多項(xiàng)式回歸模型作為子模型，不僅是基于前期的研究基礎(chǔ)[29]，還進(jìn)行了泌乳曲線擬合，各胎次泌乳曲線擬合度均高于0.95，且擬合效果較好。第2、3 胎次泌乳曲線峰值最高，同時(shí)在整個(gè)泌乳期的產(chǎn)奶量也最高；第1 胎次的擬合曲線泌乳高峰值最低，其在整個(gè)泌乳期的產(chǎn)奶量下降速度最慢，泌乳末期的日產(chǎn)奶量均高于其他胎次；第1 胎次與所有胎次的泌乳曲線很接近，其泌乳期產(chǎn)奶量達(dá)到泌乳高峰的速率，與從泌乳高峰開始下降的速率也基本一致。2 胎、3 胎泌乳性能優(yōu)于1 胎，但泌乳持續(xù)力與1 胎相比較差。同時(shí)，由于2 胎及3 胎在相同泌乳階段群體均值十分接近，使得泌乳曲線基本重合。新疆地區(qū)中國荷斯坦牛泌乳高峰出現(xiàn)在產(chǎn)犢后45～75 d，且整體泌乳持續(xù)力較好，持續(xù)力是與奶牛健康性狀及繁殖性狀相關(guān)的中等遺傳力性狀，對其進(jìn)行選擇可以有效改善群體的泌乳曲線。

最后Legendre 多項(xiàng)式回歸模型同時(shí)要作為隨機(jī)回歸效應(yīng)，用來計(jì)算每個(gè)個(gè)體泌乳5～305 d 每天的加性遺傳方差和永久環(huán)境方差，最后合成305 d 的遺傳參數(shù)[30]。新疆地區(qū)中國荷斯坦牛5～305 d 產(chǎn)奶量、乳脂量、乳蛋白量和體細(xì)胞評分的總遺傳力分別為0.18～0.22、0.12～0.15、0.14～0.19、0.06～0.10，與北京地區(qū)[24]、寧夏地區(qū)[14,31]、河南省[32]和山東省[33]利用測定日模型估計(jì)的中國荷斯坦牛產(chǎn)奶量、乳脂量、乳蛋白量和體細(xì)胞評分遺傳力的結(jié)果一致。各性狀遺傳力估計(jì)值均略低于Miglior[34]基于四階Legendre 多項(xiàng)式測定日模型對中國荷斯坦牛頭胎DHI 記錄的分析結(jié)果，產(chǎn)奶量、乳脂量、乳蛋白量和體細(xì)胞評分4 個(gè)性狀的估計(jì)遺傳力分別為0.291、0.222、0.251 和0.092。與王國龍等[10]估計(jì)的烏魯木齊種牛場荷斯坦種公牛產(chǎn)奶量遺傳力0.14、乳脂率遺傳力0.13、乳蛋白率遺傳力0.14 的結(jié)果一致。曹福存等[35]人估計(jì)了北京地區(qū)81 385 頭荷斯坦母牛的LSCC（體細(xì)胞數(shù)自然對數(shù)）遺傳力為0.12，個(gè)體測定日LSCC 的遺傳力為0.06～0.11。

種公牛的遺傳品質(zhì)直接關(guān)系到牛群遺傳改良的效果。我國種公牛站已經(jīng)處于轉(zhuǎn)型發(fā)展的關(guān)鍵階段，面對市場需求和國際化的發(fā)展趨勢，聯(lián)合育種正在創(chuàng)新發(fā)展，對于種公牛聯(lián)合后裔測定機(jī)制已經(jīng)逐步形成，為了不斷加大種公牛后裔測定力度，由種公牛站、核心育種場和育種企業(yè)聯(lián)合組成了奶牛后裔測定聯(lián)盟（北方聯(lián)盟、香山聯(lián)盟），聯(lián)合后裔測定[36]。在各省區(qū)遺傳評估的基礎(chǔ)上，積極促進(jìn)數(shù)據(jù)統(tǒng)一化，為更好的進(jìn)行聯(lián)合后裔測定夯實(shí)數(shù)據(jù)基礎(chǔ)。此外，積極進(jìn)行種公牛的遺傳評定有助于擴(kuò)大優(yōu)秀種公牛的遺傳優(yōu)勢，促進(jìn)品種的選育提高，加快群體遺傳進(jìn)展。

4 結(jié) 論

新疆地區(qū)中國荷斯坦牛305 d 產(chǎn)奶量、乳脂量、乳蛋白量和體細(xì)胞評分的遺傳力在0.07～0.22 之間，利用隨機(jī)回歸測定日模型進(jìn)行種公牛遺傳評定時(shí)，需要篩選出顯著影響性狀的效應(yīng)和擬合度較高的泌乳曲線模型才能獲得遺傳評估的最優(yōu)模型。種公牛的后裔測定工作是評定種公牛遺傳品質(zhì)的重要環(huán)節(jié)，對于牛群的遺傳改良效果影響重大，也是新疆地區(qū)中國荷斯坦牛種公牛凍精銷售的重要支撐。