基于SNP芯片的海南猪全基因组选择信号分析

冯雪燕，刁淑琪，刘玉强，徐志婷，魏 趁，袁晓龙，李加琪，张 哲*

(1.华南农业大学动物科学学院/国家生猪种业工程技术研究中心，广州 510642；2.岭南现代农业科学和技术广东省实验室茂名分中心，茂名 525000)

海南猪原产于我国海南省，该省地处我国最南端，属热带季风气候，终年高温多雨，光照充足；四面环海，饲料资源丰富。海南猪属华南型猪种，包含临高猪、屯昌猪、文昌猪和定安猪4个类群[1]。海南猪头小，耳小而薄、直立并稍向前倾。经过长期的自然选择和人工选育，海南猪形成了早熟易肥、皮薄骨细、肉质鲜美、耐热耐粗饲的特性[1]。然而，当地高温高湿的气候条件和交流封闭的地理环境在海南猪的驯化和适应性进化中的作用还不清楚，对海南猪特殊的发育规律和选择机制了解尚且不足。另外，海南猪突出的品种特质有望为猪的品种培育和种质利用工作提供优良的研究资料。在进化过程中，由于长时间的自然选择和人工选择作用在动物基因组上留下的结构特征变化，即为选择信号[2-3]。通过选择信号检测在动物群体的全基因组范围内捕获到存在的选择信号和潜在受选择位点，进一步研究种群的进化驯化历史，有助于揭示畜禽重要经济性状潜在的受选择和遗传机制，并已应用到地方猪的遗传育种工作中。目前，选择信号检测方法已经在猪[4-6]、鸡[7-8]、牛[9-10]、山羊[11]、马[12]等家养动物上得到研究。在对我国地方猪的选择信号研究上，已在通城猪[13]、金华猪[14]、莱芜猪[15-16]、皖南黑猪[17]、姜曲海猪[18]、滇南小耳花猪[19]等地方猪品种中鉴定到与猪的抗逆性状、生长和繁殖性状相关的候选基因。长片段纯合(runs of homozygosity，ROH)是二倍体生物基因中存在纯合基因型的连续片段的现象，它是由子代从亲代继承了相同单倍型而形成[20]。在家养动物育种工作中，有选择性地选种选配往往易选留存在亲缘关系的亲本个体，造成后代个体基因组上受选择区域纯合度增加，从而出现ROH片段[21-22]，对基因组上的高频ROH区段进行基因注释有利于得到受正向选择的基因[23]。利用群体ROH信息可以对动物的近交情况、群体历史等信息进行判断，有助于了解其遗传背景。目前，ROH检测已经在绵羊[24-25]、牛[26]、猪[27-28]等家养动物上得到应用。本研究关注海南猪全基因组范围内的选择信号，利用iHS方法、ROH检测和生物信息学分析的手段，旨在通过检测海南猪全基因组上的选择信号为相关研究提供参考。

1 材料与方法

1.1 试验群体

本试验以海南猪为研究对象，共68个个体，试验群体情况如表1所示。

表1 海南猪群体信息

本试验中，按照以下原则筛选待测地方猪个体：1)个体间无血缘关系；2)群体内包含更多的血统。

1.2 试验方法

1.2.1 基因组DNA的抽提 使用剪耳钳沿样本猪的耳朵边缘剪下0.5 g左右的耳样进行DNA抽提。基因组DNA抽提使用组织DNA提取试剂盒(E.Z.N.A. ©Tissue DNA Kit，D3396-02，Omega BioTek公司，美国)。基因组DNA抽提完毕后使用核酸浓度检测仪进行浓度检测，要求“DNA浓度>50 ng·μL-1，1.7

1.2.2 SNP芯片分型与数据质量控制 本试验中使用了GeneSeek Genomic Profiler(GGP)Procine SNP 80K芯片进行基因分型，共检测到57 485个常染色体SNPs位点。利用PLINK v1.90软件[29]对海南猪群体芯片数据进行质量控制，标准为：1)剔除位于性染色体和未知染色体上的位点；2)剔除SNP基因型检出率<0.90的个体；3)剔除基因型检出率<0.90的SNP位点；4)剔除品种内最小等位基因频率(minor allele frequency，MAF)<0.01的SNP位点。质控后剩余68个个体、44 578个SNPs位点用于后续分析。

