新疆褐牛是我国五大兼用牛品种之一，舍饲条件下，平均305天产奶量可达4 704 kg，其平均乳脂率达3.9%，乳蛋白率为3.5%^[1]，高于新疆地区荷斯坦奶牛乳脂率和乳蛋白率的平均水平^[2-3]，同时，因其耐寒、耐粗饲和抗逆性强等特点而深受欢迎，是新疆北疆广大农牧区黄牛改良中一直主推的品种。牛体型性状与牛健康水平^[4]、使用年限^[5]、繁殖能力^[6]和生产能力^[6]息息相关，也是目前牛群选配、新品种培育的主要信息。自1983年《新疆褐牛》品种标准发布以来，每年定期开展了新疆褐牛体型鉴定工作，然而，该方法采用的是传统的百分制(1~100分)进行评分，主观性较强。体型线性评定技术是目前国内外使用较为普遍的鉴定方法，该方法对评定的体型性状用数字化的线性尺度，来衡量其从一个极端到另一个极端的不同状态，具有客观统一的标准^[7]。2018年“新疆褐牛体型外貌线性评定方法的研究”立项，并在新疆伊犁地区、塔城地区、阿勒泰地区和乌鲁木齐市等地开展了新疆褐牛体型基础数据测量与收集，为新疆褐牛体型外貌线性评定技术规程的建立提供了数据支撑。

目前，多元线性回归被广泛应用于分析各体型性状之间的关系，研究结果表明，各体型性状之间存在一定的表型和遗传相关性。例如，王丹等^[8]研究发现，新疆褐牛体躯结构性状之间存在较强的相关性；刘建明等^[9]报道，新疆褐牛尻部性状之间具有一定的相关性；在一些其它品种牛当中也有类似的报道，比如，彭朋等^[10]发现中国荷斯坦牛后乳房高度和后乳房宽度之间具有较强的遗传相关。各体型性状间的相关性表明，性状之间存在冗余信息，避免分析大量高度相关的性状是育种学者长期以来所关心的问题。主成分分析和因子分析可以通过降低数据的维度，以最小的信息损失探索性状之间的关系^[11]。这两种方法旨在以原始变量的线性组合转化为综合变量，其中变量权重从原始数据的相关矩阵中获得^[12]。主成分解释了原始变量的全部方差，无方差损失；而因子分析除了有共性因子外，还有特殊因子，共性因子只解释了部分信息，有方差损失^[13]。鉴于这两种统计方法的重要性以及避免各体型性状之间的冗余信息，本研究旨在利用因子和主成分分析方法估计新疆褐牛体型性状的遗传参数，并进行比较，该研究结果可为新疆褐牛体型性状的遗传评估以及育种计划的制定提供理论依据。

1 材料与方法 1.1 数据来源

本试验于2018年7月至2022年7月在新疆燕犇褐牛繁育发展有限责任公司、塔城农牧科技有限公司、伊犁新褐种牛场、新疆维吾尔自治区地方国营伊犁种马场核心牛队、伊犁哈萨克自治州昭苏马场褐牛繁育中心、兰明家庭农场、尼勒克县喀拉托别乡新牧源牛羊养殖专业合作社、尼勒克县旭康农牧业有限公司、新源县牧骏农牧业发展有限公司9个新疆褐牛核心育种场进行测量和数据收集。其中新疆燕犇褐牛繁育发展有限责任公司、塔城农牧科技有限公司和伊犁新褐种牛场为全舍饲饲养，饲草料稳定，采用TMR饲喂，自由饮水、机器挤奶，每月定期进行DHI测定。剩余牛场均为半舍饲半放牧饲养，通常6~11月在草场进行自由采食，12月至次年5月进行舍饲，均无DHI测定记录。

1.2 数据处理

