尹彥鎮(zhèn)，侯黎明，劉航，陶偉，石傳宗，劉鍇月，張萍，牛培培，李強，李平華，黃瑞華

基因組分析對豬乳頭數(shù)相關(guān)數(shù)量性狀基因座鑒定

1南京農(nóng)業(yè)大學(xué)養(yǎng)豬研究所，南京 210095；2南京農(nóng)業(yè)大學(xué)淮安研究院，江蘇淮安 223001；3淮安市淮陰新淮種豬場，江蘇淮安 223322

【目的】分析乳頭數(shù)的變異，挖掘與乳頭數(shù)相關(guān)的數(shù)量性狀基因座（quantitative trait locus，QTL）和候選基因，為豬乳頭數(shù)的選育研究提供重要分子標(biāo)記?！痉椒ā繙?zhǔn)確測定了709頭蘇淮豬（335頭育肥豬和374頭種豬）的左、右和總?cè)轭^數(shù)。對蘇淮育肥豬進行80K芯片分型，并使用芯片數(shù)據(jù)計算左、右、總?cè)轭^數(shù)的遺傳力和基因組估計育種值（genomic estimated breeding value, GEBV）?；谌轭^數(shù)GEBV和表型排名，選擇前10%的個體以及后10%的個體進行群體分化指數(shù)分析（fixation Index,F）檢測高度分化的位點。接著，通過全基因組關(guān)聯(lián)分析（genome wide association analysis, GWAS）鑒定與乳頭數(shù)關(guān)聯(lián)的位點，選擇高度分化且與乳頭數(shù)顯著關(guān)聯(lián)的位點作為候選位點，選擇位于候選位點附近且功能注釋后與乳頭數(shù)相關(guān)的基因作為候選基因。最后，選擇每個染色體上最顯著的候選位點對709頭蘇淮豬進行乳頭數(shù)關(guān)聯(lián)分析，以驗證上述位點的顯著性?！窘Y(jié)果】蘇淮育肥豬左、右、總?cè)轭^數(shù)的變異系數(shù)分別為10.20%、9.26%、8.50%，遺傳力分別為0.212、0.257、0.312?；?i>F和GWAS分析，總共在7、13、16、18號染色體（Chromosome, SSC）上鑒定到20個乳頭數(shù)的候選位點，這些候選位點可解釋5.49%—8.03%的表型方差。其中，SSC7上與總?cè)轭^數(shù)關(guān)聯(lián)的位點rs80894106與文獻(xiàn)中報道的影響大白和杜洛克豬總?cè)轭^數(shù)的候選位點一致，但左乳頭的候選位點rs81444134（26.51 Mb, SSC13）和rs81233299（8.13 Mb, SSC18）均為新發(fā)現(xiàn)的與乳頭數(shù)相關(guān)的位點。左、右、總?cè)轭^的候選位點主要集中在SSC16上的6.36—10.66 Mb區(qū)間;連鎖不平衡（linkage disequilibrium, LD）分析發(fā)現(xiàn)，區(qū)間內(nèi)7.47—8.27 Mb的候選位點擬合成了一個795 kb的單倍型塊，且該單倍型塊是新發(fā)現(xiàn)的影響乳頭數(shù)的候選區(qū)域；單倍型塊內(nèi)的rs337606862（7.47 Mb）與右乳頭和總?cè)轭^最顯著關(guān)聯(lián)，單倍型塊內(nèi)的3個位點均位于cadherin 18（）基因的內(nèi)含子上，編碼Ⅱ型鈣黏附素，且鈣黏附素與發(fā)育中組織細(xì)胞的識別、分選、增殖、凋亡以及乳腺癌的發(fā)生有關(guān)。因此，可能是新的影響豬乳頭數(shù)的候選基因。再者，本研究對4個染色體上最顯著的位點：rs81444134、rs80894106、rs337606862、rs81233299在709頭蘇淮豬中基因分型，經(jīng)關(guān)聯(lián)分析后發(fā)現(xiàn)，這些位點均與乳頭數(shù)顯著相關(guān)，可以作為潛在分子標(biāo)記用于乳頭數(shù)的選育?！窘Y(jié)論】本研究通過基因組分析在蘇淮豬群體中鑒定到20個與乳頭數(shù)顯著相關(guān)的位點。其中SSC13上的26.51 Mb和SSC18上的8.13 Mb是新的乳頭數(shù)的候選QTLs，SSC16上的7.47—8.27 Mb也是新發(fā)現(xiàn)的乳頭數(shù)的候選QTL，且區(qū)間內(nèi)可能是新的影響豬乳頭形成的候選基因。

F；GWAS；數(shù)量性狀基因座（quantitative trait locus，QTL）；蘇淮豬；乳頭數(shù)；候選基因

0 引言

