萨能奶山羊体尺性状的多元逐步回归分析与生长曲线拟合

何志成,秦晓晨,吕永强,李秀金,边会龙,罗 军,李 聪*

(1.西北农林科技大学动物科技学院 陕西省动物遗传育种与繁殖重点实验室,杨凌 712100;2.仲恺农业工程学院,广州 510225;3.陇县畜牧工作站,宝鸡 721200)

萨能奶山羊是世界著名的奶山羊品种之一,且作为父本改良地方品种效果显著,参与了关中奶山羊、崂山奶山羊等品种的育成,在我国奶山羊产业发展中占有重要地位。萨能奶山羊具有乳用特征明显、体格大、产奶量高、泌乳期长、繁殖力强、抗病力强、遗传性能稳定等特点,主要分布于我国陕西省、山东省和甘肃省等地[1],目前,全国各地饲养的奶山羊中,有半数以上、程度不同地含有萨能奶山羊的血统。但在奶山羊产业快速发展的同时,也面临着良种覆盖度低、良种繁育体系不健全、种羊鉴定和选育工作开展缓慢等诸多问题[2],尤其是标准化的生产性能测定方法、详细准确的系谱与生产资料记录等的缺乏与不完善,更是制约了奶山羊遗传评估等常规育种工作的实施。在奶牛上,国内外均已积累了大量的产奶、体型、繁殖、长寿性、乳房健康和产犊难易等生产性能测定数据,建立了良种繁育技术体系,定期开展群体遗传评估,各国均建立了奶牛综合遗传性能选择指数。在我国奶牛群体中,建立了中国奶牛性能指数(China Performance Index, CPI)和中国奶牛基因组选择性能指数(Genomic China Performance Index, GCPI),这为乳用种公羊遗传品质的客观评价和奶山羊选育方法的建立提供了重要参考。体型外貌是评定种羊的重要标准之一,公羊需要满足一定体重、体尺标准才能进行留种和选配,因此开展萨能奶山羊体尺性状的分析是选种、选配和群体遗传改良的重要前提基础。

产奶性状是奶畜的重要经济性状,而体尺性状与产奶性状具有一定的相关性。研究发现,中国荷斯坦牛的体重、体高、体长、胸围和管围与305 d产奶量有显著的相关性[3],中国荷斯坦牛主要体尺性状与产奶量存在较高的表型相关与遗传相关[4]。对崂山奶山羊的研究显示,其产奶量及乳成分性状与体尺性状存在显著的相关性[5],这表明通过对体尺性状的选择可间接选择产奶性状。此外,产奶性状的获得需要母畜经历配种、妊娠、分娩等过程,基于体尺性状与产奶性状存在的显著相关性,揭示了体尺性状在一定程度上可以预测产奶性状,进而实现对产奶性状的早期选择。

目前,国内关于萨能奶山羊体尺性状的相关、回归及主成分分析尚未有报道,因此,本研究对萨能奶山羊11个体尺性状进行测定,通过主成分分析、多元逐步回归分析和模型拟合,旨在探究不同月龄萨能奶山羊不同体尺性状间的相关与回归关系,并构建其生长曲线。以期为奶山羊生产性能测定和遗传评估等常规育种提供理论基础,为规模化奶山羊场选种、选配和生产管理提供实际指导。

1 材料与方法

1.1 数据测定与整理

本试验数据来源于陕西省宝鸡市陕西宝泰良种繁育有限公司的奶山羊育种场,对场内980只不同月龄萨能奶山羊个体进行单次体尺测量,共测定11个体尺性状,将其分为一般体尺性状(体高、体斜长、胸围、管围)、反映公羊繁殖性能的体尺性状(尻斜长、尻宽、睾丸周径)和母羊乳房性状(乳房围度、乳房深度、乳头间距、乳房距地面高度)3类。由于青年母羊乳房性状难以测定,因此在数据量上相较于一般体尺性状少。剔除月龄信息缺失、测定值缺失和错误的个体记录后,最终以上3类性状分别有874、219和150只个体的数据保留,其中874条一般体尺性状记录包含219条公羊记录和655条母羊记录,各体尺性状测定条数如表1所示。

