榮昌豬初產(chǎn)繁殖性狀的全基因組關(guān)聯(lián)研究

吳平先，陳 力*，龍 熙，柴 捷，張廷煥，徐順來，郭宗義，王金勇

(1.重慶市畜牧科學(xué)院，榮昌 402460； 2.國家生豬技術(shù)創(chuàng)新中心，榮昌 402460)

榮昌豬是我國著名的地方良種豬，享有“華夏之寶”美譽，具有毛色獨特、肉質(zhì)優(yōu)良、耐粗飼、雜交配合力高及繁殖性能好等突出優(yōu)點。繁殖性狀是養(yǎng)豬生產(chǎn)中重要的經(jīng)濟性狀，榮昌豬初產(chǎn)仔數(shù)為10頭左右，經(jīng)產(chǎn)仔數(shù)在12頭左右。近年來，榮昌豬的各項生產(chǎn)性能缺乏持續(xù)的選育，在產(chǎn)仔數(shù)、瘦肉率、料肉比等生產(chǎn)性能指標(biāo)上與外種豬差距逐漸拉大。繁殖力是由微效多基因控制的數(shù)量性狀，遺傳力低，通過傳統(tǒng)的選育手段，遺傳評估準(zhǔn)確性低，遺傳進展十分緩慢[1]。因此，迫切需要通過聯(lián)合常規(guī)育種與分子育種技術(shù)開展科學(xué)選育，提升榮昌豬遺傳改良速度，振興榮昌豬種業(yè)。榮昌豬繁殖性狀的遺傳解析能為分子育種提供可靠的基因資源和分子標(biāo)記，有助于提升分子育種的準(zhǔn)確性和育種進程。但是，目前有關(guān)榮昌豬繁殖性狀特異性功能基因和位點的研究鮮有報道。

豬繁殖性狀功能基因的鑒定一直是豬育種領(lǐng)域的研究熱點，也是一個極具挑戰(zhàn)性的研究難點。隨著全基因組基因分型技術(shù)和功能基因鑒定方法的發(fā)展，全基因組關(guān)聯(lián)分析已經(jīng)成為鑒定動物疾病和數(shù)量性狀遺傳變異的重要手段。近年來，國內(nèi)外眾多研究已經(jīng)鑒定出大量影響豬繁殖性狀的重要基因和分子標(biāo)記。Rothschild等[2]最早發(fā)現(xiàn)ESR基因是影響豬產(chǎn)仔數(shù)的因果基因，BB型母豬每窩的產(chǎn)仔數(shù)比AA型母豬多1.5～2.3頭，并且發(fā)現(xiàn)ESR在不同胎次中對產(chǎn)仔數(shù)的影響明顯不同。大量研究報道證實，AHR基因與母豬繁殖性狀有緊密聯(lián)系[3-5]。Bosse等[6]通過基因組分析發(fā)現(xiàn)，AHR基因的非同義突變能顯著增加母豬的產(chǎn)仔數(shù)。研究表明，F(xiàn)SHβ基因?qū)Σ煌贩N母豬繁殖性能的影響差異較大，大白豬BB型母豬比其他基因型母豬每窩多產(chǎn)0.41～1.49頭仔豬[7]；長白豬AB型產(chǎn)仔數(shù)高于其他基因型；中國地方豬AA基因型母豬的產(chǎn)仔數(shù)高于BB基因型母豬[8]。此外，對外種豬繁殖性狀的研究已鑒定出大量與繁殖性狀緊密相關(guān)的變異位點[9]，如RFRP基因g.45859759 C>T位點的突變顯著影響母豬產(chǎn)仔數(shù)[10]，NR4A1基因g3952 A>G突變與總產(chǎn)仔數(shù)、活產(chǎn)仔數(shù)和窩重顯著相關(guān)[11]，GNB2L1基因g2373 T>C突變顯著影響仔豬的初生重[11]。根據(jù)豬PigQTL數(shù)據(jù)庫統(tǒng)計顯示，已經(jīng)鑒定到384個與總產(chǎn)仔數(shù)相關(guān)的QTLs；228個與活產(chǎn)仔數(shù)相關(guān)的QTLs；85個與窩重相關(guān)的QTLs；138個與死產(chǎn)仔數(shù)相關(guān)的QTLs。