【研究意義】豬的乳頭一般呈圓柱狀，分布在腹部中線的兩側(cè)。根據(jù)分布位置的差異，豬的乳頭劃分為左乳頭和右乳頭。發(fā)育正常的乳頭具備泌乳的功能，能夠為仔豬提供生長發(fā)育所必須的營養(yǎng)物質(zhì)。因此，擁有更多功能正常乳頭的母豬能夠哺乳更多的初生仔豬。研究表明，仔豬利用的母豬功能乳頭的數(shù)目與仔豬的成活率密切相關(guān)[1]。同時，KIM等[2]研究了商業(yè)豬種乳頭數(shù)與產(chǎn)仔數(shù)的關(guān)系，發(fā)現(xiàn)乳頭數(shù)多的母豬產(chǎn)仔數(shù)和斷奶仔豬數(shù)也更多。因此，乳頭數(shù)是生豬生產(chǎn)中非常重要的繁殖性狀。此外，隨著育種技術(shù)的應(yīng)用，以丹麥大白、長白為代表的多豬種的產(chǎn)仔數(shù)性狀得到了有效的選育提升，但其乳頭數(shù)未經(jīng)系統(tǒng)選育提升，導(dǎo)致仔豬數(shù)多于乳頭數(shù)，會使得一部分仔豬需要通過寄養(yǎng)或人工哺乳的方式存活，增加生產(chǎn)成本和管理難度。因此，乳頭數(shù)的性狀改良是當(dāng)下育種研究的重點之一?！厩叭搜芯窟M展】乳頭數(shù)性狀由多基因控制，遺傳機理復(fù)雜。僅基于表型對乳頭數(shù)性狀進行選育的遺傳進展較慢。結(jié)合分子標(biāo)記選育能有效加快性狀改良的進程[3]。SNP芯片分型技術(shù)能夠高效、準(zhǔn)確地獲得大量SNPs的基因型，推動了F[4]分析和GWAS分析在鑒定復(fù)雜性狀候選QTLs和基因上的應(yīng)用。HE等[5]按照產(chǎn)仔數(shù)EBV排名選擇前10%和后10%的二花臉豬進行F分析，鑒定到、和等3個潛在影響產(chǎn)仔數(shù)的候選基因。李開軍等[6]采取相同的選擇策略針對蘇淮豬的中性洗滌纖維表觀消化率進行F分析，鑒定到多個與纖維表觀消化率性狀相關(guān)的基因和SNPs。在商業(yè)豬種乳頭數(shù)的GWAS分析中，MOSCATELLI等[7]使用60K芯片分型了821頭意大利大白豬，在6、7號染色體上分別鑒定到顯著影響前后乳頭對稱性以及乳頭數(shù)的QTLs，并將和作為兩種性狀對應(yīng)的候選基因。LI等[8]在杜洛克豬的GWAS分析中新發(fā)現(xiàn)了影響乳頭數(shù)的和影響左右乳頭對稱性的。在國內(nèi)地方豬種的GWAS分析中，TANG等[9]在白杜洛克和二花臉的雜交群體中定位到影響乳頭數(shù)的和。ZHOU等[10]基于GWAS分析在1、13和X染色體上鑒定到多個與蘇山豬乳頭數(shù)顯著關(guān)聯(lián)的SNPs。目前，豬QTLdb數(shù)據(jù)庫中總共收錄了624個與乳頭數(shù)關(guān)聯(lián)的QTLs，每條染色體上均有分布（https://www. animalgenome.org/cgi-bin/QTLdb/index/）?！颈狙芯壳腥朦c】盡管已有多數(shù)關(guān)于乳頭數(shù)的研究報道，但是影響乳頭數(shù)的主效基因和因果突變?nèi)晕粗Ｒ延醒芯恐饕性趪馍虡I(yè)豬種，其研究結(jié)果通常存在群體異質(zhì)性，而國內(nèi)地方豬種或含有地方豬種血緣的培育豬種的研究有待系統(tǒng)開展。蘇淮豬是以本研究團隊作為技術(shù)支撐單位，歷經(jīng)12年時間，于2010年培育而成的國家級瘦肉型豬新品種，含75%的大白豬和25%的淮豬血統(tǒng)[11]。該群體在包括乳頭數(shù)在內(nèi)的性狀上存在一定的變異，是鑒別影響豬乳頭數(shù)QTL、候選基因以及解析乳頭數(shù)變異遺傳機制的良好素材?！緮M解決的關(guān)鍵問題】本研究以蘇淮豬為試驗動物，統(tǒng)計了群體內(nèi)乳頭數(shù)的分布狀況，并使用SNP芯片進行個體基因分型?；赟NP芯片數(shù)據(jù)估計了乳頭數(shù)的遺傳力，選取多乳頭數(shù)個體和少乳頭數(shù)個體進行F分析，并在整個群體進行了表型和基因型間的GWAS分析，結(jié)合F分析和GWAS分析挖掘顯著影響乳頭數(shù)的候選QTLs和基因，并對候選基因進行功能注釋。本研究的結(jié)果旨在為研究豬乳頭數(shù)的遺傳基礎(chǔ)提供理論參考，同時為蘇淮豬乳頭數(shù)的分子育種提供關(guān)鍵分子標(biāo)記。

1 材料與方法

1.1 試驗動物

試驗選用335頭蘇淮育肥豬、374頭蘇淮種豬（25頭公和349頭母豬）共709頭蘇淮豬均在江蘇省淮安市淮陰新淮種豬場內(nèi)飼養(yǎng)。在2019年內(nèi)，準(zhǔn)確記錄其左、右乳頭數(shù)及總?cè)轭^數(shù)。蘇淮育肥豬作為試驗群體用于鑒定與乳頭數(shù)顯著關(guān)聯(lián)的SNPs，全部群體用于驗證部分顯著的SNPs的顯著性。

1.2 試驗方法

1.2.1 DNA提取與基因分型 采集709頭蘇淮豬的耳組織并貯存于2 mL離心管中。離心管內(nèi)裝有2/3體積的75%酒精且外側(cè)管壁貼有膠帶以記錄耳號。取小部分耳組織用于提取DNA，剩余樣品保存至-20 ℃冰箱。DNA提取的試劑及步驟參見天根DNA提取試劑盒說明。試驗群體中，濃度和質(zhì)量符合要求的DNA進一步使用Geneseek GGP 80K芯片進行SNP分型。分型后，使用PLINK[12]軟件基于以下標(biāo)準(zhǔn)進行質(zhì)量控制：SNP在個體的檢出率及個體的SNP檢出率均大于0.9，次等位基因頻率大于0.05。總共有334個個體（其中有一個個體因未達(dá)到質(zhì)控標(biāo)準(zhǔn)被去掉）及50 504個SNPs通過了質(zhì)控并用于后續(xù)的分析。

1.2.2 主成分分析 群體內(nèi)分層可能會導(dǎo)致不同的亞群內(nèi)部個體之間的關(guān)系比實際更緊密，進而影響總體的關(guān)聯(lián)性評估。本研究使用PLINK[12]軟件執(zhí)行了主成分分析，計算了前10個主成分的特征值并以主成分1和主成分2分別作為橫、縱坐標(biāo)繪制PCA圖。

1.2.3 遺傳力評估 使用DMU[13]軟件估計試驗群體單個性狀的方差組分及個體育種值。模型如下：

y=+sex+ad+e

