應(yīng)用隨機(jī)回歸模型對(duì)北京地區(qū)荷斯坦牛產(chǎn)奶性狀的遺傳分析

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

（1.中國(guó)農(nóng)業(yè)大學(xué)動(dòng)物科學(xué)技術(shù)學(xué)院，北京 100193；2.中國(guó)奶業(yè)協(xié)會(huì)，北京 100193；3.北京奶牛中心，北京 100192）

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

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

本研究收集了北京地區(qū)1998—2018 年中國(guó)荷斯坦牛測(cè)定日數(shù)據(jù)，采用4 階Legendre 多項(xiàng)式隨機(jī)回歸模型對(duì)中國(guó)荷斯坦牛的產(chǎn)奶性狀進(jìn)行遺傳參數(shù)和育種值估計(jì)，并進(jìn)行遺傳進(jìn)展分析，為北京地區(qū)制定更加完善的育種計(jì)劃和遺傳評(píng)估提供理論參考。

1 材料與方法

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

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

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

其中，yijklm：第l 頭母牛的第m個(gè)觀察值；HTDi：第i個(gè)場(chǎng)-測(cè)定日效應(yīng)；CYMSj：第j個(gè)產(chǎn)犢年-季效應(yīng)；bkn：第k個(gè)產(chǎn)犢月齡效應(yīng)水平下第n個(gè)固定回歸系數(shù)；aln：系譜中第l 個(gè)個(gè)體加性遺傳效應(yīng)下第n個(gè)隨機(jī)回歸系數(shù)；pln：第l 個(gè)個(gè)體永久環(huán)境效應(yīng)下第n個(gè)隨機(jī)回歸系數(shù)；Xbn：Legendre 多項(xiàng)式第n個(gè)協(xié)變量；Zan（w）和Wpn（w）：加性遺傳和永久環(huán)境效應(yīng)的第n個(gè)協(xié)變量；w：標(biāo)準(zhǔn)化后的泌乳天數(shù)；eijklm：隨機(jī)殘差，并假設(shè)殘差同質(zhì)。

對(duì)應(yīng)模型的混合模型方程組MME 為：

使用估計(jì)育種值間的線性相關(guān)作為性狀間的遺傳相關(guān)，使用觀測(cè)值之間的線性相關(guān)作為性狀之間的表型相關(guān)，其公式如下：

2 結(jié)果

2.1 產(chǎn)奶性狀描述性統(tǒng)計(jì)量 對(duì)原始數(shù)據(jù)質(zhì)控后，本研究首先統(tǒng)計(jì)了5 個(gè)產(chǎn)奶性狀的表型信息，結(jié)果見表1。

表1 測(cè)定日產(chǎn)奶性狀描述性統(tǒng)計(jì)

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

隨機(jī)回歸模型可以體現(xiàn)整個(gè)泌乳期內(nèi)不同方差組分的變化趨勢(shì)（圖1）。產(chǎn)奶性狀加性遺傳方差在泌乳開始后的第一個(gè)月內(nèi)均呈大幅下降趨勢(shì)，泌乳中后期逐漸趨于平緩或小幅回升。產(chǎn)奶性狀永久環(huán)境效應(yīng)方差在泌乳開始后的第一個(gè)月內(nèi)下降至最低值，隨后保持平穩(wěn)。

圖1 產(chǎn)奶性狀加性方差（Va）和永久環(huán)境效應(yīng)方差（Vpe）隨泌乳天數(shù)的變化曲線

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

圖2 產(chǎn)奶性狀的遺傳力隨泌乳天數(shù)的變化曲線

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

表2 產(chǎn)奶性狀表型相關(guān)和遺傳相關(guān)

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

圖3 產(chǎn)量性狀（MY、FY、PY）遺傳進(jìn)展趨勢(shì)圖

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

圖4 乳脂率、乳蛋白率遺傳進(jìn)展趨勢(shì)圖

3 討 論

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

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

表3 前人相關(guān)研究整理表

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

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

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

4 結(jié) 論

本研究采用單性狀隨機(jī)回歸模型對(duì)北京地區(qū)荷斯坦牛5 個(gè)產(chǎn)奶性狀的遺傳參數(shù)和育種值進(jìn)行估計(jì)，得出以下結(jié)論：

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

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

3）本研究結(jié)果可為北京地區(qū)進(jìn)行遺傳評(píng)估和制定合理育種目標(biāo)提供理論參考。