肃南牦牛群体遗传结构、选择信号分析和ROH检测

鲍 麒，包鹏甲，马晓明，吴晓云，孟光耀，褚 敏，曹红梅，郎建英，安玉峰，梁春年，阎 萍

(1.中国农业科学院兰州畜牧与兽药研究所，兰州 730050；2.农业农村部青藏高原畜禽遗传育种重点实验室，兰州 730050；3.甘肃省牦牛繁育工程重点实验室，兰州 730050；4.肃南裕固族自治县畜牧兽医服务中心，肃南 734400)

牦牛(Bosgrunniens)在生物学分类上属于偶蹄目(Artiodactyla)牛科(Bovidae)牛亚科(Bovinae)牛属(Bos)，主要分布于青藏高原及其毗邻高山、亚高山高寒地区[1]。目前，全世界现有牦牛约2 200万头，其中95%以上分布在中国[1]。据国家畜禽遗传资源品种名录(2021年版)，中国牦牛现有18个地方品种以及2个培育品种，主要分布在中国西藏、青海、甘肃、四川等地[2]。肃南牦牛是青藏高原型牦牛的地方类群，主要分布于张掖市肃南裕固族自治县，群体数量约6万头。肃南牦牛毛色以黑褐色为主，体型外貌上带有野牦牛的特征，是肃南地区唯一在培的牦牛品种。从全基因组层面解析肃南牦牛的遗传变异，有助于进一步了解这一品种牦牛的特性以及挖掘经济性状相关的候选基因。

群体结构的鉴定是群体基因组学的研究基础，探索不同群体之间的系统发育关系、聚类分析结果和祖先成分组成，可以了解不同物种的分类地位和分化程度，在群体亲缘关系鉴定、物种起源分析、种群历史等研究中发挥作用[3-4]，为下游分析提供背景信息。选择信号(selection signal)包括选择性清除和搭车效应，是挖掘与表型相关基因的重要手段[5]。随着高通量测序技术成本的不断降低，利用选择信号解析物种之间性状差异的研究日益普遍。刘真[6]利用遗传分化指数(Fst)法对短脂尾的蒙古羊和瘦尾的藏羊群体进行选择信号分析，筛选出50个与脂类代谢相关的基因，如过氧化物酶体增生激活受体γ(PPARG)、视黄酸受体γ(RXRG)、溶质载体家族27成员2(SLC27A2)、酰基辅酶A合成酶长链家族成员6(ACSL6)等。Shen等[7]利用跨群体复合似然比(XP-CLR)、Fst和θπ比值(θπratio)方法对延边牛和非洲热带牛的重测序数据进行选择信号分析，筛选出3个与冷适应性相关的候选基因(皮质醇稳定蛋白(CORT)、成纤维细胞生长因子5(FGF5)、脂肪酸合酶(CD36))，从基因组层面解析了延边牛的耐寒机制。吴林慧等[8]在对恩施黑猪的选择信号研究中发现一系列基因分别与胴体长度、精子形成、骨骼肌分化和总产仔数相关，包括溶血磷脂酸受体2(LPAR2)、NADH 脱氢酶1α复合体亚基(NDUFA13)、肌细胞增强因子2B(MEF2B)和芳烃受体基因(AHR)。连续纯合片段(ROH)是指二倍体生物基因组中出现的连续纯合区域[9]，能够用来评估家畜近交系数、推断种群历史、辅助解析某些复杂性状的形成机制[10]。Ferenakovi等[11]利用4个家牛(瑞士褐牛、德国弗莱维赫牛、挪威红牛、蒂罗尔灰色牛)SNP数据，基于ROH长度和系谱估计分别计算近交系数FROH与FPED，结果发现ROH长度>8和16 MB时估计的FROH和FPED接近，且利用FROH评估各牛种之间近交水平的结果与FPED评估的结果一致。Wu等[12]在对20头皖南黑猪ROH分析中鉴定到1 913个ROHs，利用ROH长度计算的平均基因组近交系数为0.5234，呈现高度近交的状态，揭示了皖南黑猪目前出现繁殖性能衰退的原因。莫家远等[13]对包括东山猪在内的4个广西地方猪种进行ROH检测，注释的候选基因主要富集在与免疫、繁殖和肉品质相关的通路。Xu等[14]对金华猪群体筛选出的高频ROH区域进行注释，鉴定到与生殖相关的基因(同源框A3(HOXA3)、HOXA7、HOXA10和HOXA11)，肉质相关基因(成肌细胞决定基因1(MYOD1)、脂素家族蛋白3(LPIN3)和β-连环蛋白样蛋白1(CTNNBL1))、食欲相关基因(核连蛋白2(NUCB2))和抗病力相关基因(黏蛋白抗原4(MUC4)、MUC13、MUC20、Leishmanolysin-like肽酶(LMLN)、整合素蛋白β5(ITGB5)、心脏发育蛋白1(HEG1)、溶质载体家族12成员8(SLC12A8)和肌球蛋白轻链激酶(MYLK))。