其中，y代表左、右、總?cè)轭^數(shù)；為群體表型均值；sex為個體的性別；個體的加性效應(yīng)ad，服從正態(tài)分布ad-(0,)；殘差項e，服從正態(tài)分布e-(0,)；是由invgmatrix[14]軟件使用50 504個SNPs構(gòu)建的基因組親緣關(guān)系矩陣；是單位陣；和分別為加性方差和殘差方差?；蚪M遺傳力的計算公式為：2/(+)。

1.2.4 群體分化指數(shù)分析 研究進行了試驗群體內(nèi)多乳頭數(shù)和少乳頭數(shù)兩組間的F分析，以檢測多乳頭數(shù)個體與少乳頭數(shù)的個體之間顯著分化的位點。分組規(guī)則如下：首先，按照左、右、總?cè)轭^數(shù)的GEBV分別進行排序，各選擇排名前10%的個體和排名后10%的個體分別作為F分析的多乳頭組和少乳頭組；其次，由于GEBV是模型估計的值，可能存在一定的偏差，因此剔除了多乳頭組中表型值低于平均值的個體以及少乳頭組中表型值高于平均值的個體。最終，左、總?cè)轭^數(shù)的多乳頭組及右、總?cè)轭^數(shù)的少乳頭組均有32個個體，其余組為33個個體，并使用獨立樣本t檢驗估計極端乳頭組間表型和GEBV值的差異性。F分析的模型參照WEIR等[15]的研究模型并且使用VCFTOOLS[16]軟件計算全部SNPs的F值并取絕對值。取F值的絕對值前5%的位點作為顯著分化的位點。左、右、總?cè)轭^數(shù)F分析的閾值分別為0.1674、0.1469、0.1482，均接近WRIGHT等[4]劃分的高度分化的閾值條件（0.15＜F＜0.25）。

1.2.5 全基因組關(guān)聯(lián)分析 試驗群體乳頭數(shù)與全基因組50 504個SNPs間的關(guān)聯(lián)性檢驗通過LDAK[17]軟件的Wald檢驗實現(xiàn)。模型如下：

=++++

其中，為乳頭數(shù)表型；為群體均值；為SNP基因型矩陣，為基因型效應(yīng)值；為固定效應(yīng)矩陣，為性別效應(yīng)；為剩余多基因效應(yīng)，服從正態(tài)分布-(0,)，為校正的基因組親緣關(guān)系矩陣；為殘差，服從正態(tài)分布-(0,)；其余字母含義同遺傳力估計的模型。由于GWAS分析共執(zhí)行了50 504次假設(shè)檢驗，因此，采用BONFERRONI[18]方法計算基因組顯著的閾值為0.05/50504（9.90×10-7）?？紤]到BONFERRONI[18]方法的嚴(yán)格性，染色體顯著的閾值定義為1/50504（1.98×10-5）[19]。其次，在估計單個SNP效應(yīng)時，GWAS模型同時添加了SNP的固定效應(yīng)和多基因效應(yīng)，造成了效應(yīng)的重復(fù)校正。參考SCHMID等[20]的方法，即每次估計SNP效應(yīng)時，剔除該SNP所在染色體，使用剩余染色體的SNPs構(gòu)建矩陣。GWAS的結(jié)果包括曼哈頓圖和Q-Q圖由R包CMplot[21]可視化。最后，使用Haploview4.2[22]軟件繪制16號染色體上6.36—10.66 Mb的候選SNPs的連鎖不平衡圖，且只計算1 Mb內(nèi)的SNPs之間的LD值。

1.2.6 顯著SNPs解釋的表型方差比例 為進一步評估試驗群體的GWAS分析，計算了顯著SNPs所解釋性狀的表型方差的比例。計算公式同ZHANG等[23]：=2(1-)2/。為解釋的表型方差比例；為最小等位基因頻率；為SNP的效應(yīng)值；為表型方差，并且代表觀測值的離散程度。

1.2.7 候選SNPs所在候選基因的鑒別與功能注釋 本研究使用的80K芯片上的SNPs的物理位置是根據(jù)豬11.1參考基因組注釋。候選SNPs所在或者附近的基因作為候選基因。候選基因名稱通過Ensembl數(shù)據(jù)庫（http://www.ensembl.org/）查找。候選基因的功能使用GeneCards（https://www.genecards.org/）和Pubmed數(shù)據(jù)庫（https://pubmed.ncbi.nlm.nih.gov/）注釋。

1.2.8 候選SNPs與全部群體乳頭數(shù)關(guān)聯(lián)分析 本研究將左、右、總?cè)轭^數(shù)的F以及GWAS分析共同定位的顯著SNPs作為候選SNPs。本研究選取了每條染色體上最顯著的候選SNPs以驗證其在全部群體中的顯著性。上述SNPs的引物在NCBI的Primer-BLAST模塊（https://www.ncbi.nlm.nih.gov/tools/primer-blast/）設(shè)計，引物信息如表1，PCR擴增流程見表2，擴增產(chǎn)物通過Sanger測序判定SNPs的基因型。

表1 候選SNPs的引物信息

表2 PCR擴增流程

使用SAS 9.4（https://odamid-apse1.oda.sas.com/ SASStudio）的GLM模型對SNPs基因型和全部群體乳頭數(shù)進行關(guān)聯(lián)分析，并使用Bonferroni方法進行多重比較。模型如下：

y=+sex+g+e

其中，g為SNPs的固定效應(yīng)，其余字母含義同遺傳力估計的模型。

2 結(jié)果

2.1 蘇淮豬試驗群體乳頭數(shù)描述性統(tǒng)計及遺傳力估計

表3展示了試驗群體334頭蘇淮豬乳頭數(shù)的分布、變異情況及遺傳力。左、右、總?cè)轭^數(shù)的范圍分別是5—11、5—10、11—21，表型的分布比較分散。左、右、總?cè)轭^數(shù)的變異系數(shù)分別為10.20%、9.26%、8.50%。蘇淮豬左、右、總?cè)轭^數(shù)的遺傳力分別為0.212、0.257、0.312，屬于中高等遺傳力性狀。

