應用隨機回歸模型對北京地區荷斯坦牛產奶性狀的遺傳分析

楊文靜，王 曄，閆青霞，麻 柱，劉劍鋒*，張勝利*

（1.中國農業大學動物科學技術學院，北京 100193；2.中國奶業協會，北京 100193；3.北京奶牛中心，北京 100192）

自農業農村部實施奶牛生產性能測定（DHI）補貼項目以來，2008—2017 年十年間，我國參測牛群的平均日產奶量由2008 年的23.8 kg 提高到2017 年的29.0 kg，乳脂率由3.64% 提高到3.89%，乳蛋白率由3.19% 提高到3.35%[1]。北京作為我國最早開展奶牛生產性能測定工作的重要地區之一，其在生產性能測定規模、種公牛選育、后裔測定和育種數據庫建設等方面都處于全國領先水平。然而，雖然我國在奶牛育種方面取得了一定的進展，但由于我國常規育種工作起步較晚，基礎薄弱，與歐美等奶業發達國家相比，仍有較大差距，主要表現在優秀種公牛依賴于進口，奶牛生產效率普遍不高，牛奶總產量的增加依靠奶牛數量的增加，遺傳改良體系不夠完善等方面。由此可見，我國奶業生產水平提升還有較大空間。奶牛生產的首要目標是提高母牛單產、提升牛奶品質、獲得更高預期收益。因此，積極開展DHI測定及產奶性狀的遺傳評估，以育種手段提升奶牛群體選擇進展十分必要[2]。

目前，估計奶牛產奶性狀育種值所用的模型主要是動物模型和隨機回歸測定日模型。研究發現，不管在理論研究還是實踐應用中，隨機回歸測定日模型都具有較強的優越性[3-4]。它的優點主要表現為，可以在不預測全期泌乳量的情況下，使用所有的測定日記錄為每頭奶牛擬合不同的泌乳曲線。目前，隨機回歸模型因其靈活性已經被應用于越來越多的育種實踐中[5-8]。測定日模型擬合泌乳曲線的方法包括Wood 不完全伽馬函數模型[9]，Ali-Schaeffer 回歸法[10]，Willkmink 指數函數[11]以及多項式[12]等方法。而大量研究表明，利用高階Legendre多項式擬合泌乳曲線的方法要優于其他方法[4,7]。

本研究收集了北京地區1998—2018 年中國荷斯坦牛測定日數據，采用4 階Legendre 多項式隨機回歸模型對中國荷斯坦牛的產奶性狀進行遺傳參數和育種值估計，并進行遺傳進展分析，為北京地區制定更加完善的育種計劃和遺傳評估提供理論參考。

1 材料與方法

1.1 數據收集與整理 本研究數據由中國奶牛數據中心提供，包含了北京市1998—2018 年間131 個牛場的99 823 頭中國荷斯坦牛第一胎DHI 記錄。本研究首先對測定日數據進行了嚴格的篩選，篩選標準為：泌乳天數在5～305 d；泌乳期內記錄條數≥3；測定日產奶量1～80 kg；乳脂率1.4%～6.2%；乳蛋白率2.0%～5.0%。經篩選后，剩余的頭胎測定日記錄共計817 208 條，所分析的性狀為產奶量、乳脂量、乳蛋白量、乳脂率、乳蛋白率。

同時，對收集到的系譜作如下處理：①未知的出生日期以及父母號缺失的位置記錄為0；②對系譜個體號以及父母號進行重編；③基于DHI 記錄中提供的基本個體號信息，將有表型的母牛進行系譜追溯，向上至少追溯三代；④將系譜文件排列成DMU 軟件可識別的格式文件。整理后的系譜文件共187 695 條記錄，由四列組成：個體號、父號、母號和出生日期，個體號列包括所有在父號列和母號列出現過的個體，系譜與數據文件的個體編號一一對應。

1.2 遺傳評估模型 采用單性狀隨機回歸模型，其中，固定效應包括場-測定日效應、產犢年-季效應（產犢季節每年劃分3 個水平[13-14]：1—4 月、11 月和12 月為一季；5 月和10 月為一季；6—9 月為一季）、以及固定回歸擬合產犢月齡效應，隨機效應包括加性遺傳效應和永久環境效應。固定回歸和隨機回歸均采用4 階Legendre 多項式擬合回歸曲線。方差組分和育種值估計均采用DMU（Derivativefree Multivariate）軟件（v 6.0）的AI 模塊[15]，模型如下：

其中，yijklm：第l 頭母牛的第m個觀察值；HTDi：第i個場-測定日效應；CYMSj：第j個產犢年-季效應；bkn：第k個產犢月齡效應水平下第n個固定回歸系數；aln：系譜中第l 個個體加性遺傳效應下第n個隨機回歸系數；pln：第l 個個體永久環境效應下第n個隨機回歸系數；Xbn：Legendre 多項式第n個協變量；Zan（w）和Wpn（w）：加性遺傳和永久環境效應的第n個協變量；w：標準化后的泌乳天數；eijklm：隨機殘差，并假設殘差同質。

