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

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

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

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

1 材料與方法

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

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

1.3 關(guān)聯(lián)統(tǒng)計(jì)分析

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

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

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

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

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

2 結(jié)果

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

圖1 群體主成分分析圖

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

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

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

表2 性狀間表型及遺傳相關(guān)程度

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

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

圖3 SNPs 標(biāo)記在29 條染色體的分布情況

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

表3 單性狀GWAS 和多性狀GWAS 基因注釋結(jié)果

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

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

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

2.4 重要候選基因 對(duì)CW、CL、CD 3 個(gè)性狀分別進(jìn)行單性狀GWAS 和不同性狀之間組合的多性狀GWAS結(jié)果進(jìn)行匯總統(tǒng)計(jì)。結(jié)果得到，單性狀與多性狀GWAS共檢測(cè)到9 個(gè)與華西牛胴體性狀顯著相關(guān)的SNPs 位點(diǎn)，根據(jù)顯著SNPs 位點(diǎn)的物理位置，在這9 個(gè)顯著SNPs 位點(diǎn)上下游50 kb 的區(qū)間內(nèi)尋找距離最近的候選基因，共在8 個(gè)位點(diǎn)附近找到6 個(gè)候選基因，分別是突觸結(jié)合蛋白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 個(gè)基因在單性狀GWAS 與多性狀GWAS 中被重復(fù)檢測(cè)到。

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

圖6 單性狀GWAS 與多性狀GWAS 結(jié)果韋恩圖

3 討 論

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

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

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

4 結(jié) 論

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