張夢華,羅漢鵬,魏 趁,胥 磊,王 丹 ,張曉雪,譚世新,李紅燕,黃錫霞*,王雅春*
(1.新疆農業大學動物科學學院,新疆烏魯木齊 830052;2.中國農業大學動物科學技術學院,北京 100193;3.新疆天山畜牧生物工程股份有限公司,新疆昌吉 831100)
奶牛育種的核心是在現有優良品種中系統和持續地開展選育工作,實現群體的持續遺傳改良[1]。在奶牛群體持續性遺傳改良中,種公牛是影響牛群遺傳品質最主要的因素,據統計,種公牛對奶牛的遺傳改良貢獻率可達總遺傳進展的75% 以上,所以一個奶牛群體的生產性能水平多取決于與配公牛的遺傳水平[2-3]。培育和選擇優秀種公牛,采用冷凍精液配種,是提高奶牛群體遺傳改良的根本措施[4],種公牛選擇方法可以根據不同來源的評定信息劃分為系譜選擇法[5]、本身信息選擇法[5]、后裔測定法[5-6]和基因組選擇法[7]。后裔測定是目前評定種公牛遺傳價值最可靠的方法。
相比于其他家畜的育種技術,奶牛育種技術的研究和應用都走在了世界前列[8]。但是,由于規模化奶牛養殖成本的增加,對于奶牛生產性能的遺傳基礎要求變得越來越高,這給奶牛育種提出了更高要求[9]。而我國大部分奶牛養殖場管理人員對于奶牛育種方法及技術的認識和應用依然非常落后,尤其是對奶牛后裔測定工作缺乏深刻認識,致使積極性不高,嚴重影響了奶牛育種工作的開展和相關技術的發展。在實際育種工作中,核心是育種公司的選種和牛場的選種及選配[10]。而種公牛的后裔測定工作則是選擇優秀種公牛的主要環節。與配種公牛的好壞直接關系到新疆地區中國荷斯坦牛品種改良的遺傳進展和良種化程度的高低,因此對種公牛進行后裔測定是非常必要和緊迫的工作,對加快新疆地區中國荷斯坦牛群體遺傳改良進展具有積極作用。
國內外的奶牛群體遺傳改良歷史證明,具有系譜、生產性能等多種信息資料的牛群,其群體遺傳水平提高的速度遠高于無信息牛群[11-12]。基于新疆地區生產性能測定積累的大量生產性能測定(Dairy Herd Improvement,DHI)記錄和較為完整的系譜數據(生產性能數據屬于有重復記錄的測定類數據,所以這類性狀的遺傳評估多采用測定日模型[13-15]),本研究先通過非遺傳因素分析確定遺傳評估模型中的固定效應,再通過泌乳曲線擬合確定測定日模型中的子模型(泌乳曲線模型),最后運用測定日模型進行種公牛遺傳參數估計。
1.1 數據來源 數據來自于新疆地區2010—2015 年使用過新疆天山畜牧生物工程有限公司種公牛站凍精的千余家牛場。選取出生于2001—2010 年的1 557 頭中國荷斯坦種公牛的54 085 頭女兒,共計922 753 條DHI測定日記錄數據,建立數據集。
1.2 數據處理
1.2.1 系譜整理 數據集的系譜包括牛號、父號、母號、祖父號、祖母號、外祖父號、外祖母號。遺傳評估需要根據系譜所提供的親緣關系構建A 陣(分子親緣關系矩陣),基于新疆地區DHI 記錄中提供的基本系譜信息,有表型記錄的個體要至少向上追溯三代,剔除系譜記錄不完整個體的全部數據。初步整理后,用于方差組分估計群體及育種值估計群體各胎次記錄的系譜情況見表1。