表1 萨能奶山羊体尺性状数据条数Table 1 Number of body measurement traits data for Saanen dairy goat

11个体尺性状的测定方法如下,体高:鬐甲最高点至地面的垂直距离;体斜长:肩端前缘至坐骨端后缘的距离;胸围:肩胛骨后缘垂直体躯的周径;管围:左前肢管部下三分之一最细处的水平周径;尻斜长:腰角前缘至臀端后缘的直线距离;尻宽:髋部两侧外缘间的最大距离;睾丸周径:睾丸最宽处的周径;乳房围度:乳房最大处的周径;乳房深度:乳头基部至乳房基部的垂直距离;乳头间距:两乳头下端的直线距离;乳房距地面高度:乳头基部至地面的垂直距离。

1.2 数据分析

1.2.1 主成分分析 主成分分析常用于大样本和多变量的数据,方法基于原始变量间的相关性,在保留原始数据最大方差的情况下,降低原始变量数,达到减少变量数与简化分析的目的,探究降维后新主成分的含义可进一步了解数据的结构及其主要变异。

采用Excel 2021对数据初步整理后,分别对表1中一般体尺性状、反映公羊繁殖性能性状和母羊乳房性状等数据进行主成分分析。采用SPSS 26.0统计软件,选择相关性矩阵对数据进行分解,保留累计方差解释大于80%的主成分,计算各主成分的因子载荷,进而探究各个主成分的含义。

1.2.2 多元逐步回归分析 利用一般体尺性状分别对反映公羊繁殖性能的体尺性状和母羊乳房性状进行多元逐步回归分析,采用数据集为反映公羊繁殖性能的体尺219只、乳房性状的体尺150只和同时拥有以上两部分个体的一般体尺性状数据。采用R语言Lm和step函数,程序的运行步骤为:1)从回归模型中依次剔除1个体尺性状后计算AIC值;2)对比原模型和所有剔除单个体尺性状模型的AIC值,选择最小AIC值的模型;3)不断重复步骤“1)”、“2)”,直到从模型中剔除任意单个体尺性状均不减小AIC值为止,此时的模型为最优模型,进而得到一般体尺性状对反映公羊繁殖性能体尺性状和乳房性状的回归模型参数。

1.2.3 一般体尺性状拟合分析 采用SPSS 26.0软件的非线性回归和Bertalanffy模型拟合不同月龄公、母羊的平均一般体尺性状生长曲线,通过非线性最小二乘法得到参数估计和拟合度值,其中Bertalanffy模型公式为:Y=A×(1-B×e-kt)3,并且可计算指标:拐点体尺8A/27,拐点月龄ln(3B)/k,公式中A、B、k为模型参数,A代表极限生长量、B代表常数尺度、k代表瞬时相对生长率。

2 结 果

2.1 体尺性状描述性统计

对萨能奶山羊体尺性状的描述性统计分析如表2所示,各性状取值均在合理范围内。

表2 体尺性状描述性统计分析Table 2 Descriptive statistical analysis of body measurement traits cm

2.2 体尺性状的相关性及主成分分析

对于体尺性状间的相关性分析分别在3类体尺性状内独立进行,结果如表3所示,可见3类体尺性状内的任意2个性状间均具有极显著相关关系(P<0.01),其中一般体尺性状间的相关系数均大于0.79,反映公羊繁殖性能的性状间相关系数均大于0.8,母羊乳房性状间的相关系数处于-0.492～0.730之间。

表3 体尺性状指标间的相关性分析Table 3 Correlation analysis between body measurement traits

