郭家中，孙学良，向秋楠，肖凤敏，张红平

（四川农业大学动物科技学院，成都 611130）

我国地方畜禽品种资源丰富，但群体规模通常较小，因此需要控制近交从而避免近交衰退[1]。传统上，个体的近交系数（inbreeding coefficients，F）是利用系谱信息进行估计[2]；但这类方法主要有两方面不足：一是需假定祖先个体是非近交个体，这往往背离事实从而低估个体的近交水平；另一方面，系谱近交系数反映特定亲缘关系下等位基因同源相同的期望值，无法准确反映个体真实的近交水平[3]。此外，很多群体系谱记录不完善甚至缺乏系谱。而利用基于芯片或重测序的SNP数据估计的近交系数在全基因组水平上反映了个体真实的等位基因同源性的概率或相关性，克服了基于系谱方法的局限性，这类近交系数统称为基因组近交系数。基因组近交系数估计方法较多，但从原理上主要包括基于纯合子过度的估计值（excess of homozygosity inbreeding coefficient，FHOM）[4]、长 纯 合 片 段（runs of homozygosity，ROH）占基因组比例（FROH）[5]和基于频率校正的多位点平均纯合度[6]。在实践中，第三类方法理论上是基于基因组亲缘关系矩阵（genomic relationship matrix，GRM）而获得近交系数（FGRM）。目前，关于牛[7]、猪[8]和绵羊[9]等物种的基因组近交系数及不同方法的报道较多，但是在山羊[10]上研究较少。

全基因组选择信号研究已成为鉴定影响畜禽重要性状遗传变异位点的一种重要方法，并在山羊上已有诸多报道[10-12]。川中黑山羊（Chuanzhong Black goat，CZ）是我国优良地方山羊品种，又进一步分为金堂型和乐至型，该品种具有体型大、生长速度快和产羔数多等优点[13]。该品种在西南地区很受欢迎，被多个地区引入作为父本，改良地方黑山羊品种的生长性状；但关于川中黑山羊的遗传研究相对较少。我们前期利用9只公羊的重测序数据，对该品种进行了群体遗传分析[14]；但样本量太少。因此有必要扩大样本量对川中黑山羊近交水平、有效群体大小和选择信号等群体遗传特征进行更深入的探究。

本研究利用前期收集的41只川中黑山羊（金堂型）的重测序数据获得了高密度SNP基因分型数据，主要开展基因组近交系数、有效群体大小估计和选择信号分析，旨在为川中黑山羊遗传资源的发掘利用提供理论参考。

1 材料和方法

1.1 试验群体和高通量重测序数据

川中黑山羊（金堂型）样本均来自于成都市某公司育种场，共包括41只种羊，其中9只公羊、32只母羊。我们前期研究已经发表并公布了上述41个个体的短序列基因组重测序数据[14-15]。为进行选择信号分析，还使用了前期获得的30只建昌黑山羊的重测序数据[10]并下载了21只野山羊的重测序数据（NCBI accession number:PRJEB3136）[16]。

1.2 高通量测序数据处理和短变异检测

参考我们之前使用的分析流程[10]，进行重测序数据比对和短变异检测。使用BWA软件（v 0.7.17）[17]将短序列映射到山羊参考基因组（ARS1[18]，GCA_001704415.1），获得比对结果。使用GATK软件（v 4.0.5.2）[19]初步检测SNPs和Indels并进行硬过滤，再使用VCFtools[20]进行群体遗传学质控，仅保留最小等位基因频率大于0.05、缺失率小于10%的位点；在进行各项遗传分析前，进一步过滤掉哈迪温伯格平衡P值小于10-10的双等位基因SNPs。使用SnpEff软件（v 4.3）[21]对短变异进行功能注释。

1.3 有效群体大小和全基因组ROH分析

使用 PopLDdecay（v 3.4.1）软件[22]分析群体中SNPs之间的连锁不平衡（linkage disequilibrium，LD）。使用基于连锁不平衡理论开发的SNeP（v 1.1）[23]和默认参数估算山羊历史有效群体大小（effective population size，Ne）。使用PLINK（v 1.9）软件[24]中的“--homozyg”命令检测川中黑山羊基因组中的长纯合片段（runs of homozygosity，ROH）。主要参数值设置如下：每个ROH片段最少包含10个SNPs，最小长度为100 kb，最低SNP密度为10 kb/SNP，每个ROH内SNP的最大间隔为100 kb；滑动窗口大小为50个SNPs，每个滑动窗口中最多允许1个杂合位点和5个缺失位点，滑动窗口阈值为0.05。

