華西牛胴體性狀GWAS 分析

高 翰，葛 菲，王澤昭，李宏偉，安炳星，李海鵬，徐凌洋，朱 波，蔡文濤，張路培，高 雪，高會江，陳 燕，李俊雅

（中國農(nóng)業(yè)科學院北京畜牧獸醫(yī)研究所，北京 100193）

我國是肉牛生產(chǎn)大國，也是牛肉消費大國[1]。隨著居民生活水平的提高，消費者對牛肉的需求持續(xù)增長，供需缺口呈擴大趨勢[2]。胴體性狀作為肉牛生產(chǎn)實踐中最重要的經(jīng)濟性狀之一，是典型的由多基因調(diào)控的數(shù)量性狀，遺傳調(diào)控機制復雜[3-4]。Risch 于1996 年首次提出全基因組關聯(lián)分析（Genome-Wide Association Study，GWAS）的應用[5]，隨著測序技術的不斷發(fā)展和測序成本的降低，GWAS 被廣泛應用于研究家畜數(shù)量性狀的研究中[6-10]。利用GWAS，在肉牛中鑒定出了PLAG1[11]、IGF2[12]、MTPN[13]等影響胴體性狀的重要基因。大部分GWAS 都是基于單一性狀與SNPs 的關聯(lián)分析，即使同時分析多個相關表型，也只是逐個性狀分別進行單一性狀GWAS[14-15]。Bolormaa 等[16]研究指出，對于具有高估計育種值（Estimated Breeding Value，EBV）準確性的性狀，多性狀GWAS 的效力至少會與單性狀GWAS 的效力一樣高。如果性狀間高度相關，則多性狀分析相較于單性狀分析更具有優(yōu)勢。因為多性狀GWAS 只進行一次統(tǒng)計檢驗，同時考慮多個性狀的性狀內(nèi)和性狀間方差組分[17]，降低了多重檢驗造成的誤差[18]，提高了檢驗功效[19-20]和參數(shù)估計的精度[21]。因此，與逐個性狀與遺傳位點間的關聯(lián)分析相比，多個性狀聯(lián)合信息的多性狀GWAS 通常具有更好的關聯(lián)分析效果[22-24]。

“華西?！笔侵袊r(nóng)業(yè)科學院北京畜牧獸醫(yī)研究所歷經(jīng)40 余年潛心培育的肉牛新品種，2021 年通過國家畜禽遺傳資源委員會審定，并獲得國家畜禽新品種證書，其具有生長速度快、飼料轉化率高、凈肉率高、產(chǎn)肉和繁殖性能好、抗逆性強等特性。與國際同類型肉牛品種相比，“華西?！钡娜赵鲋?、屠宰率、凈肉率處于國際先進水平。本研究以內(nèi)蒙古烏拉蓋管理區(qū)建立的華西牛資源群體為研究對象，基于Illumina BovineHD 芯片基因型數(shù)據(jù)，對胴體重（Carcass Weight，CW）、胴體長（Carcass Length，CL）、胴體胸深（Chest Depth of Carcass，CD）3 個胴體性狀分別進行單性狀與多性狀GWAS 分析，旨在挖掘影響華西牛胴體性狀的顯著SNPs 位點及候選基因，為全基因組選擇育種提供有效的分子標記，為進一步研究華西牛胴體相關性狀的潛在遺傳機制提供有價值的參考。

1 材料與方法

1.1 實驗數(shù)據(jù)收集 實驗動物群體選自中國農(nóng)業(yè)科學院北京畜牧獸醫(yī)研究所牛遺傳育種創(chuàng)新團隊在內(nèi)蒙古錫林郭勒盟烏拉蓋管理區(qū)組建的華西牛資源群體，該群體2008 年開始建立，每年進行擴群，并且使用Illumina BovineHD 芯片對每個個體進行基因分型。胴體重、胴體長和胴體胸深3 個胴體性狀表型數(shù)據(jù)嚴格按照GB/T 27643-2011《牛胴體及鮮肉分割》[25]的要求進行測定。