1.2.3 主成分分析 为了检验本试验中的试验群体是否出现“分层现象”，对试验群体进行主成分分析(principle component analysis，PCA)。使用PLINK v1.90软件[29]筛选到22 815个相邻位点间连锁不平衡<0.5的SNPs位点(--indep-pairwise 50 5 0.5)用于主成分分析。随后使用GCTA v1.93.0软件[30]计算前3个主成分(--pca 3)，该软件中主成分的计算方法以Price等[31]提出的方法作为参照。使用R v4.0.2软件[32]对PCA分析结果进行可视化。

1.2.4 全基因组选择信号检测 本研究中，利用基于扩展单倍型纯合原理、适用于群体内选择信号检测的整合单倍型分数(integrated haplotype score，iHS)检验对海南猪群体进行全基因组选择信号检测。当iHS分数为较大负值时，表示该单倍型具有衍生新基因的可能；当iHS分数为较大正值时，表示该单倍型具有携带祖先等位基因的可能[2]。在实际应用中，基因组选择存在反复性，新衍生的等位基因和祖先等位基因均可能成为受选择的潜在作用位点，因此，iHS分数的绝对值经常被用于挖掘选择信号的实际研究中[33]。

iHS方法的计算公式如下：

其中，iHS表示单倍型积分值，iHHA表示祖先等位基因的EHH积分值，iHHD表示推断等位基因的EHH积分值。

本试验中，使用iHS软件[34]按位点计算iHS分数。首先利用fastPHASE软件[35]进行单倍型推测，随后利用R v4.0.2软件[32]将输出的单倍型整理成iHS软件[34]所需要的格式用于计算iHS分数，并按照Voight等[33]提出的校正方法对计算得到的iHS分数进行标准化，以“标准化iHS分数>1.96(P<0.05)”[36]为阈值筛选位点，得到潜在受选择位点。

1.2.5 ROH检测与近交系数计算 本研究使用PLINK v1.90软件[29]对海南猪进行全基因组ROH检测。检测参数为：1)每个滑动窗口50个SNP位点；2)窗口中纯合片段重合的比例大于0.05；3)每个ROH片段中最少有100个连续SNPs位点；4)SNP标记密度最小为每个SNP 50 kb；5)纯合片段两个SNPs之间的距离小于100 kb；6)在一个ROH片段内仅允许最多两个SNPs位点缺失和一个杂合子存在。

基于检测到的ROH计算基因组近交系数FROH：

其中，∑LROH为常染色体上ROH片段长度之和，Lauto为常染色体的物理总长度。

1.2.6 基因组注释与QTL探索 对iHS方法得到的潜在受选择位点向上、下游各自延伸200 kb[4]作为iHS方法得到的候选区域，定义其中与ROH检测得到的候选片段发生“完全重叠”的区段为潜在受选择区域，即此步骤得到的潜在受选择区域在iHS方法和ROH检测中均被检测到。

通过Ensembl数据库[37]选用11.1版本的猪基因组数据作为参考，以“基因在染色体的物理位置与潜在候选区域重叠”作为标准进行候选基因的注释，此步骤使用R v4.0.2软件[32]完成。进行QTL探索时，使用R v4.0.2软件[32]在Animal QTL数据库[38]中寻找与潜在受选择区域在物理位置上发生重叠的数量性状基因座(quantitative trait locus，QTL)，并统计该区段上每一个QTL的报导道次数。