表3 334頭蘇淮豬乳頭數(shù)描述性統(tǒng)計及遺傳力估計

LTN，左乳頭；RTN，右乳頭；TTN，總?cè)轭^；SE，標(biāo)準(zhǔn)誤；CV，變異系數(shù)。表7同

LTN, Left Teat Number; RTN, Right Teat Number; TTN, Total Teat Number; SE Standard Error; CV, Coefficient Variation. The same as table 7

2.2 蘇淮豬試驗群體分層檢驗

基于50 504個SNPs的基因型計算的主成分1和主成分2占總體特征值的百分比是相近的，并且群體內(nèi)沒有明顯離群的個體（圖1）。因此，334頭蘇淮豬的群體內(nèi)不存在更小的亞群。

2.3 蘇淮豬試驗群體乳頭數(shù)FST分析

試驗群體中用于F分析的多乳頭組的表型均值及GEBV均值均顯著高于少乳頭組（表4）。3個乳頭數(shù)性狀的F分析均鑒定到2 525個在多乳頭組和少乳頭組之間高度分化的SNPs（圖2）。左乳頭F值排名前五的SNPs位于SSC10的28.10—28.37 Mb；右乳頭F值排名前五的SNPs位于SSC8、SSC13、SSC14上；總?cè)轭^F值排名前五的SNPs位于SSC8、SSC13、SSC16上（表5）。

圖1 334頭蘇淮豬的PCA圖

表4 極端乳頭數(shù)個體GEBV及表型分布

極端乳頭組間同列不同字母表示差異顯著

Different letters in the same column between groups with extreme teats indicated significant differences

A：左乳頭；B：右乳頭；C：總?cè)轭^。X軸為1-18號染色體，Y軸為FST值，紅色虛線代表閾值線

表5 蘇淮豬乳頭數(shù)FST值排名前五的位點

2.4 蘇淮豬試驗群體乳頭數(shù)GWAS分析

通過GWAS分析以更精確地定位影響乳頭數(shù)的位點。左、右、總?cè)轭^數(shù)的GWAS分析的基因膨脹因子分別為1.07、1.179、1.035，并且Q-Q圖（圖3）中的觀測值和期望值沒有明顯的分離，說明不存在群體分層導(dǎo)致的膨脹現(xiàn)象。如表6和圖3所示，GWAS分析在SSC7、SSC13、SSC16、SSC18上鑒定到了22個與乳頭數(shù)顯著關(guān)聯(lián)的SNPs，其中16號染色體存在大量連續(xù)顯著的位點。顯著關(guān)聯(lián)的位點中，除了rs81315330 和rs319610259，剩余SNPs的F值均超過了F分析的閾值線。

將在極端乳頭數(shù)個體間高度分化且與乳頭數(shù)顯著關(guān)聯(lián)的20個位點認(rèn)定為潛在地影響乳頭數(shù)的候選位點。在候選位點中，SSC13上的rs81444134以及SSC18上的rs81233299均為左乳頭的候選位點，均未在之前的研究中被報道，且分別解釋了6.16%和6.62%的表型方差。其中，rs81444134在右、總?cè)轭^的F分析中的F值排名前五。SSC7上的rs80894106為總?cè)轭^的候選位點，解釋了8.03%的表型方差。其余的左、右、總?cè)轭^的候選位點集中位于SSC16的6.36—10.66 Mb區(qū)間內(nèi)。LD分析顯示，SSC16上7.47—8.27 Mb內(nèi)的7個候選SNPs形成了一個795 kb的單倍型塊（圖4），且其中的位點rs325271937（7.99 Mb）、rs343435106（8.04 Mb）、rs336743934（8.27 Mb）均位于cadherin 18（）基因上。單倍型塊中的位點rs337606862（7.47 Mb，SSC16）在右、總?cè)轭^的GWAS分析中均達(dá)到了全基因組顯著水平，且其右乳頭以及總?cè)轭^的PVE分別為7.73%和7.57%。目前在QTL數(shù)據(jù)庫中，SSC16上的7.47—8.27 Mb區(qū)間沒有發(fā)現(xiàn)已報道的影響乳頭數(shù)的QTL，因此該區(qū)域可能是一個新的影響乳頭數(shù)的QTL，且rs337606862是一個重要的候選SNP。

2.5 候選基因注釋

SSC16上7.47—8.27 Mb的單倍型塊內(nèi)包含了和cadherin 18（）2個編碼蛋白的基因（圖4）。編碼一類鈣黏蛋白，與組織發(fā)育中細(xì)胞的識別、分選等有關(guān)，是潛在影響乳頭數(shù)的候選基因。SSC7上的rs80894106的上游28.35 kb處為，下游5.66 kb處為，其中為已報道的影響乳頭數(shù)的候選基因。距離rs81444134（SSC13）和rs81233299（SSC18）最近的編碼蛋白的基因分別為golgi associated kinase 1A（，38.9 kb）和taste 2 receptor member 3（，2.5 kb）,暫未發(fā)現(xiàn)這兩個基因與乳頭發(fā)育有關(guān)。

A：左乳頭；B：右乳頭；C：總?cè)轭^。X軸為1-18號染色體，Y軸為-log10(P)。曼哈頓圖的實線和虛線分別代表染色體顯著閾值線及基因組顯著閾值線

方格內(nèi)的數(shù)字表示兩個SNP之間的D’值，深紅色意味著高的連鎖程度。共有兩個蛋白編碼基因位于795 kb的單倍型塊內(nèi)

表6 蘇淮豬乳頭數(shù)GWAS鑒定的顯著位點的詳細(xì)信息

2.6 SSC7、SSC13、SSC16、SSC18上最顯著位點與蘇淮豬全部群體乳頭數(shù)關(guān)聯(lián)分析

