牦牛是重要的大型哺乳动物[1]，主要分布在3 000～5 000 m的高海拔地区[2]。作为我国高海拔地区畜牧业中的优势畜种[3]，牦牛为应对高原地区的极端环境和天然草场的季节性变化形成了适应能力强和生命力顽强等的各种适应机制[4]，牧民的衣食住行都与其紧密相关[5,6]。由于牧民对牦牛品种繁殖等知识相对缺乏，饲养管理方式落后，造成了牦牛的选育水平普遍不高[7]，进而限制了牦牛的生长发育，这就使得牦牛的经济价值不能得到应有的体现。目前《国家畜禽遗传资源品种名录（2021年版）》收录的地方牦牛品种共有18个，培育牦牛品种2个。当前，随着第三次全国畜禽遗传资源普查工作的进行，会有更多的地方牦牛品种被鉴定收录。青海省位于“世界屋脊”青藏高原的东北部，省内山脉高耸，地形多样[8]，这种特殊的地理位置和多变的气候，孕育了多个牦牛遗传资源。祁连县隶属于青海省海北藏族自治州，位于祁连山中部的群山峻岭和盆地之中[9]，平均海拔3 169 m，属于高原大陆性气候，地貌多样，适合畜牧业发展[10]。祁连牦牛作为祁连县境内的牦牛类群之一，具有区别于其他牦牛类群的独特外貌特征和生产性能。祁连牦牛体质坚实、结构紧凑匀称、体躯深长、头粗而长（公牛粗短）、多有角；其肉色深红、呈明显的大理石纹。在当地独特的地理环境和气候下，当地饲养的家牦牛与祁连山型野牦牛进行了较为频繁的基因交流，从而形成了一个独特的牦牛遗传资源———祁连牦牛[11]。

全基因组重测序是在已知基因组序列的前提下，对物种群体内不同个体的基因组序列进行测序，然后与参考基因组进行比对，根据比对信息挖掘群体或个体间的变异信息[12]，进而解析群体结构，挖掘重要基因等。近年来，随着高通量测序技术的不断发展，全基因组重测序也在各领域得到了广泛的应用，并逐渐代替了其他测序手段[13,14]。目前，在牛[15,16]和羊[17,18]等各类畜禽中已经有大量使用基因组重测序的方法与研究。如Kang等[19]对6个牦牛群体进行了重测序研究，结果显示6个牦牛群体存在明显差异，并认为科才牦牛是中国牦牛的一个独立群体。Chai等[20]对我国31个牦牛群体进行了群体基因组学研究，发现牦牛最早的驯化发生在青藏高原东南部，然后向青藏高原西部和东北部扩散到四川、甘肃和青海。鲍麒等[21]利用全基因组重测序评估了肃南牦牛的遗传结构，并注释了多个与组织分化、胚胎发育和肉品质等相关的选择区域。李在文等[22]利用GBS构建了麦洼牦牛3个保种群的系谱，为麦洼牦牛的后续保种工作提供了基础。

当前，已有研究者对其他牦牛群体进行了全基因组水平的研究[23,24,25]。但是对祁连牦牛和其他牦牛资源在基因组水平上还未进行系统的研究。因此，本研究基于全基因组重测序对不同牦牛群体进行比对，同时结合群体基因组学分析构建牦牛群体基因组特征图谱，解析祁连牦牛与其他牦牛群体的差异意义重大。

1、材料与方法

1.1 样本采集

本研究选择九龙牦牛（Jiulong yak）、天祝白牦牛（Tianzhu White yak）、斯布牦牛（Sibu yak）、帕米尔牦牛（Pamir yak）、木里牦牛（Muli yak）、雪多牦牛（Xueduo yak）、祁连牦牛（Qilian yak）和肃南牦牛（Sunan yak）共8个牦牛类群的46个样本进行测序分析，样本具体信息详见表1。本研究采集所有牦牛个体静脉血10 m L，于含有EDTA的一次性负压采血容器内，放入-20℃冰箱备用。提取牦牛血液DNA时使用Easy Pure血液基因组DNA试剂盒（全式金，北京），使用Qubit荧光定量仪检测样本DNA浓度，1%琼脂凝胶电泳检测样本DNA完整性，检测合格的样本在illumina PE150平台上机测序。 

