畜牧兽医学报  2023, Vol. 54 Issue (10): 4174-4185. DOI: 10.11843/j.issn.0366-6964.2023.10.016    PDF    
利用高密度SNP芯片评估中国地方肉牛品种基因组亲缘关系
马浩然1, 张路培1, 金生云2, 宝金山3, 李红艳4, 高会江1, 徐凌洋1, 王泽昭1, 李俊雅1     
1. 中国农业科学院北京畜牧兽医研究所 牛遗传育种创新团队,北京 100193;
2. 乌拉盖管理区农牧技术推广中心,乌拉盖 026321;
3. 乌拉盖管理区供销合作社,乌拉盖 026321;
4. 内蒙古通辽市畜牧业发展中心,通辽 028000
摘要:旨在在系统对比分析亲缘关系不同计算方法基础上探索适合我国肉牛地方品种的亲缘关系评估方法。本研究主要以柴达木牛等10个肉牛地方品种为研究对象,使用重抽样方法获得地方品种的模拟数据,以此为基础使用PCA聚类结果为参照,分别利用预测误差方差、广义决定系数、预测误差相关系数及SNP与QTL的连锁一致性程度系统对比了不同评估方法对地方品种亲缘关系的分类结果,并探索了遗传力因素对不同亲缘关系评估方法的影响。通过PCA分析可知10个肉牛地方品种可分为3大类,分别为北方牛品种(柴达木牛、西藏牛、蒙古牛以及延黄牛)、南方牛品种(文山牛、南丹牛以及雷琼牛)和西南牛品种(平武牛、凉山牛以及昭通牛),该分类结果与上述品种地理分布较为一致。与PCA结果对比发现,基于LD一致性的亲缘关系评估方法的评估结果与PCA聚类结果一致,且该方法能够使用皮尔逊相关系数量化品种间亲缘关系,具有较好的准确性。PEVD法、CD法与r法3种方法与上述方法相比评估群体间亲缘关系时容易受到性状估计育种值的误差方差影响,从而造成种间亲缘关系评估结果出现误差。评估结果影响因素分析发现,PEVD法、CD法和r法与PCA和LD一致性评估法相比,易受到评估性状遗传力的影响,评估结果稳定较差。而PCA和LD一致性评估法由于仅依赖基因组数据,不受到遗传力影响,能够更稳定的量化评估品种间亲缘关系,具有更好的可靠性。因此,基于LD一致性评估方法结果可准确量化评估肉牛品种间亲缘关系且评估结果稳定性较高。
关键词多品种基因组选择    亲缘关系    肉牛地方品种    高密度SNP芯片    连锁不平衡    
Assessment of the Genomic Relationships for Chinese Indigenous Beef Cattle Using High-density SNP Chip
MA Haoran1, ZHANG Lupei1, JIN Shengyun2, BAO Jinshan3, LI Hongyan4, GAO Huijiang1, XU Lingyang1, WANG Zezhao1, LI Junya1     
1. Laboratory of Molecular Biology and Bovine Breeding, Institute of Animal Science, Chinese Academy of Agricultural Sciences, Beijing 100193, China;
2. Agricultural and Animal Husbandry Technology Extension Center of Ulagai Precinct, Ulagai 026321, China;
3. Ulagai Management District Supply and Marketing Cooperative, Ulagai 026321, China;
4. Tongliao Animal Agriculture Development Service Center, Tongliao 028000, China
Abstract: This study aimed to explore a suitable method for assessing the genetic relationships of indigenous beef breeds in China, based on a systematic comparative analysis of different computational methods for determining kinship. This study focused on 10 local beef cattle breeds, such as Chaidamu cattle, and utilized the resampling method to generate simulated data for these local breeds. Based on this, the study utilized PCA clustering results as a reference and systematically compared different evaluation methods, including predictive error variance, generalized coefficient of determination, predictive error correlation coefficient, and the degree of linkage consistency between SNPs and QTLs, to classify the genetic relationships of these local breeds. Additionally, the study explored the influence of genetic factors on different methods of assessing kinship. The PCA analysis showed that the 10 indigenous beef cattle breeds could be classified into 3 major categories, northern cattle breeds (Chaidamu cattle, Xizang cattle, Mongolian cattle and Yanhuang cattle), southern cattle breeds (Wenshan cattle, Nandan cattle and Leiqiong cattle) and southwestern cattle breeds (Pingwu cattle, Liangshan cattle and Zhaotong cattle). The classification results were consistent with the geographical distribution of the above breeds. Compared with the PCA results, the correlation of persistence of LD phase affinity assessment method was consistent with the PCA clustering results, and the method was able to quantify the relationship using Pearson correlation coefficient relationships with better accuracy. The 3 methods, PEVD, CD and r, were susceptible to the influence of the error variance of the estimated breeding values of the traits when assessing the relationship compared to the above methods, resulting in errors in the assessment results of the relationship. The analysis of factors affecting the assessment results revealed that the PEVD, CD and r methods were susceptible to the influence of the heritability of the assessed traits compared to the PCA and LD phase consistent assessment methods, and the assessment results were less stable. On the other hand, the PCA and LD phase consistent assessment methods are more stable and quantitative in assessing relatedness among breeds because they rely only on genomic data and are not affected by heritability and have better reliability. Therefore, the results based on LD phase consistent assessment method can accurately and quantitatively assess the relationship between breeds and the assessment results are more stable.
Key words: multi-breed genomic selection    relationship    indigenous beef cattle breed    high-density SNP chip    linkage disequilibrium    