本试验共测定了27个体型性状，其中有17个体型性状为测量性状，10个体型性状为评分性状(1~9分)，具体评定部位如图 1所示。共鉴定新疆褐牛泌乳牛1 816头；根据“平均值±3倍标准差”的原则剔除非正常记录，并剔除无系谱记录牛只，最终保留1 016头新疆褐牛泌乳牛体型数据进行进一步分析。其中系谱记录由186头公牛和2 067头母牛构成。

1.3 主成分分析

运用SAS软件的PRINCOMP程序进行主成分分析，该方法同样采用降维的思想将原始变量综合为一组维度少的新变量，将信息的损失降到最低，并探讨性状之间的关系。Kaiser^[14]研究表明，入选主成分的标准是选择特征值≥1的主成分。

主成分的模型表达式如下所示：

$ P C_m=\alpha_{1 m} X_1+\alpha_{2 m} X_2+\cdots \alpha_{n m} X_n $

式中，α_nm是观察变量的协方差矩阵的特征向量的元素，与第m个特征值相对应；X_n是第n个原始变量。

1.4 因子分析

运用SAS软件的FACTOR程序进行因子分析，该分析方法通过寻求一组新的变量(X₁, …, X_n)，称为共同的潜在因素，来综合包含在一组n个观察变量(y₁, …, y_n)中的信息。使用Kaiser^[14]所描述的方差最大旋转法来保持提取因子的正交性，仅保留特征值≥1的因子进行分析。通过观察各体型性状在各因子中的载荷，从生物学的角度解释分析结果。最后用标准化的评分系数计算每头牛的因子得分，因子得分有助于之后的分析中将公共因子当作分析变量，而不再分析降维前的原始变量。

根据Morrison^[15]的研究结果，因子分析的模型表达式如下所示：

$ y_n=b_{n 1} X_1+b_{n 2} X_2+\cdots b_{n p} X_p+e_n $

式中，y_n是第n个原始变量；b_np是每个变量n在各因子上的载荷；X_p是第n个变量的第p个公共因子；e_n反映第n个变量的特殊因子。

1.5 遗传参数估计

利用BLUPF90软件的AIREMLF90程序，采用平均信息约束最大似然法(average information restricted maximum likelihood，AI-REML)配合动物模型估计27个体型性状、各主成分和各因子的方差组分。模型如下所示：

$ Y=\boldsymbol{X} \beta+\boldsymbol{Z} \alpha+e \;\;\;\;\;\;\;模型1 $

式中，Y表示27个体型性状的观察值向量；β为固定效应向量，包括场效应(9个水平)、泌乳天数效应(9个水平：10~40天、41~80天、81~120天、121~160天、161~200天、201~240天、241~280天、281~320天、＞320天)、胎次效应(4个水平：1胎、2胎、3胎、≥4胎)；α为加性遗传效应向量；e为随机残差；X和Z分别为固定效应向量和随机效应向量的关联矩阵。

此外，构建双性状动物模型，将主成分和因子得分作为Y，分别估计各主成分和各因子得分之间的遗传相关。双性状动物模型的矩阵形式如下：

$ \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} g_1 \\ g_2 \end{array}\right]+\left[\begin{array}{l} e_1 \\ e_2 \end{array}\right]\;\;\;\;\;\;\;\; 模型2 $

式中，y₁和y₂分别为各因子得分和各主成分的观察值；b₁和b₂为性状1和性状2的固定效应向量，固定效应与模型1一致；g₁和g₂为性状1和性状2的加性遗传效应向量；e₁和e₂为性状1和性状2的随机残差；X₁、X₂、Z₁和Z₂分别为相应的关联矩阵。

将计算得到的方差组分估计值代入以下公式计算遗传参数^[16]：

$ \text { 遗传力 }: h^2=\frac{\sigma_a^2}{\sigma_a^2+\sigma_e^2} \text {; 遗传相关: } r_A=\frac{\operatorname{Cov}\left(a_i, a_j\right)}{\sqrt{\sigma_i^2 \sigma_j^2}} $

