几个中国地方猪群体遗传结构和产仔数性状选择信号分析

路玉洁，莫家远，綦文晶，朱思燃，杨丽丽，刘巧玲，卜亚歌，兰干球，梁 晶

(广西大学动物科学技术学院，南宁 530004)

猪与人类生活息息相关，其毛、皮和脂肪可作为工业原料，粪便可作为优良的肥料使用，最重要的是猪能为人类提供优质的蛋白质来源。虽然2020年我国猪肉产量为4 113万吨，生猪出栏52 704万头，占全年肉类产量的53.84%，位居世界第一，但目前我国猪生产水平还较低，母猪年断奶合格仔猪数(PSY)能达到25头的规模猪场并不多，与国外先进水平差距达5～6头[1]。产仔数是衡量母猪PSY高低的关键指标[2]，按照2020年生猪出栏量计算，若我国母猪的平均PSY能提高1头，预计可以减少82.52万头母猪的饲养量，这将节省大量的人力、物力和财力，大大有利于我国生猪产业的发展，因此提高产仔性能在养猪生产中尤为重要。

产仔数性状是一个复杂的数量性状，受微效多基因控制，其遗传力约为0.1，属于低遗传力性状[3-4]，实践证明基因组选择是提高此类性状的优良育种手段[5]，而性状相关候选基因的确定与基因组选择的准确性有很大关系，利用基因组序列信息对家畜品种进行选择信号分析常作为一种鉴定性状相关候选基因的方法。Bovo等[6]对欧洲猪种进行选择信号分析，不仅发现OCA2基因与毛色有关，还发现NCAPG-LCORL和CASP10基因与猪的体型有关；Reimer等[7]对汉普夏、杜洛克等猪种进行选择信号分析发现，雄激素AR基因可能与猪的小体型有关；Trenhaile等[8]对西方杂交群体进行选择信号分析发现，P2X3R基因可能与产仔数性状相关；Xu等[9]对梅山猪、金华猪等品种进行选择信号分析发现，EPAS1、CXCL2和TLR2等14个基因可能与猪传染性肺炎易感性相关，李开军等[10]对淮猪进行选择信号分析发现了7个与粗纤维消化率相关的基因。本研究利用固定指数(Fst)方法进行选择信号检测，以鉴定可能影响中国地方猪产仔数性状的候选基因，拟为中国地方猪种遗传资源的开发利用和未来精准育种的实施提供理论依据。

1 材料与方法

1.1 数据来源与分组

本研究下载了102头中国地方猪的Illumina PorcineSNP60芯片数据(http://dx.doi.org/10.5061/dryad.30tk6)[11]，包括迪庆藏猪(dq)19头、明光小耳猪(mg)16头、五指山猪(wz)16头、姜曲海猪(jq)11头、蓝塘猪(lt)20头和梅山猪(ms)20头。 根据各品种平均窝产仔数分为高、低产仔数组，其中高产仔数组包含姜曲海猪、蓝塘猪和梅山猪3个品种共51头猪，低产仔数组包含迪庆藏猪、明光小耳猪和五指山猪3个品种共51头猪(表1)。

表1 6个中国地方猪种平均窝产仔数统计

1.2 数据过滤与填充

使用plink 1.90软件对下载的芯片数据进行质控，去除最小等位基因频率(MAF)小于0.05的位点，分型成功率小于0.9的个体和位点；将质控后的数据按照品种划分为6个群体，使用Beagle 5.0软件进行基因型填充，将填充完毕的6个品种数据整合为一个vcf文件用于后续分析。

1.3 亲缘关系、主成分、群体结构和进化树分析

使用Tassel_v5软件计算6个中国猪种个体之间的亲缘关系(Kinship，K)矩阵和主成分(Principal component analysis，PCA)矩阵，并使用R 4.0.2绘图；使用plink 1.90软件对6个品种位点进行过滤，窗口大小设为100 kb，移动步数为10，R2阈值设为0.05，数据过滤后使用Structure软件进行群体结构分析，并用R 4.0.2绘制群体结构图；使用Tassel_v5软件将vcf格式转化为phy格式后，使用phylip 3.698软件以邻接法(NJ)构建进化树。

1.4 选择信号分析

使用vcftools 4.2软件对6个品种分别与高产仔数组或低产仔数组进行Fst分析，包括：迪庆藏猪与高产仔数组、明光小耳猪与高产仔数组、五指山猪与高产仔数组、姜曲海猪与低产仔数组、蓝塘猪与低产仔数组、梅山猪与低产仔数组。为进一步确定产仔数选择信号，本研究进行以下4个步骤：1) 选取6个分析组中Fst值前1%的位点作为一组中受到选择的选择位点；2) 将迪庆藏猪与高产仔数组、明光小耳猪与高产仔数组、五指山猪与高产仔数组中相同的受选择位点视为低产仔数品种与高产仔数组的选择信号；3) 将姜曲海猪与低产仔数组、蓝塘猪与低产仔数组、梅山猪与低产仔数组中相同的受选择位点视为高产仔数品种与低产仔数组的选择信号；4) 将第二和第三步筛选到的共同位点作为受到产仔数选择的选择信号。

