长寿性状是奶牛养殖中具有重要经济价值的功能性状之一，当奶牛寿命延长时，可减少后备牛的数量，降低总饲养成本，因此牛群寿命长短对牛场收益有着重要影响^[1-2]。长期以来，奶牛育种主要强调提高泌乳性能，包括产奶量以及乳成分^[3-4]。但随着产奶量等泌乳性能的提高，包括被动淘汰率高、寿命短在内的奶牛功能性状问题日益明显^[5]。目前，大多数奶业国家平均淘汰胎次小于3胎，导致利用年限低、生产寿命短，影响终身产奶量^[6-8]。随着平衡育种理念的不断加强，许多奶业发达国家已逐步将长寿性状纳入综合选择指数并日益提高权重^[9-11]。因此，研究长寿性状的遗传结构，对实施平衡育种和提高奶牛整体生产性能具有重要意义。

目前，国内外关于奶牛长寿性状遗传参数估计的报道较多^[12-19]。如，Saowaphak等^[14]发现泰国杂交荷斯坦奶牛生产寿命的遗传力为0.1；Wasana等^[15]报道韩国荷斯坦奶牛1~3胎次生产寿命的遗传力为0.06~0.13；Sewalem等^[16]提出加拿大娟姗牛和爱尔夏牛生产寿命遗传力在0.01~0.25之间。在国内，李想等^[17]使用不同动物模型估计中国荷斯坦牛的生产寿命遗传力为0.047~0.069。董明明等^[18]报道新疆地区荷斯坦牛在群寿命和繁殖寿命遗传力为0.11。Zhang等^[19]对北方地区中国荷斯坦牛的长寿性状遗传参数进行估计，发现遗传力介于0.038~0.09。上述的研究均证明，荷斯坦牛长寿性能具有一定的可遗传性，因此通过遗传育种的手段来提高牛群的长寿性能是可行的。然而，有关我国南方地区奶牛长寿性状遗传参数估计尚未见报道。

本研究基于江苏地区9个规模化牧场荷斯坦牛的长寿性及淘汰情况的记录，拟开展：1)估计以江苏为代表的南方地区规模化牧场奶牛长寿性状的遗传参数，包括方差组分、遗传力、遗传相关和表型相关等；2)研究近年来江苏地区规模化牧场长寿性状的遗传趋势。本研究可为以江苏为代表的南方地区奶牛长寿性状的选育提供参考。

1 材料与方法 1.1 数据来源与整理

本研究的数据来源于我国江苏地区9个规模化荷斯坦牛牧场，牛场存栏均在1 000头以上。原始数据从牧场管理软件中导出，主要包括2017—2023年间的离群记录以及所有离群个体出生、产犊以及干奶等繁殖事件记录。将收集的数据进行过滤整理，剔除以下数据：1)成母牛头胎产犊日龄小于600 d或大于1 280 d；2)离群原因为出售、转场等非淘汰因素。质控后，共保留19 756头成母牛用于后续分析。

1.2 长寿性状定义与说明

衡量奶牛长寿性的方法有多种，本研究中以个体在牧场中整个或某个阶段的生存天数来统计长寿性，主要包括在群寿命、生产寿命、利用胎次以及胎次寿命性状L1~L4等的7个指标(具体描述见表 1)。除利用胎次外，其他6个性状均以天为单位。其中，L1~L4为胎次寿命性状，参考Clasen等^[20]对胎次寿命统计分析方法。统计L1-L4观测天数时，为尽可能矫正因繁殖问题产生的无效泌乳天数，如久配不妊等繁殖因素导致胎次间隔过长，本研究将每个胎次最高生产寿命设置为365天。因此，L1、L2、L3和L4的最大观测值分别为365、730、1 095和1 460天。

1.3 数据分析

1.3.1 长寿性状描述性统计 　　本研究首先对淘汰成母牛不同淘汰年份以及不同离群胎次中各淘汰原因的占比进行分析，并分别统计各淘汰原因所占比例。其次，使用19 756头成母牛的淘汰、产犊、干奶等记录，对在群寿命、生产寿命、利用胎次以及胎次寿命性状分别进行描述统计。

