我国荷斯坦青年公牛基因组选择效果分析

苏丁然，彭 朋，闫青霞，陈绍祜，张胜利，李 姣，刘丑生，孙东晓*

(1. 中国农业大学动物科技学院，北京 100193； 2. 中国奶业协会，北京 100192； 3. 全国畜牧总站，北京 100125)

基因组选择是通过覆盖全基因组的高密度单核苷酸多态性(single nucleotide polymorphism, SNP)标记，并结合大量的后裔测定遗传评估结果来计算个体的基因组估计育种值(genomic estimated breeding value, GEBV)[1-2]。2001年，Meuwissen等[3]首次提出基因组选择(genomic selection, GS)的概念。2006年，Schaeffer[4]基于加拿大荷斯坦牛群体的研究结果表明，通过基因组选择可以节约92%的育种成本。畜禽全基因组高密度 SNP 芯片技术逐渐发展成熟，2007年首款商业化牛50K全基因组SNP芯片问世，由美国Illumina公司研发[5]，由此，基因组选择技术在奶牛育种工作中开始得到应用，并逐步成为重要的奶牛育种策略[6-8]，其中，美国荷斯坦协会于2009年1月首次发布荷斯坦青年公牛的基因组遗传评定结果[9]。2012年，中国农业大学奶牛育种团队及北京奶牛中心、北京首农畜牧、上海奶牛中心、全国畜牧总站和中国奶业协会等单位合作，成功建立了中国荷斯坦牛基因组选择技术平台，被农业部指定为我国荷斯坦青年公牛遗传评定的唯一平台，并开始在全国推广应用[10]。在基因组选择提出之前，后裔测定是评定奶牛种用价值最可靠的方法，但后裔测定存在世代间隔长，育种成本高等缺点。基因组选择相比后裔测定有以下优势：1) 能够捕获与目标性状有关的所有的遗传变异，具有较高的选择准确性(>70%)，高于传统系谱指数估计的可靠性[11-18]，且利用生物学先验信息可以更好地提高基因组选择的可靠性[19-21]。2) 可以不依赖表型信息进行选择[22]。3) 可以大幅度缩短世代间隔，加快遗传进展，基因组选择使公牛的世代间隔从7年下降到小于2.5年[23-24]，母牛的世代间隔从4年 左右下降到接近2.5年[25-27]；产奶性状的年遗传进展从50%～100%增加到300%，低遗传力性状的年遗传进展从3倍增加到4倍[28-31]，加快了奶牛繁殖性状、健康性状等低遗传力性状的遗传进展[32-33]。4) 大大降低了育种成本[34-36]。

我国自2012年开始，在全国所有公牛站实施荷斯坦青年公牛基因组遗传评估，截至2019年12月已有3 031头青年公牛参与基因组遗传评估，其中1 686头 公牛目前已获得了女儿至少一个完整泌乳期的生产性能测定(dairy herd improvement，DHI)和体型鉴定数据。基于此，本研究旨在利用我国荷斯坦青年公牛基因组遗传评估结果及女儿表型值对基因组选择的效果进行分析，验证基因组遗传评估的准确性，为进一步完善我国荷斯坦青年公牛基因组选择技术体系及我国奶牛群体遗传改良计划工作提供理论依据。

1 材料与方法

1.1 数据来源与整理

1.1.1 数据来源 基于中国农业大学等单位建立的中国荷斯坦牛基因组选择参考群体和技术平台，2019年共计对3 031头荷斯坦公牛进行了基因组遗传评估，包括基因组性能指数(genomic China performance index, GCPI)及产奶量、乳蛋白率、乳脂率、体型总分、泌乳系统、肢蹄和体细胞评分等7个 性状的GEBV；其中，1 686头公牛目前已经获得后裔测定成绩，而包含性能指数(China perfor-mance index，CPI)及上述7个性状的估计育种值(estimated breeding value，EBV)。

