张梦华，罗汉鹏，魏 趁，胥 磊，王 丹 ，张晓雪，谭世新，李红燕，黄锡霞*，王雅春*

（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.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 结果与分析

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 讨 论

本研究所使用的测定日模型包括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]。在各省区遗传评估的基础上，积极促进数据统一化，为更好的进行联合后裔测定夯实数据基础。此外，积极进行种公牛的遗传评定有助于扩大优秀种公牛的遗传优势，促进品种的选育提高，加快群体遗传进展。

4 结 论

新疆地区中国荷斯坦牛305 d 产奶量、乳脂量、乳蛋白量和体细胞评分的遗传力在0.07～0.22 之间，利用随机回归测定日模型进行种公牛遗传评定时，需要筛选出显著影响性状的效应和拟合度较高的泌乳曲线模型才能获得遗传评估的最优模型。种公牛的后裔测定工作是评定种公牛遗传品质的重要环节，对于牛群的遗传改良效果影响重大，也是新疆地区中国荷斯坦牛种公牛冻精销售的重要支撑。