泌乳曲线作为描述牛群产奶特性的基本方法,在奶牛生产和育种中有重要作用。利用泌乳曲线,可以将奶牛不同泌乳时间的累积产奶量统一校正到305 d产奶量,以此作为奶牛一个泌乳期的生产记录,用于奶牛遗传评估或者评价奶牛的营养和健康状况。自1927年Gains首先利用数学模型描述了产奶量与时间的函数关系以来,各国研究者对泌乳曲线建立了许多数学模型^[1],开辟了泌乳曲线的生物统计模型研究领域。1967年Wood^[2]提出的不完全伽玛模型,是当前应用最广的模型。1987年Wilmink^[3]在Wood模型的基础上提出了改进Wood模型。同年Schaeffer等^[4]推出回归模型,并与Wood模型进行了全面比较。2000年Pollot等^[5]提出了描述奶绵羊泌乳曲线的Pollot模型,其优势在于利用了2次对数曲线,2008年Pollot^[6]将其用于描述奶牛泌乳曲线。

奶牛泌乳曲线的拟合研究,国内也有报道。1990年郭智群等^[7]用北京郊区奶牛的数据对3个模型进行了拟合比较;1993年郭智群等^[8]又对模型拟合的参数进行了遗传分析;1999年王雅春等^[9]利用7个模型分别对西门塔尔牛、蒙贝利亚牛和黑白花牛进行了拟合。

近年来,随着我国荷斯坦牛群遗传素质的改良、饲养条件及生产管理技术的进步,牛群生产性能有了显著的提高,因此有必要根据最新的产奶数据拟合泌乳曲线。笔者选取6种描述泌乳曲线的数学模型,利用SAS 9.2非线性回归(NLIN)过程对中国荷斯坦牛各胎次产奶量进行拟合,筛选拟合效果最好的模型,进而根据此模型计算累积产奶量,提出产奶时间的校正系数。在此基础上,结合实际累积产奶量和相应的校正系数,将数据库中所有未达到305 d产奶量记录的奶牛校正到305 d。

本研究收集了近10年黑龙江、北京、山东、陕西、上海、四川等10省市80个牛场213 818头荷斯坦奶牛的生产数据资料,形成了包含1 084 321条测定日记录的总数据库。在此基础上根据以下原则整理数据:①删除泌乳时间大于305 d或者小于5 d的记录;②删除产奶量小于0.5 kg的记录;③删除测定间隔大于60 d的记录;④删除测定记录少于5次的泌乳期记录;⑤对各个泌乳期的产犊年龄进行限制:第1胎的产犊年龄范围为18~48月龄,第2胎为25~66月龄,第3胎为30~80月龄。

数据经过整理后,共得到71个牛场43 812头中国荷斯坦牛77 490个泌乳期的609 784条测定日记录。其中第1胎奶牛35 141头,有测定日记录280 263条;第2胎奶牛26 213头,有测定日记录204 615条;第3胎奶牛16 136头,有测定日记录124 906条。

本研究选用以下6种模型进行泌乳曲线拟合^{[2, 3, 4, 5, 6]}:

Wood不完全伽玛模型(又称Wood模型):

逆多项式模型:

Ali-Schaeffer模型:

多项式回归模型(六次):

Wilmink模型:

Pollott模型:

以上各模型中:x为泌乳时间(d),y为第x天的产奶量。拟合曲线时,曲线参数的估计采用SAS 9.2软件非线性回归分析(NLIN)模块的Gauss-Newton迭代算法^[10]。

采用拟合度(R²值)和估计误差标准差(SEE)为度量标准。前者表示回归模型中能够解释的变异占总变异的比例,其值越大越好;后者是实际数据与模型拟合数据之间差值的标准差,描述在平均数周围的分布情况,理想的曲线拟合结果估计标准误接近零。

利用拟合效果最好的泌乳曲线计算各泌乳时间的估计累积产奶量,以305 d的估计累积产奶量除以各泌乳时间的估计累积产奶量得到特定产奶时间的校正系数,即:

式中:α是t d的校正系数;t为泌乳时间(d),通常小于305 d。当t≥305时,则α=1,F(305)为305 d的估计累积产奶量,F(t)为t d的估计累积产奶量。

校正305 d产奶量由实际累积产奶量乘以对应的校正系数得到,即:

式中:₃₀₅为校正305 d产奶量;Y_t为t d的实际累积产奶量;α是t d的校正系数。Y_t的计算采用测定间隔法(test interval method,TIM)^[11]:①获得测定日的产奶量;②估计测定间隔的产奶量,为(上一次测定记录+本次测定记录)/2×测定间隔时间;③将所有测定间隔相加得出特定泌乳时间的累积产奶量。

在计算Y_t的程序中,为了更准确地估计产奶量,对第1次测定间隔进行了处理:针对第1次奶牛生产性能(DHI)测定等于5 d的记录,取实际产奶量为第1条记录进行计算;对第1次DHI测定大于5 d的记录,以F(5)即第5天的估计产奶量作为第1条记录进行计算。

通过执行SAS 9.2软件回归分析过程中的非线性(NLIN)回归过程,获得了中国荷斯坦牛第1、第2、第3胎次和总数据的6个泌乳曲线模型参数估计值、拟合度和拟合估计误差标准差(表1)。

