基于全基因組填充重測序關聯分析鑒別影響蘇山豬初生體尺與乳頭數性狀的遺傳位點

張昌政，李德森，黃 敏，方曉敏，趙為民，任守文，董煥聲 任 軍 周李生*

(1.青島農業大學動物科技學院，青島 266109； 2.華南農業大學動物科學學院，廣州 510642；3.江蘇省農業科學院畜牧研究所，南京 210014)

豬的初生體尺和乳頭數性狀是衡量豬生長和繁殖性能的重要指標，直接影響養豬業的經濟效益。已有研究表明，初生體尺作為衡量仔豬生存能力和生長發育的可靠預測指標，與后期的生長速度和成活率呈顯著正相關[1-3]；初生重較高的仔豬在達到屠宰體重后，具有更好的胴體品質[4]和更優秀的體尺特征[5]。有效乳頭數是評價母豬哺乳能力的基礎[6]，當能繁母豬的有效乳頭數無法滿足后代的哺乳需求，仔豬死亡率將顯著增加，從而造成較大的經濟損失[7-8]。同時，乳頭數與產仔數存在一定程度的表型相關和遺傳相關[9-10]，核心母豬群體的平均有效乳頭數與產仔數存在較為一致的趨勢[11]。因此，鑒別影響仔豬的初生體尺和乳頭數性狀的分子標記有助于豬生長和繁殖性能的遺傳改良。

蘇山豬是江蘇省農業科學院畜牧研究所團隊利用太湖豬、丹麥長白豬和英國約克夏豬雜交育成的優質瘦肉豬新品種。本研究采用全基因組重測序結合基因型填充策略，在蘇山豬群體中開展Fst與GWAS分析，挖掘與初生體尺和乳頭數性狀顯著關聯的位點與候選基因，為探究影響豬生長發育和乳頭數形成的遺傳機制提供借鑒，且為蘇山豬核心種群的進一步選育提高奠定理論基礎。

1 材料與方法

1.1 樣本采集及表型記錄

試驗群體為269頭蘇山豬初生仔豬，均飼養于江蘇省農業科學院六合基地實驗豬場。仔豬在出生后24 h內測定并記錄體重(body weight, BW)、體長(body length, BL)、胸圍(chest circumference, ChC)、管圍(cannon circumference, CaC)、左乳頭數(number of teats on the left side, LTN)、右乳頭數(number of teats on the right side, RTN)、總乳頭數(total teat number, TTN)表型數據，并采集耳組織樣品保存于含75%乙醇的2 mL離心管中。

體尺的測量方法：1)體長：在豬自然站立狀態下，使用軟尺自兩耳連線的額頂中心起，沿背線至尾根緊貼體表量??；2)胸圍：在豬自然站立狀態下，切于肩胛骨后緣用軟尺測量胸部垂直周徑，軟尺自然貼近體表為宜；3)管圍：在豬自然站立狀態下，使用軟尺測量左前肢管部最細處的周徑，測量時軟尺緊貼體表。

1.2 DNA提取、基因分型和質量控制

使用常規酚氯仿法從豬耳組織樣本中提取基因組DNA，使用Affymetrix Porcine SNP 55 K Array進行基因分型，并用Plink v1.9軟件對所得SNP數據進行質量控制，標準：剔除SNPs檢出率(call rate)小于90%，次等位基因頻率(minor allele frequency, MAF)小于0.01，最后得到39 032個SNPs用于后續分析。

1.3 全基因組重測序分析