与基因芯片技术相比，全基因组测序检测精度高，对基因组完整度依赖性小，能发现更多单核苷酸多态性(SNP)或稀有位点，因此能找到更短的ROH片段，有利于检测亲缘关系较远的物种之间的杂合度，避免出现ROH假阳性情况[15]。目前，关于肃南牦牛在全基因组范围内的群体遗传结构解析以及功能基因挖掘的研究鲜见报道。基于此，本研究对肃南牦牛以及巴州牦牛、九龙牦牛、斯布牦牛、天祝白牦牛4个著名牦牛品种进行10×全基因组重测序，对其进行遗传结构、ROH以及选择信号分析，比较不同牦牛之间的异同，以期为肃南牦牛的育种改良工作提供基础信息。

1 材料与方法

1.1 样本采集、DNA提取与全基因组重测序

本研究选择肃南牦牛(Sunan)、巴州牦牛(Bazhou)、九龙牦牛(Jiulong)、斯布牦牛(Sibu)和天祝白牦牛(Tianzhu)5个牦牛品种进行全基因组重测序分析。样本分别采集于甘肃省肃南裕固族自治县、新疆维吾尔自治区巴音郭楞蒙古自治州、四川省九龙县、西藏自治区墨竹工卡县和甘肃省天祝藏族自治县，样本数分别为8、8、12、10和10头。所有个体采集颈静脉血液10 mL于EDTA抗凝管中，置于-20 ℃待测。使用EasyPure血液基因组DNA提取试剂盒(北京全式金生物技术有限公司)提取血样DNA，通过测定A260 nm/A280 nm比值和1.0%琼脂糖凝胶电泳筛选，检测提取DNA的质量和完整性。用Covaris超声波DNA破碎仪将合格的基因组DNA样本随机分成片段。对片段DNA进行末端修复，在3′-端加入碱基A连接测序接头，吸附富集后进行PCR扩增形成测序文库。构建好的文库在BGISEQ-500平台进行测序。测序得到的原始图像数据通过BGISEQ-500 Base Calling软件转换为原始数据，以FASTQ文件格式存储。

1.2 全基因组序列比对及变异检测

使用Trimmomatic[16]对原始读长(raw reads)进行质控，去除适配器序列、污染和低质量reads，将高质量的reads与最新的牦牛参考基因组(GCA_005887515.1)比对得到SAM文件(Sequence Alignment Map)。使用SAMtools(v1.9)将SAM文件排序并转换为二进制格式(BAM)[17]。对上游得到的每个样本的BAM文件采用GATK v4.1.8.0进行SNP鉴定，采用MarkDuplicates、HaplotypeCaller和SelectVariants模块获取SNP的原始文件[18]。通过VariantFiltration模块，设置“QUAL<30.0、QD<2.0、FS>60.0、MQ<40.0、HaplotypeScore>13.0”的过滤条件，对上游得到的初始VCF文件执行硬过滤，最终得到含有高质量SNP的VCF文件。采用samtools flagstat[17]命令计算每个样本的统计信息，包括平均覆盖率、原始读长计数、映射读长和正确配对的读长比率。考虑到基因型数据的低等位基因频率以及性别因素可能会影响后续的分析，使用PLINK v1.9[19]剔除最小等位基因频率(MAFs)<0.05、个体检出率<0.95、SNP检出率<0.99、无染色体分配和性染色体上的SNP，参数设置为“--maf 0.05 --mind 0.05 --geno 0.01 --chr 1-29”。