1.4 基因组近交系数的估计

基于ROH的近交系数（FROH）定义为每个个体基因组的ROH总长度占常染色体基因组长度的比例（参考基因组ARS1的常染色体总长度为2 466 191 353 bp）。根据ROH长度与世代数之间的关系，依据不同长度将ROH分成4类：0.1～0.2、0.2～0.5、0.5～1.0和＞1 Mb，计算对应世代的FROH。使用PLINK 软件[24]中的“--het”命令计算FHOM。Van-Raden在2008年提出了3种GRM计算方法[25]，并被广泛使用。GMAT软件[26]实现了VanRaden提出的第一种方法，使用该软件中的“--grm agrm--fmt 0”命令获得GRM中对角线元素，估计FVR1。使用GCTA[27]（v 1.92）“--ibc”命令计算FVR2和FUNI。

1.5 选择信号鉴定与功能富集分析

本研究使用ROH岛、iHH12和XP-EHH 3种方法鉴定川中黑山羊基因组选择信号；其中，前两种方法分别基于群体内的长纯合片段和单倍型分布特征，第3种方法则是比较群体间扩展单倍型纯合度的分化程度。基于PLINK获得的ROH结果，使用R包detectRUNS[28]进行ROH岛检测，鉴定标准包含3个指标：①个体间共享的ROH出现的频率大于50%（即在41个个体中至少21个个体携带相同的ROH片段）；②ROH岛内至少包含10个SNPs；③ROH岛长度最短为1 000 bp。使用selscan软件[29]计算iHH12，主要流程如下：使用Beagle软件[30]对SNPs进行填充和定相，然后计算原始的iHH12值，再使用norm模块对原始值进行归一化处理。最后利用自编R脚本，以10 kb窗口和10 kb步长沿染色体滑动计算每个窗口内的iHH12平均值。将所有窗口平均值按从高到低进行排序，前0.5%的窗口作为选择信号的候选窗口。使用selscan软件计算XP-EHH，主要步骤类似于iHH12计算。由于在分析中，将川中黑山羊作为试验群体、野山羊作为参考群体，XP-EHH值为正值则代表某区域在川中黑山羊中受到正选择。因此，在XP-EHH值从高到低的分布中，选择前0.5%的窗口作为选择信号的候选窗口。使用BEDTools[31]对3种方法鉴定的候选区域分别进行基因注释，3种方法共享基因被定义为正选择基因。

使用R包clusterProfiler（v 4.4.1）[32]对正选择基因进行GO功能富集分析，显著富集的GO条目筛选标准为P＜0.05。在分析中，山羊全基因组范围的基因功能注释信息来自AnnotationHub（编号为“AH101444”）。

2 结果与分析

2.1 全基因组SNPs和Indels检测

在41只川中黑山羊常染色体基因组中，共检测到14 043 333个双等位基因SNPs、70 478个复等位基因SNPs和1 197 402个Indels。变异注释结果表明，位于基因间区和内含子区的SNPs的比例最高，分别为45.38%和44.15%，而外显子区的SNPs仅占0.92%。类似地，内含子区和基因间区的Indels的比例也是最高，分别为45.89%和42.95%。

2.2 连锁不平衡和有效群体大小分析

由图1a可知，当SNPs之间的物理距离为10 bp时，r2平均值（0.51）最高；随后LD迅速衰减，当SNPs之间的距离增加到1 000 bp时，r2等于0.2。图1b可知，川中黑山羊有效群体大小持续缩减，在999世代前Ne值为5 696只，而13世代前Ne值为190只。

图1 川中黑山羊连锁不平衡和有效群体大小分析Figure 1 Neand linkage disequilibrium decay in Chuanzhong Black goats

2.3 基于5种方法的川中黑山羊基因组近交系数估计值

在41只川中黑山羊基因组上共检测到47 831个ROH，ROH在1～29号染色体上均有分布；其中，1号染色体的ROH数量最多（3 287），27号染色体上ROH数量最少（742）。最长的ROH位于18号染色体（29 071 104～35 631 696 bp）长度为6.56 Mb，该区域包含CDH8、TK2和CMTM3等26个基因。