對應模型的混合模型方程組MME 為：

使用估計育種值間的線性相關作為性狀間的遺傳相關，使用觀測值之間的線性相關作為性狀之間的表型相關，其公式如下：

2 結果

2.1 產奶性狀描述性統計量 對原始數據質控后，本研究首先統計了5 個產奶性狀的表型信息，結果見表1。

表1 測定日產奶性狀描述性統計

2.2 產奶性狀方差組分及遺傳力估計 按上述遺傳評估模型，產奶量、乳脂率、乳蛋白率、乳脂量、乳蛋白量5～305 d 的總遺傳力估計值分別為：0.28、0.53、0.63、0.25、0.24。其中，產奶量、乳脂量、乳蛋白量3 個性狀的遺傳力為中等遺傳力，而乳脂率、乳蛋白率2 個性狀為較高遺傳力。

隨機回歸模型可以體現整個泌乳期內不同方差組分的變化趨勢（圖1）。產奶性狀加性遺傳方差在泌乳開始后的第一個月內均呈大幅下降趨勢，泌乳中后期逐漸趨于平緩或小幅回升。產奶性狀永久環境效應方差在泌乳開始后的第一個月內下降至最低值，隨后保持平穩。

圖1 產奶性狀加性方差（Va）和永久環境效應方差（Vpe）隨泌乳天數的變化曲線

根據奶牛產奶性狀的加性遺傳方差和永久環境效應方差隨泌乳天數發生變化的趨勢，本研究探索了遺傳力在整個泌乳期內的變化趨勢（圖2）。可以看出，產奶量、乳脂量、乳蛋白量3 個性狀的日遺傳力在泌乳初期取得最大值，泌乳開始后的第一個月內迅速下降至最低值，在泌乳中后期有小幅上升，漸趨于平緩。然而，乳脂率雖然與上述3 個性狀有相同的變化趨勢，但其日遺傳力最大值在泌乳末期取得。乳蛋白率的日遺傳力在整個泌乳期一直呈上升趨勢。

圖2 產奶性狀的遺傳力隨泌乳天數的變化曲線

2.3 各產奶性狀間的相關分析 本研究統計了各產奶性狀間的相關性估計結果（表2）。從中可以得知，無論是表型相關還是遺傳相關，產奶量、乳脂量、乳蛋白量彼此間都為高度正相關（0.71～0.94），其中產奶量和乳蛋白量的表型相關高達0.93，遺傳相關高達0.94；而產奶量與乳脂率、乳蛋白率的表型和遺傳相關都呈負相關（-0.38～-0.18）；乳脂率和乳蛋白率之間有較強正遺傳相關（0.42）；但乳脂量和乳蛋白率，以及乳脂率和乳蛋白量之間表型和遺傳相關性均較低（-0.29～0.001）。

表2 產奶性狀表型相關和遺傳相關

2.4 產奶性狀遺傳進展分析 本研究根據荷斯坦牛的出生年度，計算了母牛估計育種值的平均值，分析了1990—2016 年間北京地區荷斯坦牛的遺傳進展情況，結果列于圖3 和圖4?？梢钥闯?，產奶量、乳脂量、乳蛋白量3 個性狀的遺傳進展變化趨勢相似（圖3），呈波動型增長，主要可分成2 個階段：1990—2004 年間群體遺傳進展緩慢，改良效果不佳；2005—2016 年間遺傳進展每年加快，改良效果較理想，這可能與2005年原農業部開始實施的“奶牛良種補貼項目”和2008年實施的“全國奶牛群體遺傳改良計劃（2008—2020年）”及“奶牛生產性能測定項目”有一定關聯，這些政策促進了奶牛群體改良技術的落實與提高。1990—2016 年，產奶量、乳脂量、乳蛋白量的平均年遺傳進展分別為：10.52、0.30、0.32 kg。但2005—2016 年，3 個性狀的平均年遺傳進展分別為：18.77、0.72、0.64 kg，改良速度顯著提高。

圖3 產量性狀（MY、FY、PY）遺傳進展趨勢圖

1990—2016 年，北京地區荷斯坦奶牛乳脂率和乳蛋白率的平均年遺傳進展分別為：-0.00046%、-0.00011%。2 個性狀遺傳進展波動較大，并最終基本保持穩定（圖4）。其中，2008—2009 年以及2012—2016 年均呈上升趨勢。

圖4 乳脂率、乳蛋白率遺傳進展趨勢圖

3 討 論