在269頭蘇山豬中挑選20個能夠代表群體全部血緣的個體，采用Illumina NovaSeq platforms 的150 bp paired-end libraries，進行全基因組重測序(平均測序深度為20×)。此外，在NCBI(https://www.ncbi.nlm.nih.gov/)上下載了367頭豬重測序數據，其中包括184頭中國地方豬，13頭中國野豬，142頭歐洲商品豬，21頭歐洲野豬和7頭外群豬種。針對測序和下載的共387個個體的原始數據(raw data)進行質控，得到能滿足數據分析要求的有效數據(clean data)。

1.4 變異位點檢測

1.5 基因型填充

將269頭蘇山豬的55K SNP芯片數據作為目標數據集，387個重測序個體的SNP數據集作為參考數據集，利用Beagle軟件進行基因型填充。隨后利用等位基因準確性(allelic correct rate, CR)和基因型相關系數(correlation)對20頭同時具有芯片數據和重測序數據的蘇山豬進行基因型填充準確度的估計。最后綜合基因型準確度和SNP的分布密度，保留call rate>0.9、MAF>0.01和correlation>0.8的SNPs用于后續分析。

1.6 選擇信號分析

按照每個性狀表型值的前10%和后10%分別分成兩個亞群，基于填充后的SNP數據集，利用VCFtools v0.1.17軟件，窗口大小為50 kb，以20 kb步長對每個SNPs進行滑窗計算Fst值?；贔st的結果，使用R語言的Cmplot程序包繪制曼哈頓圖，Fst閾值為按大小排序Fst值的前1%，高于閾值的點視為顯著分化的區域。

1.7 全基因組關聯分析

基于基因型填充后的SNP數據集，利用GEMMA軟件對蘇山豬的7個性狀進行GWAS，為了校正由于親緣關系、性別以及胎次等其他因素對關聯分析造成的影響，本研究使用混合線性模型(mixed linear model, MLM)進行分析，模型如下：

y=Wα+Xβ+u+ε

y為表型向量；W為協變量(固定效應)的關聯矩陣，將體重、性別和胎次作為初生體尺性狀的協變量；α為包括截距在內的相應系數的向量；X為SNP矩陣；β為SNP的效應；u為一個m×1的隨機效應向量，其中u～MVNm(0,λτ-1K)，τ-1為殘差的方差，λ為τ-1和ε的比值，K為n×n親屬關系矩陣；ε為隨機殘差的向量，ε～MVNn(0,τ-1In)，In為單位矩陣；MVNn表示n維多元正態分布。

Bonferroni方法常被用來界定顯著性閾值：0.05/N，N為研究SNPs的總數?？紤]到Bonferroni校正對于本研究過于嚴格，則根據前人研究自定義P=1×10-5[12]為顯著閾值線。使用R語言的Cmplot程序包對GWAS結果進行可視化。

1.8 基因功能富集分析