基因组选择(genomic selection, GS)技术由Meuwissen等[1]于2001年提出,由于该技术能够实现后备牛早期选择,因此可以大幅缩短牛育种世代间隔,提升肉牛育种遗传进展。目前该技术已经成为世界肉牛主要育种技术手段之一[2]。世界范围内牛品种资源十分丰富,2009年,Harris等[3]尝试选择部分代表性牛品种组建参考群进行品种间基因组联合评估。Hayes等[4]使用荷斯坦牛、娟姗牛及其他奶牛品种进行多群体基因组选择研究,其研究证明使用奶牛混合参考群体进行基因组选择会提升目标性状的育种值估计准确性。

上述研究共同发现,多品种GS育种值估计准确性主要受到参考群体规模和验证群体间亲缘关系的影响[5]。全基因组选择准确性依赖于单核苷酸多态性(single nucleotide polymorphism, SNP)位点和数量性状基因座(quantitative trait locus, QTL)之间的连锁不平衡(linkage disequilibrium, LD)程度[6]。而亲缘关系较远的品种间由于受到等位基因频率分布以及基因间或基因与环境间的互作等影响导致QTL及标记LD一致性差异从而使多品种基因组选择的育种值估计准确性降低[7-8]。Wientjes等[9]研究发现,当品种间的遗传关系接近时,种间预测准确性较高。而遗传关系远时利用多品种的混合参考群体进行预测的准确性并不比单一品种准确性高。还有部分研究表明,不仅仅参考群体与候选个体之间的亲缘关系能够影响基因组选择的预测准确性,位于参考群中个体间的亲缘关系同样会影响选择准确性[10-11]。Lund等[12]的研究也同样指出,当远亲品种组合作为基因组选择参考群并使用基因组最佳线性无偏预测(genomic best linear unbiased prediction, GBLUP)模型评估时,育种值估计准确性并未提升,甚至有所降低。但通过使用更复杂的贝叶斯变量选择模型并结合更密集的标记集或标记的功能子集,同时对与QTL强连锁不平衡的基因组标记进行加权,则可能利用远缘品种的遗传信息来提高基因组预测准确性。因此,准确评估不同肉牛品种间的亲缘关系是实现肉牛多品种和跨品种基因组选择的第一步。

相较于系谱数据,高密度SNP芯片分型数据可以更加准确的评估品种间的亲缘关系[13]。目前估计群体间遗传关系的方法主要有以下4种:一是利用预测误差方差差异(prediction error variance of differences, PEVD)评估种间亲缘系数法[14],该方法通过计算品种之间育种值(estimated breeding value, EBV)差异的预测误差方差,来衡量品种间的遗传关联。二是基于广义决定系数(generalized coefficient of determination, CD)[15]评估品种间亲缘关系。CD定义为估计育种值比较的可靠性,即预测值差异与真实值差异间相关系数的平方。三是利用预测误差相关系数(prediction error correlation, r) 开展种间亲缘关系评估。该方法是由Lewis等[16-17]提出的基于预测误差方差(prediction error variance, PEV)的简化评估方法,即使用两品种之间的PEV相关系数来衡量亲缘关系。四是基于不同品种间SNP位点与QTL的连锁一致性程度评估品种间亲缘关系。即基因组关系的品种间LD一致性分析[18]。上述4种方法在我国肉牛种间亲缘关系评估工作中的评估性能尚无系统评价。

中国幅员辽阔,是拥有世界上牛品种最多的国家之一,共有130多个牛品种。目前支撑肉牛产业的品种主要有70多种[19]。近年来,随着种业振兴行动相关工作的持续推进,对地方品种优秀种质资源的深度挖掘和新品种培育工作均在持续进行,地方品种对基因选择的需求在不断提升。但大多数品种育种工作仍处于初期阶段,生产性能测定等常规基础育种体系建设仍不健全,育种群规模小、系谱记录不健全等客观因素导致了传统育种无法有效支撑肉牛高效育种工作。对全部品种逐一建设基因组选择参考群体既不经济且实现难度较大,因此探索肉牛多品种基因组选择方法势在必行。

目前,我国尚未开展肉牛多品种基因组选择相关研究工作。本研究基于模拟数据和高密度SNP芯片分型数据,通过使用5种品种间亲缘关系分析方法,对比分析了我国10个地方肉牛品种间的亲缘关系。旨在探索适合我国肉牛多品种间亲缘关系评估的最优策略,为肉牛多品种基因组选择技术研发奠定基础。

1 材料与方法 1.1 实际数据

本研究分析了我国10个地方牛品种基因型数据,分别是:柴达木牛(CDM, n=25)、雷琼牛(LQ, n=26)、凉山牛(LS, n=22)、蒙古牛(MG, n=21)、南丹牛(ND, n=25)、平武牛(PW, n=24)、文山牛(WS, n=24)、西藏牛(XZ, n=26)、延黄牛(YH, n=24)和昭通牛(ZT, n=23)。试验个体均静脉采血20 mL冻存,并用2 mL血液提取DNA,使用Illumina Bovine HD Bead Chip(770K,Illumina, Inc., San Diego, CA)对样本进行基因分型,该芯片由777 962个SNPs位点组成。基因分型和质量控制(quality control, QC)使用Genome Studio软件进行。除去性染色体及质粒DNA上的42 669个SNPs位点,常染色体上分布的SNPs位点共计735 293个。

