李望嬌，彭 夏，宋 徽，董文君，李新云，趙書紅，馬云龍

(華中農(nóng)業(yè)大學(xué)，農(nóng)業(yè)動(dòng)物遺傳育種與繁殖教育部重點(diǎn)實(shí)驗(yàn)室,武漢 430070)

中國(guó)地方豬種質(zhì)資源十分豐富，約占世界總數(shù)的30%，且具有抗逆性強(qiáng)、耐粗飼和性成熟早等特點(diǎn)，然而多數(shù)品種存在脂肪含量高、飼料利用率低等缺點(diǎn)[1]。隨著中國(guó)居民消費(fèi)習(xí)慣的變化，瘦肉型豬的市場(chǎng)需求日益增長(zhǎng)，在20世紀(jì)70年代至20世紀(jì)80年代，中國(guó)逐漸開始引入皮特蘭豬和杜洛克豬等瘦肉型豬種在生豬繁育體系中作為終端父本使用[2-3]。在生豬生產(chǎn)中，皮特蘭豬與杜洛克豬作為終端父本，其生長(zhǎng)速度和飼料利用率長(zhǎng)期受到高強(qiáng)度人工選擇作用，且選育方向相似。人工選擇作用在基因組上留下的特征痕跡稱之為選擇信號(hào)(selection signature)[4]，針對(duì)商業(yè)豬種進(jìn)行選擇信號(hào)檢測(cè)，可以篩選出與豬重要經(jīng)濟(jì)性狀相關(guān)的分子標(biāo)記和候選基因，為豬重要經(jīng)濟(jì)性狀遺傳改良提供一定的參考依據(jù)。目前，在豬上進(jìn)行選擇信號(hào)檢測(cè)的研究報(bào)道了一系列重要的選擇信號(hào)候選區(qū)域，在不同的繁殖力群體中ESR、AREG、PRLP等多個(gè)基因受到了不同程度的選擇作用[5-6]；BARX2基因作為肌肉生長(zhǎng)、再生和維護(hù)的重要調(diào)節(jié)基因，在約克夏豬與長(zhǎng)白豬中均受到選擇[7]；影響豬椎骨數(shù)及體長(zhǎng)的多個(gè)基因(如NR6A1、PLAG1、LCORL)在歐洲商業(yè)豬和大多數(shù)歐洲地方豬群體中受到了強(qiáng)烈選擇[8]；ADAMTS12、SIM1和NOS1基因顯示了藏豬自然選擇的特征，可能與高海拔適應(yīng)性相關(guān)，在具有帶型毛色的中國(guó)地方豬EDNRB基因位點(diǎn)上發(fā)現(xiàn)了強(qiáng)烈的定向選擇信號(hào)，表明EDNRB基因是中國(guó)地方豬白色帶型的潛在候選基因[9]。