1.2 表型處理和基因型質(zhì)控 對表型數(shù)據(jù)進行處理，包括剔除表型缺失值、三倍標準差以外的異常值。采用一般線性模型評估育肥場、屠宰場、屠宰年份、屠宰季節(jié)以及年齡對性狀影響的顯著性，將顯著的效應作為協(xié)變量加入到模型當中。利用PLINK 軟件[26]對群體進行主成分分析（Principal Components Analysis，PCA），使用R 包ggplot2 繪制群體PCA 圖。利用EIGENSTRAT計算各主成分是否具有顯著統(tǒng)計學意義，將達到顯著水平P＜0.05 的主成分作為協(xié)變量加入到模型當中。使用PLINK 軟件對相應個體的基因型數(shù)據(jù)進行質(zhì)量控制，標準如下：剔除SNP 缺失率大于10%的個體以及檢出率低于90%、最小哈德溫伯格平衡小于1.0×10-6、最小等位基因頻率小于0.05 的SNPs 位點。最后，使用R包ASReml[27]中的多性狀動物模型計算各性狀之間的表型相關和遺傳相關。

1.3 關聯(lián)統(tǒng)計分析

1.3.1 單性狀GWAS 使用GEMMA 軟件[21]進行全基因組關聯(lián)分析，采用考慮親緣關系的混合線性模型（Linear Mixed Model，LMM），計算公式如下：

式（1）中，y為所有個體表型性狀的n×1 向量；W為協(xié)變量矩陣，α為包括截距在內(nèi)對應系數(shù)的一列向量；x為標記基因型向量，β表示標記位點效應的大??；u為個體的隨機效應向量，u～N（0，KVg），K代表由SNP 標記計算出的已知的n×n遺傳關系矩陣，Vg為加性遺傳方差；ε代表n×1 隨機誤差向量，符合ε～N（0，IVe）分布，I代表n×n身份矩陣，Ve代表殘差組分。本次研究中主要利用Wald 測驗作為GWAS 分析顯著統(tǒng)計指標。

1.3.2 多性狀GWAS 使用GEMMA 軟件[21]中的多性狀混合線性模型（Multivariate Linear Mixed Model，mvLMM）進行多個性狀表型值和SNP 標記的GWAS，計算公式如下：

式（2）中，Y為n×t維表型矩陣，n是樣本數(shù)目，t是分析性狀個數(shù)；W為n×c維協(xié)變量矩陣，A為c維相應系數(shù)行向量，c是包含截距項在內(nèi)的行變量個數(shù)；x為n維當前檢驗SNP 的基因型指示變量列向量，β為n×t維不包括當前檢驗SNP 的剩余多基因效應矩陣；G～MN（0，Vg?K），K是由全基因組SNP 標記構建的基因組親緣關系矩陣，Vg是t×t維剩余多基因方差-協(xié)方差矩陣；E為n×t維誤差項矩陣，E～MN（0，Ve?I），Ve是t×t維誤差方差-協(xié)方差矩陣，I是單位矩陣。最終，全基因組上每個SNP 都能得到1 個一因多效Wald 統(tǒng)計量。應用Bonferroni 方法進行多重校正，確定染色體水平顯著性閾值為1/N，N 為用于分析的SNPs 數(shù)目。最后，用R 包ggplot2 繪制曼哈頓圖和QQ 圖。

1.4 位點版本轉化及基因注釋 Illumina BovineHD 芯片中的SNPs 位置信息為Bos taurus UMD 3.1 版本，使用在線網(wǎng)站UCSC 中的liftover 工具（http://genome.ucsc.edu/cgi-bin/hgLiftOver/）將其轉化成最新的?；蚪M版本ARS-UCD 1.2。根據(jù)得到的顯著SNPs 位點的物理位置，利用Ensemble 在線數(shù)據(jù)庫的biomart 模塊（http://www.ensembl.org/biomart/）搜索每個顯著SNP位點上下游50 kb 區(qū)域尋找距離最近的候選基因，作為與目標性狀關聯(lián)的重要功能基因。