体尺性状主成分分析结果如表4所示,除了母羊乳房性状降维成2个PC(主成分)外,一般体尺性状和反映公羊繁殖性能的体尺性状均只降维成1个新PC。一般体尺性状中,体高、体斜长、胸围和管围4个体尺性状降维成1个新的主成分变量,解释方差达到88.64%,4个原始体尺性状在新PC1变量的载荷系数均为正,且均在0.90以上。母羊乳房性状中,乳房围度、乳房深度、乳头间距和乳房距地高度4个性状降维成2个新的主成分变量,2个新的PC累计解释方差达到85.43%,其中乳房围度、乳房深度和乳头间距在PC1上表现为较高的正载荷,系数均在0.80以上;乳房距地高度在PC2上表现为较高的正载荷,系数为0.784。反映公羊繁殖性能的体尺性状中,尻斜长、尻宽和睾丸周径3个性状降维成1个新的主成分变量,解释方差达到88.89%,3个原始性状均在新PC上表现为较高的正载荷,系数均在0.90以上。

表4 三类体尺性状各主成分的因子载荷、特征值及解释方差Table 4 Factor loadings, eigenvalues, and explain variances of each principal component for three types of body measurement traits

2.3 一般体尺性状对其它体尺性状的多元逐步回归分析

基于最优模型原则,发现一般体尺性状对其它两类体尺性状均存在显著的回归关系,可构建一般体尺性状对反映公羊繁殖性能体尺性状和乳房性状的回归方程。一般体尺性状对其它两类体尺性状的多元逐步回归分析结果见表5,其中,母羊乳房性状的回归方程R2相对低,只有乳房深度的R2达到0.40,其余均在0.20以下;反映公羊繁殖性能的体尺性状R2较高,睾丸周径的R2为0.76,尻斜长和尻宽的R2均在0.80以上。

表5 一般体尺性状对其它体尺性状的多元逐步回归分析结果Table 5 Multiple stepwise regression analysis results of general body measurement traits compared to other body measurement traits

2.4 一般体尺性状的生长曲线

在公羊和母羊群体中对体高、体斜长、胸围和管围4个一般体尺性状进行拟合,拟合方程及拟合度如表6所示,除母羊管围的拟合度R2较小,其余性状的拟合度R2均大于0.75。公、母羊4个一般体尺性状的拐点月龄均小于0,表明拐点月龄和拐点体尺在本试验中没有意义,奶山羊一般体尺的增长随着月龄的增长而缓慢减小。对于不同月龄公母奶山羊体高(图1A)、体斜长(图1B)、胸围(图1C)和管围(图1D)等4个一般体尺性状的生长曲线进行拟合,发现4条曲线的形状类似,均呈现先较快上升,后逐渐趋于平缓的趋势。公羊与母羊在出生体尺上相似,随着月龄的增加,公羊体尺的增速明显高于母羊,相较于其它3个体尺性状,管围在公母羊之间的差异更大(图1D)。

图1 不同月龄公、母奶山羊一般体尺性状的生长曲线拟合Fig.1 Growth curve fitting of general body measurement traits of bucks and does with different months of age

表6 一般体尺性状的拟合方程及拟合度Table 6 Fitting equation and fit degree of general body measurement traits

3 讨 论

3.1 体尺性状的分类

研究表明,奶畜的一般体尺性状反映了个体的品种特征和育种价值;乳房性状主要与产奶性状和乳房炎相关,乳房周长和乳房深度与总产奶量相关[6],而乳头间距直接影响奶牛产奶量[7];尻部性状与产犊能力呈负相关[8],与生产寿命相关,是澳洲荷斯坦公牛的选育方向之一[9];荷斯坦公牛阴囊周径与精液密度、射精量和精子活率呈较高的正相关[10]。基于此,本研究将奶山羊的体尺性状分为了三类,包括体高、体斜长、胸围和管围等的一般体尺性状;尻高、尻斜长和睾丸周径等反映公羊繁殖性能的体尺性状;乳房围度、乳头间距、乳房深度和乳房距地高度等母羊乳房性状。因此,科学制定奶山羊体尺性状的测量标准有利于后续分析的进行并且能有效反映个体的泌乳、繁殖与生产潜力。

3.2 主成分分析

