基于高密度SNP 芯片评估“科尔沁肉牛”遗传背景

林雨浓，马万欣，包玲玲，何曙光，石顺利，张秋生，赵澈勒格日，高会江，李俊雅*，王泽昭*

（1.中国农业科学院北京畜牧兽医研究所，北京 100193；2.内蒙古自治区通辽市畜牧业发展中心，内蒙古通辽 028000）

肉牛产业是内蒙古自治区经济发展的支柱产业，是农牧民增收的重要途径。“科尔沁牛”是以西门塔尔牛为父本，蒙古牛及三河牛和蒙古牛的杂种牛为母本而培育成的乳肉兼用品种。作为通辽地区地方畜牧业的重要资源，具有适应性强，宜牧，耐粗饲，耐寒，抗病力强等特点。然而，随着肉牛市场需求的增加和养殖业的发展，“科尔沁牛”在生产性能和经济效益方面遇到了瓶颈。为了提升群体肉用性能，当地农业部门利用北美肉用型西门塔尔牛为父本，本地“科尔沁牛”为母本进行杂交改良，历经30 年，经过“杂交改良、横交固定和选育提高”三个阶段，培育出肉用型“科尔沁肉牛”群体。

近年来，群体遗传背景分析成为评估和研究动物群体遗传结构的重要工具。随着SNP 芯片和测序技术的发展，有助于人们更全面地了解不同群体之间的遗传多样性和遗传联系。Eva M.Strucken 等[1]采用777K 芯片对印度本地的15 个牛种进行遗传多样性与群体结构研究，发现印度瘤牛比印度本地普通牛有更高的遗传多样性，他们推断本地普通牛较小的有效群体导致了群体基因频率受到遗传漂变的影响更大，而低有效群体数量增加了近交衰退发生的概率，使其在选择的过程中遗传多样性降低速度快于印度瘤牛。在动物育种过程中，保持一定的遗传多样性有利于提高群体的环境适应性[2]。Alejandro Amaya 等[3]对约2.7 万头哥伦比亚西门塔尔牛（1975—2017年）的数据进行了分析，发现该种群的近交系数由5.06%下降到了2.25%，认为是得益于欧洲和美系西门塔尔牛之间的基因交流，才使该种群的遗传多样性得以提升。

我国疆域辽阔、资源丰富，但目前由于经济效益等问题致使国外种源占有率居高不下。物种作为战略性资源、是农业的基石，但是目前关于国内地方品种的研究存在着许多空白和滞后。马均等[4]研究发现，我国自主培育的专门化肉牛新品种—华西牛与蒙古牛、三河牛和其他家系西门塔尔牛出现不同的比例的群体结构组成，形成了自己独特的遗传结构特征。刘正喜等[5]在2022 年使用全基因组重测序数据对渤海黑牛的群体结构、遗传多样性、亲缘关系分析后发现，渤海黑牛与郏县红牛和鲁西牛的关系密切，该研究为逐年减少的渤海黑牛群体的种质资源保护和品种发展提供了关系的理论基础。

本研究目的是分析“科尔沁肉牛”、华西牛和美系西门塔尔牛这三个群体之间的遗传关联，评估“科尔沁肉牛”改良育种在遗传层面的影响，旨在为“科尔沁肉牛”新品种培育提供遗传背景支持。

1 材料与方法

1.1 试验群体及SNP 分型

“科尔沁肉牛”（Kerqin beefcattle，“KBC”），来自内蒙古科尔沁牛业股份有限公司的437 头核心群牛只；美系西门塔尔牛（American Simmental cattle，ASC）来自于通辽京缘种牛繁育有限责任公司的25 头种牛；华西牛（HuaXi cattle，HXC），来自内蒙古奥科斯牧业有限公司的55 头核心群牛只。

每个样本使用含有EDTA 抗凝剂的真空采血管采血5mL 送样检测，并取2mL 血样进行DNA提取。使用Illumina BovineHD BeadChip 芯片进行变异检测，该芯片含有覆盖牛基因组的777962个标记位点。