表1 數據集系譜統計
1.2.2 數據整理 根據數據結構及數據的完整性進行以下篩選:1)產奶量1.0~50 kg、乳脂率1.50%~9.00%、乳蛋白率1.50%~6.00%[16];2)刪除父號及產犢月齡缺失的記錄;3)刪除泌乳天數小于5 d 和大于305 d 的記錄;4)刪除在泌乳天數90 d 之前無記錄的個體;5)刪除連續2 次測定日間隔超過70 d 的個體;6)刪除測定次數小于6 次的個體;7)刪除后代女兒頭數低于20頭的公牛后代[17]。
數據集效應劃分:測定場按照參測牛場劃分為1 032個水平,將測定日期進行重新編號后與測定場合并為場-測定日效應;產犢年份效應根據產犢年份自2010—2015 年劃分為6 個水平,產犢季節效應根據新疆特殊的氣候條件,以3、4、5 月為春季,6、7、8 月為夏季,9、10、11 月為秋季,12、1、2 為冬季,劃分為4 個水平[18],將產犢年份和產犢季節合并為產犢年季效應;泌乳天數效應根據計算的泌乳天數篩選5~305 d 的記錄,每30 d一個水平,劃分為10 個水平;并按照1 胎、2 胎、≥3胎次分為3 個子集分別進行計算。
1.3 統計分析
1.3.1 混合動物模型 運用SAS8.1 軟件的GLM 過程,使用混合動物模型進行最小二乘分析[19],分析胎次、測定場、產犢年份、產犢季節和泌乳天數效應對日產奶量和DHI 測定指標的影響。模型為:

式中,Yjklmno為性狀觀察值;u為群體均值;Pj為第j個胎次效應;Hk為第k個測定場效應;Yl為第l個產犢年份效應;Sm為第m個產犢季節效應;Dn為第n個泌乳天數效應;ejklmno為隨機殘差。
1.3.2 泌乳曲線擬合 運用SAS 8.1 統計分析軟件對日產奶量進行描述性統計分析,使用NILN(非線性回歸模型)計算出模型參數。泌乳曲線的擬合采用Legendre多項式回歸模型[20]:

上述模型中,a、b、c、d、e、f、g均為模型的參數,t為泌乳天數,yt表示在第t天的產奶量。
1.3.3 隨機回歸模型 本研究中隨機回歸測定日模型的方差組分和育種值估計均采用DMU(Derivative Free Multivariate)軟件[21]。采用的隨機回歸模型包括固定效應、固定回歸、隨機回歸效應。固定效應包括:場-測定日和產犢年季效應;以固定回歸擬合產犢月齡效應,固定回歸采用4 階Legendre 多項式[22-23];隨機效應包括:加性遺傳和永久環境效應。
隨機回歸模型矩陣形式如下:y=Hc+X+Z+W+e

其中,c 和β分別表示場-測定日、產犢年季和產犢月齡固定效應(回歸系數)向量;a 表示加性隨機回歸系數向量;p 表示永久環境效應隨機回歸系數向量;H、X、Z 和W 是對應效應的關聯矩陣;A、P 分別是加性效應和永久環境效應回歸系數(協)方差矩陣,假設殘差同質。
1.4 計算5~305 d 總遺傳力和估計育種值 使用隨機回歸測定日模型Legendre 多項式子模型進行遺傳參數和方差組分的估計則需要經過進一步的運算才能得到遺傳力和估計育種值,求解過程如下[24]:
1)求相應階數z 矩陣(不同階數對應固定的z 矩陣);
2)求隨機回歸系數間的VCV 矩陣(DMU 結果文件給出),包括加性回歸系數的方差協方差矩陣∑a,永久環境效應回歸系數的方差協方差矩陣∑p 和殘差e;
4)個體育種值的計算,EBV=Z'C a,a 為DMU 估計得到的個體的各階次Legendre 多項式系數,即每個個體的5~305 d 加性效應Legendre 正交多項式函數。將個體一個泌乳期內5~305 d 的估計育種值進行累加,即可獲得個體某性狀在整個泌乳期內的估計育種值。對個體泌乳期估計育種值取平均值,即可獲得該個體相應泌乳期內的平均育種值。
1.5 育種值排隊 利用DMU 軟件的DMU4 模塊估計大群育種值。依據本研究中涉及的性狀,計算采用現行的中國奶牛性能指數2(CPI2,China Performance Index 2)[25]對中國荷斯坦公牛的綜合育種值進行排隊。