主成分分析是一种分解多表型变量的方法,也可用于分解多性状遗传协方差矩阵[11],分析得到的新主成分可作为畜禽育种的综合指标进行选择和管理。性状或变量间存在线性相关是主成分分析的前提,且高度线性相关的变量可通过主成分分析降维到较少的新变量中,因此在主成分分析前应检验变量间的相关性。本研究揭示:体尺性状适合进行主成分分析,且一般体尺性状和反映公羊繁殖性能的体尺性状的性状内相关系数均大于0.79,正是这两类性状均能降维成一个PC且保持解释方差大于88%的原因。

本研究中,一般体尺性状与反映公羊繁殖性能的体尺性状的主成分分析均发现,PC1在原始变量的因子载荷均表现为较高的正载荷,表明在PC1中各原始变量间存在较高的正相关,这与相关性分析结果一致,因此这两类性状的PC1可以分别解释为羊的体型外貌特征和公羊的整体繁殖外貌特征。而母羊乳房性状中PC1主要受到乳房围度、乳房深度和乳头间距的高正载荷影响,可以解释为母羊乳房的整体外貌特征,PC2主要受到乳房距地面高度的正高载荷影响,可以解释为与后躯高相关的特征,且与乳房整体的外貌无关。在对德南牛体高、体斜长、胸围、管围和十字部高进行主成分分析研究中,把第一主成分概括为体型外貌因子[12]。对新疆褐牛17个体尺性状主成分分析,概括PC1为整体结构特征、PC2为后躯特征、PC3为前躯特征、PC4为尻部特征、PC5为乳房特征[13]。对霍尔巴羊体尺性状进行主成分分析和方差最大化旋转后,得到第一主成分可解释为与羊只体型大小和重量相关的指标[14]。上述研究结果与本研究一般体尺性状的降维结果一致,表明在分析时可以对多变量的一般体尺性状和与公羊繁殖性能相关的体尺性状的数据进行变量简化,本研究中这两类性状的原始变量较少,为得到更准确的结果,后续可考虑增加样本量和原始变量数。有研究对中国荷斯坦奶牛乳房性状进行主成分分析,得到第一主成分与乳房综合表型相关[15],这与本研究母羊乳房性状中PC1的结果一致;而对奶牛乳房深度、乳头长度、悬韧带、后房宽度、后房高度、前房高度、前房附着进行主成分分析,得到3个主成分分别解释为乳房容积因子、乳房紧凑因子和乳房适宜因子[16]。对乳头位置、乳房深度、乳房悬垂程度和左右乳房分离度进行主成分分析,得到的PC1在前3个原始表型中高度相关,且由于其重复力在不同泌乳期较高,因此可用于早期遗传评估[17]。主成分分析的结果与利用不同的数据库和原始变量相关,这是本研究同以上两项研究不一致的原因。

各因子载荷除以对应PC的特征值开根号得到变换系数,再与标准化后的原始变量相乘可以得到每个PC的得分。在一般体尺性状、与公羊繁殖性能相关体尺性状和母羊乳房性状中PC1得分高分别代表羊外貌体型、繁殖体型和乳房外型均大,而母羊乳房性状中PC2得分高代表母羊的后躯高。利用聚类方法和各主成分得分对数据样本进行分类,可以得到数据在新变量上的分布,Arnal等[18]通过对泌乳曲线进行主成分分析和聚类,得到了关于泌乳曲线的独立新变量和泌乳曲线分类的多样性结果。

3.3 多元逐步回归分析