選擇4條染色體上最顯著位點與709頭蘇淮豬的左、右、總?cè)轭^數(shù)進行關(guān)聯(lián)分析。709頭蘇淮豬的左、右、總?cè)轭^數(shù)的變異系數(shù)分別為8.58%、8.22%和7.34%（表7）。關(guān)聯(lián)分析結(jié)果如表8，4個位點均與左、右、總?cè)轭^數(shù)顯著關(guān)聯(lián)（＜0.05），說明F和GWAS分析鑒定的候選位點在樣本量更大的全部群體中仍然與乳頭數(shù)具有關(guān)聯(lián)性。因此，這4個位點可以作為分子標(biāo)記用于蘇淮豬乳頭數(shù)的選育提升。

表7 709頭蘇淮豬乳頭數(shù)描述性統(tǒng)計

表8 4個位點與蘇淮豬全部群體乳頭數(shù)關(guān)聯(lián)分析

同行數(shù)字上標(biāo)不同字母表示差異顯著（<0.05）；rs81444134在80K芯片中有一個個體缺失，其有基因型的個體數(shù)為708頭

Different superscripts of numbers in the same line indicate significant differences; rs81444134 has one missing individual in 80K chip, and the number of genotyped individuals is 708

3 討論

3.1 蘇淮豬乳頭數(shù)變異統(tǒng)計及遺傳力估計

蘇淮豬是一個培育品種，群體內(nèi)各性狀可能存在一定的變異。本研究統(tǒng)計的334頭蘇淮育肥豬以及全部群體的709頭蘇淮豬的左、右、總?cè)轭^數(shù)的變異系數(shù)均超過了7%，這也證實了這一論述。周李生等[24]統(tǒng)計了中外豬種總?cè)轭^數(shù)的分布情況，丹麥、英國、美國等豬種的乳頭集中在12—16個，中國豬種的乳頭分布在10—22個之間。本研究的蘇淮豬的總?cè)轭^在11—21個之間，說明該群體是研究乳頭數(shù)變異機制的理想群體。在遺傳力估計上，多數(shù)研究估計的大白、長白、杜洛克豬的總?cè)轭^數(shù)的遺傳力均在0.2—0.5之間[8, 25-26]，與蘇淮豬遺傳力的結(jié)果相近。同時，LI等[8]估計的杜洛克豬的左乳頭（0.148）遺傳力低于右乳頭（0.177），這與蘇淮豬的結(jié)果一致。

3.2 影響蘇淮豬乳頭數(shù)QTL的鑒定

在實際的選種中，蘇淮豬育種者傾向選擇多乳頭的母豬留作種用，這種人工選擇作用會使得有利等位基因的頻率上升并在多乳頭的群體和少乳頭的群體之間產(chǎn)生基因頻率的差異。因此，期望通過F分析檢測極端乳頭數(shù)個體之間等位基因頻率差異顯著的位點。由于乳頭數(shù)的遺傳力中等，加上群體數(shù)量有限，因此結(jié)合育種值和表型值篩選極端乳頭的個體是較優(yōu)的方法。同時，相較于僅依靠系譜信息的BLUP模型，基于基因組信息的GBLUP模型估計的育種值更準(zhǔn)確[27]，且廣泛應(yīng)用于畜禽的育種中[28]。因此，本研究結(jié)合GEBV和表型值確定極端乳頭數(shù)的個體，而且統(tǒng)計發(fā)現(xiàn)本研究選擇的極端乳頭數(shù)組間的乳頭數(shù)GEBV和表型值均存在顯著差異，證明樣品選擇的合理性。F分析檢測到的高度分化的位點在18條常染色體上均有分布。由于F分析僅選擇了極端乳頭數(shù)的個體，群體量較小，檢測到的顯著分化的位點可能含有因抽樣誤差產(chǎn)生的假陽性的結(jié)果。因此，后續(xù)又開展了試驗群體的GWAS分析，并選定在F分析中高度分化且在GWAS分析中與乳頭數(shù)顯著關(guān)聯(lián)的位點為更可能的影響乳頭數(shù)的候選位點。

候選位點中，rs81444134（26.51 Mb，SSC13）和rs81233299（8.13 Mb，SSC18）均未在之前的研究中報道，是本研究新發(fā)現(xiàn)的乳頭數(shù)的候選位點。其中，rs81444134的F值在右、總?cè)轭^的F分析中排名前五。rs81444134和rs81233299附近的基因為和，且Genecards數(shù)據(jù)庫顯示這兩個基因分別與高爾基體激酶以及味道受體有關(guān)，但其影響乳頭形成的功能有待進一步研究。在GWAS中，SSC7上97.65 Mb的候選位點rs80894106與總?cè)轭^數(shù)顯著關(guān)聯(lián)，且與MOSCATELLI等[7]和LI等[8]的研究中報道的影響大白豬和杜洛克豬的總?cè)轭^數(shù)的位點相一致。左、右、總?cè)轭^的GWAS均在SSC16上出現(xiàn)連續(xù)的信號峰，且顯著信號集中在6.36—10.66 Mb。同時，3種乳頭數(shù)的F和GWAS分析在該區(qū)域均有一致的信號趨勢。本研究進一步對該區(qū)域進行LD分析，發(fā)現(xiàn)7.47—8.27 Mb區(qū)域內(nèi)的候選SNPs之間高度連鎖，形成了一個795 kb的單倍型塊，而目前未有研究報道該單倍型區(qū)域與乳頭數(shù)相關(guān)。單倍型塊內(nèi)的位點rs337606862與右乳頭和總?cè)轭^均是最顯著關(guān)聯(lián)的，值分別為2.19×10-7和2.95×10-7。ZHUANG等[29]在加拿大杜洛克群體的GWAS中鑒定到SSC16上6.12 Mb處的總?cè)轭^的QTL；BIDANEL等[30]使用梅山豬和大白豬的雜交群體鑒定的乳頭數(shù)的QTL位于SSC16的6.3—6.4 Mb處；這兩個QTLs分別距離本研究的7.47—8.27 Mb的單倍型塊約1.37 Mb和1.07 Mb遠(yuǎn)。因此，SSC16上的7.47—8.27 Mb區(qū)域應(yīng)是一個新的影響乳頭數(shù)的QTL，且位點rs337606862是QTL內(nèi)具有較大效應(yīng)的SNP。區(qū)域7.47—8.27 Mb包含了2個編碼蛋白的基因，其中的1個基因暫未有文獻(xiàn)研究，另外一個是一種鈣黏蛋白基因。生成Ⅱ型的鈣黏附素，在小鼠發(fā)育中的性腺細(xì)胞中高表達(dá)[31]。鈣黏附素在機體發(fā)育過程，包括胚胎期參與識別與分選不同種類的細(xì)胞，在組織形成中起到重要作用[32]。而豬的乳頭大概在胚胎期的23—26 d左右形成[33]。同時，鈣黏附素也可以調(diào)節(jié)細(xì)胞的增殖和凋亡活動[34]并且其功能的改變可能會導(dǎo)致乳腺癌的發(fā)生[35]。因此，推測是一個新的影響豬乳頭形成的候選基因，但該基因影響乳頭形成的機制還需要進一步研究。