隨著畜禽選擇信號(hào)研究的深入，近年來群體間平行選擇(parallel selection)信號(hào)的檢測(cè)逐漸成為群體遺傳學(xué)領(lǐng)域研究趨同遺傳基礎(chǔ)的重要方向[10]。Rehkamper等[11]研究發(fā)現(xiàn)，鳥類與哺乳動(dòng)物大腦組織在結(jié)構(gòu)和功能上具有大量相似性，主要表現(xiàn)于端腦中多模態(tài)信息融合能力方面，揭示了脊椎動(dòng)物間發(fā)生了平行進(jìn)化。此外，在人、三刺魚及秀麗隱桿線蟲等物種中也發(fā)現(xiàn)諸多關(guān)于平行選擇的證據(jù)[12-14]。在畜禽相關(guān)研究方面，在蛋雞和肉雞群體中檢測(cè)到多個(gè)與體型外貌、生產(chǎn)性能相關(guān)的平行選擇候選基因，包括WWP1、BCDO2、TSHR、AGTR2和OPG等[[15-16]；在杜洛克豬和大白豬群體間檢測(cè)到67個(gè)平行選擇區(qū)域，并挖掘出多個(gè)調(diào)節(jié)豬生長(zhǎng)發(fā)育、體型大小等性狀的候選基因[17]。

本研究擬通過對(duì)皮特蘭豬、杜洛克豬2個(gè)商業(yè)豬種進(jìn)行選擇信號(hào)檢測(cè)，揭示其作為終端父本時(shí)重要經(jīng)濟(jì)性狀在人工選擇作用下的基因組變化特征，利用檢測(cè)到的群體間平行選擇候選區(qū)域，結(jié)合生物信息學(xué)分析揭示表型趨同的功能候選基因，為進(jìn)一步解析皮特蘭豬、杜洛克豬2個(gè)品種重要經(jīng)濟(jì)性狀的遺傳基礎(chǔ)，探討表型趨同的遺傳機(jī)制以及商業(yè)豬重要經(jīng)濟(jì)性狀的遺傳改良提供一定的參考依據(jù)。

1 材料與方法

1.1 材料

本研究共選取2 165頭豬作為研究對(duì)象，樣本分別來自中國(guó)南方2個(gè)生豬企業(yè)的4個(gè)群體，包括來自A豬場(chǎng)的376頭皮特蘭群體(PP)和451頭杜洛克品系Ⅰ群體(DDⅠ)，以及來自B豬場(chǎng)的841頭杜洛克品系Ⅱ群體(DDⅡ)和497頭杜洛克品系Ⅲ群體(DDⅢ)。對(duì)每頭豬采集耳組織樣，-20 ℃保存，用于提取全基因組DNA。

1.2 方法

1.2.1 質(zhì)量控制及基因型填充 使用全基因組芯片進(jìn)行基因分型，獲得4個(gè)群體的50K SNP芯片數(shù)據(jù)。使用Plink軟件[18]對(duì)50K SNP芯片數(shù)據(jù)進(jìn)行質(zhì)量控制，條件為：①SNP檢出率(call rate)>90%；②樣本的檢出率>90%；③最小等位基因頻率(minor allele frequency，MAF)為0.01。利用Beagle軟件[19]對(duì)缺失的基因型進(jìn)行填充，對(duì)填充后的群體計(jì)算個(gè)體間的親緣關(guān)系系數(shù)(PI_HAT)，并質(zhì)控掉親緣關(guān)系較近(PI_HAT>0.5)的個(gè)體。

1.2.2 選擇信號(hào)檢測(cè) 利用綜合單倍型評(píng)分(iHS)[20]與等位基因頻率差(△AF)2種方法分別進(jìn)行群體內(nèi)及群體間選擇信號(hào)檢測(cè)，△AF定義為2個(gè)群體間等位基因頻率差的絕對(duì)值。研究利用R語(yǔ)言中的REHH包[21]及Plink軟件，按照100 kb窗口、50 kb步長(zhǎng)在全基因組范圍內(nèi)劃分窗口進(jìn)行選擇信號(hào)檢測(cè)，對(duì)窗口統(tǒng)計(jì)量平均值進(jìn)行秩排序，取其前5%的窗口作為選擇信號(hào)顯著候選區(qū)域。

1.2.3 平行選擇信號(hào)候選區(qū)域鑒定 結(jié)合iHS和△AF進(jìn)行皮特蘭豬和杜洛克豬群體間的平行選擇信號(hào)檢測(cè)。4個(gè)群體兩兩組合，以PP、DDⅠ群體為例，利用bedtools軟件將iHS、△AF 2種統(tǒng)計(jì)量所對(duì)應(yīng)的選擇信號(hào)候選區(qū)域按照窗口左右200 kb的范圍進(jìn)行合并。提取iHSPP和iHSDDⅠ重疊區(qū)域，再提取iHSPP_iHSDDⅠ與△AFPP-DDⅠ的重疊區(qū)域，所得區(qū)域即定義為PP和DDⅠ 2個(gè)群體間的平行選擇候選區(qū)域。

1.2.4 選擇信號(hào)候選區(qū)域注釋 基于定位到的平行選擇候選區(qū)域，利用Ensembl[22]平臺(tái)的BioMart(http:∥www.biomart.org)軟件進(jìn)行基因挖掘；使用DAVID[23]開展基因集進(jìn)行富集分析，顯著富集條目定義為P<0.05；基于Ensembl注釋的基因，通過NCBI數(shù)據(jù)庫(kù)對(duì)其進(jìn)行功能注釋，挖掘與豬繁殖、生長(zhǎng)和肉質(zhì)等性狀相關(guān)的候選基因。

2 結(jié) 果

2.1 數(shù)據(jù)質(zhì)量控制

皮特蘭豬和杜洛克豬4個(gè)群體的SNP數(shù)據(jù)質(zhì)量控制結(jié)果見表1。由表1可知，質(zhì)控后DDⅠ群體獲得的SNPs數(shù)目最多(44 056個(gè))，DDⅢ群體獲得的SNP數(shù)目最少(37 323個(gè))。2 165個(gè)樣本中通過樣本檢出率的個(gè)體共2 140個(gè)，對(duì)其經(jīng)親緣關(guān)系系數(shù)過濾(PI_HAT>0.5)后，最終得到1 003個(gè)親緣關(guān)系較遠(yuǎn)個(gè)體。

表1 皮特蘭豬和杜洛克豬群體質(zhì)量控制結(jié)果

2.2 群體內(nèi)選擇信號(hào)候選區(qū)域檢測(cè)

基于iHS方法群體內(nèi)選擇信號(hào)檢測(cè)發(fā)現(xiàn)，在PP群體的25 091個(gè)iHS統(tǒng)計(jì)量中有1 254個(gè)統(tǒng)計(jì)量達(dá)到顯著水平，長(zhǎng)度為119.59 Mb，約占全基因組總長(zhǎng)的4.94%；在DDⅠ群體的27 195個(gè)iHS統(tǒng)計(jì)量中有1 359個(gè)統(tǒng)計(jì)量達(dá)到顯著水平，長(zhǎng)度為129.60 Mb，約占全基因組總長(zhǎng)的5.35%；在DDⅡ群體的27 582個(gè)iHS統(tǒng)計(jì)量中有1 379個(gè)統(tǒng)計(jì)量達(dá)到顯著水平，長(zhǎng)度為131.51 Mb，約占全基因組總長(zhǎng)的5.43%；在DDⅢ群體的22 398個(gè)iHS統(tǒng)計(jì)量中有1 120個(gè)統(tǒng)計(jì)量達(dá)到顯著水平，長(zhǎng)度為106.81 Mb，約占全基因組總長(zhǎng)的4.41%(圖1)。

A～D，PP、DDⅠ、DDⅡ、DDⅢ群體A-D，PP，DDⅠ，DDⅡ and DDⅢ populations,respectively圖1 皮特蘭豬和不同品系杜洛克豬群體內(nèi)信號(hào)檢測(cè)結(jié)果Fig.1 Ingroup selection signature detection results of Pietrain pigs and different strains of Duroc pigs

2.3 群體間選擇信號(hào)候選區(qū)域鑒定

由表2可知，基于△AF方法群體間選擇信號(hào)顯著候選區(qū)域共有9 579個(gè)，總長(zhǎng)度為913.50 Mb，其中，PP-DDⅠ、PP-DDⅡ、PP-DDⅢ、DDⅠ-DDⅡ、DDⅠ-DDⅢ、DDⅡ-DDⅢ群體對(duì)分別檢測(cè)到1 648、1 541、1 467、1 689、1 616和1 618個(gè)基因組候選區(qū)域，分別約占基因組總長(zhǎng)的6.49%、6.07%、5.78%、6.65%、6.38%和6.38%。

表2 皮特蘭豬和不同品系杜洛克豬群體間選擇信號(hào)檢測(cè)結(jié)果

2.4 平行選擇候選區(qū)域鑒定

將iHS統(tǒng)計(jì)量按照窗口左右各200 kb的范圍進(jìn)行合并處理，共檢測(cè)出1 591個(gè)群體內(nèi)選擇信號(hào)候選區(qū)域，其中PP群體390個(gè)，長(zhǎng)度為99.75 Mb；DDⅠ群體460個(gè)，長(zhǎng)度為110.82 Mb；DDⅡ群體400個(gè)，長(zhǎng)度為106.86 Mb；DDⅢ群體341個(gè)，長(zhǎng)度為86.93 Mb(表3)。

表3 群體內(nèi)選擇信號(hào)候選區(qū)域

將△AF統(tǒng)計(jì)量按照窗口左右各200 kb的范圍進(jìn)行合并處理，共發(fā)現(xiàn)3 466個(gè)群體間選擇信號(hào)候選區(qū)域，其中PP-DDⅠ群體間有851個(gè)候選區(qū)域，長(zhǎng)度為136.14 Mb；PP-DDⅡ群體間有795個(gè)，長(zhǎng)度為127.28 Mb；PP-DDⅢ群體間有733個(gè)，長(zhǎng)度為121.74 Mb；DDⅠ-DDⅡ群體間有798個(gè)，長(zhǎng)度為137.66 Mb；DDⅠ-DDⅢ群體間有744個(gè)，長(zhǎng)度為136.42 Mb；DDⅡ-DDⅢ群體間的候選區(qū)域有678個(gè)，長(zhǎng)度為127.36 Mb(表4)。

表4 群體間選擇信號(hào)候選區(qū)域

利用bedtools軟件在4個(gè)群體中共挖掘到52個(gè)平行選擇候選區(qū)域(54個(gè)區(qū)域中包含2個(gè)重復(fù)區(qū)域)，其長(zhǎng)度約為4.67 Mb，值得注意的是，在DDⅠ-DDⅢ群體間檢測(cè)到的5個(gè)平行選擇候選區(qū)域中，位于11及14號(hào)染色體的2個(gè)區(qū)域分別在PP-DDⅢ、DDⅠ-DDⅡ群體間被檢測(cè)到，長(zhǎng)度均約為0.095 Mb(圖2)。

2.5 平行選擇信號(hào)候選區(qū)域基因功能注釋

利用BioMart軟件挖掘平行選擇區(qū)域的基因，注釋到88個(gè)平行選擇候選基因。DAVID數(shù)據(jù)庫(kù)中GO功能富集結(jié)果顯示，DDⅠ-DDⅡ群體間的平行候選基因LIPA、LIPM參與脂質(zhì)分解過程；PP-DD Ⅰ群體間平行選擇候選區(qū)域功能注釋結(jié)果見表5，挖掘出多個(gè)與生殖系統(tǒng)發(fā)育(HSD17B4、FAM170A、AXDND1)、調(diào)節(jié)脂質(zhì)代謝和脂肪酸代謝(HSD17B4、SOAT1、COMMD6、CDCP1、ACSL5)、參與內(nèi)質(zhì)網(wǎng)修飾(ZDHHC6)等生物過程相關(guān)的候選基因；PP-DDⅡ群體間注釋到4個(gè)平行選擇候選基因(GALNT18、FER、LDLRAD4、PLCB4)，參與調(diào)節(jié)機(jī)體生長(zhǎng)發(fā)育及細(xì)胞增殖等；PP-DDⅢ群體間挖掘出多個(gè)與脂肪沉積相關(guān)的基因(KLHDC3、MEA1、GPC6)；DDⅠ-DDⅡ群體間挖掘到33個(gè)平行選擇候選基因；DDⅠ-DDⅢ群體間挖掘到8個(gè)候選基因；DDⅡ-DDⅢ群體間挖掘到11個(gè)候選基因，主要調(diào)控豬的胴體性狀(ACOXL、ZNF638、GCP6、SUPV3L1、LIPN、LIPM、ACTA2、LIPA、FFAR4、PDE6C、SMARCD1、SNCG)以及繁殖性狀(TMEM132D、VPS26A、RBP4、GDNF、SPACA3、BMPR1A)。

A～F，PP-DDⅠ、PP-DDⅡ、PP-DDⅢ、DDⅠ-DDⅡ、DDⅠ-DDⅢ和DDⅡ-DDⅢA-F，PP-DDⅠ，PP-DDⅡ，PP-DDⅢ，DDⅠ-DDⅡ，DDⅠ-DDⅢ and DDⅡ-DDⅢ,respectively圖2 群體間平行選擇候選區(qū)域數(shù)目Fig.2 Number of candidate regions for intergroup parallel selection

表5 PP-DDⅠ群體間平行選擇候選區(qū)域和候選基因

3 討 論

為研究皮特蘭豬和杜洛克豬作為終端父本在人工選擇作用下性狀趨同的基因組變化特征，本研究基于來自A豬場(chǎng)的PP、DDⅠ群體和來自B豬場(chǎng)的DDⅡ、DDⅢ群體的50K SNP芯片數(shù)據(jù)進(jìn)行選擇信號(hào)檢測(cè)，通過整合iHS和△AF在PP、DDⅠ、DDⅡ和DDⅢ 4個(gè)群體中的選擇信號(hào)檢測(cè)結(jié)果，挑選至少在每2個(gè)群體中都存在的重疊區(qū)域，共得到52個(gè)性狀趨同的基因組區(qū)域，其中，皮特蘭豬與3個(gè)杜洛克豬群體間共檢測(cè)到21個(gè)候選區(qū)域，而3個(gè)杜洛克豬群體間共篩選出33個(gè)平行選擇的基因組候選區(qū)域，表明由于品種內(nèi)不同品系間具有更高的遺傳背景，因此性狀趨同的基因組特征一致性更高。然而，PP與DDⅠ群體間共檢測(cè)到10個(gè)候選區(qū)域，表明不同品種在相似育種方向的驅(qū)動(dòng)下仍然能夠形成較高的基因組特征相似。此外，DDⅠ和DDⅡ群體間共篩選出18個(gè)候選區(qū)域，與兩群體均屬于美系杜洛克豬密切相關(guān)；盡管DDⅠ和DDⅢ群體均屬于杜洛克豬品種，但DDⅢ群體為丹系杜洛克豬且來源于B豬場(chǎng)，在DDⅠ和DDⅢ群體間僅檢測(cè)到5個(gè)性狀趨同候選區(qū)域，表明不同品系間杜洛克豬群體的人工選擇方向可能存在一定的差異。

本研究對(duì)52個(gè)性狀趨同區(qū)域進(jìn)行基因注釋，共挖掘出88個(gè)受選擇的候選基因，主要參與調(diào)節(jié)脂肪沉積、生長(zhǎng)發(fā)育、繁殖、免疫等性狀。GO功能富集到脂質(zhì)分解過程，涉及LIPA、LIPM2個(gè)候選基因。研究表明，小鼠酸性脂肪酶基因受組織mRNA表達(dá)差異的影響，LIPA基因表現(xiàn)為廣泛的組織表達(dá)，而LIPM基因表達(dá)受限，僅在表皮組織中表達(dá)[24-25]；ZNF638、FFAR4、PDE6C基因共同參與調(diào)節(jié)脂肪細(xì)胞的發(fā)育和分化[26-28]；ACOXL、ACSL5、CDCP1、COMMD6、HSD17B4、SOAT1、PPP2R5D、KLHDC3和ACTA2基因通過調(diào)節(jié)脂質(zhì)代謝影響豬的脂肪沉積，進(jìn)而影響豬的肌內(nèi)脂肪含量[29-36]；GNMT基因通過利用和分配營(yíng)養(yǎng)素影響肉質(zhì)性狀[37]；PTPRM基因可調(diào)節(jié)肌肉生長(zhǎng)發(fā)育[38]；SLC16A12基因參與肌酸的轉(zhuǎn)運(yùn)，從而調(diào)節(jié)機(jī)體肌肉的生長(zhǎng)發(fā)育[39]。豬的肉質(zhì)與肌內(nèi)脂肪含量和肌肉密切相關(guān)，因此上述基因在4個(gè)群體中受到趨同選擇。豬的生長(zhǎng)速度與瘦肉率長(zhǎng)期以來一直受到養(yǎng)殖戶和企業(yè)的關(guān)注，SUPV3L1基因是一種重要的發(fā)育調(diào)控基因，可調(diào)節(jié)機(jī)體的生長(zhǎng)發(fā)育，SUPV3L1基因的表達(dá)始于囊胚期，在所有胎兒組織和細(xì)胞類型中廣泛表達(dá)，成年小鼠中敲除SUPV3L1基因會(huì)導(dǎo)致早衰，包括肌肉、脂肪組織的損失和嚴(yán)重的皮膚異常[40]；PLCB4基因能調(diào)節(jié)豬的生長(zhǎng)和骨架，從而影響豬體重和身體構(gòu)象性狀[41]；CUL7基因?qū)ωi四肢和蹄筋堅(jiān)固性起著重要作用[42]。