在奶畜体尺性状的多元回归分析中,常把精液品质和产奶性状作为预测的Y变量,以实现对难以测量的经济或繁殖性状的预测。有研究通过体重、体况评分、皮下脂肪厚度、体长、腰围和臀宽建立了牛一岁平均日增重(R2>0.5)、断奶至一岁平均日增重(R2<0.15)和一岁体况评分(0.15

本研究发现,萨能奶山羊一般体尺性状对与公羊繁殖性能相关的体尺性状建立了高拟合度的回归方程,表明前者对后者有较好的预测效果。在新疆褐牛中,成功构建了公牛体高、体直长、体斜长、胸围、腹围、管围和估测体重关于睾丸周径的多元回归方程(R2=0.876)[21]。本研究发现,母羊乳房性状对于一般体尺性状有一定的独立性,受到月龄和其它因素的影响。有研究发现山羊在妊娠期乳腺容积呈指数级增长[22],并且明显的乳房发育发生在妊娠期的后三分之一。第2、3胎的萨能奶山羊相比1胎羊乳房下垂,乳头更粗、长,乳头离地间隙短[23]。奶牛乳房深度、乳房长度、悬韧带及前房附着与胎次呈极显著相关关系[24]。因此,通过构建一般体尺性状对与公羊繁殖性能相关的体尺性状的多元回归方程,可为生产上预测繁殖性状提供重要参考。

3.4 一般体尺性状生长曲线的拟合

研究表明,Bertalanffy模型能较好拟合安格斯牛[25]、藏绵羊[26]和哈萨克羊[27]的生长曲线,且模型拟合度均大于0.95。基于以上研究,本研究采用Bertalanffy模型对萨能奶山羊的一般体尺性状进行拟合分析,计算得到1.5岁公羊体高、体斜长和胸围的预测值分别为80.21、85.47和95.48 cm,母羊分别为73.48、79.65和89.11 cm。结果与已有的发现中,1.5岁萨能奶山羊公羊(体高:80.2 cm;体斜长:87.3 cm;胸围:92.8 cm)和母羊(体高:69.6 cm;体斜长:80.2 cm;胸围:86.3 cm)的体尺数据[28]相似,表明本研究拟合的曲线符合萨能奶山羊品种的生长规律。

本研究中体尺性状的测定结果在月龄上分布不均匀,有月龄测定数据较少甚至仅有1条,模型的构建受到了该单个样本点影响,导致方程拟合度低于0.80。而管围性状测定数据在同一月龄下的变异程度更大,导致管围性状的曲线拟合度仅为0.60。因此,为了得到精准的群体曲线形状和拟合度,需要考虑更密集的测定月龄、更小的测定误差、更适用的拟合模型和更大的样本量。

据报道,采用Bertalanffy模型对青海藏羊体长、体高、胸围进行拟合分析,发现这3个体尺性状拐点在0.23～0.27月龄[29];采用logistic模型拟合豫东肉山羊体高和体长,发现拐点在20～50 d[30];采用三次四项式拟合波尔山羊的体高、体长和胸围,发现体尺性状的生长拐点在7～14 d[31]。以上研究表明,不同品种羊体尺性状的拐点主要在羔羊生长的早期,揭示生产中保障羔羊早期营养需求和精细管理的重要性。本研究与上述研究不一致表明了萨能奶山羊品种特有的生长规律,有研究发现萨能奶山羊的出生体高、体斜长、胸围分别占成年时的47%、37%和35%[32],表明萨能奶山羊体尺在生长早期高度发育,但这也可能由于分析模型的不同所导致。采用Gompertz模型对绒山羊体重进行拟合并估计模型参数、拐点年龄和拐点体重的遗传参数,发现三者的遗传力均低于0.15,但存在较高的遗传相关[33]。采用Brody模型对绵羊体重进行拟合并估计模型参数的遗传力为0.39,具有较高的遗传力[34-35]。以上研究均表明通过模型参数计算的拐点指标可以通过育种进行间接选择,进一步明确模型参数的意义,对于探究并选择参数作为新的育种目标是可行的方案。

4 结 论

本研究发现萨能奶山羊一般体尺性状、与公羊繁殖性能相关的体尺性状和母羊乳房性状3类性状内部间均存在显著的相关性,3类性状均可被新的主成分变量解释;一般体尺性状对与公羊繁殖性能相关的体尺性状有良好的预测效果;Bertalanffy模型可较好拟合一般体尺性状,并反映萨能奶山羊的品种特征。研究结果为生产上萨能奶山羊体尺规律的探究和体尺测定的简化提供参考,为萨能奶山羊体尺性状的遗传评估和选种选育提供数据基础。