式中，σ_a²为加性遗传方差；σ_e²为残差方差；Cov(a_i, a_j)为第i个性状和第j个性状之间的加性效应协方差。

2 结果 2.1 遗传力估计

由表 1可知，新疆褐牛体型性状的遗传力估计值范围在0.09~0.49。其中，十字部高、后腿高度和尻长属于高遗传力性状(h²＞0.35)；蹄踵深度、蹄角度和乳房平衡性为低遗传性状(h²＜0.10)；其余体型性状均为中高遗传力性状(0.35＞h²＞0.1)。

2.2 主成分分析

表 2显示了由体型性状的各主成分所解释的特征值、方差比例和累计方差占比。有8个主成分的特征值≥1，占总方差的57.69%。第一主成分(PC1)在27个体型性状中占总方差的比例最大(15.55%)。表 3为提取的主成分在各体型性状中的载荷系数，由于较低的载荷系数对变量的影响不大，因此仅显示了载荷系数≥|0.4|的结果^[17]。PC1在十字部高(0.61)、体深(0.73)、胸宽(0.54)、后腿半围(0.63)、尻长(0.51)、后乳房宽度(0.61)、前乳区长度(0.58)和后乳区长度(0.51)上的载荷系数较高，PC1反映了体躯结构、肌肉度、尻部和泌乳系统的部分信息。PC2的方差占比为10.81%，在中央悬韧带(0.74)、乳房深度(0.71)、前乳区长度(0.68)、前乳头长度(0.57)和前乳头直径(0.48)上的载荷系数较高，主要反映了泌乳系统的信息。PC3占总方差的6.99%，在十字部高(0.47)、乳房深度(0.45)和前乳头直径(0.44)的载荷系数较高，主要反映了泌乳系统的信息。PC4的方差占比为6.31%，在尻角度(0.46)和前乳区长度(0.49)的载荷系数较高，主要反映了尻部和泌乳系统的信息。PC5至PC8的方差占比分别为5.61%、4.73%、3.97%和3.72%，主要反映了肢蹄、肌肉度和泌乳系统的信息。由于剩余特征值＜1的主成分所解释的方差比例较小，因此不考虑对它们进行进一步分析。

2.3 因子分析

表 2列出了每个因子的特征值和所解释表型变异的比例，8个潜在因子的特征值＞1。与主成分相似，前8个潜在因子解释了57.67%的信息，但在每个因子的特征值、方差比例和载荷系数上存在差异。第一因子(F1)占总方差的比例最大(10.53%)，使用最大方差旋转法得到的因子载荷系数如表 4所示，仅显示因子载荷系数≥ 0.4的结果^[17]。在F1中，十字部高(0.81)、体深(0.66)、后腿半围(0.62)和尻长(0.64)具有较高的载荷系数，F1主要反映了体躯结构、肌肉度和尻部的信息，与主成分的结果类似。F2解释了10.26%的方差，在后腿高度(0.53)、蹄踵深度(0.75)、蹄角度(0.70)和中央悬韧带(0.68)有较高的载荷系数，F2主要反映了肌肉度、肢蹄和泌乳系统的信息。F3的方差占比为8.83%，在乳房深度(0.68)、前乳区长度(0.71)和后乳区长度(0.75)有较高的载荷系数，主要反映了乳房大小的信息。F4的载荷系数在前乳头位置(0.82)和后乳头位置(0.86)上较高，反映了乳头位置的信息。F5集中反映了前胸宽窄的信息，胸宽和髻甲上部宽的载荷系数分别为0.52和0.71。在F6中，前乳头长度(0.72)和前乳头直径(0.63)具有较高的载荷系数，主要反映了乳头大小的信息。F7和F8分别解释了5.01%和4.63%的方差，主要反映了肌肉度、尻部和泌乳系统的信息。各体型性状公因子方差的范围在0.38~0.78，与其它较低公因子方差的体型性状相比，公因子方差越高，说明该变量能被因子解释的程度越高。