然而，絕大部分重要基因和位點是通過外種豬群體篩選到的，對榮昌豬產(chǎn)仔性狀的特異性基因和位點的研究較少，在一定程度上限制了榮昌豬繁殖性狀分子選育的實施。本研究基于純種榮昌豬群體，統(tǒng)計總產(chǎn)仔數(shù)、活產(chǎn)仔數(shù)、死產(chǎn)仔數(shù)以及初生窩重等繁殖表型數(shù)據(jù)，通過豬50K基因芯片對群體進行基因分型，利用GWAS方法鑒定繁殖性狀相關(guān)的候選基因和位點，為榮昌豬本品種定向化持續(xù)選育以及新品種培育提供重要的分子標(biāo)記，為榮昌豬種源“卡脖子”技術(shù)攻關(guān)奠定理論基礎(chǔ)。

1 材料與方法

1.1 試驗動物與繁殖性狀

本試驗所用的429頭純種榮昌豬均來源于國家生豬核心育種場-重慶琪泰佳牧畜禽養(yǎng)殖有限公司。所有豬只均為初產(chǎn)母豬，在相同的飼養(yǎng)環(huán)境、營養(yǎng)水平和飼養(yǎng)管理等條件下飼養(yǎng)。收集2021年1月至2022年6月共429條初產(chǎn)繁殖記錄，包括總產(chǎn)仔數(shù)(total number born, TNB)、活產(chǎn)仔數(shù)(total number born alive, NBA)、死胎數(shù)(number of stillborn, NS)和初生窩重(litter weight born alive, LWB)。采集所有母豬的耳組織樣品于離心管中，加入75%乙醇，-20 ℃保存。

1.2 基因組DNA提取與質(zhì)檢

利用氯仿抽提法提取豬耳組織樣的基因組DNA，通過紫外分光光度計(Nanodrop 2000)和1%的瓊脂糖凝膠電泳進行DNA質(zhì)量檢測，當(dāng)濃度>100 ng·μL-1OD260 nm/OD280 nm介于1.8～2.0之間，條帶清晰，無拖尾降解，則DNA質(zhì)檢合格。

1.3 基因分型與基因型填充

質(zhì)檢合格的基因組DNA樣本利用“中芯一號”芯片(北京康普森農(nóng)業(yè)科技有限公司，北京)進行SNP基因分型，獲得的基因分型數(shù)據(jù)使用Plink(V1.90)軟件[12]進行質(zhì)量控制，剔除無染色體位置信息以及Y染色體上的位點、最小等位基因頻率(minor allele frequency，MAF)小于0.01的位點、個體的SNP檢出率小于90%的位點、哈迪-溫伯格平衡檢驗P<10-6的位點?；谫|(zhì)控后的基因型數(shù)據(jù)，采用Beagle(V5.0)軟件[13]對有缺失的基因型數(shù)據(jù)進行填充，填充后的數(shù)據(jù)用于后續(xù)分析。

1.4 群體結(jié)構(gòu)分析

群體分層會造成全基因組關(guān)聯(lián)分析出現(xiàn)假陽性結(jié)果[14]。因此采用Plink和GCTA軟件[15]對填充后的基因型數(shù)據(jù)進行遺傳距離分析和主成分分析，計算基于狀態(tài)同源(idengtical by state, IBS)的遺傳距離，構(gòu)建基因組親緣關(guān)系矩陣，以研究榮昌豬群體的親緣關(guān)系和群體分布。采用R軟件的ggplot2包繪制熱圖，對解釋方差最大的前3個主成分進行散點圖繪制。

1.5 全基因組關(guān)聯(lián)分析

首先，采用R軟件GenABEL包[16]將繁殖數(shù)據(jù)校正至標(biāo)準(zhǔn)正態(tài)分布。然后，采用GEMMA軟件[17]的單變量線性混合模型，將出生年、出生月作為固定因子，主成分分析的前3個主成分作為協(xié)變量，進行榮昌豬繁殖性狀的關(guān)聯(lián)分析。單變量線性混合模型如下：

