太湖流域地方豬種遺傳結構評估和體重體尺性狀差異相關選擇位點鑒定

1 材料與方法

1.1 試驗動物

采集440頭二花臉豬(EHL)，楓涇豬(FJ)，米豬(MI)，嘉興黑豬(JXB)，沙烏頭豬(SWT)，紅燈籠豬(HDL)，中型梅山豬(MMS)和小型梅山豬(SMS)的耳組織樣。二花臉豬(n=111)分別從常州市，常熟市和蘇州市的保種場采集，小型梅山豬(n=71)分別從太倉市和蘇州市的保種場采集，中型梅山豬(n=65)從昆山市保種場采集，楓涇豬(n=23)從蘇州市的保種場采集，米豬(n=59)從常州市的保種場采集，嘉興黑豬(n=40)從嘉興市的保種場采集，沙烏頭豬(n=41)從南通市的保種場采集，紅燈籠豬(n=30)從溧陽市的保種場采集。這440頭個體可以涵蓋江蘇省太湖流域地方豬保種群的全部血統。采樣的所有公豬和母豬在3代內沒有共同血親。

1.2 個體基因組信息鑒定和收集

將440頭豬耳組織樣本用常規的苯酚/氯仿方法提取DNA，并稀釋至最終濃度50 ng/μL，然后采用北京康普森公司研發的豬50 K基因型芯片(測序平臺為美國Illumina公司提供)進行基因分型，芯片總計含有51 315個SNP位點。基因型輸入文件轉換為PLINK v1.07格式文件[11]用于后續分析。

對PLINK格式文件進行SNP位點質控，具體條件為：單個個體SNP位點檢出率大于90%，每個SNP檢出率大于90%，微小等位基因頻率大于5%。所有未知位置信息和性染色體上的SNP位點都剔除。經過質控，獲得總計29 754個滿足條件的SNP位點。此外，采用上述SNP位點質控原則，在65頭中型梅山豬和71頭小型梅山豬群體內部進行質控，經過質控得到總計21 850個有效SNP位點；在111頭二花臉和59頭米豬個體進行質控并得到29 754個有效SNP位點。

1.3 群體間遺傳距離(1-Dst)計算

為了估計個體間的遺傳距離，使用PLINK v1.07軟件計算了等位基因共享比例(Dst)[11]。所有個體之間的遺傳距離計算公式為(1-Dst)[12]。基于(1-Dst)矩陣，使用MEGA v6.0軟件[13]計算兩兩群體之間平均遺傳距離的大小。使用群體間遺傳距離(1-Dst)矩陣構建聚類樹熱圖，通過內部R腳本進行可視化。

1.4 群體間遺傳分化系數(Fst)計算

使用Wright’s F-derived統計法來計算品種間遺傳分化系數Fst[14]，使用GENEPOP v4.2軟件[15]進行相關計算。不同的Fst值代表的意義為：Fst在0～0.05表示較低的遺傳分化，Fst在0.05～0.15表示中等的遺傳分化，Fst在0.15～0.25表示較高的遺傳分化，Fst高于0.25表示高度遺傳分化[16-17]。使用群體間遺傳分化系數(Fst)矩陣構建了聚類樹熱圖，并通過內部R腳本進行可視化。

1.5 群體結構分析

基于SNP基因型數據，使用ADMIXTURE v1.23軟件[18]推斷出最可能的祖先血統組成(K)，即假定有K個祖先存在的前提下，每個群體內含有的各祖先血統的比例，根據最小交叉驗證誤差(CV-error)的原則確定最大K值。對于太湖流域地方豬，使用內部R腳本將K=2～8的結果進行可視化。為了增加ADMIXTURE分析的準確性，并避免由這些群體的不同樣本量引起的潛在偏差，我們從每個群體中隨機選擇了23頭個體進行ADMIXTURE分析。

1.6 受選擇SNP位點鑒定

使用GENEPOP v4.2軟件[16]進行中型梅山豬和小型梅山豬群體間Fst分析，每個SNP位點計算得到相應遺傳分化(Fst)值。對于所有SNP位點的Fst值，采用Z檢驗，鑒定極顯著(P<0.01)偏離整體Fst均值的SNP位點[19]，并將這些SNP位點作為在梅山豬亞群間受選擇的候選SNP位點。采用相同計算原理鑒定二花臉和米豬群體間受選擇的候選SNP位點。最終，將梅山豬亞群間候選SNP位點與二花臉和米豬群體間候選SNP位點取交集，作為更有可能的受選擇位點。