如图2a所示，41只川中黑山羊FROH值范围为0.06～0.21，平均值为0.12。在群体水平上，当前世代川中黑山羊所累积的近交水平主要来自于250～500世代（FROH0.1-0.2Mb）和100～250世代（FROH0.2～0.5Mb）。相关分析表明，当前世代个体间FROH的变异与50～100世代（FROH0.5-1.0Mb，r=0.94，P＜2.2×10-16）相关性最强。如图2b所示，在41只川中黑山羊中，FVR1近交系数最小值为-0.03，最大值为0.30，平均值为0.21。FVR2近交系数最小值为-0.02，最大值为0.27，平均值为0.18。类似地，FUNI近交系数最小值为-0.03，最大值为0.26，平均值为0.18。而FHOM近交系数最小值为0.07，最大值为0.27，平均值为0.19。如表1所示，除FROH与FVR1和FVR2之间无显著相关，其余不同方法获得的川中黑山羊近交系数之间均正呈现显著性正相关（P＜0.05）。其中，FUNI和FVR1两种近交估计值的线性相关最高（r=0.983，P=2.2×10-16），而FUNI和FHOM两种估计值的相关性也较高（r=0.893，P=4.0×10-15）。

图2 基于5种方法的川中黑山羊基因组近交系数估计值Figure2 Summary of estimated genomic inbreeding coefficients in Chuanzhong Black goats using five methods

表1 基于5种方法的川中黑山羊基因组近交系数的皮尔逊相关Table 1 The Pearson’s correlations between the genomic inbreeding coefficients in Chuanzhong Black goats based on five methods

2.4 川中黑山羊选择信号的检测

如图3a所示，基于iHH12、ROH岛和XP-EHH统计量在川中黑山羊中分别鉴定到1 218（iHH12＞7.71）、88（ROH共享率大于50%）和1 227（XP-EHH＞3.10）个离群值窗口并分别注释到361、164和324个基因（以Ensemble ID为标准）。其中，NCAPG（chr6:37 858 170～37 903 004 bp）、LCORL（chr6:37 905 295～38 068 616 bp）、ESR1（chr9:76 096 964～76 376 135 bp）、KIT（chr6:70 711 312～70 794 841 bp）等67个基因被3种方法均检测到，被定义为川中黑山羊的正选择基因。

图3 川中黑山羊选择信号和正选择基因Figure 3 Summary of selection signals and positively selected genes in Chuanzhong Black goats

上述67个正选择基因显著富集在269项GO条目中（P＜0.05），其中在生物学过程、细胞组分和分子功能条目上各富集到175、50和44项条目。表2展示了显著富集的前10个GO条目，其中最显著富集的条目是上皮运输（GO:0070633，transepithelial transport）和平滑肌细胞分化调节（GO:0051150，regulation of smooth muscle cell differentiation）生物学过程；这些过程包含ABCG2、AHCYL1、KIT和MED28共4个基因。另外，NCAPG基因显著富集在有丝分裂染色体凝聚信号通路（GO:0007076，mi-totic chromosome condensation，P=0.024）。

表2 川中黑山羊正选择基因显著富集的前10个GO条目Table 2 Top ten enriched GO terms for the positively selected genes in Chuanzhong Black goats

2.5 NCAPG-LCORL选择信号区域重要变异的鉴定

由图3b可知，山羊6号染色体NCAPG-LCORL座位所在区域在全基因组范围内显示出最高的iHH12值（iHH12值=48.52），表明该区域在川中黑山羊中是一个强烈的正选择区域。如图4a所示，该区域在川中黑山羊和野山羊之间的平均加权Fst值为0.39，远高于全基因组Fst的平均值（0.16）。另外，在川中黑山羊中该区域Tajima’sD（平均值=-0.53）也远低于全基因组平均值（1.31）。在NCAPGLCORL座位内共检测到435个SNPs和63个indels；其中包括位于NCAPG基因第6外显子的1个错义突变（c.858A＞G，p.Ile286Met）和位于LCORL基因内的3个错义突变（c.4397C＞T，p.Ala1466Val；c.1433A＞G，p.Asn478Ser；c.1298A＞G，p.Tyr433Cys）。费希尔精确检验表明，c.858A＞G（P=1.36×10-5）和c.1298A＞G（P=6.08×10-5）位点在川中黑山羊和野山羊群体中的基因型频率分布存在显著性差异。但只有c.858A＞G基因型分布在川中黑山羊和建昌黑山羊群体中（参考型等位基因频率=93.33%）存在显著性差异（P=4.68×10-12）。而63个indels中包括两个移码突变位点（c.1615_1619delTTAAA，p.Leu539fs；c.828dupA，p.Ser277fs），它们均位于LCORL基因的第7外显子内。费希尔精确检验显示，移码突变位点c.1615_1619delTTAAA的基因型频率分布在川中黑山羊（参考型等位基因频率=13.16%）和野山羊群体（参考型等位基因频率=100%）存在显著性差异（P=5.89×10-11）。另外，在川中黑山羊和建昌黑山羊群体中（参考型等位基因频率=93.33%）该位点的基因型分布也存在显著性差异（P=4.3×10-12）。