得到基因分型数据后使用PLINK(V1.9)软件对每个品种常染色体位点单独进行质控,其SNPs位点保留标准为:位点检出率(call rates, CR)大于等于90%,个体检出率大于90%,最小等位基因频率(minor allele frequencies, MAF)大于0.01,哈迪-温伯格平衡(Hardy-Weinberg equilibrium, HWE)P值大于1.0×10-6。10个肉牛地方品种基因型质控样本的描述性统计结果见表 1。将通过质控后的数据按品种使用BEAGLE(V5.0)对缺失的位点进行填充,填充参数为软件默认参数。

表 1 试验样本个体及标记数量统计 Table 1 Statistics on the number of individuals and markers in the test samples
1.2 模拟数据

本研究使用了重抽样模拟方法,在10个地方品种的真实基因型数据基础上,将每个品种基因组数据模拟至1 500头。具体模拟方法及参数设置见Xu等[20]的研究。

基于模拟基因型数据,获得每个品种的表型数据。本研究中共计模拟3个不同遗传力性状表型,htrait12、htrait22和htrait32分别为0.1、0.3和0.6。表型模拟具体步骤为:一是抽取QTL集合。从单一品种全部SNPs位点集合中随机抽取位点作为QTL位点,QTL位点数量见表 2所示。二是获得QTL标记效应。每个QTL的标记效应分别从3种正态分布中随机抽取获得。本研究中模拟了3种不同效应的QTL(小效应、中效应以及大效应),服从的分布分别为:~N(0, 0.0001σg2)、~N(0, 0.001σg2)、~N(0, 0.01σg2),其中σg2为性状真实加性遗传方差参考值。三是获得模拟加性遗传效应值(true breeding value, TBV)。将所得的SNP效应值与对应位点的基因型值相乘,通过累加得到加性遗传效应。四是获得模拟表型。依据设定遗传力计算得到残差,最终将残差与加性遗传效应相加得到模拟表型数据。

表 2 表型模拟的遗传参数设置 Table 2 Genetic parameter setting of phenotypic simulation

在模拟中,残差效应直接从服从均值为0、方差为σg2(1-h2)/h2的正态分布中随机抽取,其中σg2为加性遗传方差,h2为表型的遗传力。个体的表型值直接使用育种值与残差的累加值,真实育种值为各QTL效应之和,公式如下:

$ T B V_i=\sum\limits_{j=1}^n x_{i j} a_j $

其中,xij是个体j的第i个QTL的基因型,编码为0、1、2;aji个QTL的加性效应;n是QTL的数量。

表型值为:

$ P_i=T B V_i+\sigma_{e i} $

其中,Pi为模拟表型值,TBVi为第i个个体育种值,σei为第i个个体残差。

1.3 育种值估计模型

在本研究中,使用GBLUP模型进行相应的加性遗传方差以及残差方差估计,其模型如下:

$ y=\boldsymbol{X} b+\boldsymbol{Z} g+e $

其中,y是个体的表型值;gb分别是随机加性遗传效应和固定效应,ZX分别为随机加性遗传效应和固定效应对应的关联矩阵;e是残差向量。

1.4 不同亲缘关系评估方法

本研究首先对模拟数据和实际数据进行主成分分析(principal component analysis,PCA),分别获得模拟数据和真实数据的品种间聚类结果,并将其作为后续分析结果的基准参考。

1.4.1 品种间LD一致性评估法   使用PopLDdecay软件分别计算了模拟数据和实际数据的LD衰减距离,并绘制LD衰减图。然后利用LD的r2值计算品种间的皮尔逊相关系数,并将该值作为不同品种基因组LD一致性衡量指标,评价品种间亲缘关系。r2计算公式如下:

$ r^2=\frac{D^2}{f(A) f(a) f(B) f(b)} $

公式中D=f(AB)-f(A)f(B), 其中f(AB)、f(A)、f(a)、f(B)和f(b)分别为单倍型AB的基因型频率,A、a、B和b以及等位基因的频率。

皮尔逊相关的计算公式如下:

$ r_{i j}=\frac{COV(i, j)}{\sigma_i \sigma_j} $

公式中rij为计算所得皮尔逊相关系数,COV(i, j)为品种i与品种j的协方差,σiσj分别为两个品种的标准差。

1.4.2 预测误差方差   本研究中,育种值估计模型的混合模型方程组为:

$ \left[\begin{array}{cc} \boldsymbol{X}^{\prime} \boldsymbol{X} & \boldsymbol{X}^{\prime} \boldsymbol{Z} \\ \boldsymbol{Z}^{\prime} \boldsymbol{X} & \boldsymbol{Z}^{\prime} \boldsymbol{Z}+\boldsymbol{G}^{-1} \lambda \end{array}\right]\left[\begin{array}{l} \hat{b} \\ \hat{g} \end{array}\right]=\left[\begin{array}{c} \boldsymbol{X}^{\prime} y \\ \boldsymbol{Z}^{\prime} y \end{array}\right] $

其中,ZX分别为随机加性遗传效应和固定效应对应的关联矩阵,G为个体间的加性遗传相关矩阵(G矩阵),其中,$\lambda=\frac{\sigma_e^2}{\sigma_a^2}=\frac{\left(1-h^2\right)}{h^2} $σe2为随机残差,σa2为加性方差,h2为性状的遗传力。则该混合模型方程组系数矩阵的逆矩阵为[17]

$ \boldsymbol{C}^{-1}=\left[\begin{array}{cc} \boldsymbol{X}^{\prime} \boldsymbol{X} & \boldsymbol{X}^{\prime} \boldsymbol{Z} \\ \boldsymbol{Z}^{\prime} \boldsymbol{X} & \boldsymbol{Z}^{\prime} \boldsymbol{Z}+\boldsymbol{G}^{-1} \lambda \end{array}\right]^{-1}=\left[\begin{array}{cc} \boldsymbol{C}^{11} & \boldsymbol{C}^{12} \\ \boldsymbol{C}^{21} & \boldsymbol{C}^{22} \end{array}\right] $