1.3.2 成母牛长寿性状固定效应分析 　　为估计成母牛长寿性状的遗传参数，本研究将可能影响奶牛寿命的因素作为固定效应进行方差分析，随后将显著的固定效应纳入评估模型。将牧场和出生年合并为场－出生年效应，其他因素包括出生季节、头胎产犊月龄也作为固定效应。

本研究中，将牛出生季节划分为4个水平，包括3~5月为春季、6~8月为夏季、9~11月为秋季、12~次年2月为冬季；将头胎产犊月龄划分为13个水平，包括22月龄及以下、23月龄、24月龄、25月龄、26月龄、27月龄、28月龄、29月龄、30月龄、31月龄、32月龄、33月龄和34月龄及以上。随后，采用R语言中GLM模型对影响表型的固定效应进行检测，模型如下：

$ y_{i j k m}=H Y_i+\mathrm{S}N_j+\mathrm{AF}C_k+e_{i j k m} $

式中，y_ijkm为长寿性状表型观测值，HY_i为第i个牧场-出生年的效应，SN_j为第j个出生季节的效应，AFC_k为第k个水平的头胎产犊月龄效应，e_ijkm为随机误差。

1.3.3 成母牛长寿性状方差组分估计 　　基于上述固定效应分析结果，以成母牛在群寿命、生产寿命等7个长寿性状分别作为观测值，场-出生年、出生季节、头胎产犊月龄作为固定效应，个体加性遗传效应作为随机效应，估计各方差组分。单性状混合线性模型方程如下：

$ \mathrm{y}=\boldsymbol{X} b+\boldsymbol{Z} a+\mathrm{e}, $

其中，y为7个长寿性状表型观测值向量，b为固定效应向量，X为b的关联矩阵；a和e分别为加性遗传效应向量和随机误差向量，Z为a的关联矩阵。一般地，假设a和e均服从正态分布：a~N(0, Aσ_a²)，e~N(0, Iσ_e²)，其中A为亲缘关系矩阵，I是单位矩阵，σ_a²和σ_e²分别为加性遗传方差和随机误差方差。

1.3.4 长寿性与淘汰情况遗传参数评估 　　基于上述单性状混合模型得到各性状的方差组分后，可计算获得各性状的遗传力(h²)：

$ h^2=\frac{\sigma_a^2}{\sigma_a^2+\sigma_e^2} $

此外，为了评估不同的长寿性状之间是否具有共同的遗传基础，进一步对7个长寿性状开展多性状模型分析。但由于该分析没有收敛，因此使用上述线性混合模型对7个长寿性状间分别开展双性状模型分析，来估计性状之间的遗传相关(r_g12)，模型如下：

$ \left[\begin{array}{l} y_1 \\ y_2 \end{array}\right]=\left[\begin{array}{cc} \boldsymbol{X}_1 & 0 \\ 0 & \boldsymbol{X}_2 \end{array}\right]\left[\begin{array}{l} b_1 \\ b_2 \end{array}\right]+\left[\begin{array}{cc} \boldsymbol{Z}_1 & 0 \\ 0 & \boldsymbol{Z}_2 \end{array}\right]\left[\begin{array}{l} a_1 \\ a_2 \end{array}\right]+\left[\begin{array}{l} e_1 \\ e_2 \end{array}\right], $

其中，y_n为第n个长寿性状表型观测值向量，b_n为第n个长寿性状的固定效应向量，a_n为第n个长寿性状加性遗传效应向量，e_n为第n个长寿性状随机误差向量，X_n、Z_n为对应的关联矩阵。随后，遗传相关可计算为:

$ r_{g_{12}}=\frac{\sigma_{a 12}}{\sqrt{\sigma_{a_1}^2 \sigma_{a_2}^2}}, $

式中，σ_a12为性状1和性状2之间的加性遗传协方差，σ_a1²和σ_a2²分别为性状1和性状2的加性遗传方差。

1.3.5 长寿性状育种值估计与遗传趋势分析 　　对江苏地区规模化牧场荷斯坦成母牛在群寿命、生产寿命等7个长寿性状的育种值进行估计，评估模型与上述单性状混合线性模型一致。随后，对成母牛估计育种值可靠性进行质控，剔除估计育种值可靠性低于0.15的个体，并对保留个体进行遗传趋势分析。个体育种值的可靠性计算为：Rel=1-(SE_m)²/σ_a²，其中，SE_m为个体m的估计育种值的标准误。