皮特蘭豬和杜洛克豬作為終端父本，盡管其在總產(chǎn)仔數(shù)性狀上選擇強(qiáng)度始終較低，然而精液品質(zhì)始終是重點(diǎn)關(guān)注的性狀類型。研究發(fā)現(xiàn)，F(xiàn)AM170A基因可能通過調(diào)節(jié)其他基因的表達(dá)來介導(dǎo)其對(duì)精子細(xì)胞頭形和精子形成的影響[43]；AXDND1基因僅在人和小鼠的圓形和細(xì)長(zhǎng)的精子細(xì)胞中表達(dá)，可能通過調(diào)節(jié)精子管的動(dòng)力學(xué)、精子頭的形狀和精子鞭毛的組裝對(duì)精子發(fā)生和雄性生育起著重要作用[44]；KLHDC3與MEA1基因在人與小鼠中的同源基因均被證明與精子生成相關(guān)[45]；GDNF基因被證明與精原干細(xì)胞的維持密切相關(guān)[46]；SPACA3基因與精子頂體相關(guān)，直接影響精子和卵子結(jié)合過程，從而影響母畜受胎率[47]；SNCAIP和VPS26A基因調(diào)節(jié)胚胎發(fā)育，可間接影響產(chǎn)仔數(shù)；RBP4基因可影響豬繁殖性狀，相較于純合子個(gè)體，在雜合子個(gè)體中表現(xiàn)出較高的產(chǎn)仔數(shù)[48]。因此，在生產(chǎn)過程中，終端父本的繁殖性狀同第一母本一樣一直受到長(zhǎng)時(shí)間、高強(qiáng)度的人工選擇，這些基因的發(fā)現(xiàn)可為繁殖性狀的遺傳標(biāo)記提供一定參考。