由拟合结果可以发现:6种泌乳曲线在各胎次的R²值均较高,达到0.94以上。3个胎次间相比,第1胎的拟合度最高,第2胎次之,第3胎最低。在各胎次内,6种模型的拟合效果总体上非常接近,但也存在一定差异。六次多项式回归模型在3个胎次中R²值都是最大,Ali-Schaeffer模型和Wilmink模型在前2个胎次中R²也是最大的,其次是Wood不完全伽玛模型,再次是逆多项式模型,而Pollott模型由于其参数过少,所以在3个胎次中都是最小。六次多项式回归模型的估计误差标准差也最小,其次为Ali-Schaeffer模型,然后为Wilmink模型、Wood不完全伽玛模型和逆多项式模型,最差的仍是Pollott模型。由此得出6个模型在本研究中从优到劣依次为:六次多项式回归模型、Ali-Schaeffer模型、Wilmink模型、Wood不完全伽玛模型、逆多项式模型、Pollott模型。因此,本研究选择六次多项式回归模型作为泌乳曲线模型来拟合泌乳曲线,计算估计累积产奶量,进而计算产奶时间的校正系数。

选用六次多项式回归模型,绘制第1、第2、第3胎次和总数据的泌乳曲线(图1)。从图中可以看出,各条泌乳曲线波形趋势一致,都是在40~50 d间上升到一个峰值,然后再缓慢下降。第2胎和第3胎曲线十分接近,出现峰值的时间早于第1胎,同时峰值后的下降速度也明显快于第1胎,即第1胎的持续力高于后面的胎次。总体而言,第2、第3胎曲线高于第1胎曲线,第3胎曲线略高于第2胎曲线。

	图 1 中国荷斯坦牛各胎次的泌乳曲线 Fig. 1 Lactation curves for different parities in Chinese Holsteins
图选项

根据拟合效果最好的模型,即六次多项式回归模型拟合泌乳曲线,利用公式计算出中国荷斯坦牛各胎次90~305 d产奶时间的校正系数,结果见表2。在实际应用中,用90 d或更早的实际累积产奶量去估计305 d累积产奶量的偏差非常大,而泌乳时间在200 d以上的产奶量记录校正具有高可靠性。

根据TIM方法计算出实际积累产奶量,泌乳时间大于等于305 d的取实际积累产奶量,泌乳时间不足305 d的,用本研究所得到的校正系数,根据其所在的胎次将其校正到305 d,即由实际积累产奶量乘以相应的校正系数(公式(8))得到。从中国荷斯坦牛各胎次的平均305 d产奶量(表3)可以看出,第1胎次的平均产奶量与第2、第3胎次存在极显著差异,而第2胎和第3胎的平均产奶量差异没有达到显著水平。

本研究收集了80个牛场213 818头中国荷斯坦牛132 609个泌乳期的1 084 321条测定日记录,经过筛选后剩余71个牛场43 812头中国荷斯坦牛77 490个泌乳期的609 784条测定日记录,删除率达到43.76%。泌乳胎次大于3胎和1头牛在一个泌乳期内的测定条数过少是重要的原因。结合实际生产中奶样的采集和DHI数据的测定工作,反映出我国DHI测定工作开展的艰辛,当前奶牛DHI测定的系统性和完整性还有待提高,所以在数据分析前需要对异常记录进行识别和剔除,才能保证后续的研究。

国内对非线性拟合大都使用常规的线性化方法,本研究对6种模型均采用非线性最小二乘法求解。按照拟合度(R²)值和估计误差标准差(SEE)对6种模型从优到劣进行排序如下:模型4(0.943,7.708)、模型3(0.941,7.752)、 模型5(0.941,7.753)、模型1(0.941,7.756)、 模型2(0.940,7.775)、 模型6(0.940,7.788)。这与Santos等^[12]给出的:模型4 、 模型3、模型5 、模型1排序一致,与王雅春^[9]给出的几种模型按照R²排名的顺序也一致,说明六次多项式回归是对中国荷斯坦牛的产奶量性状拟合度效果最好的泌乳曲线模型。用6个模型对数据拟合,发现所有模型的R²值都较高,接近国外文献报道的0.89~0.95^[13],高于国内滕晓红^[11]报道的0.6~0.8,这主要原因是本研究数据量大,且数据都来源于大型高产牛场,牛群一致性较好。

胎次、饲养管理水平和产犊季节对泌乳曲线的形状均有影响^[14]。对不同胎次的泌乳时间校正系数的比较可以看出,校正系数随着时间的增加而下降,相同时间时,第1胎的校正系数总是最大,其次是第2胎,再次是第3胎;3个胎次间的差异随着时间的增加逐渐缩小。这也印证了Togashi等^[13]的报道,第1胎的持续力大于后面的胎次,第2胎的持续力大于第3胎次,即奶牛泌乳持续力随着胎次的增加而递减。不难看出,校正系数受胎次的影响很大,为了准确地计算校正305 d产奶量,对各胎次应分别使用一套校正系数。本研究中使用的数据都来自于大型高产牛场,牛群饲养管理水平比较一致。数据量较大,各产犊季节数据分布均匀。因此没有考虑不同牛场的饲养管理差异和产犊季节对泌乳曲线的影响,所得结果为大型高产牧场各产犊季节的平均校正系数。

本研究用校正系数计算获得中国荷斯坦牛各胎次的305 d产奶量从低到高依次是,1胎(7 888.84 kg)、3胎(8 563.05 kg)、2胎(8 568.05 kg),说明第1胎次的泌乳性能低于以后的胎次。但是第2胎次与第3胎次之间的关系和国内外其他文献报道不同^{[11, 15]},本研究中第2胎次略高于第3胎次。其原因一方面是第3胎次数据量少于第2胎次,另一方面可能是我国荷斯坦奶牛随着胎次的增加,奶牛乳房炎的发病率有所增加^[16],造成第3胎次牛的记录中低产牛增加。

综上,本研究利用中国荷斯坦牛大规模DHI测定日记录数据,进行泌乳曲线拟合研究,获得了拟合效果最好的泌乳曲线及其参数,并据此得到了305 d产奶量的校正系数,为奶牛产奶性能的早期预测和遗传评估提供了可靠的方法。