图4 川中黑山羊基因组的NCAPG-LCORL选择信号深入分析Figure 4 In-depth analysis of the selection signal NCAPG-LCORL locus in Chuanzhong Black goats

3 讨论

基因组水平的畜禽遗传研究依赖于全基因组范围的分子标记，本研究利用重测序技术获得了川中黑山羊群体的高密度SNP图谱。与低密度的SNP芯片（例如，山羊52 k SNP芯片[33]）技术相比，利用短序列高通量测序在山羊[10]、绵羊[34]、牛[35]、猪[36]和鸡[37]等物种上均可获得千万级的SNP基因型数据，从而显著提高了全基因组关联研究等多种遗传分析的成功率和精确性。与其他物种的注释结果相似，在川中黑山羊群体中绝大部分SNPs位于基因间区和非编码区，而外显子区的变异占比极低。

有效群体大小本质上反映了群体的遗传多样性丰富程度。本研究发现，在最近1 000世代内川中黑山羊Ne随着时间持续减少，这和我国其他山羊品种的变化特征一致[10,38-39]。导致上述变化的主要原因是，我国大部分地方品种生长性能不突出、养殖效益低，导致养殖户的不断退出，群体规模逐渐减少。另一方面，由于研究样本往往仅来自保种群或育种群，这些群体初始世代包含的种羊数量偏少。另外，川中黑山羊最近世代的Ne高于我国其他地方品种。在未来的川中黑山羊遗传改良工作中，应保持甚至增加家系数量，防止遗传多样性的丧失。

尽管相比于系谱近交系数，基因组近交系数更能反映个体真实的近交水平；究竟哪种方法更适应于有效群体含量较小的畜禽品种尚无统一答案[6-7,40]。因此，本研究采用了5种常用方法对川中黑山羊近交系数进行估计和比较。总体上，基于ROH（＞100 kb）估计的川中黑山羊的FROH值与瑞士山羊的近交程度接近[41]，而低于建昌黑山羊的近交水平[10]。F.Bertolini等[42]则利用SNP芯片获得的长ROH（＞1 Mb）估计了全世界117个山羊群体的近交程度，并依据FROH将近交水平划分为低（FROH＜0.1）、中（0.1＜FROH＜0.2）和高（FROH＞0.2）3类。据此，总体上川中黑山羊处于中等程度的近交水平；但如果排除掉川中黑山羊基因组中短ROH（＜1 Mb），则川中黑山羊的近交水平较低。与建昌黑山羊的结果类似，川中黑山羊个体近交系数之间的差异主要归因于最近50～100世代的近交。另外，川中黑山羊FROH和FHOM之间存在较高的相关性，这与其他家畜上的结果[7,43]一致。而利用GRM矩阵所估计川中黑山羊的FUNI、FVR1和FVR2值之间相关较高，主要是因为这些方法的基本原理相同，都是均基于频率校正的多位点平均纯合度。

川中黑山羊具有体型大、生长速度快和产羔数多等优点，但关于上述性状的遗传基础鲜有报道。本研究表明，6号染色体的NCAPG-LCORL座位是川中黑山羊中最强烈的选择信号。NCAPG基因编码非SMC凝聚素Ⅰ复合亚基G，该基因除了在有丝分裂和减数分裂中调节染色体的稳定和压缩外，还在肿瘤发生中扮演重要作用[44]。LCORL基因编码配体依赖性核受体共抑制因子样蛋白，该基因最初被认为是一个在精子细胞中表达的转录因子[45]。NCAPG和LCORL在多个动物基因组中均是彼此相邻，故被合称为NCAPG-LCORL座位。基于多个群体研究表明，NCAPG-LCORL座位与牛生长（例如，采食量）、体型（例如，体重）和繁殖性状均呈现显著性关联，具有一因多效性[46]。综合马[47]、犬[48]和猪[49]等其他物种的研究，NCAPG-LCORL被公认为是影响动物体型大小（体高、体长等）的一个重要遗传座位。基于3个群体，我们发现山羊NCAPG和LCORL基因均存在大效应突变位点，但其他物种上LCORL基因内的变异位点与动物体型性状的关联效应更大[46]。具体哪个基因更可能是影响山羊体型的候选基因，以及上述突变位点的具体效应和作用机制值得扩大样本量进一步探究。