PEVD计算公式如下[14]

$ \begin{array}{c} PEVD\left(\hat{\mu}_i-\hat{\mu}_j\right)=\left[PEV\left(\hat{\mu}_i\right)+PEV\left(\hat{\mu}_j\right)-\right. \\ \left.2 PEC\left(\hat{\mu}_i, \hat{\mu}_j\right)\right]=\left(\boldsymbol{C}_{i i}^{22}+\boldsymbol{C}_{j j}^{22}-2 \boldsymbol{C}_{i j}^{22}\right) \sigma_e^2 \end{array} $

其中$ PEV\left(\hat{\mu}_i\right)$$ P E V\left(\hat{\mu}_j\right)$分别为第ij个品种的个体加性效应的预测误差方差,$PEC\left(\hat{\mu}_i, \hat{\mu}_j\right) $为第ij个品种的个体加性效应的预测误差协方差,σe2为残差方差。C22为系数矩阵的逆矩阵的子矩阵。

1.4.3 广义决定系数   广义决定系数计算公式如下[17]

$ C D=1-\lambda \frac{\boldsymbol{C}_{i i}^{22}+\boldsymbol{C}_{i i}^{22}-2 \boldsymbol{C}_{i j}^{22}}{\boldsymbol{K}_{i i}+\boldsymbol{K}_{j j}-2 \boldsymbol{K}_{i j}} $

$\lambda=\frac{\sigma_e^2}{\sigma_a^2}=\frac{\left(1-h^2\right)}{h^2} $σe2为随机残差,σa2为加性方差,h2为表型的遗传力。C22为系数矩阵的逆矩阵的子矩阵, K矩阵为基因组关系矩阵。

1.4.4 预测误差相关系数   预测误差相关系数计算公式如下:

$ r_{i j}=\frac{PEC\left(\hat{a}_i, \hat{a}_j\right)}{\sqrt{PEV\left(\hat{a}_i\right) PEV\left(\hat{a}_j\right)}} $

其中$PEV\left(\hat{a}_i\right) $$PEV\left(\hat{a}_j\right) $分别为第ij个品种的个体加性效应的预测误差方差,$PEC\left(\hat{a}_i, \hat{a}_j\right) $为第ij个品种的个体加性效应的预测误差协方差。

本研究中群体间PEVD、CD、r、PEV以及PEC的计算使用了不同群体的所有个体两两配对的均值。其中CD值和r值与品种间亲缘关系成正比(数值越大关系越近),而PEVD值则相反(数值越小关系越近)[13]

2 结果 2.1 不同地方牛品种PCA分析结果

图 1展示的是10个不同地方品种模拟数据和实际数据PCA聚类结果。图 1a中,PC1及PC2为第一主成分和第二主成分,分别解释了48.76%以及16.98%的变异。由图可见,10个地方品种模拟数据明显聚集为3大类,其中延黄牛、蒙古牛、西藏牛以及柴达木牛聚为一类,可初步判断上述4个品种间亲缘关系较近。平武牛、昭通牛和凉山牛聚为一类,文山牛、雷琼牛以及南丹牛之间虽然聚类显示一定距离,但相较于其他品种仍然可以被归为一类,品种间亲缘关系较远。图 1b展示的是真实基因型数据聚类结果。如图所示,PC1及PC2两个主成分分别解释了56.71%以及4.92%的变异,且在真实群体中10个地方品种整体上也聚类为3大类,每大类包含的品种聚类结果与模拟数据一致。

图 1 10个不同地方牛品种PCA分析结果图 Fig. 1 Principal component analysis of 10 different indigenous cattle breeds
2.2 不同地方牛品种K-means聚类结果

图 2展示的是10个不同地方品种模拟数据的K-Means聚类结果,设定10个初始化聚类中心进行聚类,取最终的聚类中心进行绘图。由图可见,10个地方品种模拟数据明显聚集为3大类,首先,延黄牛、蒙古牛、西藏牛以及柴达木牛聚为一类,其中该类中蒙古牛与延黄牛聚为一类,柴达木牛与西藏牛聚为一类,其次,平武牛、昭通牛和凉山牛聚为一类,其中平武牛与另外两个品种亲缘关系较远,最后,文山牛、雷琼牛以及南丹牛分为一类,但同样,其中文山牛与另外两个品种亲缘关系较远,该结果得到的亲缘关系与PCA聚类分析所展示的结果图一致,在一定程度上证实了PCA分析的可靠性。

图 2 10个不同地方牛品种聚类结果图 Fig. 2 Clustering results of 10 different indigenous cattle breeds
2.3 种间LD一致性评估法

2.3.1 不同地方牛品种LD衰减结果   图 3分别展示了模拟数据和实际数据基因组r2的衰减趋势。模拟数据结果显示,延黄牛、蒙古牛、柴达木牛与西藏牛LD衰减趋势一致,其LD衰减距离分别为82.29、85.41、89.04和87.13 kb,其品种间衰减距离较为相似。南丹牛及雷琼牛2个品种LD衰减距离分别为120.12和123.79 kb。剩余4个品种LD衰减距离分别为105.78、101.26、109.47和102.58 kb,其衰减距离值较为相似。

图 3 10个不同地方牛品种LD衰减结果图 Fig. 3 LD decay of 10 different indigenous cattle breeds

