榮昌豬初產繁殖性狀的全基因組關聯研究

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

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

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

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

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

1 材料與方法

1.1 試驗動物與繁殖性狀

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

1.2 基因組DNA提取與質檢

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

1.3 基因分型與基因型填充

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

1.4 群體結構分析

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

1.5 全基因組關聯分析

首先，采用R軟件GenABEL包[16]將繁殖數據校正至標準正態分布。然后，采用GEMMA軟件[17]的單變量線性混合模型，將出生年、出生月作為固定因子，主成分分析的前3個主成分作為協變量，進行榮昌豬繁殖性狀的關聯分析。單變量線性混合模型如下：

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

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

1.6 候選基因篩選

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

2 結 果

2.1 表型數據分析

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

表1 榮昌豬繁殖性狀表型數據統計分析Table 1 Descriptive statistics of reproductive traits in Rongchang pig population

2.2 質量控制

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

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

2.2 群體結構分析

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

橫坐標和縱坐標都代表樣本個體編號；圖中每個小方格代表其所連接的兩個個體之間的IBS遺傳距離，遺傳距離越遠越接近紅色，反之亦然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 親緣關系矩陣熱圖Fig.2 Kinship matrix heatmap

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

圖3 PCA群體結構圖Fig.3 Principal components analysis of population structure

2.3 全基因組關聯分析結果

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

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

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

初生窩重關聯分析結果顯示，在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位點與多個性狀存在顯著關聯，其中SNP位點SSC1: 279 214 647 bp為總產仔數和活產仔數的共享位點；SNP位點SSC17: 57 315 180 bp同時與總產仔數、活產仔數、初生窩重存在潛在關聯(表2)，暗示該位點存在一因多效性，可能同時影響榮昌豬總產仔數、活產仔數和初生窩重等繁殖性狀。

3 討 論

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

3.1 群體分層與胎次效應

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

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

3.2 候選基因功能分析

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

左圖為繁殖性狀全基因組關聯分析的曼哈頓圖，圖中紅色實線代表全基因組水平顯著閾值(1.43×10-6)，藍色實線代表全基因組水平潛在顯著閾值(2.85×10-5)；右圖為繁殖性狀全基因組關聯分析的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 繁殖性狀全基因組關聯分析的曼哈頓圖和Q-Q圖Fig.4 Manhattan and Q-Q plots of genome-wide association study for reproductive traits in Rongchang pig population

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

4 結 論

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