y=γCov+Xβ+Zα+Wμ+e

多重檢驗后，采用Bonferroni校正法來設(shè)定顯著閾值，基因組水平顯著閾值設(shè)置為0.05/N=0.05/35 046=1.43×10-6；基因組潛在顯著閾值設(shè)置為1/N=1/35 046=2.85×10-5，其中N表示參與分析的SNPs數(shù)量。使用R軟件的qqman包[18]繪制GWAS結(jié)果的曼哈頓圖和Q-Q圖。

1.6 候選基因篩選

本研究以顯著位點上、下游0.5 Mb范圍內(nèi)的基因組區(qū)域作為QTLs區(qū)域，通過Ensembl網(wǎng)站(http://ensembl.org/Sus_scrofa/Info/Index)檢索QTLs區(qū)域內(nèi)的候選基因，通過NCBI(https://www.ncbi.nlm.nih.gov/)和Genecards在線網(wǎng)站(https://www.genecards.org/)以及已報道文獻(xiàn)進一步查找基因的功能，利用Animal QTL Database下載影響豬產(chǎn)仔數(shù)、死胎數(shù)和初生窩重的QTLs信息，與本研究鑒定出的QTLs進行比對，綜合分析所有候選基因的生物學(xué)功能和QTLs信息，初步篩選出影響榮昌豬繁殖性狀的重要候選基因。

2 結(jié) 果

2.1 表型數(shù)據(jù)分析

表1展示了榮昌豬初產(chǎn)繁殖性狀的平均值、標(biāo)準(zhǔn)差和變異系數(shù)等信息。表1結(jié)果顯示，總產(chǎn)仔數(shù)的平均值和變異系數(shù)分別是10.18頭和25.57%；活產(chǎn)仔數(shù)的平均值和變異系數(shù)分別是8.89頭和28.71%；死胎數(shù)的平均值和變異系數(shù)分別是0.38頭和264.49%；初生窩重的平均值和變異系數(shù)分別是8.20 kg和28.77%。該結(jié)果表明榮昌豬繁殖性狀的變異系數(shù)較高，具有較大的遺傳改良潛力。

表1 榮昌豬繁殖性狀表型數(shù)據(jù)統(tǒng)計分析Table 1 Descriptive statistics of reproductive traits in Rongchang pig population

2.2 質(zhì)量控制

利用中芯一號50K SNP芯片對429頭純種榮昌母豬進行基因分型檢測，質(zhì)量控制后，得到35 046個有效的SNPs，其中1號染色體上的位點數(shù)最多，有3 815個；18號染色體上的位點數(shù)最少，有879個。質(zhì)控后的SNPs分布密度圖見圖1，同一條染色體上相鄰SNPs間的平均距離最大為0.089 Mb，位于11號染色體上；最小為0.055 Mb，位于X染色體上。

圖1 SNP分布密度圖Fig.1 SNP distribution density map

2.2 群體結(jié)構(gòu)分析

圖2為429頭榮昌母豬基于IBS的親緣關(guān)系矩陣熱圖，親緣關(guān)系矩陣熱圖的顏色越紅表示個體間遺傳距離越遠(yuǎn)。結(jié)果顯示，榮昌母豬群體的IBS距離值為0.000 1～0.343，平均遺傳距離為0.242 7，絕大部分個體的遺傳距離都大于0.1，表明榮昌豬群體個體間的平均遺傳距離較遠(yuǎn)。

橫坐標(biāo)和縱坐標(biāo)都代表樣本個體編號；圖中每個小方格代表其所連接的兩個個體之間的IBS遺傳距離，遺傳距離越遠(yuǎn)越接近紅色，反之亦然The x-axis and y-axis indicate the sample number. Each small square shows the IBS genetic distance value between tow connected individuals, the farther genetic distance, the color is closer to red, and vice versa圖2 親緣關(guān)系矩陣熱圖Fig.2 Kinship matrix heatmap

圖3是榮昌豬群體結(jié)構(gòu)分布情況，從圖中可以看出試驗群體分成了3個部分，說明該群體在人工選育過程中引起群體結(jié)構(gòu)差異，但是各部分之間存在一定的聯(lián)系，說明各部分間存在基因交流。因此，在后續(xù)的全基因組關(guān)聯(lián)分析中將群體分層情況作為協(xié)變量，對群體結(jié)構(gòu)進行校正，以降低群體分層帶來的影響。

圖3 PCA群體結(jié)構(gòu)圖Fig.3 Principal components analysis of population structure

2.3 全基因組關(guān)聯(lián)分析結(jié)果

本研究采用單標(biāo)記關(guān)聯(lián)分析對榮昌豬初產(chǎn)的總產(chǎn)仔數(shù)、活產(chǎn)仔數(shù)、死胎數(shù)和初生窩重進行全基因組關(guān)聯(lián)分析，獲得與目標(biāo)性狀顯著關(guān)聯(lián)的位點和候選基因，結(jié)果見圖4和表2?？偖a(chǎn)仔數(shù)的關(guān)聯(lián)分析結(jié)果顯示，共檢測到5個潛在關(guān)聯(lián)的SNPs(P<2.85×10-5)，分別位于1、8、9、14、17號染色體上，其中最顯著的SNPs位于SSC1: 279 214 647 bp(P=5.94×10-6)，但是在其上、下游500 kb范圍內(nèi)沒有發(fā)現(xiàn)候選基因。

活產(chǎn)仔數(shù)關(guān)聯(lián)分析結(jié)果顯示，共檢測出3個潛在關(guān)聯(lián)的SNPs(P<2.85×10-5)，分別位于1、5、17號染色體上，其中最顯著的SNPs SSC17: 57 315 180 bp(P= 6.02×10-6)位于BMP7、TFAP2C基因附近。

死胎數(shù)關(guān)聯(lián)分析結(jié)果顯示，在1、2、12、13、15號染色體上共檢測到11個顯著的SNPs(P<2.85×10-5)，其中有1個SNPs(SSC2: 89 287 506 bp，P=8.72×10-8)達(dá)基因組顯著水平，并且在其上、下游500 kb內(nèi)檢測到一個候選基因MSH3。在13號染色體上檢測到3個連續(xù)的顯著SNPs，位于CBLB基因內(nèi)。

初生窩重關(guān)聯(lián)分析結(jié)果顯示，在5和17號染色體上共檢測到2個基因組水平顯著的SNPs(P<1.43×10-6)。其中最顯著的SNPs(SSC17: 57 315 180 bp，P=5.61×10-8)位于BMP7、TFAP2C、RBM38基因附近。

2.4 顯著的共享位點

本研究中鑒定到2個共享SNP位點與多個性狀存在顯著關(guān)聯(lián)，其中SNP位點SSC1: 279 214 647 bp為總產(chǎn)仔數(shù)和活產(chǎn)仔數(shù)的共享位點；SNP位點SSC17: 57 315 180 bp同時與總產(chǎn)仔數(shù)、活產(chǎn)仔數(shù)、初生窩重存在潛在關(guān)聯(lián)(表2)，暗示該位點存在一因多效性，可能同時影響榮昌豬總產(chǎn)仔數(shù)、活產(chǎn)仔數(shù)和初生窩重等繁殖性狀。

3 討 論

總產(chǎn)仔數(shù)、活產(chǎn)仔數(shù)和初生窩重是衡量母豬繁殖力的重要指標(biāo)，死胎數(shù)是直接衡量母豬繁殖力損失的重要特征。榮昌豬作為地方豬品種之一，由于缺乏長期的持續(xù)選育，繁殖性能的改良仍在原地踏步。為了提高榮昌豬的生產(chǎn)效率，育種研究者正在積極探索針對繁殖性能的育種策略。隨著分子育種技術(shù)的發(fā)展，利用基因組選擇等分子技術(shù)有望快速提高榮昌豬繁殖性能的遺傳改良速度。目前，與榮昌豬繁殖性能相關(guān)的特征性遺傳變異位點和基因的挖掘鑒定仍然十分缺乏，迫切需要加快繁殖性狀特征遺傳變異位點的鑒定與利用。全基因組關(guān)聯(lián)研究是解析表型變異的強有力方法，為數(shù)量性狀的遺傳解析提供了許多關(guān)鍵的遺傳信息，在畜禽育種中做出了重大貢獻(xiàn)。但是，榮昌豬繁殖性狀的全基因組關(guān)聯(lián)研究鮮有研究報道。因此，本研究基于50K芯片數(shù)據(jù)對榮昌豬第一胎次的總產(chǎn)仔數(shù)、活產(chǎn)仔數(shù)、初生窩重和死胎數(shù)進行全基因組關(guān)聯(lián)研究，以期鑒定出影響榮昌豬繁殖性狀的重要變異位點和基因，為榮昌豬多個繁殖性狀的基因組協(xié)同選擇提供基礎(chǔ)數(shù)據(jù)。

3.1 群體分層與胎次效應(yīng)

群體分層是直接導(dǎo)致假陽性關(guān)聯(lián)分析結(jié)果的主要因素，也是衡量關(guān)聯(lián)分析結(jié)果可靠性的決定因素。本研究中，主成分分析結(jié)果說明榮昌豬群體存在一定的群體分離情況。為了減少因群體分層帶來的假陽性結(jié)果，本研究在線性混合模型中將群體主成分作為協(xié)變量，出生年、出生月作為固定因子，以消除群體分層造成的影響。理想模型中，當(dāng)基因組膨脹系數(shù)為1.00時，說明群體分層對關(guān)聯(lián)分析沒有影響。在實際的研究中，當(dāng)基因組膨脹系數(shù)小于1.05時，說明不存在群體分層的影響。本研究中Q-Q圖分析結(jié)果說明該關(guān)聯(lián)分析很好的控制了群體分層帶來的影響。

豬的繁殖性狀是數(shù)量性狀，受微效多基因的綜合調(diào)控，并且因母豬個體生理、飼養(yǎng)環(huán)境等因素的影響，不同胎次間同一基因的表達(dá)模式不同，直接造成不同胎次間產(chǎn)仔性能的差異[19]。不同胎次的繁殖性狀之間遺傳相關(guān)較低[20]，說明不同胎次的遺傳機制存在明顯差異，母豬胎次效應(yīng)可能是影響繁殖性能的重要因素，不同胎次可定義為不同的性狀進行分析[19]。為了避免胎次效應(yīng)帶來的假陽性或假陰性結(jié)果，提升繁殖性狀候選基因鑒定的準(zhǔn)確性，本研究僅對榮昌豬第一胎次的繁殖性狀進行了關(guān)聯(lián)分析，挖掘出一系列影響榮昌豬初產(chǎn)繁殖性能的遺傳變異。

3.2 候選基因功能分析

本研究發(fā)現(xiàn)17號染色體上的位點CNC10171044同時與總產(chǎn)仔數(shù)、活產(chǎn)仔數(shù)和初生窩重顯著相關(guān)，并且在該位點附近篩選出BMP7基因。BMP7基因是骨形態(tài)發(fā)生蛋白，屬于BMP亞家族，通過調(diào)節(jié)激素生成、顆粒細(xì)胞狀態(tài)和卵泡發(fā)育來影響雌性動物的生殖力，并且在生長發(fā)育過程中發(fā)揮著關(guān)鍵作用[21]。有研究表明，母羊生殖器官中BMP7基因表達(dá)豐度與繁殖力呈正相關(guān)，繁殖力高的個體卵巢中BMP7基因表達(dá)豐度比繁殖力低的個體高[22]，并且發(fā)現(xiàn)BMP7基因g.58171886位點變異直接影響母羊第一胎次的產(chǎn)仔數(shù)[23]。Li等[24]基于全基因組重測序，通過選擇信號分析等方法鑒定了高、低繁殖力大白母豬基因組差異，結(jié)果發(fā)現(xiàn)，BMP7是影響母豬排卵率和窩產(chǎn)仔數(shù)性狀的潛在候選基因。此外，有研究表明，BMP7基因的c.2256 G>C、c.1569 A>G、g.35161 T>C都與母豬總產(chǎn)仔數(shù)、死胎數(shù)以及初生窩重顯著相關(guān)[25-26]，是影響母豬繁殖力的重要候選基因。BMP7基因可能通過調(diào)節(jié)母豬激素生成和卵泡發(fā)育進而影響母豬的繁殖力。因此，該基因可能是影響榮昌母豬產(chǎn)仔數(shù)、活產(chǎn)仔數(shù)和初生窩重的重要候選基因。

左圖為繁殖性狀全基因組關(guān)聯(lián)分析的曼哈頓圖，圖中紅色實線代表全基因組水平顯著閾值(1.43×10-6)，藍(lán)色實線代表全基因組水平潛在顯著閾值(2.85×10-5)；右圖為繁殖性狀全基因組關(guān)聯(lián)分析的Q-Q圖The figures at left represent Manhattan plots of genome-wide association study for reproductive traits, the red line in each figure shows genome-wide significant threshold: 1.43×10-6, and the blue line shows the suggestive significant threshold: 2.85×10-5. The figures at right represent quantile-quantile (Q-Q) plots of genome-wide association study for reproductive traits圖4 繁殖性狀全基因組關(guān)聯(lián)分析的曼哈頓圖和Q-Q圖Fig.4 Manhattan and Q-Q plots of genome-wide association study for reproductive traits in Rongchang pig population

與榮昌母豬死胎數(shù)顯著相關(guān)的位點位于2號染色體89 287 506 bp，其鄰近基因為MSH3基因。通過與豬QTL數(shù)據(jù)庫比對發(fā)現(xiàn)，該位點接近于已報道的大白豬死胎數(shù)相關(guān)的QTLs(87.4～87.4 Mb)[27]。MSH3基因位于染色體5q11～q13，編碼1 137個氨基酸殘基的蛋白質(zhì)，在DNA錯配識別、修復(fù)以及維持基因組穩(wěn)定中發(fā)揮著關(guān)鍵作用。MSH3基因表達(dá)降低和缺失都會引起生殖相關(guān)疾病的發(fā)生。有研究表明，MSH3基因的表達(dá)與精子活力顯著相關(guān)，并且與卵巢癌、乳腺癌等繁殖相關(guān)的疾病直接相關(guān)，進而影響生殖力[28-29]。有研究發(fā)現(xiàn)，DNA損傷修復(fù)基因表達(dá)在胚胎發(fā)育過程中發(fā)揮著關(guān)鍵作用[30]，暗示MSH3基因可能通過調(diào)節(jié)DNA錯配識別和修復(fù)進而影響胚胎的發(fā)育。母體妊娠期疾病(如：腎病、糖尿病以及心血管疾病等)是妊娠期胎兒死亡的重要原因，有研究發(fā)現(xiàn)，與母體腎病、糖尿病以及心血管疾病相關(guān)的基因可能直接引發(fā)妊娠期胎兒致死[31]。本研究在13號染色體上檢測到3個連續(xù)的SNPs影響榮昌母豬死胎數(shù)，并且在其附近發(fā)現(xiàn)CBLB基因。CBLB基因與多種疾病相關(guān)，有研究表明，CBLB基因突變是先天免疫性疾病和Ⅰ型糖尿病發(fā)生的重要原因[32]。CBLB基因SNP位點的錯義突變可能通過改變蛋白質(zhì)的結(jié)構(gòu)進而影響牛炎癥相關(guān)疾病的發(fā)生[33]。綜上，MSH3和CBLB基因可能通過誘發(fā)妊娠期母體疾病的發(fā)生，進而影響胎兒的正常生長發(fā)育。

4 結(jié) 論

本研究基于50K基因芯片對榮昌豬總產(chǎn)仔數(shù)、活產(chǎn)仔數(shù)和初生窩重進行全基因組關(guān)聯(lián)分析，挖掘影響榮昌豬繁殖性狀的重要遺傳變異信息。BMP7基因可能是影響榮昌豬總產(chǎn)仔數(shù)、活產(chǎn)仔數(shù)和初生窩重的重要功能位點；MSH3和CBLB基因可能是影響榮昌豬死胎數(shù)性狀的重要候選基因。研究結(jié)果為榮昌豬繁殖性狀提供了重要的遺傳變異位點和候選基因，也為榮昌豬繁殖性狀的基因組選擇提供了重要的理論基礎(chǔ)。