此外，在集約化飼養(yǎng)條件下瘦肉型豬種的抗性實(shí)質(zhì)上長(zhǎng)期處于高強(qiáng)度選擇壓力之下，本研究發(fā)現(xiàn)，在豬的選育過程中CAMK4、PTCRA和IFIT5基因受到強(qiáng)烈的選擇作用，其中，CAMK4基因是一種多功能絲氨酸/蘇氨酸激酶，其通過激活包括T細(xì)胞和抗原呈遞細(xì)胞的各種免疫細(xì)胞中的轉(zhuǎn)錄因子來調(diào)節(jié)基因表達(dá)[49]；PTCRA基因通過編碼前T細(xì)胞抗原受體α，以在分化的不同階段調(diào)節(jié)早期T細(xì)胞發(fā)育[50]；IFIT5基因是先天免疫反應(yīng)的重要增強(qiáng)因子，在IFIT5基因存在的情況下，宿主細(xì)胞的抗病毒反應(yīng)顯著增加[51]。

4 結(jié) 論

皮特蘭豬和杜洛克豬各群體內(nèi)選擇信號(hào)顯著候選區(qū)域共有5 112個(gè)，群體間的選擇信號(hào)候選區(qū)域共有9 579個(gè)，在群體之間存在平行選擇的候選區(qū)域共52個(gè)，包含88個(gè)受到平行選擇的基因；基因富集及功能注釋發(fā)現(xiàn)這些基因主要與豬的肉質(zhì)、生長(zhǎng)發(fā)育及繁殖性狀相關(guān)，表明高強(qiáng)度的人工選擇作用是造成繁殖力、胴體、肉質(zhì)等重要經(jīng)濟(jì)性狀表型趨同的重要因素。