应用随机回归模型对北京地区荷斯坦牛产奶性状的遗传分析

杨文静，王 晔，闫青霞，麻 柱，刘剑锋*，张胜利*

（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）本研究结果可为北京地区进行遗传评估和制定合理育种目标提供理论参考。