基于SNP芯片的海南豬全基因組選擇信號分析

馮雪燕，刁淑琪，劉玉強，徐志婷，魏 趁，袁曉龍，李加琪，張 哲*

(1.華南農業(yè)大學動物科學學院/國家生豬種業(yè)工程技術研究中心，廣州 510642；2.嶺南現(xiàn)代農業(yè)科學和技術廣東省實驗室茂名分中心，茂名 525000)

海南豬原產于我國海南省，該省地處我國最南端，屬熱帶季風氣候，終年高溫多雨，光照充足；四面環(huán)海，飼料資源豐富。海南豬屬華南型豬種，包含臨高豬、屯昌豬、文昌豬和定安豬4個類群[1]。海南豬頭小，耳小而薄、直立并稍向前傾。經過長期的自然選擇和人工選育，海南豬形成了早熟易肥、皮薄骨細、肉質鮮美、耐熱耐粗飼的特性[1]。然而，當?shù)馗邷馗邼竦臍夂驐l件和交流封閉的地理環(huán)境在海南豬的馴化和適應性進化中的作用還不清楚，對海南豬特殊的發(fā)育規(guī)律和選擇機制了解尚且不足。另外，海南豬突出的品種特質有望為豬的品種培育和種質利用工作提供優(yōu)良的研究資料。在進化過程中，由于長時間的自然選擇和人工選擇作用在動物基因組上留下的結構特征變化，即為選擇信號[2-3]。通過選擇信號檢測在動物群體的全基因組范圍內捕獲到存在的選擇信號和潛在受選擇位點，進一步研究種群的進化馴化歷史，有助于揭示畜禽重要經濟性狀潛在的受選擇和遺傳機制，并已應用到地方豬的遺傳育種工作中。目前，選擇信號檢測方法已經在豬[4-6]、雞[7-8]、牛[9-10]、山羊[11]、馬[12]等家養(yǎng)動物上得到研究。在對我國地方豬的選擇信號研究上，已在通城豬[13]、金華豬[14]、萊蕪豬[15-16]、皖南黑豬[17]、姜曲海豬[18]、滇南小耳花豬[19]等地方豬品種中鑒定到與豬的抗逆性狀、生長和繁殖性狀相關的候選基因。長片段純合(runs of homozygosity，ROH)是二倍體生物基因中存在純合基因型的連續(xù)片段的現(xiàn)象，它是由子代從親代繼承了相同單倍型而形成[20]。在家養(yǎng)動物育種工作中，有選擇性地選種選配往往易選留存在親緣關系的親本個體，造成后代個體基因組上受選擇區(qū)域純合度增加，從而出現(xiàn)ROH片段[21-22]，對基因組上的高頻ROH區(qū)段進行基因注釋有利于得到受正向選擇的基因[23]。利用群體ROH信息可以對動物的近交情況、群體歷史等信息進行判斷，有助于了解其遺傳背景。目前，ROH檢測已經在綿羊[24-25]、牛[26]、豬[27-28]等家養(yǎng)動物上得到應用。本研究關注海南豬全基因組范圍內的選擇信號，利用iHS方法、ROH檢測和生物信息學分析的手段，旨在通過檢測海南豬全基因組上的選擇信號為相關研究提供參考。

1 材料與方法

1.1 試驗群體

本試驗以海南豬為研究對象，共68個個體，試驗群體情況如表1所示。

表1 海南豬群體信息

本試驗中，按照以下原則篩選待測地方豬個體：1)個體間無血緣關系；2)群體內包含更多的血統(tǒng)。

1.2 試驗方法

1.2.1 基因組DNA的抽提 使用剪耳鉗沿樣本豬的耳朵邊緣剪下0.5 g左右的耳樣進行DNA抽提?；蚪MDNA抽提使用組織DNA提取試劑盒(E.Z.N.A. ?Tissue DNA Kit，D3396-02，Omega BioTek公司，美國)。基因組DNA抽提完畢后使用核酸濃度檢測儀進行濃度檢測，要求“DNA濃度>50 ng·μL-1，1.7