2018 年全國DHI 參測奶牛的產奶性狀平均值分別為：產奶量29.60 kg，乳脂率3.76%，乳蛋白率3.23%，乳脂量1.11 kg，乳蛋白量0.95 kg。整體而言，北京地區中國荷斯坦牛產奶性狀均值高于同期全國DHI 測定群體的產奶性狀均值，其主要原因可能是北京地區的生產管理水平較高，并且較早開展生產性能測定，這使得牛群的遺傳潛能得到充分發揮，遺傳改良效果較好。

本研究將產奶性狀的遺傳力估計結果與前人的研究進行了比較（表3）。周云鵬等[16]2007 年對北京、黑龍江和天津地區的中國荷斯坦牛采用隨機回歸模型配合5 階固定效應和加性效應、7 階永久環境效應Legendre多項式對產奶量、乳脂量、乳蛋白量的遺傳參數進行估計。結果表明：產奶量、乳脂量、乳蛋白量的305 d 遺傳力分別為0.22、0.14、0.19。Miglior 等[17]研究中國荷斯坦牛第1 胎305 d 遺傳力分別是：產奶量為0.29，乳脂量0.22，乳蛋白量0.25，李想[18]使用隨機回歸模型估計北京地區中國荷斯坦牛產奶性狀第一胎305 d遺傳力分別是：產奶量0.26，乳脂量0.17，乳蛋白量0.25。2016 年加拿大奶牛工作網（https://www.cdn.ca/home.php）公布的結果為：產奶量遺傳力0.28，乳脂量遺傳力0.26，乳蛋白量遺傳力0.26，本研究遺傳力估計結果整體高于周云鵬[16]、李想[18]的研究結果，與Miglior[17]、加拿大估計遺傳力結果較為相近，其原因可能與不同研究所使用的數據量大小以及群體結構差異有關。另外，使用的Legendre 多項式的階數不同，也可能對結果產生一定的影響。

表3 前人相關研究整理表

本研究對測定日遺傳力趨勢探究表明，在泌乳期前1 個月內，產奶量性狀遺傳力呈下降趨勢，并且遺傳力在泌乳期前10 d 最高，這與Jamrozik[19]、Miglior 等[17]研究結果類似。而在泌乳期第2 個月內取得最低值，說明該時段內環境因素對產奶表型產生很大的影響，推測該時間段可能是泌乳性能相關基因差異表達的關鍵時間點，可通過調控營養、管理等環境條件改善產奶性狀的表現。同時，不同泌乳階段遺傳力的差異也佐證了測定日模型的有效性，即將測定日記錄簡單看成重復度量值是不合理的。

本研究各產奶性狀間的表型相關和遺傳相關的結果與前人研究結果類似[17-18]，產奶量、乳脂量、乳蛋白量之間無論表型相關還是遺傳相關均為高度正相關，而產奶量與乳脂率、乳蛋白率間存在中等負遺傳相關，乳脂率和乳蛋白率之間有較強正遺傳相關，在做育種規劃時，對產量性狀的選擇可以充分利用性狀間的遺傳相關，在對乳脂量和乳蛋白量選擇的同時產奶量也能得到提高，可以降低育種成本。

產奶性狀的遺傳進展分析結果反映了我國奶牛育種工作的階段性成果。遺傳進展趨勢在前人的研究中也有所體現[20]。自2008 年以來，產奶性狀的遺傳進展明顯加快，這是因為2008 年農業部發布了中國奶牛群體遺傳改良計劃[21]，并實施了“奶牛生產性能測定補貼項目”，說明奶牛遺傳改良工作卓有成效，北京地區的奶牛遺傳素質有所提高。乳脂率、乳蛋白率的平均年度遺傳進展分別為：-0.00046%、-0.00011%，基本保持穩定，其原因可能與中國荷斯坦奶牛的階段性育種目標有關。產奶量一直以來是中國荷斯坦牛的主要育種目標，奶牛場選配時主要追求產奶量，對其他性狀的重視不夠，因此造成了乳脂率和乳蛋白率遺傳進展波動較大且進展緩慢。自2008 年以來，乳脂率、乳蛋白率的平均育種值有所提高，說明北京地區奶牛育種在保持高產奶量的同時也加強了對乳成分的選育。

4 結 論

本研究采用單性狀隨機回歸模型對北京地區荷斯坦牛5 個產奶性狀的遺傳參數和育種值進行估計，得出以下結論：

1）產奶量、乳脂率、乳蛋白率、乳脂量、乳蛋白量的305 d 泌乳期遺傳力分別為0.28、0.53、0.63、0.25、0.24。

2）1990—2016 年間產奶量、乳脂率、乳蛋白率、乳脂量、乳蛋白量的年度平均遺傳進展分別為10.52 kg、-0.00046%、-0.00011%、0.30 kg、0.32 kg。

3）本研究結果可為北京地區進行遺傳評估和制定合理育種目標提供理論參考。