2.4 主成分和因子的遗传力和遗传相关估计

8个主成分和因子得分的遗传力和遗传相关如图 2和图 3所示。两种方法的遗传力估计值具有类似的结果，但主成分的平均遗传力估计值(0.26)略高于因子得分(0.23)。基于主成分分析估计遗传相关的范围在-0.76(PC3-PC6)~0.74(PC3-PC4)，因子分析估计遗传相关的范围在-0.75(F3-F6)~0.86(F1-F6)，由结果可知，第三主成分和第六主成分及第三因子和第六因子均呈现强的负遗传相关。

3 讨论 3.1 遗传力

本研究发现，新疆褐牛体躯结构性状(十字部高、体深和胸宽)均属于中高遗传力性状，但遗传力估计值略低于Rendena牛^[18]和德系西门塔尔牛^[19]的遗传力估计结果，与瓦尔多斯塔纳牛^[20]的研究结果相似。肌肉度性状(髻甲上部宽、后腿半围、后腿高度和中部肋骨)均属于中高遗传力性状，遗传力估计值高于新疆乳肉兼用型西门塔尔牛^[21]肌肉度性状的遗传力，与德系西门塔尔牛^[19]研究结果类似。新疆褐牛尻部性状(尻长、尻宽和尻角度)的遗传力同样均属于中高遗传力性状，均高于德系弗莱维赫牛^[19]和山东省荷斯坦奶牛^[22]的遗传力估计结果，与意大利褐牛^[23]和Rendena牛^[18]研究结果相似。肢蹄性状(蹄踵深度、蹄角度、后肢侧视、飞节质地和后肢后视)中，蹄角度和蹄踵深度属于低遗传力性状，国内外多项研究^[24-25]报道了大多数肢蹄性状均属于低遗传力性状。乳房性状中除乳房平衡性外，其余乳房性状均属于中高遗传力性状，先前有研究^[26]报道了乳房平衡性属于高遗传力性状，但该研究对表型数据的定义是测量前乳头和后乳头离地面的垂直距离的平均值，与本研究乳房平衡性评分表型存在差异。通常实际测量值属于数量性状，而评分性状属于分类性状，因此体型表型数据的收集方式对于遗传力的估计起到至关重要的作用，此外，体型鉴定员和系谱的完整性也是影响遗传力估计结果的主要因素。

3.2 主成分和因子分析

本研究主成分分析解释的总方差高于Gutiérrez和Goyache^[27]报道的49.7%(10个性状选取前2个主成分)，低于刘建明等^[9]报道的62.91%(17个性状选取前5个主成分)，然而刘建明等^[9]分析的体型性状中未考虑肢蹄性状和大部分乳房性状，本研究所分析的性状多于前人对新疆褐牛体型研究的报道，在一定程度上增加了数据的维度。对于因子分析，Mazza等^[28]在Rendena牛和奥斯塔红皮牛的研究中，选择了6个特征值大于1的因子，分别占两个兼用牛品种总方差的63%(20个性状)和58%(22个性状)，Olasege等^[29]对中国荷斯坦牛23个体型性状的研究中选择7个特征值＞1的因子，占总方差的60.37%。上述研究中的总方差均高于本研究的总方差，然而，本研究的结果高于Chu和Shi^[30]的报道，其选择了2个因子(15个性状)，占总方差的49.1%。所选取的主成分和因子的数量及累计方差的大小，可能是由于不同的品种和体型性状的数量以及样本量所造成的。此外，提取的主成分和潜在因子所解释的方差比例不同，可能是由于因子分析估计了原始变量的共性方差，而主成分分析并没有将原始变量的方差进行剖分，将公因子方差的估计值设定为1^[31]。主成分分析的方法是把给定的一组相关变量通过线性变换转换成另一组不相关的变量，适用于各体型性状权重的赋予；因子分析的方法则是对因子载荷矩阵进行正交旋转，更有利于实际意义的解释。