1 686头荷斯坦公牛的女儿DHI数据来自全国27个省市自治区2012—2019年期间2 018个牛场的416 086头女儿的720 393条DHI以及1 216头荷斯坦公牛的56 902头女儿的113 710条体型鉴定记录。DHI数据包括305 d产奶量、305 d乳脂率、305 d乳蛋白率与305 d体细胞评分；体型鉴定数据包括体型总分、肢蹄评分与泌乳系统评分。其中3个 305 d产奶性状和305 d体细胞评分的表型值是基于DHI测定日记录，根据泌乳曲线拟合的数学模型计算得到的；体型总分、肢蹄评分和泌乳系统评分为各部位中所包含的体型性状按各自权重计算得到。全部表型数据由中国奶牛数据中心提供。

1.1.2 数据整理 通过对荷斯坦青年公牛女儿生产性能记录原始数据的分布进行分析，结合奶牛的生理规律，并参照行业标准《中国荷斯坦牛生产性能测定技术规范》，对生产性能数据按以下标准进行筛选：1) 牛号编码正确；2) 第1、2、3泌乳期产奶性状记录完整；3) 无重复记录；4) 305 d产奶量为1 500～18 000 kg；5) 305 d乳脂率为2.0%～6.0%；6) 305 d乳蛋白率为2.0%～5.0%；7) 305 d体细胞评分为0～9；8) 系谱完整。

按照国家标准《中国荷斯坦牛体型鉴定规程》，对体型鉴定数据按以下质控标准进行筛选：1) 母牛的牛号正确；2) 第1泌乳期体型测定记录完整；3) 牛 只体型测定记录完整；4) 各部位评分和体型总分介于0～100分之间；5) 无重复记录；6) 系谱完整。

1.2 统计分析方法

1.2.1 荷斯坦公牛的基因组遗传评估结果与后裔测定结果相关性分析 利用R软件，对基因组遗传评估结果与后裔测定结果进行相关性分析，分别计算以下指标间的Spearman相关系数：1)荷斯坦青年公牛GCPI与后裔CPI；2)305 d产奶量、乳蛋白率、乳脂率、体细胞评分和体型总分、泌乳系统评分、肢蹄评分等7个性状各自的GEBV与EBV；3)对 荷斯坦青年公牛的GCPI、CPI分别由高到低进行排名，计算GCPI、CPI排名均在前 10%、前 20%、前 30%、前 40%、前 50%的公牛以上7个性状GEBV 与 EBV。

1.2.2 基于女儿性能数据分析基因组选择效果 根据305 d产奶量、乳蛋白率、乳脂率、体细胞评分和体型总分、泌乳系统、肢蹄评分等7个性状的GEBV分别对公牛由高到低进行排名并划分为2种分组，分组1划分标准见表1，分组2划分标准见表2。

表1 荷斯坦青年公牛GEBV划分标准

表2 荷斯坦青年公牛高低组GEBV划分标准

通过R软件对各组公牛女儿表型值进行方差分析，并利用Bonferroni法进行多重比较。

产奶性状分析模型：y=μ+F+P+G+e，体型性状分析模型：y=μ+F+G+e，式中，y为女儿表型值，μ为群体均值，F为场-固定效应，P为胎次效应-固定效应，G为公牛GEBV分组-固定效应，e为随机残差。

2 结 果

2.1 产奶与体型性状的基本统计量

由表3可知，产奶量、乳脂率和乳蛋白率的均值分别为8 382.99 kg、3.85%和3.31%，变异系数均小于0.30；体细胞评分的均值为4.02，变异系数为0.36；体型总分、肢蹄和泌乳系统评分的均值分别为85.68、85.46和85.91，变异系数均小于0.10。

表3 产奶与体型性状的基本统计量

图1为产奶性状和体型性状的表型值分布直方图，可以看出产奶量、乳脂率、乳蛋白率和SCS基本上都服从正态分布，而体型总分、肢蹄评分和泌乳系统评分均偏离正态分布，正态分布检验也显示，产奶性状的4个性状服从正态分布，而体型性状的3个性状偏离正态分布。

图1 产奶性状和体型性状表型值分布直方图

2.2 公牛的基因组遗传评估与后裔测定成绩相关性分析