2 結果

2.1 主成分分析 由圖1 可知，實驗群體存在分層現(xiàn)象。利用EIGENSTRAT 計算發(fā)現(xiàn)，有2 個達到P＜0.05 顯著水平的主成分，因此在后續(xù)分析中將前2 個主成分作為協(xié)變量加入到模型中。

圖1 群體主成分分析圖

2.2 表型及基因型描述性統(tǒng)計

2.2.1 表型值的描述性統(tǒng)計 分別對840 頭華西牛公牛的CW、CL 和CD 3 個性狀開展單性狀GWAS 與多性狀GWAS 分析。表1 為各性狀的原始表型描述性統(tǒng)計分析結果。3 個性狀的表型頻率分布圖見圖2，可以看出，各性狀的表型值基本都符合正態(tài)分布。此外，對各表型值進行了表型相關分析與遺傳相關分析（表2）。結果顯示，胴體重、胴體長、胴體胸深3 個性狀間均具有較高的表型相關（r＞0.6）；胴體長與胴體胸深具有中等遺傳相關性（r=0.339 2），胴體重與胴體長、胴體胸深具有強遺傳相關性（r＞0.4）。

表1 胴體重、胴體長、胴體胸深的描述性統(tǒng)計

表2 性狀間表型及遺傳相關程度

圖2 CW、CL 和CD 的頻率分布

2.2.2 基因型的描述性統(tǒng)計 按篩選原則對基因型數(shù)據(jù)進行質(zhì)控后，共剩余607 198 個SNPs 標記。由圖3 可知，各染色體上的SNP 位點分布較均勻，可以用于華西牛群體的GWAS 分析。

圖3 SNPs 標記在29 條染色體的分布情況

2.3 GWAS 分析 使用GEMMA 軟件中考慮親緣關系的混合線性模型（LMM）對3 個性狀分別做單性狀GWAS 分析。曼哈頓和QQ 圖如圖4 所示，在CW 和CL 性狀中都沒有檢測到顯著的遺傳變異位點，僅在CD 性狀中發(fā)現(xiàn)了6 個顯著的SNPs 位點，分別位于4、7、10、23、24、25 號染色體上，共注釋到了5 個基因（表3）。

表3 單性狀GWAS 和多性狀GWAS 基因注釋結果

圖4 單性狀GWAS 曼哈頓和QQ 圖

使用mvLMM 分別進行了CW-CL、CW-CD、CLCD、CW-CL-CD 4 個組合的多性狀GWAS 分析（圖5）。在CW-CL 組合的多性狀GWAS 中，發(fā)現(xiàn)1 個超過染色體顯著水平的SNP 位點，位于1 號染色體上，注釋到了1 個基因；在CW-CD 組合的多性狀GWAS 中，發(fā)現(xiàn)5 個超過染色體顯著水平的SNPs 位點，分別位于4、7、23、25 號染色體上，共注釋到3 個基因；在CL-CD組合的多性狀GWAS 中，發(fā)現(xiàn)2 個超過染色體顯著水平的SNPs 位點，分別位于4、25 號染色體上，注釋到了1 個基因；在CW-CL-CD 組合的多性狀GWAS 中，發(fā)現(xiàn)3 個超過染色體顯著水平的SNPs 位點，分別位于4、7、23 號染色體上，共注釋到3 個基因。合并多性狀GWAS 的4 個組合的結果發(fā)現(xiàn)，有4 個SNPs 位點在多性狀GWAS 中被重復檢測到，其中Hapmap36353-SCAFFOLD29708_3468 在CW-CD、CL-CD、CW-CLCD 3 個組合中被重復檢測到，BovineHD0700018227和BovineHD2300004986均在CW-CD、CW-CL-CD 2個組合中被重復檢測到，BovineHD2500006960 在CLCD、CW-CD 2 個組合中被重復檢測到。去除重復位點，多性狀GWAS 共找到6 個一因多效的SNPs 位點。

圖5 多性狀GWAS 曼哈頓和QQ 圖