1.2.2 SNP芯片分型與數(shù)據質量控制 本試驗中使用了GeneSeek Genomic Profiler(GGP)Procine SNP 80K芯片進行基因分型，共檢測到57 485個常染色體SNPs位點。利用PLINK v1.90軟件[29]對海南豬群體芯片數(shù)據進行質量控制，標準為：1)剔除位于性染色體和未知染色體上的位點；2)剔除SNP基因型檢出率<0.90的個體；3)剔除基因型檢出率<0.90的SNP位點；4)剔除品種內最小等位基因頻率(minor allele frequency，MAF)<0.01的SNP位點。質控后剩余68個個體、44 578個SNPs位點用于后續(xù)分析。

1.2.3 主成分分析 為了檢驗本試驗中的試驗群體是否出現(xiàn)“分層現(xiàn)象”，對試驗群體進行主成分分析(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 全基因組選擇信號檢測 本研究中，利用基于擴展單倍型純合原理、適用于群體內選擇信號檢測的整合單倍型分數(shù)(integrated haplotype score，iHS)檢驗對海南豬群體進行全基因組選擇信號檢測。當iHS分數(shù)為較大負值時，表示該單倍型具有衍生新基因的可能；當iHS分數(shù)為較大正值時，表示該單倍型具有攜帶祖先等位基因的可能[2]。在實際應用中，基因組選擇存在反復性，新衍生的等位基因和祖先等位基因均可能成為受選擇的潛在作用位點，因此，iHS分數(shù)的絕對值經常被用于挖掘選擇信號的實際研究中[33]。

iHS方法的計算公式如下：

其中，iHS表示單倍型積分值，iHHA表示祖先等位基因的EHH積分值，iHHD表示推斷等位基因的EHH積分值。

本試驗中，使用iHS軟件[34]按位點計算iHS分數(shù)。首先利用fastPHASE軟件[35]進行單倍型推測，隨后利用R v4.0.2軟件[32]將輸出的單倍型整理成iHS軟件[34]所需要的格式用于計算iHS分數(shù)，并按照Voight等[33]提出的校正方法對計算得到的iHS分數(shù)進行標準化，以“標準化iHS分數(shù)>1.96(P<0.05)”[36]為閾值篩選位點，得到潛在受選擇位點。

1.2.5 ROH檢測與近交系數(shù)計算 本研究使用PLINK v1.90軟件[29]對海南豬進行全基因組ROH檢測。檢測參數(shù)為：1)每個滑動窗口50個SNP位點；2)窗口中純合片段重合的比例大于0.05；3)每個ROH片段中最少有100個連續(xù)SNPs位點；4)SNP標記密度最小為每個SNP 50 kb；5)純合片段兩個SNPs之間的距離小于100 kb；6)在一個ROH片段內僅允許最多兩個SNPs位點缺失和一個雜合子存在。

基于檢測到的ROH計算基因組近交系數(shù)FROH：

其中，∑LROH為常染色體上ROH片段長度之和，Lauto為常染色體的物理總長度。

1.2.6 基因組注釋與QTL探索 對iHS方法得到的潛在受選擇位點向上、下游各自延伸200 kb[4]作為iHS方法得到的候選區(qū)域，定義其中與ROH檢測得到的候選片段發(fā)生“完全重疊”的區(qū)段為潛在受選擇區(qū)域，即此步驟得到的潛在受選擇區(qū)域在iHS方法和ROH檢測中均被檢測到。

通過Ensembl數(shù)據庫[37]選用11.1版本的豬基因組數(shù)據作為參考，以“基因在染色體的物理位置與潛在候選區(qū)域重疊”作為標準進行候選基因的注釋，此步驟使用R v4.0.2軟件[32]完成。進行QTL探索時，使用R v4.0.2軟件[32]在Animal QTL數(shù)據庫[38]中尋找與潛在受選擇區(qū)域在物理位置上發(fā)生重疊的數(shù)量性狀基因座(quantitative trait locus，QTL)，并統(tǒng)計該區(qū)段上每一個QTL的報導道次數(shù)。