本研究中方差及育种值估计分析均使用BLUPF90家族软件进行^[21]。遗传方差分析采用平均信息限制最大似然(AI-REML)估计，遗传效应预测采用最佳线性无偏预测(BLUP)方法。

2 结果 2.1 长寿性状与淘汰情况描述性统计分析

研究首先统计了成母牛的淘汰类型情况(图 1)以及7个长寿性状的表型信息(表 2)。从图 1可知，江苏地区规模牧场成母牛主要由于产能过低、繁殖障碍以及乳房健康等因素被淘汰。从表 2可知，江苏地区规模牧场荷斯坦成母牛平均在群寿命为1 640.87 d，约为4.50年；平均生产寿命为917.15 d，约为2.51年，平均生产寿命仅为在群寿命的55.89%；平均利用胎次为2.60胎。

2.2 长寿性状方差组分及遗传力估计

江苏地区规模化牧场荷斯坦成母牛长寿性状方差组分及遗传力结果如表 3所示。遗传力估计值变化范围为0.07~0.15，表明长寿性状属低遗传力性状。在胎次寿命性状L1-L4中，加性遗传方差和遗传力随胎次的增加呈现上升趋势。

2.3 长寿性状间的相关性分析

本研究统计了7个长寿性状间的表型和遗传相关，如表 4所示。无论是表型相关还是遗传相关，在群寿命、生产寿命与利用胎次两两之间均呈现高度正相关(0.88~0.99)。其中，在群寿命与生产寿命的表型和遗传相关性趋近于1。在群寿命、生产寿命、利用胎次与L1-L4之间具有较高的遗传相关性(0.79~0.99)，且随着胎次的增加，相关性越高。在群寿命、生产寿命、利用胎次与L1-L4之间的表型呈中等程度相关(0.27~0.99)，且相关性也随着胎次的增加而增加。

2.4 长寿性状遗传趋势分析

为进一步研究近年来成母牛长寿性状的遗传趋势，本研究计算了不同年份出生的成母牛估计育种值的平均值，分析了江苏地区2010—2020年间出生牛只长寿性的遗传趋势，详细结果见图 2和图 3。从图 2和图 3可以看出，在群寿命、生产寿命、利用胎次以及L1-L4的遗传进展变化趋势基本类似，且波动明显，具体可分成2个阶段：2010~2014年间出生牛只群体遗传进展整体呈上升趋势；2015年及以后出生牛只群体遗传进展整体呈下降趋势。

3 讨论 3.1 长寿性状与淘汰情况分析

一般认为，荷斯坦牛在第三胎之后才会实现较好的盈利，并且在第4~7个泌乳期能达到最大的产奶潜能^[22-24]。然而，目前大多数国家奶牛的平均利用胎次仅为3胎左右^[25-28]，因此绝大部分奶牛在尚未达到最大产奶量潜能前就被淘汰了，给牧场带来了严重的经济损失。以荷兰为代表以终生产奶量为育种目标的奶业国家，其牛群的生产寿命较长，平均为3.43年^[12]。本研究中，江苏地区荷斯坦牛成母牛平均利用胎次约为2.60胎，尚未达到最大产奶潜能的利用胎次。但该结果与周福振等^[29]报道江苏地区牛只离群胎次2.47胎，彭朋等^[30]报道北京、上海地区奶牛平均利用2.71胎次，以及董明明等^[18]报道的新疆地区中国荷斯坦牛平均利用2.78胎次接近。随着世界各奶业发达国家越来越重视奶牛的生产寿命，将长寿性能纳入选择指数，我国南方地区的牧场也应重视奶牛的长寿性。

3.2 长寿性状与淘汰情况方差组分及遗传力分析

本研究估计了长寿性状的加性遗传方差，发现胎次寿命性状中，加性遗传方差从L1到L4呈递增趋势，这一现象与Clasen等^[20]的研究结果相似。其中，L1、L2的加性遗传方差较小的原因可能是由于大多数奶牛在前两个胎次中存活下来，而每个胎次的生产寿命最高被设置为365 d。因此，不排除很大一部分个体的L1和L2观测值被分别设置为最大限度365和730 d，造成了这两个性状观测值的变异较小。