3.3 主成分和因子的遗传力和遗传相关

一般体躯深而高，乳房宽而大的牛只在产奶性能方面也会有较好的表现^[32-34]，PCA1和F1载荷系数较高的体型性状通常与牛只的产奶性能相关，这与Kern等^[34]的研究结果相似，因此，育种目标旨在提升产奶性能的情况下，可以将PC1或F1纳入选择指数。Manafiazar等^[35]研究表明，十字部高、胸宽和后乳房宽度与剩余采食量有较强的遗传相关，而剩余采食量是衡量饲料效率的指标。本研究发现，PC1和F1均属于高遗传力性状，因此，通过表型选择将会获得较快的遗传进展。

产奶量较高的个体通常会面临乳房疾病和肢蹄疾病^[36-37]，从而造成牛只过早的被动淘汰，F2和PC5中载荷系数较高的是肢蹄性状，因此，选择F2和PC5得分较高间接意味着选择了蹄角度较大、蹄踵较深、肢蹄较直、飞节质地较为细致的牛只，一定程度上将会降低肢蹄疾病的发生率。给牧场造成最大的经济损失和最常见的疾病之一就是临床乳房炎^[38]，而体细胞数目前被广泛应用于指示牛只是否患有乳房炎^[39]，多项研究表明^[40-42]，体细胞评分与前乳房附着、前乳头位置、乳房深度有较强的正遗传相关，本研究中，PC2、F3、F4和F8载荷系数较高的体型性状与体细胞数相关，且从因子分析的结果来看，F3反映了乳房大小的信息，F4反映了乳头位置的信息，F8则反映了乳房形态的信息，且F3与F8之间存在较强的负遗传相关。DeGroot等^[43]研究表明，附着性较强的乳房以及较近的乳头位置有利于降低体细胞数；此外，Dube等^[44]研究也表明，较浅的乳房和较近的乳头位置与较低的体细胞数有强的相关性。此外，PC3与PC6以及F3与F6均呈现强的负遗传相关，主成分和因子载荷系数表明，以上主成分和因子均在乳房深度、前乳头长度和前乳头直径处有较高的载荷系数。Xue等^[45]对中国荷斯坦牛的研究表明，乳房深度与前乳头长度呈现强的负遗传相关，与本研究结果一致；同时Sinha等^[46]的研究指出乳房越深，乳头越长的牛只更容易感染乳房炎。

尻长、尻宽和尻角度一般与牛只的繁殖性能相关，例如，Eaglen等^[47]研究表明尻宽和尻角度与产犊难易度存在中等的正遗传相关，尻宽与母体的妊娠期长度存在强的正遗传相关，Wall等^[48]研究表明过窄的臀部和过低的尻角度通常被认为是不利于产犊的表现。本研究中，PC1、PC4和F7均解释了尻部的信息，且PC1和PC4存在中等的正遗传相关，Olasege等^[29]的研究结果显示PC1与PC4存在中等的正遗传相关，与本研究结果一致。因此，在制定选择指数时，可以考虑将PC1、PC4和F7纳入育种目标性状当中，从而达到改善群体繁殖性能的目的。此外，大多数繁殖性状属于低遗传力性状，通过表型选择，其遗传进展缓慢，而大部分体型性状、主成分和因子表型属于中高遗传力性状，通过体型或主成分和因子表型进行间接选择也是一种提高牛只繁殖力的方法。

4 结论

本研究发现，新疆褐牛体型性状中除蹄踵深度、蹄角度和乳房平衡性外，其它性状均属于中高遗传力性状。各主成分和潜在因子的遗传力除F3外，均属于中高遗传力，从而表明将其纳入选择指数的可行性。在体型性状的评估中可以考虑使用主成分和潜在因子，将多维度的体型性状简化为新的变量，这种方法可以避免分析大量相关性较高的性状，从而降低分析大量数据而造成的计算负担。