為了進一步驗證試驗群體鑒定的乳頭數(shù)的候選位點的可靠性，本研究進一步選擇了包含試驗群體在內(nèi)的共709頭蘇淮豬群體并進行了最顯著的SNPs：rs80894106（SSC7）、rs81444134（SSC13）、rs337606862（SSC16）、rs81233299（SSC18）的基因型與蘇淮豬群體乳頭數(shù)的關(guān)聯(lián)分析。結(jié)果表明，上述位點仍然與709頭蘇淮豬的乳頭數(shù)顯著關(guān)聯(lián)。盡管試驗群體使用了混合線性模型，而整個群體的關(guān)聯(lián)分析由于沒有足夠的親緣信息，沒有考慮隨機加性效應(yīng)，但上述結(jié)果一定程度上驗證了4個最顯著位點與蘇淮豬乳頭數(shù)關(guān)聯(lián)的可靠性，可作為乳頭數(shù)選育的潛在分子標(biāo)記。

4 結(jié)論

本研究發(fā)現(xiàn)蘇淮豬群體內(nèi)乳頭數(shù)存在一定的變異，通過基因組F分析和GWAS分析，在蘇淮豬中鑒定到20個與乳頭數(shù)顯著相關(guān)的位點。其中，SSC13的26.51 Mb和SSC18的8.13 Mb是新的乳頭數(shù)的候選QTLs；SSC16的7.47—8.27 Mb區(qū)間存在連續(xù)與乳頭數(shù)顯著關(guān)聯(lián)的位點，是新發(fā)現(xiàn)的影響乳頭數(shù)的候選QTL區(qū)域，該區(qū)域內(nèi)可能是新的影響豬乳頭形成的候選基因。

[1] ALEXOPOULOS J G, LINES D S, HALLETT S, PLUSH K J. A review of success factors for piglet fostering in lactation. Animals: an Open Access Journal from MDPI, 2018, 8(3): 38.

[2] KIM J S, JIN D I, LEE J H, SON D S, LEE S H, YI Y J, PARK C S. Effects of teat number on litter size in gilts. Animal Reproduction Science, 2005, 90(1/2): 111-116.

[3] ROHRER G A, NONNEMAN D J. Genetic analysis of teat number in pigs reveals some developmental pathways independent of vertebra number and several loci which only affect a specific side. Genetics, Selection, Evolution: GSE, 2017, 49(1): 4.

[4] WRIGHT S. The genetical structure of populations. Annals of Eugenics, 1951, 15(4): 323-354. doi:10.1111/j.1469-1809.1949. tb02451.x.

[5] HE L C, LI P H, MA X, SUI S P, GAO S, KIM S W, GU Y Q, HUANG Y, DING N S, HUANG R H. Identification of new single nucleotide polymorphisms affecting total number born and candidate genes related to ovulation rate in Chinese Erhualian pigs. Animal Genetics, 2017, 48(1): 48-54.

[6] 李開軍, 侯黎明, 蒲廣, 劉航, 劉根盛, 石傳宗, 金通, 周娟, 李平華, 黃瑞華. 基于全基因組Fst和nSL分析鑒別蘇淮豬中性洗滌纖維表觀消化率相關(guān)候選基因位點. 畜牧獸醫(yī)學(xué)報, 2021, 52(7): 1809-1819.

LI K J, HOU L M, PU G, LIU H, LIU G S, SHI C Z, JIN T, ZHOU J, LI P H, HUANG R H. Identification of candidate gene loci related to apparentdigestibility of suhuai pigs based on genome-wide fst and nSLAnalyses. Acta Veterinaria et Zootechnica Sinica, 2021, 52(7): 1809-1819. (in Chinese)

[7] MOSCATELLI G, DALL'OLIO S, BOVO S, SCHIAVO G, KAZEMI H, RIBANI A, ZAMBONELLI P, TINARELLI S, GALLO M, BERTOLINI F, FONTANESI L. Genome-wide association studies for the number of teats and teat asymmetry patterns in Large White pigs. Animal Genetics, 2020, 51(4): 595-600.

[8] LI Y, PU L, SHI L Y, GAO H D, ZHANG P F, WANG L X, ZHAO F P. Revealing new candidate genes for teat number relevant traits in duroc pigs using genome-wide association studies. Animals: an Open Access Journal from MDPI, 2021, 11(3): 806.

[9] TANG J H, ZHANG Z Y, YANG B, GUO Y M, AI H S, LONG Y, SU Y, CUI L L, ZHOU L Y, WANG X P, ZHANG H, WANG C B, REN J, HUANG L S, DING N S. Identification of loci affecting teat number by genome-wide association studies on three pig populations. Asian-Australasian Journal of Animal Sciences, 2017, 30(1): 1-7.

[10] ZHOU L, ZHAO W, FU Y, FANG X, REN S, REN J. Genome-wide detection of genetic loci and candidate genes for teat number and body conformation traits at birth in Chinese Sushan pigs. Animal Genetics, 2019, 50(6): 753-756.

[11] WANG, LI, ZHOU, GAO, LIU, LI, NIU, ZHANG, LI, ZHOU, HUANG. Association of twelve candidate gene polymorphisms with the intramuscular fat content and average backfat thickness of Chinese suhuai pigs. Animals, 2019, 9(11): 858.