1.3 群体基因组学分析

基于SNP信息，采用PLINK进行主成分分析(PCA)以确定群体间的遗传结构。PCA的可视化基于R包ggplot2完成。使用ADMIXTURE软件[20]中最大似然估计方法进行祖先成分分析，最佳K值的确定参考Evanno等[21]的方法。采用PLINK软件过滤强连锁SNP位点(--indep-pairwise 50 10 0.3)，随后构建遗传距离矩阵，使用脚本转换为MEGA格式后采用邻接法构建系统进化树。

1.4 亲缘关系系数计算

为了消除亲缘关系较近的样品对下游分析的影响，采用KING软件[22]计算两两个体之间的亲缘关系系数，剔除存在三代以内亲缘关系的个体，即关系判断阈值>0.0442的样本。亲缘关系矩阵可视化使用R包corrplot完成。

1.5 连锁不平衡分析

采用PopLDdecay软件[23]计算每个品种的连锁不平衡(LD)程度，并使用软件自带的Plot_MultiPop.pl脚本绘制LD衰减曲线图。

1.6 ROH分析及近交系数计算

每个个体的ROH使用PLINK软件计算，参数设置为“--homozyg --homozyg-density 50 --homozyg-gap 100 --homozyg-kb 300 --homozyg-snp 50 --homozyg-window-het 3 --homozyg-window-snp 50 --homozyg-window-threshold 0.05”。对得到的ROH数量和长度进行统计和分类，使用GraphPad Prism 8.0软件进行展示。参考Mcquillan等[24]方法计算近交系数(FROH)，计算公式为：FROH=常染色体上ROH片段长度之和/常染色体总长度之和。将各自群体内检测到ROH的SNP进行频次统计，以频率Top 5%的SNP定义为高频SNP，将其构成的ROH区域作为高频ROH区域。使用snpEff和BEDtools对各品种高频ROH区域进行基因注释，并对定位到的基因及重叠基因分别进行GO功能和KEGG信号通路富集分析。

1.7 选择信号分析和候选基因定位

ROH的形成受到自然选择和人工选择的影响，检测ROH片段中位点的选择强度可能会鉴定到与表型相关的关键基因[12]。采用rehh包[25]分别计算各高频ROH区域内SNP的综合单倍型评分(iHS)，将iHS评分Top 1%的位点定义为强受选择位点。针对iHS法得到的受选择位点，使用1.6的方法对其进行注释筛选出候选基因。为了缩短计算时间，使用R包parallel进行并行计算。

2 结 果

2.1 5个牦牛群体变异检测结果

本研究对巴州牦牛、九龙牦牛、斯布牦牛、肃南牦牛以及天祝白牦牛共48个样本进行重测序，共获得816 G高质量数据，平均每个样本获得了17 G的数据，平均测序深度为7.14×。将高质量测序数据比对到参考基因组之后，平均比对率为98.864%。为了进行后续群体遗传分析，对5个群体一同进行SNP变异检测，经过滤后得到15 092 883个SNP位点。此外，5个群体SNP的注释结果显示，各牦牛群体中的SNP位点分布相似，主要都位于内含子区域，其次是基因间区、基因下游区域、基因上游区域、外显子区域，平均占比分别为51.47%、37.96%、4.57%、4.49%、1.09%(表1)。