1.2.7 GO功能及KEGG通路富集分析 以“P< 0.05”作為篩選標準，采用DAVID v6.8數(shù)據庫(https://david.ncifcrf.gov/home.jsp)對潛在受選擇區(qū)域內的基因進行GO功能及KEGG通路富集分析，其中，GO功能富集分析涵蓋細胞組分(cellular component)、分子功能(molecular function)、生物學過程(biological process)3類條目。

2 結 果

2.1 主成分分析結果

對海南豬群體進行主成分分析，結果如圖1所示，前3個主成分共解釋了19.19%的遺傳變異。分析結果顯示，海南豬兩個群體分層不明顯，僅有部分個體出現(xiàn)分層，由于解釋的遺傳變異較小，因此在后續(xù)的分析中可作為一個整體進行全基因組選擇信號檢測分析。

2.2 海南豬群體全基因組選擇信號檢測

用iHS軟件[34]進行選擇信號檢測，共得到17 686個SNPs位點的iHS分數(shù)，分布情況見圖2A，為了便于分析和比較，將計算得到的iHS分數(shù)進行標準化，標準化iHS分數(shù)近似服從標準正態(tài)分布(圖2B)；觀察到大多數(shù)位點的標準化iHS分數(shù)介于-3～3 之間。

本試驗中，以“標準化iHS分數(shù)>1.96(P<0.05)”作為選擇信號閾值，在海南豬全基因組上共篩選到395個潛在受選擇位點。潛在受選擇位點不均勻地分布在1～18號染色體上，其中SSC9(SusScrofachromosome 9)、SSC12和SSC14上最多，SSC9上篩選到81個潛在受選擇位點，占總潛在受選擇位點數(shù)的20.51%。海南豬全基因組范圍內的選擇信號iHS分數(shù)分布情況如圖2C所示，最高的iHS分數(shù)的受選擇位點位于9號染色體上，其標準化iHS分數(shù)為4.104。本研究已將分析獲得的全部潛在受選擇位點上傳至Figshare公共數(shù)據庫(網址：https://figshare.com/s/dbaf15b23b9e07e0f991)。

2.3 海南豬全基因組ROH檢測及近交系數(shù)計算

在海南豬全基因組上共檢測到172個ROH(圖3)片段，ROH片段在1～18號染色體上均有分布，位于SSC12上ROH片段數(shù)量最多，出現(xiàn)頻率最高的ROH片段位于SSC13(200.93～201.51 Mb)上，最長ROH片段長度為12.39 Mb，位于SSC16上。ROH檢測分析完整結果已上傳至Figshare數(shù)據庫(https://figshare.com/s/0a84906c1e9 da3be7a02)。根據ROH檢測結果計算近交系數(shù)，定安豬FROH=0.039 94，屯昌豬FROH=0.014 17，海南豬FROH=0.028 53。

2.4 潛在受選擇區(qū)域基因組注釋結果與潛在受選擇區(qū)域重疊的QTL

將iHS方法得到的潛在受選擇位點向上、下游各自延伸200 kb[4]作為iHS方法得到的候選區(qū)域，定義其中與ROH檢測得到的候選片段發(fā)生“完全重疊”的片段為潛在受選擇區(qū)域。對得到的136個潛在受選擇區(qū)域進行基因注釋，共注釋到469個候選基因。其中，前10個潛在受選擇區(qū)域共注釋到41個候選基因(表2)，全部潛在受選擇區(qū)域的注釋情況已上傳至Figshare公共數(shù)據庫(網址：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分數(shù)頻數(shù)分布；B.標準化iHS分數(shù)頻數(shù)分布；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

對本試驗檢測到的潛在受選擇區(qū)域與Animal QTL數(shù)據庫[38]進行比對，共有382個QTLs與檢測到的潛在受選擇區(qū)域發(fā)生重疊。其中，與胸膜肺炎放線桿菌易感性相關的QTL在93處潛在受選擇區(qū)域被報道。同時，對比得到的QTL大多與影響海南豬的平均日增重(131處、424次報道)、平均背膘厚(109處、232次報道)、滴水損失(94處、364次報道)、乳頭數(shù)量(125處、208次報道)等重要經濟性狀的QTL重疊。此外，本研究結果顯示影響海南豬初情期啟動日齡的QTL主要和海南豬SSC12上的潛在受選擇區(qū)域重疊。

2.5 GO功能和KEGG通路分析

利用DAVID數(shù)據庫，對注釋得到的469個基因進行GO功能和KEGG通路富集分析，共有246個基因被識別。按照“P<0.05”為顯著條件篩選到21項條目，在生物學過程、細胞組分、分子功能條目上各富集到5、2、4項條目，在KEGG通路上富集到10項條目。表3展示了前10項富集條目，本研究富集分析的完整結果已上傳至Figshare數(shù)據庫(https://figshare.com/s/aa6f4435ce4ca4 db6413)?；蚋患Y果顯示，共有91個注釋基因共同富集在生物學過程、細胞組分及分子功能3大類GO類別和KEGG通路類別上。經分析發(fā)現(xiàn)，海南豬全基因組的潛在受選擇區(qū)域上的基因功能主要集中在生長代謝和免疫抗病相關的通路上，其中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)上發(fā)生富集。

3 討 論

本研究利用iHS方法對海南豬種群體(定安豬和屯昌豬2個類群)進行群體內全基因組選擇信號檢測，檢測到395個潛在受選擇位點，并向上、下游各擴充200 kb得到候選區(qū)域。同時進行ROH檢測得到172個片段，對iHS方法和ROH檢測得到的片段以發(fā)生“完全重疊”為標準得到136個潛在受選擇區(qū)域，對潛在選擇區(qū)域進行基因注釋和QTL探索分析，共注釋到469個候選基因，QTL探索分析結果顯示，與潛在受選擇區(qū)域發(fā)生重疊的QTL大多與豬的平均日增重、眼肌面積、平均背膘厚度及其他肉質性狀相關。

本研究中，與耳朵直立相關的QTL與潛在受選擇區(qū)域有39處發(fā)生重疊(報道43次)，與耳朵重量、大小、區(qū)域性狀相關的QTL與潛在受選擇區(qū)域共有56處發(fā)生重疊(報道92次)，此結果與海南豬耳小而薄、直立并稍向前傾[39]的外貌特征一致。海南豬種性成熟早，母豬3～4月齡達初情期，7～8月齡 可配種[39]。本研究中，影響初情期啟動日齡的QTL與SSC12上的潛在受選擇區(qū)域重疊，分別在15個潛在區(qū)域中被注釋到，總報道次數(shù)為17次，研究結果符合海南豬性成熟早的特點。本研究中，在SSC17的35.92～36.24 Mb上注釋到NOL4L基因，該基因已在巴馬香豬的17號染色體上被鑒定為與繁殖性狀相關的基因[40]。另外，有多個有關肉質指標的QTL在潛在受選擇區(qū)域中被報道多次，如影響滴水損失、肉色、肌肉含水量、大理石紋和肌間脂肪含量的QTL等，分析結果與海南豬具有的肉質優(yōu)良特性相一致；其中，影響滴水損失的QTL主要與SSC1上的潛在受選擇區(qū)域重疊，在全基因組范圍內共在94個潛在受選擇區(qū)域上發(fā)生364次報道。

表3 海南豬候選基因富集分析結果(前10條通路)

本研究中，SSC1上的潛在候選基因ESR1基因在催乳素信號通路(ssc04917:Prolactin signaling pathway，P=0.498 2)和雌激素信號通路(ssc04915:Estrogen signaling pathway，P=0.032 99)中發(fā)生顯著富集，該基因同樣在Wang等[16]的研究中被鑒定到，已有相關研究發(fā)現(xiàn)該基因突變與豬的產仔數(shù)性狀相關[41]。iHS方法檢測到具有最高iHS分數(shù)的位點位于SSC14上，在該位點所在區(qū)段上注釋到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)通路上發(fā)生顯著富集，Lu等[42]的研究表明該基因與乳腺分支形態(tài)發(fā)生過程相關。

本試驗中，除了對iHS方法中根據“標準化iHS分數(shù)>1.96(P<0.05)”得到的潛在受選擇位點所在的候選區(qū)域與ROH檢測得到的片段發(fā)生“完全重疊”的區(qū)段作為潛在受選擇區(qū)域，并進行QTL探索、基因注釋及下游富集分析外，單獨對iHS方法中得到的潛在受選擇位點所在區(qū)域(395個)進行了QTL探索、基因注釋及下游富集分析，對應的QTL探索及基因注釋結果已經上傳至Figshare數(shù)據庫(https://figshare.com/s/67f725a3c9887b71fcf6)，對應的基因富集分析結果已上傳至Figshare數(shù)據庫(https://figshare.com/s/ad7ff315257b9a523058)。

基因注釋和通路富集結果顯示，在SSC17上36.44～36.84 Mb區(qū)域注釋到的候選基因有BPIFB2、BPIFB6、BPIFB4、BPIFB3、BPIFA2、BPIFA3、BPIFA1和BPIFB1，與王亞楠[43]在萊蕪豬基因組上進行選擇信號分析過程中發(fā)現(xiàn)的多個參與脂類代謝過程的基因存在重合，其中包括BPIFB2、BPIFA1、BPIFA2和BPIFB4基因，因此推測該區(qū)域內的候選基因可能與海南豬的高脂肪含量有關。此外，與豬產仔數(shù)性狀相關的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)上發(fā)生顯著富集。

在SSC9的57.1～57.5 Mb片段上注釋到的ADAMTS8和ADAMTS15基因，均屬于ADAMTS家族，該家族基因多與腫瘤疾病有關，有研究表明這兩個基因對惡性腫瘤有抑制作用[44]；此外，對候選基因進行GO功能和KEGG通路富集分析結果顯示，趨化因子CCL2在與美洲錐蟲病發(fā)生(ssc05142:Chagas disease (American trypanosomiasis)，P=0.028 63)相關的條目上發(fā)生顯著富集，本研究在SSC12上40.33～41.00 Mb區(qū)段上多次注釋到的趨化因子CCL8、CCL11、CCL2、CCL1基因，屬巨噬細胞趨化因子、T淋巴細胞趨化因子，在多個研究中被證實參與巨噬細胞炎癥發(fā)生、腫瘤發(fā)生和復發(fā)等疾病通路的調節(jié)[45-48]。以上結果均表明，海南豬在其進化馴化過程中曾在抗逆抗病性狀上受到過強力選擇。

4 結 論

本研究以海南豬GeneSeek Genomic Profiler Procine SNP 80K芯片基因組數(shù)據為對象，利用iHS方法對海南豬群體的常染色體進行選擇信號檢測，共檢測到395個潛在受選擇位點，ROH檢測得到172個片段；在發(fā)生“完全重疊”的136個潛在受選擇區(qū)域上注釋到469個候選基因。對潛在受選擇區(qū)域的QTL探索結果揭示，海南豬的選擇信號多與肉質、生長和抗病性狀相關，與海南豬的初情期啟動日齡相關的QTL主要與SSC12上的潛在受選擇區(qū)域重疊；對候選基因進行GO功能和KEGG通路富集分析的結果顯示，有91個候選基因在21條GO功能和KEGG通路條目上發(fā)生顯著富集，主要集中在生長代謝、免疫應答和雌激素信號通路相關的條目上。本研究揭示了海南豬群體在其進化馴化歷史上可能受到的選擇情況。研究結果表明，在海南豬全基因組范圍內進行選擇信號檢測，有助于進一步了解海南豬的進化歷史及其重要經濟性狀的遺傳機制，在一定程度上為華南型地方豬種質資源的保存利用和相關研究提供參考。