其中,Milk 表示產奶量;Fatpct 表示乳脂率;Propct 表示乳蛋白率;SCS 表示體細胞評分;δmilk表示產奶量的標準差;δfatpct表示乳脂率的標準差;δpropct表示乳蛋白率的標準差;δSCS表示體細胞評分的標準。
2.1 描述性統計 通過數據篩選,最后用于計算的樣本量為77 6759 條測定日記錄。如表2 所示,新疆地區中國荷斯坦牛平均日產奶量為25.74 kg,乳脂量為0.897 2 kg,乳蛋白量為0.823 7 kg,體細胞評分為3.75。

表2 日產奶量和DHI 測定指標的描述性統計
2.2 最小二乘分析結果 由表3 可知,胎次、測定場、產犢年份、產犢季節和泌乳天數均對產奶量、乳脂量、乳蛋白量和體細胞評分有極顯著影響,所以,遺傳參數估計模型中需要考慮以上效應。

表3 日產奶量、乳脂量、乳蛋白量和體細胞評分的方差分析表
2.3 Legendre 多項式回歸模型擬合泌乳曲線的效果 利用Legendre 多項式回歸模型擬合新疆地區中國荷斯坦牛泌乳曲線,通過計算,一胎、二胎、三胎和所有胎次的泌乳曲線擬合度分別為0.975 8、0.980 4、0.969 1、0.940 1(表4)。圖1 為利用Legendre 多項式回歸模型擬合中國荷斯坦牛不同胎次的泌乳曲線,從圖中可以看出,中國荷斯坦牛從開始泌乳,產奶量迅速增加,在泌乳的第40~60 d 出現產奶高峰,一直維持到第100 天左右,之后開始緩慢下降,這一泌乳特征與該品種的實際泌乳規律一致。所以,利用Legendre 多項式回歸模型作為遺傳參數估計模型的子模型是可行的。

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

表4 Legendre 多項式回歸模型各胎次泌乳曲線擬合參數及擬合度
2.4 最優隨機回歸模型遺傳參數結果 由表5 可以看出,新疆地區中國荷斯坦牛產奶量5~305 d 總遺傳力在0.18~0.22 之間,乳脂量的5~305 d 總遺傳力在0.12~0.15之間,乳蛋白量的5~305 d 總遺傳力在0.14~0.19 之間,體細胞評分的5~305 d 總遺傳力在0.06~0.10 之間。

表5 不同性狀各胎次5~305 d 方差組分及遺傳參數
2.5 公牛育種值排隊 表6 所示為新疆地區排名前10位的后裔測定種公牛,冠軍公牛是65109051,CPI2 為4 439.86,國別為中國。其父親HOUSAM17058140 的GLPI(Genomic Lifetime Performance Index)為1 807(R2為96%),母親為60635835,為來自美國的優秀荷斯坦牛。中國奶牛數據中心對該公牛的官方描述為:進口美國優秀胚胎移植,父親對其女兒的體軀結構、泌乳系統、肢蹄等均有明顯的改良效果,其女兒的平均日產奶量為29 kg,平均乳脂率為3.87%,平均乳蛋白率為3.23%,平均體細胞數為18.2 萬個/mL。