2.3.2 LD一致性评估亲缘关系   根据品种间r2计算得到的10个肉牛地方品种间亲缘关系结果如图 4所示。如图 4a所示,在模拟数据中,蒙古牛与延黄牛亲缘关系最高,品种间r2相关系数为0.64。其次,南丹牛与雷琼牛、延黄牛与西藏牛以及昭通牛与凉山牛3组品种间r2相关系数均为0.63。雷琼牛与延黄牛品种间r2相关系数为0.22,表明两个品种间亲缘关系较远。上述结果与PCA分析展示的聚类结果一致。但其中昭通牛与文山牛的品种间r2相关系数为0.59,且LD衰减趋势较为一致,LD一致性评价结果显示两品种间亲缘关系较为密切,但该发现与PCA结果存在差异。在真实群体的分析结果中(图 4b),延黄牛与蒙古牛、昭通牛与凉山牛以及南丹牛与雷琼牛品种间r2相关系数最高(r2=0.74),LD一致性评价结果显示上述品种组合间亲缘关系较为密切。其次,凉山牛与平武牛品种间r2相关系数为0.72,昭通牛与文山牛品种间r2相关系数为0.71,表明上述品种对间存在较高遗传联系。但昭通牛与文山牛LD一致性评价结果与PCA结果存在差异。第三是延黄牛与雷琼牛品种间r2相关系数仅为0.25,表明品种间亲缘关系较远。

图 4 10个不同地方牛品种亲缘关系 Fig. 4 Genomic relationship of 10 different indigenous cattle breeds
2.4 预测误差方差法

表 3展示的是基于预测误差方差法计算的群体间亲缘关系结果,表中数字表示预测误差方差,数值越低代表亲缘关系越高。由于实际数据部分表型值缺失,因此本分析方法仅对模拟数据使用。由表可知,10个肉牛地方品种PEVD值范围在0.80~0.87之间。延黄牛与蒙古牛、西藏牛以及柴达木牛PEVD值范围在0.80~0.81之间,初步表明上述4个品种间亲缘关系较近,与PCA聚类分析结果一致。平武牛与昭通牛和凉山牛间聚为一类,其PEVD值范围在0.80~0.82之间。文山牛与雷琼牛以及南丹牛之间PEVD数值为0.83~0.84,说明上述3个品种可聚为一类。但根据表 3展示结果发现,品种间PEVD值较为集中,与上述方法相比种间亲缘关系分层情况不明显。

表 3 预测误差方差法计算亲缘关系 Table 3 The genomic relationship calculated by PEVD
2.5 广义决定系数法

表 4展示的是基于广义决定系数法计算的群体间亲缘关系结果,数值越高代表亲缘关系越高。同样由于实际数据部分表型值缺失,本分析方法仅对模拟数据使用。由表可知,10个肉牛地方品种CD值范围在0.72~0.79之间。其中,延黄牛与蒙古牛、西藏牛以及柴达木牛间CD值范围在0.78~0.79之间,表明上述4个品种间亲缘关系较近,与PCA聚类分析结果一致。平武牛与昭通牛和凉山牛聚为一类,其CD值范围在0.77~0.78之间。文山牛与雷琼牛以及南丹牛之间CD数值为0.76,说明上述3个品种可聚为一类。但根据表 4展示结果发现,品种间CD值相差不大,与上述方法相比种间亲缘关系分层情况不明显。

表 4 广义决定系数法计算亲缘关系 Table 4 The genomic relationship calculated by CD
2.6 预测误差相关系数法

表 5展示的是基于预测误差相关系数法计算的群体间亲缘关系结果,数值越高代表亲缘关系越近。同样由于实际数据部分表型值缺失,本方法仅对模拟数据使用。由表可知,10个地方品种r值范围在0.000 7~0.001 3之间。延黄牛与蒙古牛、西藏牛以及柴达木牛间r值范围为0.001 2~0.001 3,初步表明上述4个品种间亲缘关系较近,与PCA聚类分析结果一致。平武牛与昭通牛和凉山牛可聚为一类,其r值均为0.001 2。文山牛与雷琼牛以及南丹牛之间r数值为0.000 9~0.001 0,说明上述3个品种可聚为一类。

表 5 预测误差相关系数法计算亲缘关系 Table 5 The genomic relationship calculated by r
2.7 遗传力大小对亲缘关系评估的影响

本研究基于10个地方牛品种的实际数据对表型数据进行了模拟,通过设置不同大小遗传力(htrait12、htrait22和htrait32分别为0.1、0.3和0.6),来探究遗传力对3种基于GEBV评估结果对群体间遗传关系进行估算方法结果的影响。研究中使用了蒙古牛与延黄牛、凉山牛与延黄牛以及雷琼牛与延黄牛3组品种组合(3种组合代表了前述研究中的高、中、低3个不同梯度的亲缘关系),如表 6所示,以蒙古牛和延黄牛为例,随着模拟性状遗传力上升PEVD法评估结果从1.01降低到0.80,下降了19.8%。表明利用PEVD方法评估的蒙古牛和延黄牛种间亲缘关系受到了性状遗传力的影响。CD法和r法的评估结果分别从0.26上升至0.79,从0.000 3上升至0.001 2。

表 6 不同遗传力对亲缘关系评估的影响 Table 6 The effect of different heritabilities on the assessment of relationship
3 讨论