1.7 候選基因注釋

對SNP位點上下游各50 kb[20]的區域內進行基因注釋，使用豬11.1版本參考基因組(Sscrofa 11.1)，注釋網站為ensemble官網(http://www.ensembl.org)。

1.8 受選擇SNP位點間連鎖不平衡分析

在人為或自然選擇的壓力下，優勢等位基因型頻率會不斷上升，并且由于“搭車效應”形成選擇性清除區域，選擇性清除區域附近的SNP位點連鎖不平衡程度加深[21]。為了進一步確定候選SNP位點受選擇的可能性，我們使用PLINK v1.07軟件[13]來檢測位于同一染色體且位置距離接近的多個受選擇SNP位點間的連鎖不平衡情況(r2)。執行命令“Plink-file-r2”，通過SNP位點間的相關系數r2值來估算連鎖不平衡(LD)程度，并以r2值0.3作為閾值來衡量位點間是否存在連鎖不平衡[12,22]。

1.9 受選擇SNP位點與豬體重體尺等性狀QTLs注釋

對于1、3和6號染色體受到選擇的SNP位點形成的受選擇區域(相應染色體彼此位置距離接近且強連鎖的多個受選擇SNP位點，以位置距離最小的SNP位點作為起點，位置距離最大的SNP位點作為終點，形成的受選擇區域)進行豬體重體尺等相關性狀QTLs注釋，注釋網站為PigQTLdb官網(https://www.animalgenome.org/cgi-bin/QTLdb/SS/index)。

2 結果與分析

2.1 太湖流域地方豬群體內部結構

2.1.1 太湖流域地方豬各品種及梅山豬品種內兩個亞群間遺傳關系

由表1可知，在太湖流域地方豬各品種及梅山豬品種內兩個亞群間的遺傳距離(1-Dst)分析中，群體間遺傳距離最為接近的是二花臉和米豬，群體間遺傳距離為0.100，其次較為接近的是中型梅山豬與小型梅山豬，群體間遺傳距離為0.153。二花臉和米豬間的遺傳距離比梅山豬兩個亞群更為接近。

表1 太湖流域8個地方豬群體間遺傳距離(1-Dst)

8個地方豬在群體間遺傳分化方面的分析結果如表2所示。中型梅山豬與小型梅山豬群體間的遺傳分化程度最小，遺傳分化系數為0.264，其次為楓涇和嘉興黑間的遺傳分化系數，數值為0.278，排在第三的是二花臉和米豬間的遺傳分化系數，數值為0.282。所有群體間的遺傳分化系數都超過0.250。

表2 太湖流域8個地方豬群體間遺傳分化系數(Fst)

2.1.2 太湖流域地方豬各品種及梅山豬品種內兩個亞群間ADMIXTURE分析

為了判定太湖流域地方豬各品種及梅山豬品種內兩個亞群間的進化路線，進行了ADMIXTURE分析，結果如圖1所示。圖1展示了K值為2～8的結果，在K值為8時，ADMIXTURE分析擁有最小的交叉檢驗誤差(CV-error)，梅山豬兩個亞群的K值為2～5，這兩個群體都保持著一致的祖先血統組成，K值為6時，分化成兩個不同的群體。二花臉和米豬的K值為2～7時，這兩個群體都保持著一致的祖先血統組成，直到K為8時才最終分化成兩個群體。二花臉和米豬擁有最為一致的進化路線，這與二花臉和米豬間遺傳距離最小的結果一致。此外，楓涇和嘉興黑間的祖先血統組成在K值為2～5時保持一致。在K值為8(最小CV-error值)時，這8個群體都可以區分成不同的群體。

注：EHL代表二花臉豬；MMS代表中型梅山豬；SMS代表小型梅山豬；FJ代表楓涇豬；HDL代表紅燈籠豬；JXB代表嘉興黑豬；MI代表米豬；SWT代表沙烏頭豬；每種顏色代表一個祖先群，群體間均以虛線分隔圖1 太湖流域地方豬8個群體的祖先血統組成

2.2 太湖流域地方豬受選擇SNP位點鑒定

2.2.1 梅山豬亞群間以及二花臉豬和米豬群體間共同分化SNP位點鑒定及基因注釋結果

對梅山豬亞群間、以及二花臉和米豬群體間進行兩組Fst分析，鑒定這兩組共同分化的SNP位點，并作為太湖流域地方豬受選擇的候選SNP位點。將閾值線以上的點作為群體間顯著分化的SNP位點(Z檢驗，P<0.01)，結果如圖2所示。一共發現24個SNP位點在這兩組Fst分析中都顯著分化，其中，有7個位于1號染色體，6個位于3號染色體，6個位于6號染色體，彼此間距離接近，其余分布在2、4、5和16號染色體。

對這24個候選SNP位點進行基因注釋，得到20個可以編碼蛋白的基因，其中有8個基因與體重體尺性狀相關，具體信息見表3。這些基因主要有如下功能：體重體尺(DENND1A、LHX2、NR6A1、FANCL、KCNK3、CDH11、CDH5和AGBL4)，繁殖性狀(STRBP、FANCL、MRGPRF和PSME4)，神經發育(MVB12B、ARHGAP39、CNTN1、PLCG2和CTNND2)，抗病力(NEK6和SCAI)，嗅覺(NR2C1)，運動能力(VPS13D)。

表3 受選擇SNP位點基因注釋結果

2.2.2 受選擇SNP位點間連鎖不平衡分析

針對1、3和6號染色體中受選擇候選SNP位點，多數位置距離彼此接近，因此我們對這3條染色體上的位置距離彼此接近的候選SNP位點進行連鎖不平衡分析，結果見表4。上述3條染色體內的候選SNP位點都存在較強的連鎖不平衡，r2值遠大于0.3，其中1號染色體多個SNP位點達到完全連鎖，r2值為1。并且1號染色體彼此位置距離接近的7個SNP位點，分別在中型梅山豬和二花臉群體中形成了一種主要的單倍型“TCCGAAA”和“TCCGAAA”，同時在小型梅山豬和米豬群體內變異較低；對于3號染色體彼此位置距離接近的5個SNP位點，分別在中型梅山豬，小型梅山豬，二花臉和米豬群體中形成了一種主要的單倍型“GGCAG”，“AATGA”，“GGCAG”和“AATGA”；對于6號染色體彼此位置距離接近的3個SNP位點，分別在中型梅山豬和米豬群體中形成了一種主要的單倍型“GTG”和“ACA”。可見1、3和6號染色體在群體間經歷了不同方向的選擇，導致群體間分化增大，并使得群體內變異程度減少，逐漸形成單一的單倍型。

a.中型梅山和小型梅山群體間Fst分析；b.米豬和二花臉群體間Fst分析；x軸為18條染色體，y軸為每個位點Fst值；不同的染色體用不同顏色標注，紅色實線為閾值線，閾值(Z檢驗，P<0.01)上的點代表受選擇位點，候選基因用紅色圓框標注，與體重體尺相關的候選基因用紅色字體標注圖2 梅山豬亞群間以及二花臉豬和米豬群體間共同分化SNP位點鑒定

表4 豬1號、3號和6號染色體受選擇SNP位點間的連鎖不平衡分析(r2)

2.2.3 受選擇SNP位點與豬體重體尺性狀相關QTLs注釋

對SNP位點所形成的選擇區域進行豬體重體尺等相關性狀QTLs注釋，結果見表5。1、3和6號染色體相應區域注釋到大量的QTLs與豬的體重體尺性狀相關，其中主要包括豬日增重性狀，體重性狀，椎骨數性狀，料重比性狀和采食量性狀，表明這些區域很可能對于豬的體重體尺等表型的形成起到巨大作用。

表5 受選擇SNP位點所在區域QTLs注釋結果

3 討論

3.1 太湖流域地方豬種內部結構研究

通過對太湖流域地方豬各品種，以及梅山豬品種內兩個亞群進行多種群體遺傳學分析發現，梅山豬亞群間遺傳距離接近，遺傳分化程度最低，進化路線相對一致，而二花臉與米豬間也有著最近的遺傳距離，較低的分化程度以及最一致的進化路線，這與他人的研究結果一致[6,24]。試驗發現，二花臉豬與米豬之間的遺傳距離小于梅山豬兩個亞群間的遺傳距離，并且擁有最為一致的祖先血統組成，表明二花臉豬與米豬之間的遺傳關系最為接近，其接近程度甚至超過了梅山豬的兩個亞群間的遺傳關系，該結果與Wang等[6]研究一致。此外還發現，所有太湖流域地方豬群體間的遺傳分化都超過0.25，達到高度分化水平[18-19]，并在ADMIXTURE分析中，K值為8時(最小的CV-error)，各自形成了不同的祖先血統組成情況，可見包括梅山豬兩個亞群在內的8個群體彼此間都存在一定的分化，因此在保種時應分開保種。

3.2 太湖流域地方豬種體重體尺性狀差異的候選SNP位點鑒定

對于梅山豬亞群、二花臉和米豬這兩個群體組合分別進行遺傳分化Fst的分析，并尋找在兩組分析中共同分化的SNP位點作為候選SNP位點。采用這種遺傳結構十分接近、但體重體尺性狀差異顯著的群體進行比對分析，可以減少由于群體間其他性狀差異或者本身遺傳結構差異導致SNP位點顯著分化的可能性。為了減少顯著分化SNP位點是假陽性的概率，通過取兩個分析結果交集的方法來進一步縮小差異SNP位點的范圍。本試驗共發現24個候選SNP位點，它們在這兩組分析中都顯著分化，這些候選SNP位點分布集中，主要分布在1、3、6號染色體區域內，表明1、3、6號染色體在這4個群體中更可能受到強的選擇。

將受選擇SNP位點進行基因注釋，注釋到20個可以編碼蛋白的基因，8個基因被報道與體重體尺等性狀相關。DENND1A基因的rs10818854位點與人的肥胖顯著相關[25]；LHX2基因通過調節破骨細胞中的RANKL信號傳導調節小鼠的骨重塑[26]；NR6A1基因作為豬肋骨數的主效基因，被運用于多個豬種的分子標記輔助選育改良過程中[27-28]；FANCL基因可能是導致人肥胖的候選基因[29]；KCNK3基因可以調節脂肪生熱和肥胖[30]；CDH11基因是調節小鼠骨骼結構和股骨形態的主要候選基因[31]；髖臼唇細胞中的CDH5可能參與了骨關節炎的發病機制[32]；AGBL4基因被發現與牛的生長等數量性狀顯著相關[33]。

對上述3個染色體的候選SNP位點進行連鎖不平衡分析結果發現，每條染色體上候選SNP間的r2值都遠大于0.3，達到強連鎖水平[12,22]，并且這些SNP在各個群體中都展現出較低多態性，形成了主要的單倍型。在人為或自然選擇的壓力下，優勢等位基因型頻率會不斷上升，并且由于連鎖不平衡，使得優勢基因型附近形成大片段的純合單倍型[19]，而1、3、6號染色體內候選SNP位點間有很強的連鎖不平衡，并且候選SNP位點多態性降低，表明這些候選SNP位點所在的區域受到較強的選擇。

將受選擇SNP位點與豬QTLs數據庫進行比對，結果注釋到大量與體重體尺性狀相關的QTLs。Christine等[9]利用梅山豬和皮特蘭雜交產生的F2代群體鑒定到1號染色體226 764～265 324 kp處存在影響日增重的QTL；HUANG等[33]利用杜洛克和二花臉雜交產生的F2代群體鑒定到1號染色體266 307～287 364 kp處存在影響胸椎數的QTL；Guo等[34]利用大白豬和梅山豬雜交產生的F2代群體鑒定到3號染色體4 657～140 466 kp處存在一個大的QTL影響豬體重性狀；利用皮特蘭和梅山豬雜交產生的F2代群體鑒定到6號染色體24 513～24 582 kp處存在一個影響料重比的QTL[10]。這些結果都與本試驗發現的候選SNP所在區域重合。以上證據表明，在兩個Fst分析中共同受選擇，且基因注釋結果與體重和體尺性狀相關的區域，可能與太湖流域地方豬種群體間體重體尺性狀差異有關，這些區域將作為候選區域進行進一步研究。

致謝：常州市焦溪二花臉豬專業合作社及其工作人員顧岳清、黃媛和李順等，常熟市牧工商有限公司、國家級二花臉豬保種中心及其工作人員陸志強等，蘇州蘇太企業有限公司及其工作人員黃雪根等，浙江青蓮食品股份有限公司及其工作人員趙善林等，江蘇省溧陽市乾豐養殖有限公司及其工作人員程躍強等，昆山市梅山豬保種有限公司及其工作人員沈永杰、時建新、周騰冰等，太倉市種豬場及其工作人員石曉峰、施擎天等，江蘇興旺農牧科技發展有限公司、國家級沙烏頭豬保種場及其工作人員唐慧娟等，常州市金壇米豬原種場及其工作人員田國興等，對于豬樣本采集工作予以支持和幫助，在此一并感謝。