表6 公牛后裔測定排名前10 位的種公牛
本研究所使用的測定日模型包括3 個部分:固定效應、固定回歸效應和隨機回歸效應。首先,需要利用混合動物模型進行最小二乘分析,從而確定模型中需要剖分出哪些顯著或者極顯著的固定效應,或者通過劃分數據集、分別計算的方式避免或降低固定效應的影響,以提高模型的準確性。本研究中胎次、測定場、產犢年份、產犢季節和泌乳天數均對產奶量、乳脂量、乳蛋白量和體細胞評分有極顯著影響,所以在模型中考慮場-測定日和產犢年季效應作為固定效應的同時,按照不同的胎次劃分數據集。此外,在奶牛群體初次進行遺傳分析時,常存在系譜信息記錄不完整的問題,進一步完善系譜深度可以挖掘群體內的遺傳聯系,有望提高遺傳力的估計值,并增加育種值估計的可靠性[24]。
其次,測定日模型中的子模型有很多(Legendre多項式回歸模型[14]、Wood 模型[26]、Wilmink 模型[27]、Nelder 逆多項式模型[28]、Ali-Schaeffer 模型[28]等),其作用是利用子模型來計算群體泌乳曲線,那么就需要根據當前數據集進行群體泌乳曲線擬合后選定模型擬合度高于0.95 的作為測定日模型的子模型較為適宜。本研究選用了Legendre 多項式回歸模型作為子模型,不僅是基于前期的研究基礎[29],還進行了泌乳曲線擬合,各胎次泌乳曲線擬合度均高于0.95,且擬合效果較好。第2、3 胎次泌乳曲線峰值最高,同時在整個泌乳期的產奶量也最高;第1 胎次的擬合曲線泌乳高峰值最低,其在整個泌乳期的產奶量下降速度最慢,泌乳末期的日產奶量均高于其他胎次;第1 胎次與所有胎次的泌乳曲線很接近,其泌乳期產奶量達到泌乳高峰的速率,與從泌乳高峰開始下降的速率也基本一致。2 胎、3 胎泌乳性能優于1 胎,但泌乳持續力與1 胎相比較差。同時,由于2 胎及3 胎在相同泌乳階段群體均值十分接近,使得泌乳曲線基本重合。新疆地區中國荷斯坦牛泌乳高峰出現在產犢后45~75 d,且整體泌乳持續力較好,持續力是與奶牛健康性狀及繁殖性狀相關的中等遺傳力性狀,對其進行選擇可以有效改善群體的泌乳曲線。
最后Legendre 多項式回歸模型同時要作為隨機回歸效應,用來計算每個個體泌乳5~305 d 每天的加性遺傳方差和永久環境方差,最后合成305 d 的遺傳參數[30]。新疆地區中國荷斯坦牛5~305 d 產奶量、乳脂量、乳蛋白量和體細胞評分的總遺傳力分別為0.18~0.22、0.12~0.15、0.14~0.19、0.06~0.10,與北京地區[24]、寧夏地區[14,31]、河南省[32]和山東省[33]利用測定日模型估計的中國荷斯坦牛產奶量、乳脂量、乳蛋白量和體細胞評分遺傳力的結果一致。各性狀遺傳力估計值均略低于Miglior[34]基于四階Legendre 多項式測定日模型對中國荷斯坦牛頭胎DHI 記錄的分析結果,產奶量、乳脂量、乳蛋白量和體細胞評分4 個性狀的估計遺傳力分別為0.291、0.222、0.251 和0.092。與王國龍等[10]估計的烏魯木齊種牛場荷斯坦種公牛產奶量遺傳力0.14、乳脂率遺傳力0.13、乳蛋白率遺傳力0.14 的結果一致。曹福存等[35]人估計了北京地區81 385 頭荷斯坦母牛的LSCC(體細胞數自然對數)遺傳力為0.12,個體測定日LSCC 的遺傳力為0.06~0.11。
種公牛的遺傳品質直接關系到牛群遺傳改良的效果。我國種公牛站已經處于轉型發展的關鍵階段,面對市場需求和國際化的發展趨勢,聯合育種正在創新發展,對于種公牛聯合后裔測定機制已經逐步形成,為了不斷加大種公牛后裔測定力度,由種公牛站、核心育種場和育種企業聯合組成了奶牛后裔測定聯盟(北方聯盟、香山聯盟),聯合后裔測定[36]。在各省區遺傳評估的基礎上,積極促進數據統一化,為更好的進行聯合后裔測定夯實數據基礎。此外,積極進行種公牛的遺傳評定有助于擴大優秀種公牛的遺傳優勢,促進品種的選育提高,加快群體遺傳進展。
新疆地區中國荷斯坦牛305 d 產奶量、乳脂量、乳蛋白量和體細胞評分的遺傳力在0.07~0.22 之間,利用隨機回歸測定日模型進行種公牛遺傳評定時,需要篩選出顯著影響性狀的效應和擬合度較高的泌乳曲線模型才能獲得遺傳評估的最優模型。種公牛的后裔測定工作是評定種公牛遺傳品質的重要環節,對于牛群的遺傳改良效果影響重大,也是新疆地區中國荷斯坦牛種公牛凍精銷售的重要支撐。