1.2 数据质控

使用Genome Studio 软件进行基因分型。通过PLINK v1.9[6]软件提取出常染色体后进行数据质量控制（Quality control，QC）。QC 顺序如下：①删除等位基因检出率（Call rate）小于90%的位点；②删除最小等位基因频率（Minor allele frequency，MAF）小于0.01 的位点；③删除哈代温伯格平衡检验P 值小于1×10 箒6 的SNP 位点；④删除SNP 缺失率大于10%的样本个体。通过QC 后数据用于后续分析。单独使用PLINK 中的--indep-pairwise 参数对所有个体的连锁位点进行过滤，剩余位点进行群体结构分析。过滤后的个体数见表1。

表1 试验样本数量

1.3 遗传多样性分析

对QC 后数据统计群体期望杂合度（Expected heterozygosity，He）、群体观察杂合度（Observed heterozygosity，Ho）、MAF 分布情况。 MAF 分布情况是按照不同MAF 大小将SNP 标记分为5 组，分别是（MAFSNP<0.1）、（0.1≤MAFSNP＜0.2）、（0.2≤MAFSNP＜0.3）、（0.3≤MAFSNP＜0.4）和（0.4≤MAFSNP），统计3 个群体SNP 标记的分布情况。

1.4 群体分析

使用GCAT 软件进行了群体主成分（Principal component analysis，PCA）计算，输出数据导入到R 中，使用ggplot2 包[7]进行可视化。随后使用邻接法（Neighbor-Joining，NJ）构建进化树[8]：从成对FST 值的距离矩阵生成邻居网络树，并在Splitstree 4.14.5[9]中实现建树。使用Admixture v.1.3 软件的默认设置进行群体结构分析。设置了10 倍交叉验证选择最合适的K 值[10]。使用PopLDdecay 软件[11]计算了所有数据的LD 衰减距离，并绘制了LD 衰减图。

2 结果

2.1 群体遗传多样性

由表2 可知，经过质控后的“KBC”群体具有最高的平均MAF 值（0.28），ASC 群体的平均MAF 值最低（0.23）。平均He 的范围取值从0.34（HXC）到0.37（“KBC”）。平 均Ho 最高为“KBC”（0.37），最低的是ASC（0.36），HXC和ASC 的平均Ho 略大于He，“KBC”的Ho 和He 几乎相等。

表2 群体遗传多样性参数

由图1 可知，“KBC”群体被过滤掉了最多的遗传标记。“KBC”群体过滤掉了81778 个MAF 小于0.01 的标记。QC 后，ASC 群体拥有最多的低频标记（MAF＜0.1）。随着MAF 的增大，HXC 和“KBC”群体内中高MAF 的SNPs 多于ASC 群体。在三个群体中，“KBC”具有最多的MAF≥0.3 的SNP 个数，具体数量为（413400），其次是HXC（363989），最少的是ASC（344174）。

图1 不同群体的MAF 分布

2.2 PCA 和进化树

由图2（a）可知，三个群体主要聚集成为两个簇，将“KBC”群体和HXC 群体分开。此外，“KBC”群体与ASC 群体中部分个体出现明显的聚集和重叠。根据群体在图中的空间位置，可以判断出“KBC”和HXC 的遗传距离较远，和ASC 之间遗传距离较近。

图2 三个牛群体的主成分分析结果和邻居连接系统发育树

为了观察三个群体在进化距离上的相对位置，本研究构建了环形进化树（图2b）。“KBC”群体和HXC 群体分别位于进化树的起点位置和末端位置，表明这两个群体的进化距离相对较远。ASC 群体和HXC 群体处于两个比较靠近的分支，暗示这两个群体经历了不同的进化路径，但是其分化路径相对较近。ASC 群体处于“KBC”群体的一个分支下，表明这两个群体可能具有共同的祖先，使得在基因层面上具备一定的相似性。