表1 样本基本信息  

1.2 全基因组序列比对与遗传变异检测

因为在血液样本DNA提取、建库和测序之后生成的数据会存在误差，并且误差的存在会对后续分析产生影响。因此，为了保证后续分析的准确性，对上述步骤产生的数据误差进行了质量评估和过滤。使用Trimmomatic[26]软件进行质控，GATK软件进行SNP鉴定，并用软件自带参数（Select Variants和Variant Filtration）对SNP进行过滤[27]。

1.3 群体基因组学分析

1.3.1 主成分分析和系统进化树构建

为了将个体按照主成分进行分类，了解研究群体的聚类和差异，采用主成分分析进行聚类，使用PLINK[28]软件对常染色体SNP数据进行分析。使用bcftools提取vcf文件子集，将vcf转换为plink格式，得到3个以.map、.nosex、.ped结尾的文件，基于.ped生成bed二进制文件，然后得到2个以.bim、.bed结尾的文件，之后进行PCA分析，获得2个.eigenval、.eigenvec结尾的文件，其中.eigenvec记录特征向量，用于坐标轴的绘制。最后利用R语言的ggplot2包可视化绘制PCA图。为了解群体间个体的谱系关系，使用PLINK软件计算遗传距离矩阵，采用Neighbor-Joining法构建系统进化树，使用i TOL在线软件绘制系统进化树。

1.3.2 杂合度和近交情况分析

为了评估牦牛的种质资源现状，使用ml Rho[29]软件计算牦牛个体的杂合度以及近交程度。

1.3.3 群体分化指数分析

为了进一步衡量群体间的分化程度，用VCFtools[30]软件采用滑动窗口法进行Fst分析，最后使用R语言中的corrplot安装包对结果进行可视化。

1.3.4 连续纯合子区域（ROH）分析

在动物进化过程中，ROH会在动物基因组上呈现出不同的分布模式，因此动物基因组上ROH的数目、长度和分布可以揭示动物种群的进化历史等。本研究采用PLINK[28]软件对ROH进行检测。

2、结果与分析

2.1 不同牦牛类群的聚类分化情况

主成分分析作为群体遗传学中常用的分析方法，可以将许多相关性很高的变量转化成彼此相互独立或不相关的变量，进而通过结果了解群体的聚类情况。由图1A可知，本研究牦牛群体能够得到明显的区分，说明这几个牦牛群体之间很少有杂交现象，并且群体内样本相对较聚集，群体间较分散。其中祁连牦牛与肃南牦牛和雪多牦牛2个群体的遗传距离较近，与帕米尔牦牛、木里牦牛、九龙牦牛和天祝白牦牛4个牦牛群体的遗传距离相对较远。为了解不同牦牛的谱系关系，本研究计算了个体间的遗传距离矩阵，并采用邻接法构建系统进化树。图1B显示，所有牦牛群体单独聚在一起，未出现牦牛个体混杂现象，这与前面的主成分分析结果相一致。其中斯布牦牛与帕米尔牦牛聚为一类，这也与两个牦牛群体的地理位置相符合；同处在四川省的木里牦牛和九龙牦牛聚为一类；位于甘肃省境内的天祝白牦牛和肃南牦牛聚为一类，然后与祁连牦牛聚为一大类。

图1 不同牦牛群体的聚类分化情况   

2.2 不同牦牛类群的杂合度和近交系数分析

杂合度和近交系数分析结果见表2，几个牦牛群体的理论杂合度在0.354 4～0.360 0之间。肃南牦牛的观测杂合度最大，木里牦牛的观测杂合度最小，分别为0.346 3和0.311 2。所有牦牛群体的观测杂合度均小于理论杂合度。通过对牦牛群体的近交系数计算发现，8个牦牛群体的近交系数在0.022 9～0.136 0之间，其中木里牦牛最大，肃南牦牛最小。

2.3 不同牦牛类群的群体分化指数