[12] PURCELL S, NEALE B, TODD-BROWN K, THOMAS L, FERREIRA M A R, BENDER D, MALLER J, SKLAR P, DE BAKKER P I W, DALY M J, SHAM P C. PLINK: a tool set for whole-genome association and population-based linkage analyses. American Journal of Human Genetics, 2007, 81(3): 559-575.

[13] MADSEN P, JENSEN J. A user’s guide to DMU. A package for analysing multivariate mixed models version, 2013, 6: 1-33.

[14] VANRADEN P M. Efficient methods to compute genomic predictions. Journal of Dairy Science, 2008, 91(11): 4414-4423.

[15] WEIR B S, COCKERHAM C C. Estimating f-statistics for the analysis of population structure. Evolution; International Journal of Organic Evolution, 1984, 38(6): 1358-1370.

[16] DANECEK P, AUTON A, ABECASIS G, ALBERS C A, BANKS E, DEPRISTO M A, HANDSAKER R E, LUNTER G, MARTH G T, SHERRY S T, MCVEAN G, DURBIN R, GENOMES PROJECT ANALYSIS GROUP 1 0 0 0. The variant call format and VCFtools. Bioinformatics (Oxford, England), 2011, 27(15): 2156-2158.

[17] SPEED D, HEMANI G, JOHNSON M R, BALDING D J. Improved heritability estimation from genome-wide SNPs. The American Journal of Human Genetics, 2012, 91(6): 1011-1021.

[18] YANG Q, CUI J, CHAZARO I, CUPPLES L A, DEMISSIE S. Power and type I error rate of false discovery rate approaches in genome- wide association studies. BMC Genetics, 2005, 6(Suppl. 1): S134.

[19] YE S P, CHEN Z T, ZHENG R R, DIAO S Q, TENG J Y, YUAN X L, ZHANG H, CHEN Z M, ZHANG X Q, LI J Q, ZHANG Z. New insights from imputed whole-genome sequence-based genome-wide association analysis and transcriptome analysis: the genetic mechanisms underlying residual feed intake in chickens. Frontiers in Genetics, 2020, 11: 243.

[20] SCHMID M, MAUSHAMMER M, PREU? S, BENNEWITZ J. Mapping QTL for production traits in segregating Piétrain pig populations using genome-wide association study results of F2crosses. Animal Genetics, 2018, 49(4): 317-320.

[21] YIN L L, ZHANG H H, TANG Z S, XU J Y, YIN D, ZHANG Z W, YUAN X H, ZHU M J, ZHAO S H, LI X Y, LIU X L. rMVP: a memory-efficient, visualization-enhanced, and parallel-accelerated tool for genome-wide association study. Genomics, Proteomics & Bioinformatics, 2021, 19(4): 619-628.

[22] BARRETT J C, FRY B, MALLER J, DALY M J. Haploview: analysis and visualization of LD and haplotype maps. Bioinformatics, 2004, 21(2): 263-265.

[23] ZHANG T Y, GAO H D, SAHANA G, ZAN Y J, FAN H Y, LIU J X, SHI L Y, WANG H W, DU L X, WANG L X, ZHAO F P. Genome-wide association studies revealed candidate genes for tail fat deposition and body size in the Hulun Buir sheep. Journal of Animal Breeding and Genetics = Zeitschrift Fur Tierzuchtung Und Zuchtungsbiologie, 2019, 136(5): 362-370.

[24] 周李生, 趙為民, 涂楓, 吳云鶴, 任守文, 方曉敏. 豬乳頭性狀生理學(xué)和遺傳學(xué)研究進展. 遺傳, 2019, 41(5): 384-390.

ZHOU L S, ZHAO W M, TU F, WU Y H, REN S W, FANG X M. Physiology and genetics research progress of teat traits in pigs. Hereditas, 2019, 41(5): 384-390. (in Chinese)

[25] SON M V, LOPES M S, MARTELL H J, DERKS M F L, GANGSEI L E, KONGSRO J, WASS M N, GRINDFLEK E H, HARLIZIUS B. A QTL for number of teats shows breed specific effects on number of vertebrae in pigs: bridging the gap between molecular and quantitative genetics. Frontiers in Genetics, 2019, 10: 272.

[26] TAN C, WU Z F, REN J L, HUANG Z L, LIU D W, HE X Y, PRAKAPENKA D, ZHANG R, LI N, DA Y, HU X X. Genome-wide association study and accuracy of genomic prediction for teat number in Duroc pigs using genotyping-by-sequencing. Genetics, Selection, Evolution: GSE, 2017, 49(1): 35.

[27] 李晶, 王杰, 康慧敏, 劉冉冉, 李華, 趙桂蘋. 基于BLUP和GBLUP方法估計北京油雞胴體和肉質(zhì)性狀遺傳參數(shù)的差異. 畜牧獸醫(yī)學(xué)報, 2020, 51(1): 35-42.

LI J, WANG J, KANG H M, LIU R R, LI H, ZHAO G P. The difference of genetic parameters for carcass and meat quality traits by BLUP and GBLUP methods in Beijing You chicken. Acta Veterinaria et Zootechnica Sinica, 2020, 51(1): 35-42. (in Chinese)

[28] 朱墨, 鄭麥青, 崔煥先, 趙桂蘋, 劉楊. 基于GBLUP和BayesB方法對肉雞屠宰性狀基因組預(yù)測準(zhǔn)確性的比較. 中國農(nóng)業(yè)科學(xué), 2021, 54(23): 5125-5131. doi: 10.3864/j.issn.0578-1752.2021.23.016.

ZHU M, ZHENG M Q, CUI H X, ZHAO G P, LIU Y. Comparison of genomic prediction accuracy for meat type chicken carcass traits based on GBLUP and BayesB method. Scientia Agricultura Sinica, 2021, 54(23): 5125-5131. doi: 10.3864/j.issn.0578-1752.2021.23.016. (in Chinese)