2.4 重要候選基因 對CW、CL、CD 3 個性狀分別進行單性狀GWAS 和不同性狀之間組合的多性狀GWAS結果進行匯總統(tǒng)計。結果得到，單性狀與多性狀GWAS共檢測到9 個與華西牛胴體性狀顯著相關的SNPs 位點，根據(jù)顯著SNPs 位點的物理位置，在這9 個顯著SNPs 位點上下游50 kb 的區(qū)間內(nèi)尋找距離最近的候選基因，共在8 個位點附近找到6 個候選基因，分別是突觸結合蛋白5L（Syntaxin Binding Protein 5 Like，STXBP5L）、磷酸二酯酶1C（Phosphodiesterase 1C，PDE1C）、過氧化物酶體增殖物激活受體γ輔激活因子-1β（Peroxisome Proliferator-Activated Receptor Gamma Coactivator 1-Beta，PPARGC1B）、（Metaxin 3，MTX3）、鈣調(diào)磷酸酶調(diào)節(jié)因子2（Regulator Of Calcineurin 2，RCAN2）、?？康鞍?（Docking Protein 6，DOK6）。其中，PDE1C、PPARGC1B、RCAN23 個基因在單性狀GWAS 與多性狀GWAS 中被重復檢測到。

由圖6 可知，單性狀與多性狀GWAS 重復檢測到的3 個顯著SNPs 位點分別是Hapmap36353-SCAFFOLD29708_3468、BovineHD2300004986 和Bovine HD2500006960。其中，單性狀GWAS 與多性狀GWAS中CW-CD、CL-CD 組合共同檢測到的BovineHD25000 06960 未注釋到相關基因。此外，存在不同策略檢測到的多個SNPs 位點注釋到同一個基因的情況，具體為：在CD 單性狀GWAS 檢測到的BovineHD0700018210 與CW-CD、CW-CL-CD 2 個多性狀GWAS 共同檢測到的BovineHD0700018227，同時注釋到基因PPARGC1B；CW-CD 多性狀GWAS 檢測到的BovineHD2300004985與單性狀GWAS 和CW-CD、CW-CL-CD 2 個多性狀GWAS 共同檢測到的BovineHD2300004986，同時注釋到基因RCAN2。PPARGC1B和RCAN2將作為后續(xù)重點關注的與胴體性狀關聯(lián)的重要候選基因。

圖6 單性狀GWAS 與多性狀GWAS 結果韋恩圖

3 討 論

在肉牛產(chǎn)業(yè)中，胴體性狀被認為是影響肉牛生產(chǎn)的重要經(jīng)濟性狀，對于提高生產(chǎn)效率和效益起著至關重要的作用[28-29]。然而，肉牛胴體性狀的測定只能在個體屠宰后進行，這就導致胴體性狀的相關表型數(shù)據(jù)收集困難且成本高昂[30]。因此，闡明胴體性狀的遺傳調(diào)控機制，提供可靠的分子遺傳標記，可以對肉牛胴體相關性狀進行早期遺傳評估，從而降低育種成本，提高生產(chǎn)效益。