2.3 种群结构分析

通过假设群体结构数的先验值K，推断群体内亚群个数。通过对比不同祖先群体遗传成分比例观察“KBC”、ASC 和HXC 群体间遗传相似性。由图3 可知，在K=2 时，“KBC”的祖先比例和HXC 相似，少数个体祖先比例与ASC 相似；当K=3 时，“KBC”群体和ASC 群体部分个体祖先比例相似，此时HXC 群体出现与“KBC”群体和ASC 群体截然不同的祖先成分比例；随着K 值不断增大，“KBC”群体和HXC群体差异的趋势逐渐增大；当K=5 时，“KBC”出现了更为复杂的遗传比例；K=2～4 时，ASC群体的祖先组成比例变化不明显，“KBC”的祖先比例变化复杂。

图3 当K=2～5 时，基于LD 过滤的三种牛种群的SNPs 的结构分析

2.4 连锁不平衡

如图4 所示，三个群体的LD-decay 曲线中“KBC”群体衰减速度最快，ASC 的LD 衰减速度较慢，HXC 群体的LD 衰减速度位于两者之间。三个群体的远端SNP 之间的连锁水平出现差异，在≥200Kb 之后并未衰减到相近的r2范围。随着距离逐渐增大，三个群体的LD 衰减曲线逐渐平缓，r2结果始终保持是r2ASC≥r2HXC≥r2KBC。

图4 三个群体的连锁不平衡衰减图

3 讨论

3.1 遗传多样性参数

遗传多样性是生物多样性的核心维度之一，反映了群体内部基因的多样性水平。它在维持种群的健康、适应性和进化潜力方面起着关键作用[12,13]。目前，通过基因芯片和测序技术辅助分析种群遗传进展成为热点研究内容[14]，通过评估种群遗传多样性对于制定未来的育种目标至关重要[2]。本研究对“KBC”群体、HXC 群体和ASC群体的基因型数据进行分析，通过对群体遗传多样性的统计，以期发现“KBC”在育种过程中的遗传变化。研究发现，“KBC”的MAF、He 和Ho 分别为0.28、0.37 和0.37，在三个群体中均为最高，并且“KBC”群体He 和Ho 的一致性在三个群体中为最高。这一结果反映出“KBC”群体内的遗传多样性比ASC 和HXC 群体较高，同时暗示着近亲交配的水平较低[15]。种群遗传多样性越高，意味着进化潜力越大，对环境变化的反应能力越强[16]。此外，人工选择也会极大的影响群体的多样性[17]。“KBC”是以“科尔沁牛”（三河牛×西门塔尔牛）为母本，以引进的美系西门塔尔牛为父本杂交培育而成的群体，与ASC 相比，其培育期间的选择强度相对较低，使其还保持着一定的遗传改良潜力。Maiorano AM 等[18]研究发现，亲本自交系杂交产生的杂种优势赋予了后代更高的遗传多样性，并加强了杂种个体优势性状表现[19,20]。“KBC”群体在遗传改良的过程中，保留了“科尔沁牛”群体的适应性强的同时，借助ASC 群体的基因交流提升了肉用性能。在哈迪-温伯格平衡下，“KBC”的Ho 大于He，说明杂交系统仍然是扩大品种遗传变异的重要方法。

3.2 PCA 和进化树

PCA 是利用解释群体最大遗传方差的特征值向量对群体进行分离[21]。本结果表明，“KBC”和HXC 群体各自聚集成了两个簇并保持了一定的距离。同时，“KBC”群体中部分个体与ASC群体出现明显的聚集，这些“重叠”，标志着个体之间的基因型具有一定的相似性。这种相似性是由于群体之间存在相似的遗传背景或近亲关系导致了其共享某些遗传特征或具有相似的基因型组合。