基于GWAS和Fst的結果，篩選出共有的顯著SNPs，以這些SNPs為中心，利用R語言的GALLO程序包，參考基因組為11.1版本的豬基因組(Sus scrofa 11.1)，在其上、下游500 kb范圍內搜尋候選基因。為更好的了解基因的功能，精準地鑒別出與目標性狀相關的候選基因，本研究利用Metascape(http://metascape.org)數據庫對搜尋到的基因進行GO(gene ontology)分析和KEGG(kyoto encyclopedia of genes and genomes)通路分析，搜尋與目標性狀相關的通路。利用在線富集分析可視化網站(http://www.bioinformatics.com.cn)對顯著富集的基因繪制柱狀圖和氣泡圖。

2 結 果

2.1 表型描述統計

本研究統計269頭蘇山豬的7種表型的平均值、標準差、最大值、最小值和變異系數，結果如表1所示。結果顯示，體重、體長、胸圍、管圍、左乳頭數、右乳頭數和總乳頭數性狀的變異系數為26.87%、8.97%、10.31%、9.54%、12.28%、12.69%和11.63%，其中體重的變異系數最大，高達26.87%，表明體重數據中存在較大的變異。基于標準差可以看出，右乳頭數的變異程度大于左乳頭數。

表1 蘇山豬初生體尺和乳頭數性狀表型值的統計結果Table 1 Statistical results of phenotypic values of body conformation at birth and teat number traits in Sushan pigs

2.2 基因型填充準確度

經過利用Beagle軟件對269頭蘇山豬的55K芯片數據進行基因型填充，共獲得31 123 761個SNPs，剔除MAF<0.01的SNPs后，得到8 547 616個SNPs。利用CR和correlation估算每條染色體的基因型填充準確度(圖1)，結果顯示，平均CR為0.71，最大值和最小值分別為0.83(SSC X)和0.67(SSC 5)；平均correlation為0.58，最大值和最小值分別為0.75(SSC X)和0.52(SSC 5)，基因型填充后，顯著增加了55K芯片的密度(圖2)。通過call rate<0.1、MAF<0.01和correlation<0.8過濾掉不符合分析要求的SNPs后，最終有2 676 280個SNPs用于后續數據分析，其平均相關性為0.93。

橫坐標以蘇山豬染色體繪制，縱坐標表示填充準確度，correct rate和correlation為過濾前的等位基因準確性和相關性，filter為過濾后的相關性The x-axis represents the chromosomes of Sushan pigs, and the y-axis represents the imputation accuracy, correct rate and correlation represent CR and correlation before filtering, filter is the correlation after filtering圖1 蘇山豬染色體填充準確度Fig.1 Imputation accuracy across whole chromosomes in Sushan pig population

上圖和下圖分別為染色體填充前和填充后的SNP密度分布，橫坐標表示染色體大小，縱坐標表示染色體，刻度顏色表示1 Mb窗口內的SNP數目Above and below are figures respectively represent the SNP density distribution before and after imputation, the x-axis denotes the size of chromosomes, the y-axis denotes the chromosomes, color scale denotes the number of SNP within 1 Mb window size圖2 蘇山豬染色體填充前后SNP密度分布比對Fig.2 Comparison of SNP density distribution before and after imputation of chromosomes in Sushan pig population

2.3 選擇信號分析和全基因組關聯分析

按照各性狀Fst值的前1%為標準，可得體重、體長、胸圍、管圍、左乳頭數、右乳頭數和總乳頭數性狀的Fst閾值分別為0.16、0.20、0.20、0.17、0.14、0.14和0.15，與Wright[13]定義的高度分化的閾值范圍接近(0.15

初生體尺和乳頭數性狀的群體分化指數和全基因組關聯分析曼哈頓圖進行組合。群體分化指數曼哈頓圖中虛線為Top1%的閾值線，橫坐標以蘇山豬染色體繪制，縱坐標表示群體分化指數；全基因組關聯分析曼哈頓圖中橫坐標以染色體繪制，縱坐標表示 -lg(P)，虛線表示顯著閾值線Population differentiation index and genome-wide association analysis Manhattan plots were combined for body conformation at birth and teat number traits. The dotted lines indicate the threshold line of top1% in the population differentiation Manhattan plots, the x-axis represents the chromosomes of Sushan pigs, and the y-axis represents the Fst value; genome-wide association analysis Manhattan plots are plotted with chromosomes on the abscissa, the y-axis represents the -lg(P), and the dashed lines represents significance thresholds圖3 蘇山豬初生體尺和乳頭數性狀群體分化指數和全基因組關聯分析曼哈頓圖Fig.3 Manhattan plots of fixation index and genome-wide association analysis for body conformation at birth and teat number traits in Sushan pig population

蘇山豬群體初生體尺和乳頭數性狀的GWAS分析結果顯示，與體尺性狀顯著關聯的SNP有32個，乳頭數性狀顯著關聯的SNP有90個。其中，與體重顯著關聯的SNP有1個，位于X染色體；與體長顯著關聯的位點有13個，全部位于14號染色體；與胸圍顯著關聯的位點有9個，分別位于14和X染色體；與管圍顯著關聯的位點有9個，位于X染色體；與左乳頭數顯著關聯的位點有6個，分別位于1、13和16號染色體；與右乳頭數顯著關聯位點有53個，分別位于1和13號染色體；與總乳頭數顯著關聯的位點有31個，分別位于13、16和X染色體。QQ-plot分析結果見圖4，圖中膨脹系數(λ)介于0.884～1.062之間，表明本研究試驗群體中不存在明顯的種群分層現象。

λ為基因組膨脹系數λ represents coefficient of expansion圖4 蘇山豬初生體尺和乳頭數性狀分位數-分位數圖Fig.4 The quantitle-quantitle plots for body conformation at birth and teat number traits in Sushan pig population

通過對蘇山豬群體初生體尺和乳頭數性狀進行Fst和GWAS分析，篩選出與體長、胸圍、管圍和總乳頭數性狀顯著關聯的4個候選位點，定位在13、14和X染色體上(表2，圖3)。與體長性狀顯著相關的候選位點位于14號染色體，ACTA2基因為位置功能候選基因。與管圍和胸圍性狀都顯著相關的候選位點位于X染色體上，COL4A5和COL4A6基因為位置功能候選基因?？側轭^數性狀顯著相關的候選位點分別在13和X染色體上，ATP1B4、UBE2A、NKRF、SEPTIN6、NKAP、ZBTB33、CRTAP、TRIM71、FBXL2基因為位置功能候選基因。

表2 蘇山豬初生體尺和乳頭數性狀顯著相關SNPsTable 2 The SNPs significantly associated with body conformation at birth and teat number traits in Sushan pigs

2.4 富集分析

本研究基于Metascape在線網站中的小鼠數據庫開展了 GO和KEGG功能富集分析。共富集到了46條顯著GO條目(圖5)，其中參與生物過程的GO條目為34條，參與形成細胞成分的GO 條目為11條，能解釋分子功能的GO條目為1條。圖6顯示了基因富集GO 條目的顯著性。候選基因的KEGG通路富集分析中發現4條通路與蛋白質消化、松弛素信號通路、溶酶體生物發生以及細胞自噬功能有關(表3)。

橫坐標表示GO條目；縱坐標表示GO條目富集的基因數The x-axis represents GO term; the y-axis represents the number of genes enriched by GO term圖5 候選基因GO富集分析柱狀圖Fig.5 Histogram of GO enrichment analysis of the candidate genes

橫坐標表示rich factor，rich factor 為位于該通路的差異表達基因個數/位于該通路的所有注釋基因個數；縱坐標表示GO 條目The x-axis represents the rich factor, rich factor indicate the number of differentially expressed genes in this pathway/the number of all annotated genes in this pathway; the y-axis represents the GO terms圖6 候選基因GO富集分析氣泡圖Fig.6 Bubble plot of GO enrichment analysis of the candidate genes

表3 蘇山豬初生體尺和乳頭數性狀KEGG分析顯著富集Table 3 Significant kyoto encyclopedia of genes and genomes (KEGG) pathway associated with body conformation at birth and teat number traits in Sushan pig population

3 討 論

本研究利用387頭豬的重測序數據對269頭豬已有的SNP芯片數據進行填充，得到更高質量和更高密度的基因數據集，以CR和相關性作為基因型填充準確性的評估指標[14]。據研究表明，影響填充準確性的因素主要包括參考群體大小[15]、群體結構[16]、最小等位基因頻率[17]、填充群體之間的遺傳距離[18]和填充方法[19]。Ye等[16]在對雞的研究中使用60K數據對基因組序列進行填充，填充的準確性為0.62。Van Binsbergen等[20]在對奶牛的研究中使用基因型填充手段，填充的準確性范圍分別為0.77～0.83、0.37和0.46。本研究與先前兩個研究的填充準確性較為接近，表明本研究填充后的數據可用于后續分析。

先前對豬乳頭數的數量基因座(QTL)研究過程中，除了3、4、5、18和19號染色體上沒有發現與乳頭數相關的QTL外，其余的染色體上均有QTL存在，有趣的是，在豬QTL庫中發現和乳頭數相關的QTL與胸椎數的QTN存在一些重疊[21-23]。Duijvesteijn等[24]發現，乳頭數的QTL包含了與椎體發育相關的基因，從而解釋了乳頭數和脊椎數之間存在相關的生物機制[25]，乳頭數和脊椎數呈正相關，并且豬體長性狀與胸椎數和脊椎數密切相關，體長性狀相關聯的基因很可能會影響豬乳頭數[26-27]。rs331034247位于ACTA2基因附近，與體長性狀顯著關聯。根據先前報道，Weigand等[28]利用正常女性和患有乳腺癌女性的乳腺細胞，分別進行培養并開展一系列表達分析，發現ACTA2與脂肪干細胞(adipose derived stem cells,ADSC)和乳腺上皮細胞(mammary epithelial cells,MES)之間存在正相關[29]。黃京書和熊遠著[30]利用RT-PCR探究了ACTA2基因在多個品種豬的骨骼肌表達分析，結果顯示處于不同發育階段的骨骼肌中基因的表達量存在較大差異，ACTA2基因很可能與骨骼肌發育存在一定的相關性。由上可知，ACTA2基因很可能是體長和乳頭數的候選基因。

rs1108457554位于CRTAP基因的蛋白編碼區域以及FBXL2基因的下游，與TTN性狀顯著相關。Tang等[31]、Grol等[32]和Tonelli等[33]在人類、小鼠和斑馬魚的研究中發現，CRTAP基因中存在導致隱性成骨發育不全癥的致病雜合致病突變，由此推斷CRTAP基因也可能會導致豬椎體發育異常，從而影響乳頭數量。FBXL2基因在調節骨骼成肌細胞增殖和成肌分化中發揮著重要作用，該基因的沉默會導致骨骼肌中的成肌細胞增殖[34-35]，骨骼肌的發育異常很可能會引起動物椎骨畸形和體尺特征較差，進而導致乳頭數減少和死亡率增加。在rs1108457554的下游搜尋到TRIM71 基因，Connacher和Goldstrohm[36]發現TRIM71 基因促進干細胞增殖，誘導成纖維細胞重編程為多功能干細胞，也有研究表明[37]該基因與骨骼肌的發育有關。

rs332938314與TTN性狀顯著相關，位于SEPTIN6、UBE2A、NKRF和ATP1B4基因的附近。Li等[38]研究發現，SEPTIN6 基因是影響氨基酸介導的奶牛乳腺上皮細胞生長和酪蛋白合成的關鍵積極調控因子，SEPTIN6 基因參與細胞的生理過程，包括胞質分裂、囊泡運輸、細胞形態、細胞運動和細胞周期[39]。Fingar等[40]研究發現，氨基酸可以影響細胞增殖和牛奶的生成，由此推斷該基因很可能與豬的乳頭數以及泌乳能力存在一定關聯。在對人的研究中發現，UBE2A和NKRF基因發生突變或者缺失會導致乳頭間距增大、智力障礙、全身多毛癥和泌尿生殖系統發育異常等癥狀[41-43]，這兩個基因突變和缺失是否會導致豬出現相同的癥狀，還需要開展相關研究。Pestov等[44]在對小鼠的研究中發現，ATP1B4基因表達量會對表皮的生長產生影響，在其他研究中發現ATP1B4基因編碼的蛋白質βm主要在骨骼肌中高水平表達[45-47]。

rs324391601位于COL4A5基因內和COL4A6基因上游附近，與ChC顯著關聯。Jiang 等[48]在對湖羊的研究中發現，COL4A6基因是ChC的候選基因，本研究中對蘇山豬胸圍性狀候選基因搜尋中也得出相同的結果，因此COL4A6基因可以作為ChC的候選基因。在對日本和牛的研究中發現，COL4A5基因在肌內脂肪組織中高特異性表達，與日本和牛的肌內脂肪形成有關，可以將該基因作為豬肌內脂肪研究的相關基因。

富集分析結果表明，ACTA2、CRTAP、TRIM71、COL4A5和COL4A6基因主要顯著富集在跨膜受體蛋白酪氨酸激酶信號通路，該通路可以影響乳腺癌細胞增殖和凋亡[49]。ACTA2、COL4A5和COL4A6基因同時顯著富集在松弛素信號通路，通路異?？梢杂绊懶∈蟮娜轭^嚴重發育不良，造成無法哺乳后代[50]。結果表明，這些候選基因可能在豬胚胎發育過程中對體形態和乳頭數性狀起重要作用。

4 結 論

本研究基于SNP數據集進行了Fst和GWAS分析，對顯著性SNPs進行注釋和富集分析之后，發現了1個與體長性狀相關的候選基因ACTA2；發現了1個與胸圍性狀相關的候選基因COL4A6；發現了3個與乳頭數性狀相關的候選基因ACTA2、CRTAP和SEPTIN6，為蘇山豬體尺性狀和乳頭數性狀的遺傳改良提供了新的視角。