Brøndum等[21]在研究中证明,将与验证群有亲缘关系的品种加入到参考群中,选择的准确性得到了提升。因此,为了更好的进行相关的多品种基因组选择的研究,评估各个品种中的亲缘关系尤为重要。目前,绝大部分研究仅分析了群体间遗传结构,并未针对品种或群体间的亲缘关系开展研究[22-26],而一些研究中仅使用了芯片或系谱数据进行了近交系数的计算,但其仅能代表个体之间的相关,并不能作为品种间的亲缘关系[27-28]。因此,本研究系统对比分析了5种亲缘关系评估方法。PCA聚类结果显示,延黄牛、蒙古牛、西藏牛以及柴达木牛主成分位置较近,平武牛、昭通牛和凉山牛主成分位置较近,文山牛、雷琼牛以及南丹牛主成分较为接近。依据PCA结果将10个肉牛地方品种可划分为3个类群,结果与上述品种地理分布相符。根据真实品种的地理位置分布分类:柴达木牛、西藏牛、蒙古牛以及延黄牛集中分布在我国北方或西北地区,属于北方牛品种。文山牛、南丹牛以及雷琼牛多分布在我国南方地区,属于南方牛品种。平武牛、凉山牛以及昭通牛主要集中分布在我国西南地区,属于西南牛品种。

3.1 种间LD一致性评估法评估结果比较

Ma等[18]的研究指出,种间LD一致性可以反映品种间亲缘关系。杨祎挺等[29]认为,不同地方猪品种的LD衰减差异大,代表了其种间遗传结构差异性大。因此,若两个品种的LD衰减速度较为一致,一定程度上代表了品种的亲缘关系。本研究中无论在分析模拟数据还是实际数据,种间LD一致性评估法评估结果均与PCA结果保持了高度一致。但该方法的优势在于可以量化品种间亲缘关系,为肉牛多品种基因组选择提供更为准确的参考。因此,基于本研究分析结果,种间LD一致性评估法是一种较为适合评估肉牛地方品种种间亲缘关系的方法。

3.2 品种间遗传关联度量方法比较

Foulley等[30]与Laloë等[31]认为遗传关联性是一个可预测的衡量标准,基于此,Kennedy和Trus[14]、Laloë[15]以及Lewis等[16]分别提出了使用PEVD、CD以及r值衡量群体间的遗传联系,但这3种方法均受到性状的遗传结构、QTL数目、群体大小和结构等因素的影响从而使预测得到的品种间遗传联系产生偏差[32]。在本研究中,由于地方品种的真实表型数据缺失,无法使用实际数据进行评估,但周子文等[13]的研究表明模拟数据与实际数据的结果存在差异,其原因可能是模拟数据中仅考虑了加性效应,并不能很好的反映真实群体性状的遗传结构等,从而使预测产生了偏差。此外,由于3种方法依赖于误差方差的估计,不同的评估模型、性状遗传力也会对结果产生影响[33],同时,评估的准确性也会影响评估种间亲缘关系,研究中使用了多品种基因选择预测模型进行计算,根据Xu等[34]的研究表明,多品种基因组选择的准确性低于传统的基因组选择,进而使误差方差估计出现了偏差,这可能是导致3种方法估计出现偏差的主要原因。同时,Kuehn等[35]研究表明,引入不具有血缘关系的个体会导致预测误差方差水平的降低,从而降低群体的遗传联系。周子文等[13]及Zhang等[36]的研究指出,基于CD值评估品种间亲缘关系时,即使系谱中不存在亲缘关系的个体,也会估计得到较高的遗传关联,从而会过高估计品种间的亲缘程度;基于G矩阵计算的r值较低,无法区分群体间的遗传差异,不能准确反映群体间的实际群体关联。在对品种间的亲缘关系进行评估时,其得到的结果不应与个体的表型产生关联,但本研究中的3种方法均需依赖表型进行相应的计算,因此,如何减少表型对亲缘关系预测的影响仍需进一步探究。

3.2.1 预测误差方差法评估果   PEVD法在评估10个肉牛地方品种亲缘关系时,评估结果与PCA结果较为一致。与周子文等[13]基于高密度SNP芯片估计猪群体间遗传关系的研究结果相似,本研究中10个地方品种间PEVD值范围在0.80~0.87之间,PEVD值较为集中,表明群体间的遗传关联没有显著差异,与LD法相比种间亲缘关系分层情况不明显,很难直观判断品种间的亲缘关系。基于相关结果,本研究同样认为PEVD方法不是理想的度量群体间遗传关系的方法。

3.2.2 广义决定系数法结果   本研究发现,10个地方品种CD值范围在0.72~0.79之间,表明上述品种间均存在较高的群体遗传联系,这与PCA分析结果存在较大差异。因此,该方法同样不是理想的度量不同品种间肉牛遗传关系的方法。

3.2.3 预测误差相关系数法结果   与CD值类似,r值取值范围在0~1之间。本研究发现,10个地方品种r值范围在0.000 7~0.001 3之间。不同品种间,r值聚集程度较为紧密,且均接近于0,表明依据r值结果,10个肉牛地方品种间不存在遗传关联,这与PCA结果和真实情况差距较大,因此,不能真实有效反映品种间的亲缘关系。

3.3 不同遗传力对计算方法的影响

在本研究中模拟了3种遗传力的性状(对应实际应用中的低、中、高遗传力性状),来评估遗传力对不同计算亲缘关系方法的影响。基于LD一致性的亲缘关系计算方法仅依赖于基因组数据, 因此,性状的遗传力高低对其没有任何影响。而在PEVD、CD以及r的计算中,需要依赖预测所得的估计育种值计算相关参数, 其育种值估计准确性会受到遗传力的显著影响。如使用高遗传力性状评估品种间亲缘关系时,会提升育种值估计的准确性,降低了预测中的误差,进一步提升了预测得到的品种间亲缘关系(表 6)。但在选取较低遗传力性状时,育种值估计中的误差较大,导致了种间亲缘关系评估结果与实际不符。所以,使用低、中遗传力的性状计算得到的PEVD、CD以及r值并不能够准确描述品种间的亲缘关系。