通过构建进化树，发现“KBC”群体和ASC群体进化距离较近，与HXC 群体距离较远。需要注意的是，环形进化树是无根树，仅提供了群体之间遗传和进化距离的相对关系。群体在进化树上的位置不代表起源和进化方向，而是反映了在进化距离上的相对位置。最初“KBC”群体在培育中使用了ASC 改良其肉用性能，在基因组上留下了选择足迹（Selective Footprint）[22]。经过长期的人工选择，增加了“KBC”群体和ASC 群体的遗传相似性，因此其在进化树上的位置处于同一分支。后续研究如果要全面地理解群体间的关系和演化过程，需要分析表型特征信息或分子钟估计[23]等进行更深入的研究。

3.3 群体结构分析

群体遗传结构分析基于遗传标记信息将种群分为若干亚群，对于研究群体和个体间的遗传关联具有重要意义[24]。揭示了群体间的遗传交流程度、迁移历史以及潜在的遗传壁垒[25]。研究发现，在K=2～3 时“KBC”群体内部分个体与ASC 祖先组成比例十分相似，在K=2～4 时，ASC 群体的组成比例比较稳定，而“KBC”群体祖先比例变化明显。Ocampo 等[26]发现不同品种之间的分化模式和遗传关系与每个品种的进化历史高度一致。“KBC”群体和HXC 群体在杂交改良阶段都引入过西门塔尔的血统，因此，当K值比较小的时候，更容易地捕捉到群体之间的主要遗传差距。在“KBC”群体内，不同个体之间的祖先组成比例也在不断变化，一方面是其原始群体本身的多态性所引起；另一方面，经过不断的杂交改良，致使“KBC”群体存在多祖先的基因渗入。

本研究中，在假设K 值范围内并未出现交叉验证（Cross-Validation，CV）错误率的拐点，可能是研究中ASC 群体和HXC 群体使用的样本数量较小，捕捉三个群体间的遗传分化效果较差。当K 大于4 时，三个群体祖先比例出现明显不同，且个体祖先比例变得相对混乱。这是因为，随着K 值增大，模型捕捉群体之间细微的遗传差异的灵敏度也随之提高。

3.4 LD 衰减

通过对比群体之间的连锁不平衡水平，对于合理利用种质资源至关重要。三个群体中，“KBC”群体衰减的最快，ASC 群体的LD-decay衰减的最慢。而人工选择的效果可以反映在每个群体的基因组连锁不平衡水平上[27]。在物种遗传改良过程中，在高强度的选择下，家畜群体中有利的染色体片段频率迅速增加[28]，并促进群体中LD 水平的增加[29,30]。当群体中LD 水平增加，提高后续选择的准确性，可以更快地改善目标性状。HXC 群体已经于2021 年12月通过国家畜禽遗传资源委员会审定，被认定为肉牛新品种。而“KBC”群体目前还处于新品种培育的过程中。因此，“KBC”群体内远端标记的连锁效果比其他两个群体较差。连锁不平衡分析可以反映相应的独特育种选择历史和种群结构[27,31,32]，同时，也对后续基因组信息挖掘提供指示。值得注意的是，在研究中随着标记之间的距离逐渐增大，三个群体的LD 衰减逐渐平缓，但是一直存在差距。后续如要进行GWAS 分析时，为了保证准确性，“KBC”群体可能会需要比HXC 和ASC 更多的基因组信息。

4 结论

本研究使用遗传信息数据评估了“KBC”、HXC 和ASC 群体的遗传多样性和种群结构。“KBC”群体内中高频（MAF≥0.3）的标记数量最多，且具有最高的遗传多样性和最低的连锁不平衡衰减速度。在进化方向上，“KBC”群体与HXC 群体的距离更远，但部分个体的种群组成与ASC 群体分化较差。通过进行群体遗传背景分析，发现“科尔沁肉牛”经过杂交改良、横交固定、选育提高三个阶段的持续选育，在遗传方面（分子水平）与华西牛和美系西门塔尔牛品种存在明显差异，已具备形成新培育品种的先决条件。