2.2 牦牛群体结构分析

对5个地方牦牛群体进行主成分分析、亲缘关系分析、群体祖先成分分析和系统进化树分析的结果见图1。由图1可知，5个牦牛群体按照第一主成分(2.5397%)和第二主成分可以区分开(1.8990%)，说明5个群体间没有杂交现象，与品种所处的地理位置相符合。除肃南牦牛群体内较为聚集外，其他4个牦牛品种群体内较分散(图1A)。亲缘关系结果显示，5个牦牛群体内亲缘关系较低，关系判断阈值均<0，不存在三代以内亲缘关系的个体，满足后续分析的要求(<0.0442)(图1B)。祖先成分分析显示，最佳分群数为4群；K=2时，巴州牦牛首先被鉴定出来，肃南牦牛与天祝白牦牛祖先成分相似；K=3时，巴州牦牛、九龙牦牛、天祝白牦牛可以明显区分；K=4时，肃南牦牛拥有其他4个牦牛祖先成分，且天祝白牦牛祖先成分所占比例最高(图1C)。系统进化树显示，巴州牦牛、九龙牦牛、斯布牦牛、肃南牦牛和天祝白牦牛群体各自聚集在一个分支上，肃南牦牛聚集在进化树一端，其他4个牦牛品种集中在进化树另一端(图1D)。

A，5个牦牛群体的主成分分析；B，5个牦牛群体的亲缘关系鉴定；C，祖先成分分析；D，系统进化树A,Principal component analysis of five yak populations;B,Genetic relationship identification of five yak populations;C,Ancestral component analysis;D,Phylogenetic tree图1 群体基因组学分析Fig.1 Population genomics analysis

2.3 LD分析

为了评估不同牦牛群体的LD情况，采用PopLDdecay分别计算各个群体的成对r2值，用于衡量LD水平(图2)。LD分析显示，5个牦牛LD衰减的顺序为：九龙牦牛>斯布牦牛>天祝白牦牛>巴州牦牛>肃南牦牛。

图2 LD分析Fig.2 Linkage disequilibrium analysis

2.4 ROH基本统计及基因组近交系数结果

在5个牦牛品种中共检测到8 426个ROHs，不同群体ROH数量和长度高度相关，介于0.98～0.99之间(图3A)。按照ROH长度进行分类，可将ROH分为3类(0～0.5 Mb；0.5～1 Mb；>1 Mb)，ROH数目在0～0.5 Mb大小的占绝大部分，其次是0.5～1 Mb的ROH。肃南牦牛、天祝白牦牛、巴州牦牛、九龙牦牛和斯布牦牛分别存在26、4、3、1和1个>1 Mb的ROH(图3B)。每条染色体都有ROH分布，高频ROH主要位于22、16、4、2和9号染色体(图3C)。所有品种中ROH平均长度最长为肃南牦牛(133.33 MB)，最短为天祝白牦牛(52.02 Mb)。ROH数量和高频ROH数量排序为：肃南牦牛>九龙牦牛>天祝白牦牛>巴州牦牛>斯布牦牛；平均近交系数排序为：肃南牦牛>巴州/九龙牦牛>天祝白牦牛>斯布牦牛(表2)。

A,群体内每头个体ROH总数目与ROH总长度的散点图；B,短、中和长片段ROH数量直方图；C,不同染色体ROH分布直方图A,Scatter plot of the total number of ROH and the total length of each individual in the population;B,Histogram of the number of short,medium and long ROH;C,Histogram of the number of ROH on different chromosomes图3 ROH分析Fig.3 ROH analysis

表2 ROH和近交系数统计

2.5 高频ROH区域基因注释与富集分析