本研究发现，长寿性状的遗传力较低，范围为0.07~0.15，这一发现与先前的研究遗传力估计值范围基本相一致^[31-33]。Saowaphak等^[14]以泰国杂交荷斯坦奶牛为研究对象，发现生产寿命的遗传力为0.1。Wasana等^[15]分别评估韩国荷斯坦奶牛1~3胎次生产寿命的遗传参数，发现遗传力介于0.06和0.13之间。李想等^[17]使用单性状动物模型、单性状公畜模型以及多性状动物模型等估计北京地区中国荷斯坦牛长寿性状遗传力为0.047~0.069。董明明等^[18]对新疆地区中国荷斯坦牛长寿性状遗传参数进行评估，发现在群寿命和生产寿命的遗传力均为0.11。Hu等^[33]研究宁夏地区中国荷斯坦牛的长寿性状，发现遗传力为0.009~0.025。Zhang等^[19]进一步扩大研究群体范围，对北京、河北、内蒙古等北方地区的中国荷斯坦牛的长寿性状遗传参数进行评估，发现遗传力为0.038~0.09。本研究中长寿性状的统计方法与Zhang等^[19]和Hu等^[33]类似，但遗传力估计结果整体高于上述两项研究结果。其原因可能与使用的数据量大小以及群体结构差异有关。此外，完整、准确的系谱记录也是影响遗传参数分析的关键因素。

3.3 长寿性状遗传相关分析

本研究结果表明，在群寿命、生产寿命以及利用胎次之间无论是表型相关还是遗传相关均呈高度正相关关系(0.88~0.99)，这与以往研究结果一致^{[18, 31-37]}。生产寿命是在群寿命减去第一次产犊日龄，因此，这两个性状本质上是相同的长寿性状指标，可研究其中任何一个开展奶牛长寿性的育种。在群寿命、生产寿命和利用胎次与L1-L4等胎次寿命性状具有中等到较高的遗传相关和表型相关，且随着泌乳次数的增加，遗传相关性越高，这一发现与Zhang等^[19]和Hu等^[33]的研究结果相似。这是由于随着泌乳次数的增加，胎次寿命性状更接近生产寿命。因此，在一定程度上，生产寿命与L4的表型观测值是一致的。例如，个体在第四胎次时被淘汰，生产寿命与L4的观测值均计算为淘汰与头胎产犊之间的天数。研究发现，L1-L4等阶段性寿命性状与在群天数或生产寿命的遗传相关性较高(0.79~0.99)，这将有利于开展长寿性状的早期选择。

3.4 长寿性状遗传趋势分析

本研究对江苏地区中国荷斯坦牛在群天数、生产寿命、利用胎次和L-L4等7个长寿性状遗传趋势进行分析，结果表明在群寿命和生产寿命等性状的育种值随年份的增加整体呈下降的趋势，这可能与我国牧场建立的育种目标有关。长期以来，我国多数牧场主要追求产奶量、乳脂率和乳蛋白率等泌乳性状，兼顾体型外貌、肢蹄等性状，对长寿性状的重视不够或者未纳入育种计划，因此造成了长寿性状遗传进展缓慢甚至出现下降的趋势。此外，也有研究发现生产寿命与牛奶产量之间存在负相关关系^{[6, 38-39]}，而这种负相关可能与高产奶牛繁殖性能低有关。国外多项研究发现高产奶牛的繁殖能力偏低，导致淘汰提前或增加^[40-43]。因此，在追求产奶性能的同时，也应重视长寿、繁殖等功能性状的选育。

4 结论

本研究对江苏地区荷斯坦牛7个长寿性状的遗传参数和育种值进行评估，得出以下结论：1)在群寿命、生产寿命等长寿性状遗传力较低，介于0.07至0.15之间；2)在群寿命、生产寿命以及利用胎次之间存在较强的遗传相关，在群寿命、生产寿命、利用胎次与L1-L4等胎次寿命性状具有中等到较高的遗传相关和表型相关；3)近年来，江苏地区牧场长寿性状遗传进展缓慢，且整体呈下降趋势。本研究结果可为开展以江苏为代表的南方地区奶牛重要经济性状遗传参数评估以及制定长寿性状育种规划提供理论参考。