荷斯坦公牛GCPI与后裔CPI的Spearman相关系数为0.47。各性状的GEBV与EBV的相关性分析结果如表4所示，产奶量和体细胞评分具有较强的正相关(0.4

表4 荷斯坦青年公牛各性状GEBV与EBV的相关性

本研究进一步比较分析了GCPI、CPI排名均在前10%、前20%、前30%、前40%、前50%的公牛各性状GEBV和EBV的相关性，结果如表5所示，排名前10%的公牛各性状育种值Spearman相关系数除乳蛋白率和乳脂率外都较高(rs>0.7)，排名前20%、前30%、前40%、前50%的公牛及总样本的各性状育种值Spearman相关系数均低于前10%。

表5 荷斯坦青年公牛各性状GEBV与EBV的相关系数

2.3 基于女儿生产性能测定(DHI)数据分析公牛基因组选择效果

2.3.1 公牛GEBV连续分组1 公牛女儿表型值和公牛GEBV的变化趋势如图2所示。利用全国的生产性能数据分析基因组选择效果，发现公牛产奶量、乳蛋白率GEBV与公牛女儿表型值的变化趋势一致，随着表型值下降而下降；乳脂率中除低组外，GEBV与表型值变化趋势基本一致；体细胞评分高组、中高组与中等组的GEBV与表型值变化趋势一致，中低组与低组的表型值呈上升趋势，与GEBV变化趋势不一致。

图2 各公牛组GEBV及其女儿表型值变化趋势

通过对北京和上海地区及其他省市(地区)生产性状的基因组选择效果进行比较分析发现，在产奶量中，北京及上海地区公牛GEBV与公牛女儿表型值的变化幅度较其他地区的变化幅度更为一致；在乳蛋白率中，其他地区公牛女儿表型值呈先上升后下降的趋势，与公牛GEBV的下降趋势不一致；乳脂率中除高组外，北京及上海地区公牛GEBV与公牛女儿表型值变化趋势较为一致均呈下降趋势，其他地区各组公牛GEBV呈下降趋势，而与公牛女儿表型值先下降后上升的趋势不一致；体细胞评分中，其他地区公牛女儿表型值在低组呈上升趋势，与公牛GEBV的下降趋势不一致。

生产性状各GEBV分组的最小二乘均值与多重比较结果如表6所示，利用全国生产性能数据对产奶性状基因组选择效果进行分析，结果显示，公牛各生产性状GEBV从高到低各组间女儿相应性状的表型值大部分存在极显著差异(P<0.01)。其中，公牛乳蛋白率GEBV低组与中低组组间差异不显著(P>0.05)；乳脂率GEBV排名后5%的公牛女儿的乳脂率异常偏高，且和中等组、中高组相比，差异不显著(P>0.05)；体细胞评分中等组、中低组、低组的组间差异不显著(P>0.05)。

表6 产奶性状表型值及各GEBV组的比较(最小二乘均值±标准误)

利用北京及上海地区生产性能数据对产奶性状基因组选择效果进行分析，结果显示，产奶量和乳蛋白率各组之间均存在极显著差异(P<0.01)；在乳脂率中，高组低于中高组，低组略高于中低组，差异均不显著(P>0.05)；在体细胞评分中，除中组与中低组之间差异不显著(P>0.05)，其他各组均差异显著且由高组到低组呈下降趋势。

利用其他地区生产性能数据对产奶性状基因组选择效果进行分析，结果显示，产奶量中各组由高到低呈下降趋势，且各组间差异极显著(P<0.01)；在乳蛋白率中，高组乳蛋白率水平异常偏低，低于中高组和中等组，且差异极显著(P<0.01)；在乳脂率中，高组低于中高组，但差异不显著(P>0.05)，而低组高于中低组，且差异极显著(P<0.01)。

2.3.2 公牛GEBV两尾分组2 对各性状GEBV排名高低组公牛在女儿表型值上的差异进行对比分析，结果如图3所示。利用全国的生产性能数据高低组进行分析，结果显示，产奶量、乳脂率和体细胞评分的前5%和后5%差值分别为999 kg、0.11%和0.34分，前25%与后25%差值分别为452 kg、0.02%和0.21分，前5%与后5%的差值均大于前25%与后25%的差值；前5%和后5%公牛的女儿在乳蛋白率上仅相差0.04%，与前25%和后25%的差值(0.04%)相同。

图3 公牛女儿高组与低组表型值对比图

北京及上海地区产奶量的前5%与后5%、前25%与后25%的差值分别为566 kg和419 kg，其他省市(地区)的差值分别为1 008 kg和473 kg，其他省市(地区)的高低组差值要高于北京及上海地区；在乳蛋白率中，北京及上海地区的高组与低组差值分别为0.10%和0.04%，均高于其他省市(地区)的0.01%和-0.01%；北京及上海地区的乳脂率与体细胞评分的前5%与后5%差值分别为0.05%和0.48，均高于其他省市(地区)的0.02%和0.28，而前25%与后25%的差值要低于其他省市(地区)。

2.4 基于女儿体型鉴定数据分析公牛因组选择效果

2.4.1 公牛GEBV连续分组1 公牛GEBV及其女儿表型值变化趋势如图4所示。利用全国体型鉴定数据分析基因组选择效果，结果发现，公牛女儿的体型总分、泌乳系统表型值与公牛GEBV变化趋势相反，肢蹄评分表型值与公牛GEBV趋势相同，呈下降趋势。

图4 公牛GEBV及其女儿表型值变化趋势

进一步比较分析北京和上海地区及其他省市(地区)的体型性状基因组选择效果发现，在北京和上海地区，3种体型性状公牛GEBV与公牛女儿表型值变化趋势基本一致，均呈下降趋势；在其他地区中，公牛女儿体型总分与泌乳系统评分呈上升趋势，与公牛GEBV的下降趋势相反，公牛女儿的肢蹄评分为先上升后下降，也与公牛GEBV的下降趋势不一致。

体型性状各GEBV分组的最小二乘均值与多重比较结果如表7所示。利用全国体型鉴定数据分析基因组选择效果，结果表明，在体型总分和泌乳系统表型值中，高组比低组分别低1.04分和1.09分，且差异均极显著(P<0.01)；在肢蹄评分表型值中，高组、中高组和中等组公牛间女儿肢蹄评分差异不显著(P>0.05)，而高组比低组高0.81分，差异极显著(P<0.01)。

表7 体型性状表型值及各GEBV组的比较(最小二乘均值±标准误)

利用北京及上海地区的体型鉴定数据分析基因组选择效果，结果表明，在体型总分中，除中组外其他各组均呈下降趋势，高组与低组相差0.70分，差异极显著(P<0.01)；泌乳系统各组均呈下降趋势，高组与低组相差1.80分，差异极显著(P<0.01)；肢蹄评分中，各组呈下降趋势，高组与低组相差1.20分，差异极显著(P<0.01)。

利用其他省市(地区)的体型鉴定数据分析基因组选择效果，结果表明，体型总分与泌乳系统中各组呈上升趋势，高组比低组分别低0.55分与1.10分，且差异极显著(P<0.01)；肢蹄评分各组也呈上升趋势，但高组与低组之间差异不显著(P>0.05)。

2.4.2 公牛GEBV双尾分组2 对各体型性状GEBV排名高低组的公牛在女儿表型值上的差异进行对比分析，结果如图5所示。在体型总分与泌乳系统中，前5%均低于后5%且分别相差-1.04和-1.09，而前25%与后25%分别相差0.30和0.46；在肢蹄评分中，高低组分别相差0.81和0.50，前5%与后5%的差值要大于前25%与后25%的差值。

图5 公牛女儿高组与低组表型值对比图

北京及上海地区的体型总分、肢蹄评分和泌乳系统评分前5%与后5%的差值分别为0.70、1.20和0.80，前25%与后25%的差值分别为0.30、0.82和0.46。其他省市(地区)的前5%与后5%的差值分别为-0.55、0.13和-1.10，前25%与后25%的差值分别为-2.42、0.02和-4.44，体型总分与泌乳系统评分的高组均低于低组。

3 讨 论

3.1 产奶与体型性状表型分析

本研究中，产奶性状的表型数据均服从正态分布，而体型性状的3个性状有些偏离正态分布。体型性状是由体型鉴定员按照9分制(1～9的整数)对头胎产犊后30～180 d内荷斯坦母牛20个性状进行线性评分，再按照各个性状所占权重计算得到部位评分，然后各部位评分加权得到体型总分。因此，体型性状表型值不同于数量性状的连续分布，使其偏离正态分布。

3.2 基因组遗传评估结果与后裔测定成绩的相关性

整体上看，我国荷斯坦青年公牛基因组遗传评估结果与后裔测定结果一致性较好。公牛GCPI与CPI的相关系数为0.47，呈正相关，与国内其他研究一致[37-39]。本研究中，产奶量GEBV与EBV的相关系数为0.69，要低于Wiggans等[12]报道的产奶量育种值的可靠性(0.70～0.80)，而体细胞评分的相关系数为0.46，也低于新西兰家畜遗传改良公司(Livestock Improvement Corporation Ltd.，LIC)公布的体细胞评分GEBV估计可靠性(0.50～0.67)[39]，主要原因可能为我国参考群体规模较小，且奶牛的表型数据准确性也低于发达国家[40-41]；而乳蛋白率、乳脂率的相关系数分别为0.14、0.18，要低于产奶量和体细胞评分的相关性，原因可能是乳脂率与乳蛋白率的测定易受DHI奶样采集和样品测定环节中操作不规范的影响。体型性状(体型总分、泌乳系统和肢蹄评分)GEBV与EBV的相关性均较低，分别为0.19、0.24和0.26，体型性状评分是由体型鉴定员通过观察来对每个部位进行打分，在鉴定过程中易受鉴定员的主观因素和鉴定环境的影响，而导致体型鉴定数据的准确性降低；相较于产奶性状，体型性状的数据量较少，也可能是导致体型性状准确率较产奶性状低的原因。有研究表明，通过扩大奶牛参考群体规模，可以提高公牛GEBV的准确性[42]，而我国的奶牛参考群体相较于奶业发达国家数量较少，这可能是导致我国基因组选择准确性低的原因。

3.3 基于女儿表型值分析荷斯坦青年公牛基因组评估效果

产奶性状及体细胞评分的GEBV与表型值的变化趋势基本一致，这与德国育种数据处理中心(Vereinigte Informations systeme Tierhaltung，VIT)公布的结果相一致，平均基因组育种值随着表型值升高而升高，且生产性状的GEBV可靠性达75%；加拿大奶业信息网(Canada Dairy Network，CDN)公布的基因组育种值考虑系谱信息与直接基因组评估值(direct genomic value，DGV)，平均可靠性达60%。基于全国数据的乳脂率GEBV排名后5%的公牛，其女儿乳脂率反而偏高，除北京和上海以外的其他省市(地区)排名后25%的公牛其女儿的乳脂率水平也异常偏高，其原因可能是有些省市(地区)的部分牛场在DHI奶样采集和测定环节的规范性存在不足，导致乳脂率、乳蛋白率数据的可靠性偏低，不能很好地反映出基因组选择的效果。另外体型总分、泌乳系统评分女儿表型值变化趋势同公牛GEBV变化趋势相反。荷斯坦奶牛的体型鉴定是将20个线性性状，按生物学特性的变异范围，并由各省市(地区)的鉴定员按照1～9的整数来进行线性评定。因此，体型数据易受鉴定员因素影响，使得不同省市(地区)的体型鉴定数据存在一定的误差，由此，全国的体型鉴定数据不能很好地反映公牛基因组育种值的高低。

北京及上海地区的乳蛋白率、乳脂率和体细胞评分的高组与低组差值均大于其他省市(地区)，体型性状表型值高组均显著高于低组，而其他省市(地区)的高组却低于低组。北京及上海地区相对于其他省市(地区)更能反映公牛的基因组选择效果。主要原因可能是北京及上海地区的奶牛育种历史悠久，育种水平居全国之首，牛场DHI和体型鉴定较规范，产奶、体细胞评分和体型性状的表型数据准确性较高；另外，我国的奶牛基因组选择参考群牛只主要来自北京及上海地区的规范牛场，这两个地区的荷斯坦青年公牛和参考群的遗传联系更高，因此，北京和上海地区的数据可以更好地反映基因组选择效果。

4 结 论

基于1 686头荷斯坦公牛基因组选择及其女儿表型数据的分析结果表明，公牛GCPI与CPI呈正相关(rs>0.3)，产奶量和体细胞评分的基因组育种值(GEBV)与估计育种值(EBV)呈较强的正相关(0.4