对高频区域内SNP进行注释，其中巴州牦牛、九龙牦牛、斯布牦牛、肃南牦牛和天祝白牦牛分别定位到58、73、52、465和35个基因，其中包含核内不均一性核糖核蛋白K(HNRNPK)、驱动蛋白27(KIF27)、G激酶锚定蛋白1(GKAP1)等3个重叠基因，均位于7号染色体。对肃南牦牛筛选出的候选基因进行富集分析，以Q-value<0.05作为筛选条件，结果显示肃南牦牛富集到17条GO条目(RNA聚合酶对转录的负调控Ⅱ、ATP结合、胚胎骨骼系统形态发生、脂肪氧合酶通路等)和2条KEGG信号通路(Notch信号通路和Hippo信号通路)(表3)。此外，对3个重叠基因的富集分析结果显示，这些基因参与了16个GO条目和3条KEGG信号通路，包括基因表达调控、细胞连接、胰岛素受体信号通路的正调控、体元投射通路等(图4)。

表3 肃南牦牛候选基因的富集分析结果

续表

图4 重叠基因富集分析结果Fig.4 Enrichment analysis results of overlap gene

2.6 选择信号分析

为了解高频ROH区域受到的选择作用，采用基于单体型(haplotype)的iHS方法对5个牦牛群体共有的高频ROH区域进行选择信号检测，结果显示，候选基因主要定位在16号染色体3个高频ROH区域上(Chr16：64 780 362-65 091 013；Chr16：75 041 712-75 388 110；Chr16：68 520 840-68 824 137)。N-α乙酰转移酶25(NAA25)在每个品种基因组上都筛选出来(图5)。巴州牦牛基因组上筛选到内质网分子伴侣29(ERP29)、跨膜蛋白16(TMEM116)、TRAF型锌指结构域蛋白1(TRAFD1)(图5A)，斯布和肃南牦牛基因组上筛选到ERP29、TMEM116(图5C、5D)。天祝白牦牛基因组上除了NAA25和TRAFD1外，还筛选到HECT结构域E3泛素连接酶(HECTD4)(图5E)。

A～E，分别代表巴州牦牛、九龙牦牛、斯布牦牛、肃南牦牛和天祝白牦牛共享高频ROH区域iHS分析A-E，The iHS analysis in shared ROH islands of Bazhou,Jiulong,Sibu,Sunan and Tianzhu White yak,respectively图5 5个牦牛群体共有高频ROH区域iHS分析Fig.5 iHS analysis in shared ROH islands of five yak populations

3 讨 论

3.1 5个牦牛品种群体结构分析

肃南牦牛是肃南地区在培的牦牛品种，其与国内其他牦牛品种之间的分化程度、分类地位以及进化历史尚不清楚，群体基因组学研究有助于在基因组水平解决这个问题。相比于其他牦牛群体，肃南牦牛群体内部存在最小的变异幅度，可能是因为肃南牦牛由于种质资源开发时间较晚与地理环境原因，因此其变异幅度较小。肃南牦牛与其他牦牛群体在基因组水平上差异较大，其在系统发育树分支长度远长于其他4个牦牛品种。祖先成分分析能够推断每个材料的基因组变异来源于假设亚群的可能性[20]。肃南牦牛包含了其他4种牦牛群体祖先血统，这可能与肃南牦牛培育及改良过程中引入其他牦牛品种优良个体有关，这也与前期文献报道的一致[26]。在各个成分中，其包含天祝白牦牛祖先成分的比例最大，表明肃南牦牛群体与天祝白牦牛群体存在广泛交流，这也与二者所处的地理位置较近相符合[26-27]。根据变异程度、进化距离和祖先成分推断，肃南牦牛可能是在长期外部交流下产生的一个独特类群，与地理距离最近的天祝白牦牛交流广泛。

3.2 ROH分析