在前面的分析中，本研究使用的8个牦牛群体可以得到明显的区分，两个分析结果相一致。但为了进一步对群体分化程度数字化，更加清晰地了解群体间的分化程度，基于SNP数据计算了各牦牛群体间的群体分化指数。目前认为遗传分化指数≥0.04的牦牛群体，可以拆分为新的遗传品种资源，结果见图2。由图2可知，九龙牦牛（JL,0.08）、木里牦牛（ML,0.08）、帕米尔牦牛（PMR,0.08）、祁连牦牛（QL,0.05）、斯布牦牛（SB,0.07）、肃南牦牛（SN,0.07）、天祝白牦牛（TZ,0.06）和雪多牦牛（XD,0.07）这8个牦牛群体的分化指数均大于0.04，这与上述分析结果相一致，说明祁连牦牛可以作为一个单独的遗传资源，这为祁连牦牛的科学分类和发展提供了科学依据。

表2 各牦牛群体的杂合度和近交系数  

图2 不同牦牛群体的遗传分化指数   

2.4 不同牦牛类群的连续纯合子片段分析

在本研究的8个牦牛群体中共检测到ROH 15 704个，并发现ROH的数量和长度呈正相关，介于0.91～1.00之间，结果见图3A。其中帕米尔牦牛群体检测到的ROH数量最大，肃南牦牛最少，分别为4 566个和559个。将ROH的长度进行分类（<1 Mb;1～2 Mb;>2 Mb),ROH数目<1 Mb大小的占绝大部分，其次为1～2 Mb大小的；8个牦牛群体（九龙牦牛，天祝白牦牛，斯布牦牛，帕米尔牦牛，木里牦牛，雪多牦牛，祁连牦牛，肃南牦牛）ROH>2 Mb类别的数目分别为0、1、3、3、4、0、0和0个，结果见图3B。检测到ROH的数量排序为：帕米尔牦牛>木里牦牛>天祝白牦牛>祁连牦牛>斯布牦牛>九龙牦牛>雪多牦牛>肃南牦牛。

图3 不同牦牛群体的ROH分析  

3、讨论

种群群体遗传结构的确定对于家养动物与野生动物来说都很重要，这关系到他们的育种和种群保护[31]。有大量的研究人员使用微卫星标记研究了牦牛群体的遗传结构[32,33,34]。本研究通过群体基因组学构建牦牛群体特征图谱发现：所有牦牛群体能够得到区分，群体之间较为分散，没有出现混杂情况。系统进化树与PCA分析结果相一致。通过这两个分析可以将牦牛群体进行初步的区分。遗传多样性也称为基因多样性，是某一生物品种或群体内个体之间遗传变异的总和，是评估种质资源现状的重要指标[35]，群体的遗传多样性越高，对环境的适应能力就越强[36]。同样杂合度和近交系数也可用来评估群体的遗传多样性，研究中的理论杂合度是通过计算所得，而观测杂合度则是两个样本等位基因不相同的概率。近交系数是指基因纯化的程度。相关研究发现，群体的杂合度与生长速率和繁殖力等呈正相关[37]。本研究结果表明，8个牦牛群体的观测杂合度略低于期望杂合度，但相差并不大，说明虽然群体内受到了一定程度的选择压力，但尚未发生遗传漂变等。其中肃南牦牛的观测杂合度最高（0.346 3），但近交系数最低（0.022 9）；而木里牦牛的观测杂合度最低（0.311 2），近交系数最高（0.136 0）；各牦牛群体近交系数顺序为：木里牦牛>帕米尔牦牛>祁连牦牛>斯布牦牛>九龙牦牛>天祝白牦牛>雪多牦牛>肃南牦牛，说明木里牦牛群体的近交程度最大，肃南牦牛群体的近交程度最小。本研究中8个牦牛群体的杂合度与王佟等[23]的研究结果基本相符，但低于李铎等[33]的研究结果。推测牦牛群体的遗传多样性是波动的，也可能与研究所用的位点数目、样本数量及采样地区等多因素相关，这仍需进一步的深入研究。