[29] ZHUANG Z W, DING R R, PENG L L, WU J, YE Y, ZHOU S P, WANG X W, QUAN J P, ZHENG E Q, CAI G Y, HUANG W, YANG J, WU Z F. Genome-wide association analyses identify known and novel loci for teat number in Duroc pigs using single-locus and multi-locus models. BMC Genomics, 2020, 21: 344.

[30] BIDANEL J P, ROSENDO A, IANNUCCELLI N, RIQUET J, GILBERT H, CARITEZ J C, BILLON Y, AMIGUES Y, PRUNIER A, MILAN D. Detection of quantitative trait loci for teat number and female reproductive traits in Meishan × Large White F2pigs. Animal, 2008, 2(6): 813-820.

[31] PIPREK R P, KOLASA M, PODKOWA D, KLOC M, KUBIAK J Z. Cell adhesion molecules expression pattern indicates that somatic cells arbitrate gonadal sex of differentiating bipotential fetal mouse gonad. Mechanisms of Development, 2017, 147: 17-27.

[32] HALBLEIB J M, NELSON W J. Cadherins in development: cell adhesion, sorting, and tissue morphogenesis. Genes & Development, 2006, 20(23): 3199-3214.

[33] ROBINSON G W. Cooperation of signalling pathways in embryonic mammary gland development. Nature Reviews Genetics, 2007, 8(12): 963-972.

[34] PIPREK R P, KLOC M, MIZIA P, KUBIAK J Z. The central role of cadherins in gonad development, reproduction, and fertility. International Journal of Molecular Sciences, 2020, 21(21): 8264.

[35] ASHAIE M A, CHOWDHURY E H. Cadherins: the superfamily critically involved in breast cancer. Current Pharmaceutical Design, 2016, 22(5): 616-638.

Identifying Quantitative Trait Loci Associated with Teat Number of Pig by Genomic Analysis

1Institute of Swine Science, Nanjing Agricultural University, Nanjing 210095;2Huaian Academy, Nanjing Agricultural University, Huai’an 223001, Jiangsu;3Huaiyin Xinhuai Pig Breeding Farm of Huai’an City, Huai’an 223322, Jiangsu

【Objective】The purposes of this study were to analyze the variation of teat number, to explore the quantitative trait locus (QTL) and candidate genes related to teat number, and to provide important molecular markers for the breeding of pig teat number.【Method】This study accurately measured left, right, total teat number of 709 Suhuai pigs (335 fattening pigs and 374 breeding pigs). Fattening pigs were selected for 80K chip genotyping and the heritability and genomic estimated breeding value (GEBV) of left, right and total teat number were calculated by chip data. Based on the rank of GEBV and phenotype of teat number, the top 10% individuals and the bottom 10% individuals were selected for Fixation Index (F) analysis to detect highly differentiated loci. Then, the loci associated with teat number were identified by genome wide association analysis (GWAS) and loci which were highly differentiated and significantly associated with teat number were selected as candidate loci.Genes located near candidate loci and related to teat number after functional annotation were selected as candidate genes.Finally,the association analyses between the most significant candidate loci on each chromosome and teat number of 709 Suhuai pigs were performed to verify the significance of the above loci.【Result】The variation coefficients of left, right and total teat number of Suhuai fattening pigs were 10.20%, 9.26% and 8.50%, respectively, and the heritability were 0.212, 0.257 and 0.312, respectively. Based onFand GWAS analyses, a total of 20 candidate loci onchromosomes (SSC) 7, 13, 16 and 18 for teat number were identified and these candidate loci could explain 5.49%-8.03% of the phenotypic variance. Among them, locus rs80894106 on SSC7 associated with total teat number was consistent with the reported candidate locus of total teat number based on Large white and Duroc pig populations, but candidate loci rs81444134 (26.51 Mb, SSC13) and rs81233299 (8.13 Mb, SSC18) of left teat number were newly discovered loci related to teat number. Interestingly, candidate loci of left, right and total teat number were mainly concentrated in the 6.36-10.66 Mb interval on SSC16;Linkage disequilibrium (LD) analysis found that candidateloci in 7.47-8.27 Mb interval fit into a 795 kb haplotype block, andthis haplotype block was a newly discovered candidate area that affected teat number;rs337606862 (7.47 Mb) in the haplotype block was the most significantly SNP associated with the left and total teat number, and three loci in the haplotype block were all located on the intron of cadherin 18 () gene;gene encoded type II cadherin, and cadherin was related to the identification, sorting, proliferation, apoptosis of cells in developing tissue and the occurrence of breast cancer. Thus,might be a new candidate gene that affected pig teat number.In addition, the most significant loci rs81444134, rs80894106, rs337606862 and rs81233299 on 4 chromosomes were genotyped in 709 Suhuai pigs in this study. After association analysis, these loci were significantly associated with teat number, and could be used as potential molecular markers for the selection of teat number.【Conclusion】In this study, 20 loci significantly related to teat number were identified in Suhuai pig population by genome analysis. Among them, 26.51 Mb on SSC13 and 8.13 Mb on SSC18 were new candidate QTLs for teat number. The 7.47-8.27 Mb on SSC16 was also a newly discovered candidate QTL for teat number, andgene in this interval might be a new candidate gene that affected the formation of pig teat.

F; GWAS; QTL; Suhuai pig; teat number; candidate genes

10.3864/j.issn.0578-1752.2023.10.014

2021-12-17；

2022-01-19

國家自然科學(xué)基金-河南聯(lián)合基金項目（U1904115）、江蘇現(xiàn)代農(nóng)業(yè)（生豬）產(chǎn)業(yè)技術(shù)體系建設(shè)項目（JATS[2020]399）、農(nóng)業(yè)農(nóng)村部農(nóng)業(yè)品種改良提升專項（19210387）、江蘇省農(nóng)業(yè)新品種創(chuàng)制重點項目（PZCZ201732）

尹彥鎮(zhèn)，Tel：13851711578；E-mail：2019105013@njau.edu.cn。通信作者李平華，E-mail：lipinghua718@njau.edu.cn。通信作者黃瑞華，E-mail：rhhuang@njau.edu.cn

（責(zé)任編輯 林鑒非）