家畜经过较强人工选育后，会导致特定区间纯和度的增加，使得群体内ROH富集的频率增加，基因组多样性也随之降低[28]。对5个品种的ROH分类显示，牦牛ROH片段普遍较短(<2 Mb)，统计发现0～0.5 Mb的ROH数目最多，其次是0.5～1 Mb片段，说明小片段ROH对整个基因组的近交贡献最大。相比于其他4个牦牛品种，肃南牦牛的ROH片段总长度和数目最多，其原因可能是肃南牦牛经过高度的近交选育，核苷酸多样性最低。近交会增加纯和等位基因频率，从而增加ROH的长度和数量，因此ROH也常被用来评估家畜群体的近交配程度[10]。本研究计算结果显示肃南牦牛平均近交系数要高于其余4个牦牛群体，其原因可能与其育种模式有关，即使用有限数量的高品质公牛作为父本进行人工授精，造成群体的近交水平的显著增加。由于LD状态受到家畜近交水平的影响，驯化程度越高，选择强度越大的群体，LD衰减越慢[10]。肃南牦牛和巴州牦牛相对其他3种牦牛具有较高的连锁程度，这可能与其受到较强的选择强度有关的，这同样解释了肃南牦牛群体具有较高近交水平的结果。

肃南牦牛高频ROH区域注释到465个候选基因，富集到22条显著通路，其中胚胎骨骼系统形态发生、前后模式规范、甲状腺发育通路与机体发育和大脑形态形成有关，且参与3条通路的主要为同源异形盒(Homeobox，HOX)基因家族成员(HOXA3、HOXA5、HOXD3)。HOX蛋白是一种转录因子，在胚胎发育过程参与组织分化过程[29]。此外研究还发现HOX家族成员与毛囊形成和毛发生长密切相关，参与调控角蛋白和角蛋白关联蛋白的表达，维持毛囊正常形态等[30]。甲状腺发育通路包含FGF8基因，该基因与中脑发育有关，在肢体发育的诱导、启动和维持中起作用[31-32]，且有研究显示FGF8在小鼠表皮中过表达会抑制毛囊的发育[33]。HOXA3、HOXA5、HOXD3和FGF8这4个基因可能参与调控肃南牦牛毛发生长过程，值得进一步研究。脂肪氧合酶通路包含4个基因，其中3个为脂氧合酶(LOX)基因家族成员编码基因(ALOX12B、ALOX15B、ALOXE3)，该家族主要参与氧化不饱和脂肪酸如花生四烯酸等特定的过氧化物[34]。研究显示添加花生四烯酸会影响动物肉类的适口性[35]，该类基因可能与肃南牦牛肉品质有关。参与显著通路的基因主要以基因家族的形式出现，这一特点可能与肃南牦牛品种培育过程中特定功能的强化有关。另外还富集到一条Hippo通路，该通路最初在果蝇中被发现，被认为是一种保守的信号通路，在控制器官大小、限制细胞增殖和诱导细胞凋亡方面发挥着重要作用[36]，可能与肃南牦牛机体发育有关。

对各高频ROH区域筛选的候选基因取交集，发现HNRNPK、KIF27、GKAP1 3个基因在5个牦牛品种中是共享的，均位于7号染色体1个高频ROH区域中(Chr7：12 661 870-13 045 935)。对共享基因的富集分析结果显示，3个基因参与了胰岛素受体信号通路的正调控通路。Bao等[37]研究表明胰岛素相关通路是牦牛能够适应高原环境的主要调节机制之一，胰岛素的动态调节有助于高原动物更好地利用能量，抵御恶劣的自然环境。参与脑室系统开发通路的KIF27是哺乳动物中COS2的同源基因之一，其功能是作为Hedgehog信号通路的调节因子，参与胚胎期发育、组织修复、器官形成等过程[38]，也有小鼠上的研究证明该通路在低氧条件下肺血管重建过程中发挥作用[39]。HNRNPK是一种对细胞周期有广泛影响的RNA结合蛋白。小鼠研究中发现，血红素氧合酶的破坏可能增强hnRNPK介导的蛋白翻译抑制，进而削弱β-catenin/hnRNPK调控的基因表达，这是协同肺修复和再生所必需的[40]。相比平原地区的牛，牦牛生活在高原低氧极端环境中，其肺部结构发生了一系列适应性变化，如牦牛的囊状期较短，肺泡期较长，肺脏发育时间较早等[41]，推测hnRNPK基因与牦牛低氧气浓度影响下的肺功能修复过程有关。本研究鉴定的3个共享基因可能通过胰岛素代谢、肺血管重建和肺功能修复途径参与高原适应性的调节，这一结果有助于解析牦牛高原适应性进化的分子机制和遗传原理。