群体分化指数的范围在0～1之间，数值越大，说明群体间的分化程度就越大。本研究中祁连牦牛的群体分化指数为0.05，与其他牦牛群体间存在明显的遗传分歧，可作为潜在的牦牛群体遗传资源，为今后大力发展该牦牛群体的资源挖掘和开发保护工作提供理论基础。相关研究表明，ROH在基因组中的大小和分布受到近亲交配和人工选择等因素的影响[38]，目前ROH分析常用来计算评估近交系数，并且这种方法具有较高的准确度[39]。本研究对8个牦牛群体的ROH分析结果显示，帕米尔牦牛的ROH数目最多，可能是因为帕米尔牦牛经过了近交选育，这与其群体的近交系数较高相一致。畜禽种质资源是种业振兴的基础，对牦牛群体的科学保护和优质资源的开发利用，可以有效地保护种群内的基因多样性，促进品种的改良和创新。总之，本研究为牦牛群体的资源保护提供了科学依据。

4、结论

综上，8个牦牛群体之间没有过多的基因交流，其中祁连牦牛是在长期的进化发育过程中产生的一个独特类群。本研究从全基因组水平揭示了祁连牦牛的系谱及遗传结构，为后续祁连牦牛遗传资源的保护及开发利用提供了理论依据。

参考文献:

[4]丁维芹,孙永刚,刘亚倩,等.冷季不同饲养方式下牦牛皮下脂肪转录组分析[J/OL].甘肃农业大学学报,2023:1-15[2023-08-06].

[6]都吉,和嘉华,和云凤,等.“牧繁农育”模式在中甸牦牛生产中应用的可行性分析[J].当代畜禽养殖业,2023,43(2):34-36.

[7]韩雄国,周彪.牦牛高效养殖现状及对策建议[J].四川畜牧兽医,2023,50(7):11-12.

[8]李广祯.青海牦牛全基因组遗传多样性与选择信号研究[D].西宁:青海大学,2023.

[9]杨振菲,马志杰,陈生梅,等.青海省祁连牦牛的母系遗传多样性及群体遗传结构分析[J].青海畜牧兽医杂志,2022,52(1):1-6.

[10]杨涛,罗增海.以草地为载体的种畜经济发展模式探讨———以海北州祁连县为例[J].中国畜禽种业,2018,14(10):36-37.

[11]洛藏旦增,韩学平,郭继军,等.祁连牦牛遗传资源调查与屠宰测定报告[J].青海畜牧兽医杂志,2022,52(1):44-45,49.

[12]李晓龙.基于全基因组重测序数据挖掘绵羊产奶和产羔性状的关键基因[D].兰州:甘肃农业大学,2022:4.

[13]宋明坤.基于全基因组重测序鉴定影响家兔生长性状的候选基因及验证[D].郑州:河南农业大学,2022:8.

[18]刘家鑫,魏霞,邓天宇,等.绵羊全基因组ROH检测及候选基因鉴定[J].畜牧兽医学报,2019,50(8):1554-1566.

[21]鲍麒,包鹏甲,马晓明,等.肃南牦牛群体遗传结构､选择信号分析和ROH检测[J].中国畜牧兽医,2022,49(8):2992-3005.

[22]李在文,李响,李小伟,等.基于GBS简化基因组测序数据重建麦洼牦牛保种群系谱[J].畜牧兽医学报,2023,54(9):3710-3721.

[23]王佟,马晓明,李婧,等.基于全基因组重测序对7个牦牛类群的群体结构分析[J].中国草食动物科学,2023,43(3):9-13,26.

[24]李广祯,马志杰,陈生梅,等.野牦牛及青海地方牦牛品种全线粒体基因组母系遗传多样性､分化及系统发育分析[J].畜牧兽医学报,2022,53(5):1420-1430.

基金资助:祁连牦牛､白藏羊遗传资源申报项目;中央高校基本科研业务费专项资金(31920210160-01);现代肉牛牦牛产业技术体系(CARS-37);甘肃省基础研究创新群体项目(20JR5RA580);

文章来源:王佟,马晓明,余道宁等.全基因组重测序解析不同牦牛类群的系谱及遗传结构[J].中国草食动物科学,2023,43(06):1-6.