1.5 候选基因确定

使用vcftools 4.2软件计算迪庆藏猪、明光小耳猪、五指山猪、姜曲海猪、蓝塘猪和梅山猪各品种的连锁不平衡衰减率与6个品种总的连锁不平衡衰减率，并用R 4.0.2绘制连锁不平衡(linkage disequilibrium，LD)衰减图。将筛选到的受选择SNP位点的上、下游12.68 kb作为候选区域，使用Biomart(https://may2017.archive.ensembl.org/biomart/martview/f80b3777e8381c63c577 d5980549139c)进行候选基因定位。对定位到的候选基因在David网站(https://david.ncifcrf.gov/conversion.jsp?VFROM=NA)进行GO和KEGG分析，使用R 4.0.2绘图。通过文献查阅基因的生物学功能，综合分析确定与中国地方猪产仔数性状相关的候选基因。

2 结 果

2.1 平均窝产仔数统计

对6个中国地方猪种群的平均窝产仔数进行统计发现，迪庆藏猪、明光小耳猪和五指山猪构建的低产仔数组的平均窝产仔数低于姜曲海猪、蓝塘猪和梅山猪构成的高产仔数组(表2)。

2.2 SNP质控

本研究下载的102头中国地方猪数据包含61 773个 SNPs位点，去除非常染色体和性染色体上的位点、MAF小于0.05和分型成功率小于0.9的位点，最终共有33 285个位点用于后续分析。

2.3 亲缘关系及主成分分析

使用Tassel_v5软件分析6个中国地方猪种的亲缘关系和主成分关系，得到亲缘关系矩阵与主成分矩阵，并用R 4.0.2绘制亲缘关系图和三维PCA图。结果显示,6个品种间的亲缘关系较远，品种内亲缘关系较近，在梅山猪的群体内存在3个亲缘关系较高的小群体(图1)；PC1、PC2和PC3分别解释了变异的12.99%、9.74%和5.42%，迪庆藏猪、姜曲海猪、梅山猪和五指山猪4个猪种内个体分别聚集在一起，明光小耳猪群体内个体分布较分散，且与迪庆藏猪距离近，蓝塘猪群体中有部分个体单独聚集(图2A)。

表2 高、低产仔数组平均窝产仔数统计

2.4 LD衰减计算

使用vcftools 4.2软件计算我国6个猪品种的LD衰减率，对102头中国地方猪整体LD进行计算发现，在距离为12.68 kb时整体的R2只剩下0.17，而在126.6 kb以后的距离处观察到R2小于0.10(表3)。其中迪庆藏猪和五指山猪的LD衰减很接近，在6个中国地方品种中衰减最快，在250 kb之后与整体的衰减速度基本一致；6个品种的LD衰减速率依次为梅山猪<蓝塘猪<姜曲海猪<明光小耳猪<五指山猪或迪庆藏猪(图2B)。

表3 102头中国地方猪LD衰减率

2.5 群体遗传结构分析

Plink 1.90软件质控之后剩余932个SNPs位点，使用Structure软件进行群体遗传结构分析。结果显示，当k=2时，迪庆藏猪、明光小耳猪和五指山猪低产仔数组与姜曲海猪、蓝塘猪和梅山猪高产仔数组很明显的分离出来；当k=3时，迪庆藏猪和梅山猪分离出来；当k=4时，明光小耳猪和五指山猪，姜曲海猪和蓝塘猪仍然混在一起没有分离开；当k=5时，五指山猪分离出来；当k=6时，蓝塘猪分离出来(图3A)。使用phylip 3.698绘制群体进化树，结果显示，姜曲海猪、蓝塘猪、五指山猪、迪庆藏猪和明光小耳猪群体各自在一个分支上，梅山猪群体在1个分支上存在3个小分支，明光小耳猪群体中有1个个体并未和群体中的其他个体划分在同一分支上(图3B)。

2.6 选择信号分析

分别对迪庆藏猪、明光小耳猪和五指山猪与高产仔数组之间进行Fst分析，保留其各自前1%的SNP位点，分别获得331、332和331个位点，其中在3个对比组中都出现的SNP位点共有176个；分别对姜曲海猪、蓝塘猪和梅山猪与低产仔数组之间进行Fst分析，保留其各自1%的SNP位点，分别获得331、331和332个位点，其中在3个对比组中都出现的SNP位点共有176个；高、低组共同筛选到的受选择位点为176个(表4)。

2.7 候选基因定位

将筛选到的176个SNPs位点上、下游12.68 kb距离作为选择区域，用Biomart进行基因定位，共获得66个候选基因，并使用David网站进行GO和KEGG富集分析。GO富集分析结果显示，66个基因在启动子特异性染色质结合(GO:1990841)、肌动蛋白细胞骨架组织的调节(GO:0032956)、RNA聚合酶Ⅱ调控区序列特异性DNA结合(GO:0000977)、Ras蛋白信号转导(GO:0007265)和骨化结合(GO:0001503)通路中富集；KEGG富集分析显示，66个基因只在谷氨酸能突触传递(ssc04724)通路富集。通过查阅文献发现，有CHD7、FBXO43、RUNX1、STRBP和TEDDM1共5个基因与繁殖性状相关，其中RUNX1基因参与RNA聚合酶Ⅱ调控区序列特异性DNA结合和骨化结合通路，CHD7参与启动子特异性染色质结合通路。

dq表示迪庆藏猪，mg表示明光小耳猪，wz表示五指山猪，jq表示姜曲海猪，lt表示蓝塘猪，ms表示梅山猪；下同dq means Diqing Tibetan pig, mg means Mingguang Xiaoer pig, wz menas Wuzhishan pig, jq means Jiangquhai pig, lt means Lantang pig, ms means Meishan pig. The same as below图1 6个猪种亲缘关系图Fig.1 Genetic relationship of 6 pig breeds

A. 三维PCA图；B. LD衰减图；图中total表示所有个体A. Three-dimensional PCA map; B. LD attenuation diagram; The total means all of the pigs in this research图2 6个中国地方猪品种的PCA图和LD衰减图Fig.2 The PCA map and LD attenuation map in 6 Chinese indigenous pig breeds

A.群体结构图；B.进化树A. Population structure diagram; B. Phylogenetic tree图3 群体遗传结构和进化树Fig.3 Population structure and phylogenetic tree

表4 Fst分析受选择位点

3 讨 论

亲缘关系分析可以使用亲缘系数直观地反映群体内和群体间的相关性高低，PCA分析可对群体进行聚类分析，群体遗传结构分析可以反映群体间是否存在基因交流，进化树可以直观地反映群体间的进化关系及亲缘关系。对群体进行亲缘关系、PCA、群体遗传结构分析和进化树构建，研究各群体之间的分化情况，可以更进一步了解不同品种的遗传距离和进化关系等，为深入了解中国地方猪种的进化提供理论参考[18-19]。本研究中的亲缘关系分析发现，梅山猪群体内有3个亲缘关系较近的亚群，这可能是因为采集的样品主要由3个家系构成导致的；PCA分析发现，五指山猪、姜曲海猪、蓝塘猪和梅山猪4个群体分离度高，在3个PC的情况下完整分离出来，迪庆藏猪和明光小耳猪两个群体有部分个体聚集在一起，可能是因为这2个品种都来自云南地区，基因间存在一定的交流。Structure分析发现，当k=2时，迪庆藏猪、明光小耳猪、五指山猪低产仔数组和姜曲海猪、蓝塘猪、梅山猪高产仔数组出现明显分离，但品种之间无法分离，这可能与两个分组的产仔数差异和体型差异等有关；当k=3时，迪庆藏猪和梅山猪首先被分离出来，可能因为迪庆藏猪与梅山猪是本研究中产仔数差异最大的两个品种，且迪庆藏猪属于生长在高海拔的小型猪种，而梅山猪属于生长在平原的较大型地方猪品种。进化树结果显示，姜曲海猪、蓝塘猪、五指山猪、迪庆藏猪和明光小耳猪群体各自在同一个分支上，梅山猪群体在1个分支上有3个小分支，与亲缘关系结果一致；明光小耳猪群体中有1个个体出现离群，并和群体中的其他个体在同一个分支上，这可能与样品来源不纯等因素有关。从进化树上看，群体间的距离基本上与各群体所在的地理位置距离一致，迪庆藏猪与明光小耳猪均位于云南，在进化树上是距离最近的两个品种；五指山猪与蓝塘猪分别位于海南与广东，其距离也较近；姜曲海猪和梅山猪均位于江苏地区,聚集在一个大分支上，且与五指山猪进化距离较远，与Ai等[20]利用全基因组测序方法对中国地方猪进化关系研究的结果一致。

对102头中国地方猪进行统计分析发现，高产仔数组平均产仔数比低产仔数组高7.31头，为低产仔数组的2.28倍，说明不同中国本地猪种在各地驯化后由于长时间的自然选择和人工选择在很多性状上出现了极大的变异，是研究基因组选择信号非常好的材料。Ai等[21]利用猪60K SNP芯片研究中西方猪的连锁不平衡差异发现，西方猪的平均R2为0.3时距离为125 kb，而中国地方猪仅为10.5 kb。本研究计算6个中国地方品种的LD衰减速度也有类似发现，它们的衰减速度都非常快，在12.68 kb就已经衰减到了0.17，与广西地方猪种[19]和太湖地区6个中国地方猪品种[22]的群体遗传结构研究结果一致。而西方商业猪种在500 kb时R2仍超过0.3[23]，说明中国本地猪种在近代人工选育程度远不如西方商业猪种，也说明中国本地猪种的诸多优良性状还有待进一步开发利用。在本研究中，6个群体在12.68 kb距离时其连锁不平衡程度已经较低，因此针对最终筛选到的176个受选择位点，使用其上、下游12.68 kb作为候选基因定位的范围。

前人研究发现，多个基因与中国地方猪的繁殖性状相关，如孙浩[24]利用撒坝猪基因组信息检测到PRL、ADRA1B和BMP7等基因影响撒坝猪的繁殖力，管丽飞[25]利用测序技术和AS-PCR技术检测到精子鞭毛2(SPEF2)基因与香猪产仔数之间有关联，张倩[26]利用分子标记技术鉴别到微管蛋白δ1(TUBD1)基因是影响二花脸猪产仔数的候选基因等。本研究对176个受选择位点进行基因定位，发现66个蛋白编码基因位于受选择区间内，其中CHD7是编码染色体结构域解旋酶DNA结合蛋白7的基因，其在常染色体显性遗传(AD)中发现的变体与生长发育迟缓和生殖器异常相关，并且携带CHD7的家庭在其他AD中的突变表现出了先天性性腺机能减退[27]。Schrimpf等[28]利用11匹马进行下一代全基因组测序，对数据检索到的稀有高影响力变体进行基因分型，结果发现了FBXO43等9个与种马繁殖力相关的基因。有研究发现，FBXO43在减数分裂过程中起作用[29]，还影响呈圆周运动的精子数量[30]。RUNX1是一种核转录因子，在小鼠排卵前卵泡中被黄体生成素诱导，并促进孕激素的产生[31]，还可以刺激山羊颗粒细胞增殖和孕酮分泌[32]。Zhong等[33]在对RUNX1在受H3K27me3调节的猪颗粒细胞中增殖作用的研究发现，RUNX1通过调节与类固醇生成相关基因的表达而控制雌激素、雄激素和前列腺素的分泌，并且充当媒介，增强类固醇生成相关基因对人绒毛膜促性腺激素(HCG)的反应，促进核成熟和卵泡破裂而刺激排卵[34]。STRBP是精子核周RNA结合蛋白[35]，Zou等[36]对雷州黑鸭卵泡发育过程中产卵量差异进行研究，发现了STRBP等25种与鸭卵巢发育有关的基因，它们可能影响卵子的产生。Juma等[37]对小鼠精子研究发现，敲除TEDDM1的小鼠生育力下降[38]。此外，本研究依据产仔数划分的高、低组猪的体型差异也正好与产仔数性状差异分组相一致，但有研究证实，体型与产仔数之间并无较大的相关性[39-40]，因此CHD7、FBXO43、RUNX1、STRBP和TEDDM1基因可能是影响中国地方猪产仔数性状的候选基因。

4 结 论

本研究发现，6个中国地方猪种基因组的LD衰减速度均较快，顺序依次为梅山猪<蓝塘猪<姜曲海猪<明光小耳猪<五指山猪或迪庆藏猪；梅山猪群体内存在3个亲缘关系较高的小群体，明光小耳猪有部分个体与迪庆藏猪聚集在一起，各个猪种间的进化距离与地理距离较一致。通过选择信号分析发现，CHD7、FBXO43、RUNX1、STRBP和TEDDM1基因可能作为中国地方猪产仔数性状的候选基因。