1.2.7 GO功能及KEGG通路富集分析 以“P< 0.05”作为筛选标准，采用DAVID v6.8数据库(https://david.ncifcrf.gov/home.jsp)对潜在受选择区域内的基因进行GO功能及KEGG通路富集分析，其中，GO功能富集分析涵盖细胞组分(cellular component)、分子功能(molecular function)、生物学过程(biological process)3类条目。

2 结 果

2.1 主成分分析结果

对海南猪群体进行主成分分析，结果如图1所示，前3个主成分共解释了19.19%的遗传变异。分析结果显示，海南猪两个群体分层不明显，仅有部分个体出现分层，由于解释的遗传变异较小，因此在后续的分析中可作为一个整体进行全基因组选择信号检测分析。

2.2 海南猪群体全基因组选择信号检测

用iHS软件[34]进行选择信号检测，共得到17 686个SNPs位点的iHS分数，分布情况见图2A，为了便于分析和比较，将计算得到的iHS分数进行标准化，标准化iHS分数近似服从标准正态分布(图2B)；观察到大多数位点的标准化iHS分数介于-3～3 之间。

本试验中，以“标准化iHS分数>1.96(P<0.05)”作为选择信号阈值，在海南猪全基因组上共筛选到395个潜在受选择位点。潜在受选择位点不均匀地分布在1～18号染色体上，其中SSC9(SusScrofachromosome 9)、SSC12和SSC14上最多，SSC9上筛选到81个潜在受选择位点，占总潜在受选择位点数的20.51%。海南猪全基因组范围内的选择信号iHS分数分布情况如图2C所示，最高的iHS分数的受选择位点位于9号染色体上，其标准化iHS分数为4.104。本研究已将分析获得的全部潜在受选择位点上传至Figshare公共数据库(网址：https://figshare.com/s/dbaf15b23b9e07e0f991)。

2.3 海南猪全基因组ROH检测及近交系数计算

在海南猪全基因组上共检测到172个ROH(图3)片段，ROH片段在1～18号染色体上均有分布，位于SSC12上ROH片段数量最多，出现频率最高的ROH片段位于SSC13(200.93～201.51 Mb)上，最长ROH片段长度为12.39 Mb，位于SSC16上。ROH检测分析完整结果已上传至Figshare数据库(https://figshare.com/s/0a84906c1e9 da3be7a02)。根据ROH检测结果计算近交系数，定安猪FROH=0.039 94，屯昌猪FROH=0.014 17，海南猪FROH=0.028 53。

2.4 潜在受选择区域基因组注释结果与潜在受选择区域重叠的QTL

将iHS方法得到的潜在受选择位点向上、下游各自延伸200 kb[4]作为iHS方法得到的候选区域，定义其中与ROH检测得到的候选片段发生“完全重叠”的片段为潜在受选择区域。对得到的136个潜在受选择区域进行基因注释，共注释到469个候选基因。其中，前10个潜在受选择区域共注释到41个候选基因(表2)，全部潜在受选择区域的注释情况已上传至Figshare公共数据库(网址：https://figshare.com/s/1cfa3f0fcad37faaedc0)。

A.主成分1-主成分2; B.主成分1-主成分3; C.主成分2-主成分3; D.前3个主成分三维示意图A. PC1-PC2; B. PC1-PC3; C. PC2-PC3; D. The three-dimensional diagram of the first three PCs图1 海南猪群体主成分分析(PCA)Fig.1 Principal component analysis (PCA) of Hainan pigs population

A.原始iHS分数频数分布；B.标准化iHS分数频数分布；C.全基因组iHS选择信号分析(常染色体)：横线表示95%置信水平, 横线上方的位点为检测到的潜在受选择位点A. The distribution of statistics frequency for raw iHS scores; B. The distribution of statistics frequency for standard iHS scores; C. Genome-wide analysis of selection signatures detected by iHS (autosomes): Horizontal line in figureC displays the threshold levels of 5%, the points upon the horizontal line display the potential selected SNPs detected图2 海南猪全基因组选择信号分布Fig.2 Distribution of selection signature identified on the whole genome of Hainan pigs

图3 海南猪常染色体ROH分布Fig.3 Distribution of ROH identified on chromosomes of Hainan pigs

对本试验检测到的潜在受选择区域与Animal QTL数据库[38]进行比对，共有382个QTLs与检测到的潜在受选择区域发生重叠。其中，与胸膜肺炎放线杆菌易感性相关的QTL在93处潜在受选择区域被报道。同时，对比得到的QTL大多与影响海南猪的平均日增重(131处、424次报道)、平均背膘厚(109处、232次报道)、滴水损失(94处、364次报道)、乳头数量(125处、208次报道)等重要经济性状的QTL重叠。此外，本研究结果显示影响海南猪初情期启动日龄的QTL主要和海南猪SSC12上的潜在受选择区域重叠。

2.5 GO功能和KEGG通路分析

利用DAVID数据库，对注释得到的469个基因进行GO功能和KEGG通路富集分析，共有246个基因被识别。按照“P<0.05”为显著条件筛选到21项条目，在生物学过程、细胞组分、分子功能条目上各富集到5、2、4项条目，在KEGG通路上富集到10项条目。表3展示了前10项富集条目，本研究富集分析的完整结果已上传至Figshare数据库(https://figshare.com/s/aa6f4435ce4ca4 db6413)。基因富集结果显示，共有91个注释基因共同富集在生物学过程、细胞组分及分子功能3大类GO类别和KEGG通路类别上。经分析发现，海南猪全基因组的潜在受选择区域上的基因功能主要集中在生长代谢和免疫抗病相关的通路上，其中FGFR2基因富集在磷脂酰肌醇3激酶(PI3K)/AKT(蛋白激酶B)信号通路(P=0.007 46)和ATP结合通路(P=0.014 90)等5条通路上；另外，在富集得到的显著通路中，基因ESR1、PIK3CG、MAPK3共同在繁殖性状相关的雌激素信号通路(P=0.032 28)和催乳素信号通路(P=0.049 82)上发生富集。

3 讨 论

本研究利用iHS方法对海南猪种群体(定安猪和屯昌猪2个类群)进行群体内全基因组选择信号检测，检测到395个潜在受选择位点，并向上、下游各扩充200 kb得到候选区域。同时进行ROH检测得到172个片段，对iHS方法和ROH检测得到的片段以发生“完全重叠”为标准得到136个潜在受选择区域，对潜在选择区域进行基因注释和QTL探索分析，共注释到469个候选基因，QTL探索分析结果显示，与潜在受选择区域发生重叠的QTL大多与猪的平均日增重、眼肌面积、平均背膘厚度及其他肉质性状相关。

本研究中，与耳朵直立相关的QTL与潜在受选择区域有39处发生重叠(报道43次)，与耳朵重量、大小、区域性状相关的QTL与潜在受选择区域共有56处发生重叠(报道92次)，此结果与海南猪耳小而薄、直立并稍向前倾[39]的外貌特征一致。海南猪种性成熟早，母猪3～4月龄达初情期，7～8月龄 可配种[39]。本研究中，影响初情期启动日龄的QTL与SSC12上的潜在受选择区域重叠，分别在15个潜在区域中被注释到，总报道次数为17次，研究结果符合海南猪性成熟早的特点。本研究中，在SSC17的35.92～36.24 Mb上注释到NOL4L基因，该基因已在巴马香猪的17号染色体上被鉴定为与繁殖性状相关的基因[40]。另外，有多个有关肉质指标的QTL在潜在受选择区域中被报道多次，如影响滴水损失、肉色、肌肉含水量、大理石纹和肌间脂肪含量的QTL等，分析结果与海南猪具有的肉质优良特性相一致；其中，影响滴水损失的QTL主要与SSC1上的潜在受选择区域重叠，在全基因组范围内共在94个潜在受选择区域上发生364次报道。

表3 海南猪候选基因富集分析结果(前10条通路)

本研究中，SSC1上的潜在候选基因ESR1基因在催乳素信号通路(ssc04917:Prolactin signaling pathway，P=0.498 2)和雌激素信号通路(ssc04915:Estrogen signaling pathway，P=0.032 99)中发生显著富集，该基因同样在Wang等[16]的研究中被鉴定到，已有相关研究发现该基因突变与猪的产仔数性状相关[41]。iHS方法检测到具有最高iHS分数的位点位于SSC14上，在该位点所在区段上注释到FGFR2基因，该基因在GO和KEGG通路富集分析中在GO:0005524～ATP binding(P=0.014 90)、ssc04151:PI3K-Akt signaling pathway(P=0.007 46)、ssc05215:Prostate cancer(P=0.021 82)、ssc04550:Signaling pathways regulating pluripotency of stem cells(P=0.023 49)和ssc04810:Regulation of actin cytoskeleton(P=0.036 33)通路上发生显著富集，Lu等[42]的研究表明该基因与乳腺分支形态发生过程相关。

本试验中，除了对iHS方法中根据“标准化iHS分数>1.96(P<0.05)”得到的潜在受选择位点所在的候选区域与ROH检测得到的片段发生“完全重叠”的区段作为潜在受选择区域，并进行QTL探索、基因注释及下游富集分析外，单独对iHS方法中得到的潜在受选择位点所在区域(395个)进行了QTL探索、基因注释及下游富集分析，对应的QTL探索及基因注释结果已经上传至Figshare数据库(https://figshare.com/s/67f725a3c9887b71fcf6)，对应的基因富集分析结果已上传至Figshare数据库(https://figshare.com/s/ad7ff315257b9a523058)。

基因注释和通路富集结果显示，在SSC17上36.44～36.84 Mb区域注释到的候选基因有BPIFB2、BPIFB6、BPIFB4、BPIFB3、BPIFA2、BPIFA3、BPIFA1和BPIFB1，与王亚楠[43]在莱芜猪基因组上进行选择信号分析过程中发现的多个参与脂类代谢过程的基因存在重合，其中包括BPIFB2、BPIFA1、BPIFA2和BPIFB4基因，因此推测该区域内的候选基因可能与海南猪的高脂肪含量有关。此外，与猪产仔数性状相关的ESR1基因在ssc04961:Endocrine and other factor-regulated calcium reabsorption(P=0.005 17)、GO:0003677～DNA binding(P=0.025 97)、ssc04919:Thyroid hormone signaling pathway(P=0.036 98)和ssc04915:Estrogen signaling pathway(P=0.049 81)上发生显著富集。

在SSC9的57.1～57.5 Mb片段上注释到的ADAMTS8和ADAMTS15基因，均属于ADAMTS家族，该家族基因多与肿瘤疾病有关，有研究表明这两个基因对恶性肿瘤有抑制作用[44]；此外，对候选基因进行GO功能和KEGG通路富集分析结果显示，趋化因子CCL2在与美洲锥虫病发生(ssc05142:Chagas disease (American trypanosomiasis)，P=0.028 63)相关的条目上发生显著富集，本研究在SSC12上40.33～41.00 Mb区段上多次注释到的趋化因子CCL8、CCL11、CCL2、CCL1基因，属巨噬细胞趋化因子、T淋巴细胞趋化因子，在多个研究中被证实参与巨噬细胞炎症发生、肿瘤发生和复发等疾病通路的调节[45-48]。以上结果均表明，海南猪在其进化驯化过程中曾在抗逆抗病性状上受到过强力选择。

4 结 论

本研究以海南猪GeneSeek Genomic Profiler Procine SNP 80K芯片基因组数据为对象，利用iHS方法对海南猪群体的常染色体进行选择信号检测，共检测到395个潜在受选择位点，ROH检测得到172个片段；在发生“完全重叠”的136个潜在受选择区域上注释到469个候选基因。对潜在受选择区域的QTL探索结果揭示，海南猪的选择信号多与肉质、生长和抗病性状相关，与海南猪的初情期启动日龄相关的QTL主要与SSC12上的潜在受选择区域重叠；对候选基因进行GO功能和KEGG通路富集分析的结果显示，有91个候选基因在21条GO功能和KEGG通路条目上发生显着富集，主要集中在生长代谢、免疫应答和雌激素信号通路相关的条目上。本研究揭示了海南猪群体在其进化驯化历史上可能受到的选择情况。研究结果表明，在海南猪全基因组范围内进行选择信号检测，有助于进一步了解海南猪的进化历史及其重要经济性状的遗传机制，在一定程度上为华南型地方猪种质资源的保存利用和相关研究提供参考。