3.3 受选择基因分析

iHS统计量可以灵敏地检测正选择信号及近期选择，可用来揭示选择信号背后的功能基因[42]。本研究用iHS方法对共有高频ROH区域内的位点进行选择信号分析，发现5个牦牛16号染色体上3个高频ROH区域(Chr16：64 780 362-65 091 013；Chr16：75 041 712-75 388 110；Chr16：68 520 840-68 824 137)内部分SNP受到强烈选择，这3个相邻区段可能是5个牦牛品种在染色体序列发生不断重组的进化过程中保留下来的受选择热点区域。注释结果显示，该区段得到5个不同的基因(NAA25、ERP29、TMEM116、TRAFD1、HECTD4)，这些基因部分或全部牦牛品种中检测到，可能与滑动窗口和步长的大小设置有关。NAA25在每个品种中都鉴定出，该基因编码N-端乙酰转移酶B复合体(NatB)的辅助亚基，能够乙酰化蛋氨酸残基，可能在正常细胞周期进程中起作用[43]。除九龙和天祝白牦牛外，ERP29在其他牦牛都鉴定出，该基因编码蛋白是内质网腔内的一种驻留蛋白，参与内质网分泌蛋白的加工过程[44]。TMEM116在巴州牦牛、斯布牦牛和肃南牦牛上鉴定出，该基因编码蛋白属于跨膜蛋白(TMEM)家族，与细胞间、细胞内的信号传导、免疫相关疾病相关[45]。TRAFD1在巴州牦牛和天祝白牦牛中鉴定出来，该基因负向调控人类和哺乳动物的Toll样受体(TLR)和RIG-Ⅰ样受体(RLR)信号转导[46]。TLR介导的信号通路与免疫反应中病原体的消除、获得性免疫体系的建立有关[47]。RLR是识别胞浆中病毒RNA的主要受体，在抵抗病毒感染的过程起关键作用[48]。TMEM116、TRAFD1受到强烈选择可能与牦牛品种良好的抗病性有关。在天祝白牦牛中鉴定到HECTD4，该基因在人上可能编码E3泛素蛋白连接酶[49]。其在牦牛中的作用尚不清楚，需要进一步的验证。本研究受选择分析中显著位点集中在16号染色体3个相邻区段上，推测是由于牦牛之间分化程度较低，牦牛共有高频ROH区段内位点受到较为一致的选择压力。筛选到受选择基因主要与抗病性、内质网分泌蛋白加工、以及细胞周期调节相关，这一结果有助于进一步了解不同牦牛类群整体的进化方向。

4 结 论

研究结果表明，肃南牦牛是在长期外部交流下产生的一个独特类群，受天祝白牦牛血统影响较大。相比于其他4个成熟牦牛品种，具有较高的近交水平，可能与品种培育过程受到的强化选择有关。在肃南牦牛高频ROH区域中鉴定到与胚胎发育与组织分化过程(HOXA3、HOXA5、HOXD3)，肉品质(ALOX12B、ALOX15B、ALOXE3)和中脑发育(FGF8)相关的基因。此外，与其他4个牦牛品种共享3个与高原适应性相关的基因hnRNPK、KIF27、GKAP1。在共有高频ROH区域鉴定到的受选择基因主要与抗病性、内质网分泌蛋白加工及细胞周期调节相关，包括NAA25、ERP29、TMEM116、TRAFD1、HECTD4基因。本研究结果有助于了解肃南牦牛的品种特性，为正在进行的基因组选择育种和下一步育种计划的制定提供参考。