4 结论

本研究对比了5种不同计算品种间亲缘关系的方法,其中以PCA聚类结果为参照,基于LD一致性的亲缘关系评估方法的评估结果与PCA聚类结果一致,且该方法能够使用皮尔逊相关系数量化品种间亲缘关系,具有较好的准确性。PEVD法、CD法与r法3种方法与上述方法相比评估群体间亲缘关系时容易受到性状估计育种值的误差方差影响,从而造成种间亲缘关系评估结果出现误差。因此,基于LD一致性的亲缘关系评估方法是一种较为适合评估肉牛地方品种种间亲缘关系的方法。

参考文献
[1]
MEUWISSEN T H E, HAYES B J, GODDARD M E. Prediction of total genetic value using Genome-Wide dense marker maps[J]. Genetics, 2001, 157(4): 1819-1829. DOI:10.1093/genetics/157.4.1819
[2]
MEUWISSEN T, HAYES B, GODDARD M. Genomic selection: A paradigm shift in animal breeding[J]. Anim Front, 2016, 6(1): 6-14. DOI:10.2527/af.2016-0002
[3]
HARRIS B L, JOHNSON D L, SPELMAN R J. Genomic selection in New Zealand and the implications for national genetic evaluation[C]//Identification, Breeding, Production, Health and Recording of Farm Animals. Proceedings of the 36th ICAR Biennial Session. Niagara Falls, USA: ICAR Technical Series, 2009: 325-330.
[4]
HAYES B J, BOWMAN P J, CHAMBERLAIN A C, et al. Accuracy of genomic breeding values in multi-breed dairy cattle populations[J]. Genet Sel Evol, 2009, 41(1): 51. DOI:10.1186/1297-9686-41-51
[5]
GAO H, SU G, JANSS L, et al. Model comparison on genomic predictions using high-density markers for different groups of bulls in the Nordic Holstein population[J]. J Dairy Sci, 2013, 96(7): 4678-4687. DOI:10.3168/jds.2012-6406
[6]
成海建, 姜富贵, 张清峰, 等. 全基因组选择技术在肉牛育种中的应用[J]. 中国牛业科学, 2018, 44(6): 68-72.
CHENG H J, JIANG F G, ZHANG Q F, et al. Application of genomic selection in beef cattle[J]. China Cattle Science, 2018, 44(6): 68-72. (in Chinese)
[7]
DUENK P, BIJMA P, CALUS M P L, et al. The impact of non-additive effects on the genetic correlation between populations[J]. G3 (Bethesda), 2020, 10(2): 783-795. DOI:10.1534/g3.119.400663
[8]
LEGARRA A, GARCIA-BACCINO C A, WIENTJES Y C J, et al. The correlation of substitution effects across populations and generations in the presence of nonadditive functional gene action[J]. Genetics, 2021, 219(4): iyab138. DOI:10.1093/genetics/iyab138
[9]
WIENTJES Y C J, VEERKAMP R F, CALUS M P L. The effect of linkage disequilibrium and family relationships on the reliability of genomic prediction[J]. Genetics, 2013, 193(2): 621-631. DOI:10.1534/genetics.112.146290
[10]
CLARK S A, HICKEY J M, DAETWYLER H D, et al. The importance of information on relatives for the prediction of genomic breeding values and the implications for the makeup of reference data sets in livestock breeding schemes[J]. Genet Sel Evol, 2012, 44(1): 4. DOI:10.1186/1297-9686-44-4
[11]
PSZCZOLA M, STRABEL T, MULDER H A, et al. Reliability of direct genomic values for animals with different relationships within and to the reference population[J]. J Dairy Sci, 2012, 95(1): 389-400. DOI:10.3168/jds.2011-4338
[12]
LUND M S, SU G S, JANSS L, et al. Genomic evaluation of cattle in a multi-breed context[J]. Livest Sci, 2014, 166: 101-110. DOI:10.1016/j.livsci.2014.05.008
[13]
周子文, 王雪, 丁向东. 基于高密度SNP标记估计群体间遗传关联[J]. 遗传, 2021, 43(4): 340-349.
ZHOU Z W, WANG X, DING X D. Measuring genetic connectedness between herds based on high density SNP markers[J]. Hereditas, 2021, 43(4): 340-349. DOI:10.16288/j.yczz.20-351 (in Chinese)
[14]
KENNEDY B W, TRUS D. Considerations on genetic connectedness between management units under an animal model[J]. J Anim Sci, 1993, 71(9): 2341-2352. DOI:10.2527/1993.7192341x
[15]
LALOË D. Precision and information in linear models of genetic evaluation[J]. Genet Sel Evol, 1993, 25(6): 557. DOI:10.1186/1297-9686-25-6-557
[16]
LEWIS R M, CRUMP R E, SIMM G, et al. Assessing connectedness in across-flock genetic evaluations[J]. Proc Br Soc Anim Sci, 1999, 1999: 121. DOI:10.1017/S1752756200002763
[17]
YU H P, SPANGLER M L, LEWIS R M, et al. Genomic relatedness strengthens genetic connectedness across management units[J]. G3 (Bethesda), 2017, 7(10): 3543-3556. DOI:10.1534/g3.117.300151
[18]
MA P P, HUANG J, GONG W J, et al. The impact of genomic relatedness between populations on the genomic estimated breeding values[J]. J Anim Sci Biotechnol, 2018, 9(1): 64. DOI:10.1186/s40104-018-0279-4
[19]
昝林森, 梅楚刚, 王洪程. 中国黄牛选育改良及肉牛种业发展建议[J]. 中国牛业科学, 2016, 42(6): 1-4.
ZAN L S, MEI C G, WANG H C. Breeding and improvement of Chinese yellow cattle and suggestions on beef cattle industry[J]. China Cattle Science, 2016, 42(6): 1-4. (in Chinese)
[20]
XU L, ZHU B, WANG Z Z, et al. Evaluation of linkage disequilibrium, effective population size and haplotype block structure in Chinese cattle[J]. Animals (Basel), 2019, 9(3): 83.
[21]
BRØNDUM R F, RIUS-VILARRASA E, STRANDÉN I, et al. Reliabilities of genomic prediction using combined reference data of the Nordic Red dairy cattle populations[J]. J Dairy Sci, 2011, 94(9): 4700-4707. DOI:10.3168/jds.2010-3765
[22]
高超群, 曹然然, 杜文苹, 等. 基于全基因组SNP标记分析中国地方鸡品种的遗传多样性和种群结构[J]. 畜牧兽医学报, 2023, 54(2): 554-562.
GAO C Q, CAO R R, DU W P, et al. Genetic diversity and population structure analysis of Chinese native chicken breeds using genome-wide SNPs[J]. Acta Veterinaria et Zootechnica Sinica, 2023, 54(2): 554-562. (in Chinese)
[23]
李隐侠, 牙生江·那斯尔, 赛里克·都曼, 等. SNP芯片评估柯尔克孜羊群体遗传多样性和遗传结构[J]. 畜牧兽医学报, 2023, 54(2): 572-583.
LI Y X, NASIER Y S J, DUMAN S L K, et al. Evaluation of genetic diversity and genetic structure in Kirgiz sheep population based on SNPs chip[J]. Acta Veterinaria et Zootechnica Sinica, 2023, 54(2): 572-583. (in Chinese)
[24]
DE A K, SAWHNEY S, MUTHIYAN R, et al. Legacies of domestication, Neolithic diffusion and trade between Indian subcontinent and Island Southeast Asia shape maternal genetic diversity of Andaman cattle[J]. PLoS One, 2022, 17(12): e0278681. DOI:10.1371/journal.pone.0278681
[25]
LIU B, TAO W K, FENG D H, et al. Revealing genetic diversity and population structure of endangered Altay white-headed cattle population using 100 k SNP markers[J]. Animals (Basel), 2022, 12(22): 3214.
[26]
刘晨龙, 卢丹, 周泉勇, 等. 利用高密度SNP芯片分析杭猪的群体遗传结构[J]. 畜牧兽医学报, 2022, 53(8): 2502-2513.
LIU C L, LU D, ZHOU Q Y, et al. Analysis of population genetic structure of hang pigs by high density SNP chip[J]. Acta Veterinaria et Zootechnica Sinica, 2022, 53(8): 2502-2513. (in Chinese)
[27]
KAWAGUCHI F, NAKAMURA M, KOBAYASHI E, et al. Comprehensive assessment of genetic diversity, structure, and relationship in four Japanese cattle breeds by Illumina 50 K SNP array analysis[J]. Anim Sci J, 2022, 93(1): e13770. DOI:10.1111/asj.13770
[28]
UEMOTO Y, SUZUKI K, YASUDA J, et al. Evaluation of inbreeding and genetic diversity in Japanese Shorthorn cattle by pedigree analysis[J]. Anim Sci J, 2021, 92(1): e13643. DOI:10.1111/asj.13643
[29]
杨祎挺, 甘麦邻, 刘杨, 等. 基于SNP芯片挖掘伍隍猪新遗传材料[J]. 中国畜牧杂志, 2022, 58(8): 178-188.
YANG Y T, GAN M L, LIU Y, et al. Mining new genetic materials for Wuhuang pigs based on SNP chips[J]. Chinese Journal of Animal Science, 2022, 58(8): 178-188. (in Chinese)
[30]
FOULLEY J L, HANOCQ E, BOICHARD D. A criterion for measuring the degree of connectedness in linear models of genetic evaluation[J]. Genet Sel Evol, 1992, 24(4): 315. DOI:10.1186/1297-9686-24-4-315
[31]
LALOË D, PHOCAS F, MÉNISSIER F. Considerations on measures of precision and connectedness in mixed linear models of genetic evaluation[J]. Genet Sel Evol, 1996, 28(4): 359. DOI:10.1186/1297-9686-28-4-359
[32]
YU H P, SPANGLER M L, LEWIS R M, et al. Do stronger measures of genomic connectedness enhance prediction accuracies across management units?[J]. J Anim Sci, 2018, 96(11): 4490-4500.
[33]
TARRÉS J, FINA M, PIEDRAFITA J. Connectedness among herds of beef cattle bred under natural service[J]. Genet Sel Evol, 2010, 42(1): 6.
[34]
XU L, WANG Z Z, ZHU B, et al. Theoretical evaluation of Multi-Breed genomic prediction in chinese indigenous cattle[J]. Animals (Basel), 2019, 9(10): 789.
[35]
KUEHN L A, NOTTER D R, NIEUWHOF G J, et al. Changes in connectedness over time in alternative sheep sire referencing schemes[J]. J Anim Sci, 2008, 86(3): 536-544.
[36]
ZHANG S Y, OLASEGE B S, LIU D Y, et al. The genetic connectedness calculated from genomic information and its effect on the accuracy of genomic prediction[J]. PLoS One, 2018, 13(7): e0201400.

(编辑   郭云雁)