GWAS 方法已經(jīng)廣泛應用于動植物復雜性狀的遺傳解析[31-33]。本研究基于Illumina BovineHD 芯片基因分型，對華西牛CW、CL、CD 3 個胴體性狀進行了關聯(lián)分析，綜合考慮單性狀和多性狀GWAS 分析得到的結果，共得到9 個顯著的SNPs 位點和6 個相關基因，分別是STXBP5L、PDE1C、PPARGC1B、MTX3、RCAN2和DOK6。其中，多性狀與單性狀策略共同鑒定到3 個基因，分別為PDE1C、PPARGC1B和RCAN2。PDE1C基因在小鼠胚胎早期發(fā)育過程中中樞神經(jīng)系統(tǒng)內(nèi)遷移的神經(jīng)元細胞中表達[34]。在人類疾病的研究中，PDE1C基因的異常會導致兒童發(fā)育遲緩[35]。PPARGC1B基因是一種調(diào)控基因轉錄的共激活劑，對轉錄因子和核受體活動、脂肪氧化、糖酵解以及能量代謝等多種生物學過程具有調(diào)控作用[36-39]。Jiang 等[40]研究了秦川牛犢牛和成年發(fā)育階段牛脂肪組織中circRNA 的表達，發(fā)現(xiàn)circFUT10 在脂肪細胞增殖和分化中有重要作用，circFUT10 間接促進PPARGC1B的轉錄，然后調(diào)節(jié)細胞的增殖和分化，且PPARGC1B在成年牛脂肪組織中表達水平較高，說明該基因對胴體性狀的表型差異可能存在重要影響。在小鼠模型中的相關研究表明，RCAN2基因對體重調(diào)節(jié)有重要作用[41-43]。RCAN2在成纖維細胞、心臟、腦、肝臟和骨骼肌中表達，能夠抑制鈣調(diào)磷酸酶活性，在骨骼發(fā)育和維護中有重要作用[44]。由此可見，該基因可能是影響華西牛胴體性狀的重要候選基因。MTX3與DOK6是單性狀GWAS 分析中鑒定到的與CD 性狀相關的2 個特有基因。MTX3在蛋白質(zhì)轉運到線粒體這一生物過程中發(fā)揮重要作用[45-46]。DOK6能促進RET 介導的神經(jīng)突生長，可能在大腦發(fā)育或維護中發(fā)揮作用[47]。Jiao 等[48]對杜洛克豬群體進行GWAS 分析，提出DOK6基因為平均日采食量、育肥期日增重、活體背膘厚3 個重要經(jīng)濟性狀最有可能的候選基因之一。STXBP5L基因為CW-CL 組合的多性狀GWAS 中單獨鑒定到的基因，它在維持葡萄糖穩(wěn)態(tài)[49]、囊泡運輸和胞吐作用抑制[50]等過程中發(fā)揮著重要作用。

對單性狀GWAS 與多性狀GWAS 結果進行比較，結果發(fā)現(xiàn)，在相同的顯著性閾值下，兩者檢測到的基因數(shù)目無明顯差別。但是單一性狀GWAS 分析中，CW、CL 性狀都沒有檢測到顯著SNPs 位點及相關基因。通過將CW 與CL、CD 兩兩組合進行多性狀GWAS，在CW-CL、CW-CD、CL-CD、CW-CL-CD 4 個不同組合中均檢測到了不同數(shù)目的顯著SNPs 位點和相關基因。同時，在單性狀GWAS 以及4 組多性狀GWAS 中，Hapmap36353-SCAFFOLD29708_3468、BovineHD070 0018227、BovineHD2300004986、BovineHD2500006960共4 個SNPs 位點，以及PDE1C、PPARGC1B、RCAN2共3 個基因被重復檢測到，推測認為這些位點及基因極大可能是潛在的一因多效位點和一因多效基因?？傮w而言，多性狀GWAS 可以提高對位點的檢測效率和準確性，這與高進等[51]、Bolormaa 等[16]的研究結果相一致。在實際育種當中，CW、CL 和CD 的表型通常是同時被選擇的，由于多性狀GWAS 分析檢測到的SNPs 位點具有一因多效的優(yōu)越性，因此綜合考慮多性狀GWAS 分析在生物學和實際應用層面更有意義。

4 結 論

本研究基于Illumina BovineHD 高密度基因分型結果，利用840 頭華西牛胴體重、胴體長和胴體胸深表型數(shù)據(jù)，采用單變量和多變量線性混合模型方法進行單性狀和多性狀的GWAS 分析，共檢測到9 個顯著關聯(lián)的SNPs 位點，注釋到STXBP5L、PDE1C、PPARGC1B、MTX3、RCAN2、DOK6共6個候選基因，其中PPARGC1B、RCAN2被不同方法檢測到的多個SNPs 位點同時注釋到，且參與體重調(diào)節(jié)、脂肪細胞增殖和分化等過程，推測其可能是影響胴體性狀的重要候選基因。