李柯安寧,杜麗麗,安炳星,鄧天宇,2,梁 忙,曹 晟,3,杜悅瑩,4,徐凌洋,高 雪,張路培,李俊雅,高會江*

(1.中國農(nóng)業(yè)科學(xué)院北京畜牧獸醫(yī)研究所,北京 100193;2.西北農(nóng)林科技大學(xué)動物科技學(xué)院,楊凌 712100;3.天津農(nóng)學(xué)院計算機與信息工程學(xué)院,天津 300384;4.青島農(nóng)業(yè)大學(xué)動物科技學(xué)院,青島 266109)

中國是牛肉產(chǎn)出和消費大國,但長期以來牛肉進出口處于貿(mào)易逆差[1]。不同部位牛肉口感差異較明顯,牛肉消費也因此細化至不同分割部位[2]。不同部位牛肉價值不同,胴體重量以及不同分割部位重量與經(jīng)濟效益直接相關(guān),美國、英國、意大利等國家此前均有對分割部位重量以及不同零售部位牛肉分級的相關(guān)研究[3-7]。而國內(nèi)肉牛領(lǐng)域相關(guān)研究多關(guān)注于肉質(zhì)性狀、生長性狀以及胴體重、屠宰率等部分屠宰性狀[8-9],對分割肉塊重量的研究較少。肉牛各零售部位分割肉塊的重量表型記錄多集中于由背最長肌分割而出的上腦、眼肉、外脊重量以及里脊重量[5,10-11],軀干外的分割肉塊諸如米龍、牛霖等部位的記錄較少,因而國內(nèi)外對分割肉塊重量這一類性狀的研究較少。Judge等[12]的研究表明,針對性的育種計劃能夠在不增加胴體重的情況下改變單個分割肉塊的重量。對分割肉塊的遺傳參數(shù)估計能夠幫助制定更精細的育種方案,優(yōu)化肉牛產(chǎn)業(yè)經(jīng)濟效益的同時,也盡力滿足消費者的需求。

現(xiàn)有研究根據(jù)自然變異現(xiàn)象鑒定到了與分割肉塊重量有關(guān)的功能基因,如肌肉生長抑制素基因[13](myostatin gene,MSTN)以及多肋基因[14](vertnin gene,VRTN)等。MSTN的自然純合突變以及基因編輯敲除個體均表現(xiàn)出雙肌肉表型[15],其個體的肩胛部、背部以及臀部的肌肉較普通個體更為發(fā)達[16]。VRTN純合個體較野生型個體的胸椎數(shù)目和肋骨數(shù)目多[14],在豬上表現(xiàn)尤為明顯,胸椎數(shù)目和肋骨數(shù)目的改變會引起體長變化以及分割肉塊“肋排”重量占胴體重比例的改變[17-18],從而對胴體的經(jīng)濟價值產(chǎn)生影響。

本研究以華西牛資源群體為研究素材,基于Illumina BovineHD SNP 770K芯片進行基因分型,對胴體性狀和原始分割肉塊性狀進行遺傳參數(shù)估計和全基因組關(guān)聯(lián)分析(genome-wide associate study, GWAS),以期獲得影響華西牛分割肉塊重量的基因及顯著位點,為華西牛后續(xù)育種方案的優(yōu)化以及進一步解析胴體相關(guān)性狀遺傳機制提供參考。

1 材料與方法

1.1 試驗群體

本研究選取的1 585頭18月齡華西牛來源于中國農(nóng)業(yè)科學(xué)院北京畜牧獸醫(yī)研究所牛遺傳育種科技創(chuàng)新團隊組建的華西牛資源群體。該資源群體于2008年創(chuàng)建于內(nèi)蒙古錫林郭勒盟烏拉蓋管理區(qū)并開始進行數(shù)據(jù)記錄。將斷奶后犢牛送至育肥場按照統(tǒng)一的飼養(yǎng)管理方案進行育肥。宰前按照統(tǒng)一管理要求斷食斷水。群體所有牛只按照國家標(biāo)準(zhǔn)GB/T 27643—2011進行數(shù)據(jù)測定、屠宰和胴體分割[19]。

1.2 基因型數(shù)據(jù)獲得與質(zhì)控

試驗群體所有牛只在育肥期間采集血樣,抗凝處理后將達到標(biāo)準(zhǔn)的DNA樣品送至芯片服務(wù)商,利用Illumina BovineHD SNP芯片進行基因型檢測。使用PLINK(V1.07)[20]軟件對芯片數(shù)據(jù)做如下處理:1)剔除最小等位基因頻率<1%的SNPs位點;2)剔除SNPs檢出率<95%的位點;3)剔除哈代-溫伯格平衡檢驗P<10-6的SNPs位點;4)剔除SNPs缺失率>5%的個體。最終保留597 082個SNPs位點。

1.3 影響因素統(tǒng)計與顯著性檢驗

表型值校正使不同條件下的個體在較一致的基礎(chǔ)上進行比較,能夠提高遺傳參數(shù)估計的準(zhǔn)確性。本研究選取的群體中,出生年相同的牛只均在同一育肥場育肥,因此固定效應(yīng)考慮場效應(yīng)與性別效應(yīng),按育肥場不同劃分為7個效應(yīng),并以進場重和育肥天數(shù)作為協(xié)變量。

利用R程序中的線性模型(linear model, LM)與R包car中的Anova函數(shù)進行固定效應(yīng)顯著性分析,采取的線性模型如下所示:

yij=μ+genderi+farmj+enterweight+
fattendays+eij

式中,yij為性狀表型值,μ為群體均值,genderi為性別固定效應(yīng),farmj為育肥場固定效應(yīng),enterweight為進場重協(xié)變量,fattendays為育肥天數(shù)協(xié)變量,eij為殘差。利用car包中Anova對模型各參數(shù)顯著性進行分析時發(fā)現(xiàn),當(dāng)場效應(yīng)與年效應(yīng)共存時,回歸模型中存在多重共線性,這表明場效應(yīng)與年效應(yīng)高度相關(guān)或一致,這與出生年相同的牛只在同一育肥場育肥的實際情況相符,因此最終模型呈現(xiàn)如上,以場效應(yīng)、性別效應(yīng)作為固定效應(yīng),進場重與育肥天數(shù)作為協(xié)變量。

1.4 遺傳參數(shù)估計模型

遺傳參數(shù)估計模型如下:

y=Xβ+Wu+e

使用GCTA軟件[21]完成W矩陣的構(gòu)建以及方差組分和遺傳參數(shù)的估計。

1.5 全基因組關(guān)聯(lián)分析

使用rMVP軟件包[22]的FarmCPU模型[23]進行單性狀GWAS分析,使用PLINK(V1.9)[24]完成主成分分析。該模型交替使用固定效應(yīng)和隨機效應(yīng)來提高位點檢出率,并使用偽數(shù)量性狀核苷酸(pseudo quantitative traits nucleotides)作為協(xié)變量,其中固定模型如下:

yi=Xb+a1PC1+…+a3PC3+b1Gi1+…+
btGit+djSij+ei

式中,yi為第i個個體的觀測值,b為固定效應(yīng)與協(xié)方差,與“1.4”中的β一致,X為對應(yīng)的關(guān)聯(lián)矩陣,PC為主成分效應(yīng),本研究中再加入前3個主成分作為協(xié)變量,前3個主成分能夠解釋68.56%的方差比例,a為主成分效應(yīng)對應(yīng)的回歸系數(shù),Git為第t個加入至模型中的偽QTN的基因型,第一次迭代為空集,bt為對應(yīng)的效應(yīng)值,Sij為第i個個體的第j個遺傳標(biāo)記基因型,dj為對應(yīng)的效應(yīng)值,ei為殘差。

隨機效應(yīng)模型如下:

yi=ui+ei

式中,yi為第i個個體的觀測值,ui為第i個個體的總遺傳效應(yīng)。

根據(jù)Bonferroni校正多重檢驗確定顯著性閾值,全基因組顯著水平閾值為8.37×10-8,全基因組水平潛在顯著閾值為1.67×10-6。

1.6 基因功能注釋

使用Ensembl數(shù)據(jù)庫對顯著位點所在位置上、下游50 kb范圍內(nèi)的基因進行檢索,收集該范圍內(nèi)的基因信息,使用NCBI網(wǎng)站對基因進行查詢以及功能的注釋,對于能夠翻譯蛋白質(zhì)的基因,利用基因名對應(yīng)的蛋白名稱在uniprot數(shù)據(jù)庫中進行蛋白質(zhì)的功能注釋與生物功能相關(guān)信息的檢索。

2 結(jié) 果

2.1 研究性狀的基本統(tǒng)計

本研究對1 585頭華西牛的胴體性狀與原始分割肉塊重量性狀進行了表型測定,描述性統(tǒng)計結(jié)果見表1。屠宰率與凈肉率性狀的表現(xiàn)相對穩(wěn)定,而不同部位的胴體原始分割肉塊重量性狀表現(xiàn)的浮動較大,變異系數(shù)的范圍在17.71%～33.85%。

2.2 固定效應(yīng)及協(xié)變量顯著性檢驗

各固定效應(yīng)水平的觀測值見表2,表3列出了性別、育肥場、進場重、育肥天數(shù)等固定效應(yīng)和協(xié)變量的顯著性檢驗結(jié)果。其中,性別效應(yīng)對胴體重、里脊、背最長肌以及臀肉重量性狀的檢驗結(jié)果為不顯著,其余效應(yīng)對各性狀的檢驗結(jié)果均為顯著。

表2 固定效應(yīng)不同水平的觀察個數(shù)Table 2 The number of observations at different levels of fixed effets

表3 不同性狀影響因素的顯著性檢驗結(jié)果Table 3 Results of the significance test of factors for different traits

2.3 遺傳力估計

使用GCTA軟件通過GBLUP模型計算所得方差組分與遺傳力結(jié)果如表4所示,分割肉塊部分中四肢及臀部區(qū)域分割肉塊重屬于中高遺傳力性狀,遺傳力范圍在0.41～0.57之間;軀干部分分割肉塊重以及出欄重、胴體重、屠宰率、凈肉率屬于中等遺傳力性狀,遺傳力范圍在0.19～0.39之間。

表4 胴體性狀及原始分割肉塊重量性狀的遺傳力Table 4 Heritability of carcass traits and primal cuts weight traits

2.4 表型相關(guān)與遺傳相關(guān)

圖1為各性狀的表型相關(guān)及個體的散點分布圖,紅線代表擬合趨勢,由圖可見屠宰率和凈肉率與其他性狀的表型相關(guān)程度較低,范圍在0.074～0.46之間;出欄重、胴體重與其余分割肉塊重量之間的表型相關(guān)均較高,范圍在0.57～0.79之間;各性狀的樣本個體分布符合正態(tài)分布。圖2為各個性狀的遺傳相關(guān),四肢及臀部區(qū)域各分割肉塊之間的遺傳相關(guān)較高,軀干部分各分割肉塊之間的遺傳相關(guān)較高,兩部分區(qū)域相互之間遺傳相關(guān)與前二者相比較低。腹肉部分與眼肉部分直接相連,在遺傳相關(guān)性上表現(xiàn)的較同區(qū)域上腦、外脊、里脊等部分更高,屠宰率、凈肉率與其余所有性狀的遺傳相關(guān)范圍在0.22～0.65之間。表5為遺傳相關(guān)和表型相關(guān)的結(jié)果匯總,上三角為表型相關(guān),下三角為遺傳相關(guān)。

上三角為表型相關(guān)數(shù)值,數(shù)值與字體越大,相關(guān)性越強,下三角為性狀兩兩之間相關(guān)性擬合分布,中間為各性狀表型值分布The value of phenotypic correlations is given above the diagonal, the larger the value and the font, the stronger the correlation, and the correlation fitting distribution between traits two by two below the diagonal, the distribution of phenotypic values of each trait is in the middle圖1 華西牛胴體性狀與原始分割肉塊重量性狀表型相關(guān)熱圖Fig.1 Heatmap of phenotypic correlation for carcass traits and primal cuts weight traits in Huaxi cattle

上三角為遺傳相關(guān)數(shù)值,下三角圓越小表示相關(guān)程度越低The value of genetic correlations are given above the diagonal, and the smaller the circle below the diagonal, the lower the correlation degree圖2 華西牛胴體性狀與原始分割肉塊重量性狀遺傳相關(guān)熱圖Fig.2 Heatmap of genetic correlation for carcass traits and primal cuts weight traits in Huaxi cattle

表5 胴體性狀與原始分割肉塊重量性狀間遺傳與表型相關(guān)Table 5 Genetic and phenotypic correlation among carcass traits and primal cuts weight traits

2.5 全基因組關(guān)聯(lián)分析

全基因組關(guān)聯(lián)分析結(jié)果的QQ圖及曼哈頓圖見圖3,QQ圖顯示,模型對各性狀影響因素校正良好。曼哈頓圖顯示本研究共檢測到26個顯著SNPs位點,通過數(shù)據(jù)庫檢索共篩選得到24個候選基因。11號染色體上的BovineHD4100009281位點以及19號染色體上的BovineHD4100014011位點均由兩個性狀同時檢出,分別為胴體重與出欄重以及胴體重與后腱子重。在兩個位點上、下游50 kb的區(qū)間內(nèi)各自檢測到5個候選基因,分別是位于11號染色體的PNPLA7、NSMF、NOXA、EXD3、ENTPD8以及19號染色體的ADAP2、RNF135、RHOT1、TEFM、ATAD5。本研究中的最小P值2.61×10-7在上腦重量性狀中檢出,該SNP位點位于2號染色體,定位到候選基因RF00001,檢出的所有位點及定位的候選基因詳細結(jié)果見表6。

A.出欄重、胴體重、屠宰率、凈肉率QQ圖與曼哈頓圖;B.軀體部分分割肉塊上腦、眼肉、外脊、背最長肌、腹肉QQ圖與曼哈頓圖;C.四肢及臀部區(qū)域分割肉塊米龍、牛霖、臀肉、后腱子QQ圖與曼哈頓圖;A. Quantile-quantile (Q-Q) plot and Manhattan plot of marketing weight, carcass weight, dressing percentage, meat percentage; B. Quantile-quantile (Q-Q) plot and Manhattan plot of the primal cuts in trunk: high rib, ribeye, striploin, longissimus dorsi, thin flank; C. Quantile-quantile (Q-Q) plot and Manhattan plot of the primal cuts in extremities and hip area: topside, knuckle, rump, shank圖3 華西牛胴體性狀與原始分割肉塊性狀關(guān)聯(lián)分析結(jié)果的QQ圖及Manhattan圖Fig.3 Manhattan and QQ plots of the GWAS for carcass traits and primal cuts weight traits in Huaxi cattle

表6 華西牛胴體性狀與原始分割肉塊重量性狀顯著位點及候選基因Table 6 Significant loci and candidate genes for carcass traits and primal cuts weight traits in Huaxi cattle

3 討 論

我國牛肉消費需求不斷增長,而國內(nèi)的產(chǎn)能不足以滿足市場供給,做好本土肉牛專門化品種資源開發(fā)與利用勢在必行[25]。華西牛是中國農(nóng)業(yè)科學(xué)院北京畜牧獸醫(yī)研究所經(jīng)過40余年不斷改良的專門化肉牛新品種。對華西牛的屠宰性狀進行遺傳參數(shù)估計以及全基因組關(guān)聯(lián)分析,能夠為后續(xù)育種策略的細化提供參考。

3.1 遺傳參數(shù)估計

本研究利用770K芯片數(shù)據(jù)與GBLUP模型計算得到華西牛出欄重、胴體重、屠宰率與凈肉率的遺傳力分別為0.29、0.32、0.21、0.19,屬于中等遺傳力性狀。牛紅等[9]利用BLUP模型計算得到的出欄重、胴體重、屠宰率、凈肉率性狀遺傳力分別為0.43、0.38、0.31、0.39;Utrera和Van Vleck[26]綜述得到的胴體重遺傳力以及Pritchard等[27]利用大數(shù)據(jù)計算14個品種的胴體重遺傳力均值分別為0.40以及0.13～0.45。對上腦、眼肉、外脊、背最長肌的遺傳力計算結(jié)果為0.22～0.37。Naserkheil等[10,28]對韓牛的上腦、眼肉、外脊及其之和的遺傳力估計結(jié)果為0.21～0.42,其中韓牛的屠宰分割對于上腦性狀的分割起始位置較國內(nèi)的更靠前,這部分分割肉塊與國內(nèi)的頸肉部位有重合,因此影響了計算的遺傳力結(jié)果。Zhu等[29]對中國肉用西門塔爾牛經(jīng)濟性狀的一項研究表明,上腦、眼肉、外脊性狀的遺傳力為0.24～0.27。以上報道與本研究中的相關(guān)結(jié)果0.22～0.37較為一致。Berry等[30]對腹肉與里脊部位重量性狀的遺傳力計算結(jié)果分別為0.14及0.37,Pabiou等[5]的研究結(jié)果為0.25與0.29。本研究對腹肉及里脊部位重量的遺傳力計算結(jié)果為0.25與0.39。此外,本研究對前腱子、后腱子、米龍、臀肉、黃瓜條、牛霖等四肢及臀部區(qū)域分割部位重量性狀進行計算得到遺傳力為0.41～0.57,結(jié)果表明這些性狀屬于高遺傳力性狀,這與此前不同肉牛品種研究中Zhu等[29]的結(jié)果0.40～0.62,Berry等[30]的結(jié)果0.45～0.75及Naserkheil等[10]的結(jié)果0.50～0.52較為一致。

四肢及臀部區(qū)域分割肉塊之間的遺傳相關(guān)性較高,軀干部分分割肉塊之間的遺傳相關(guān)性也較高,而二者相互之間的遺傳相關(guān)程度較之則稍低,這與不同區(qū)域分割肉塊的生理分布以及表型相關(guān)程度相一致,總體上各屠宰性狀之間為高遺傳相關(guān),Judge等[11]的研究也印證了這一點。本研究計算所得四肢及臀部區(qū)域遺傳相關(guān)為0.81～0.97,軀干區(qū)域分割肉塊之間遺傳相關(guān)為0.49～0.92,這與此前韓牛(0.48～0.90、0.49～0.68)以及愛爾蘭地區(qū)肉牛(0.63～0.87、0.70～0.83)的研究結(jié)果較為一致,相關(guān)研究均不涉及黃瓜條以及牛霖。在本研究中,黃瓜條與臀肉相關(guān)達到0.97,考慮是因其均屬于臀部肌肉,且為避免過度校正,未將胴體重作為協(xié)變量加入模型中,胴體重加入模型能夠使計算結(jié)果數(shù)值減小[11,30]。軀干部分腹肉與里脊的遺傳相關(guān)為0.49,腹肉僅與眼肉的遺傳相關(guān)程度較高,達到0.90,其余相關(guān)均為中等水平,這可能與在胴體分割標(biāo)準(zhǔn)指導(dǎo)下,商業(yè)分割模式下腹肉與眼肉區(qū)域之間直接相連有關(guān)。出欄重及胴體重與其余分割肉塊重量性狀的遺傳相關(guān)為0.53～0.90,屠宰率、凈肉率與其余性狀的遺傳相關(guān)較之稍低(0.29～0.65),國內(nèi)外關(guān)于屠宰率、凈肉率遺傳相關(guān)的研究較少,牛紅等[9]計算屠宰率、凈肉率與宰前活重、胴體重的遺傳相關(guān)為0.23～0.52,Pariacote等[31]計算美國短角牛屠宰率與胴體重表型相關(guān)與遺傳相關(guān)為0.41±0.03及0.65±0.19,以上報道與本研究結(jié)果較為一致。

3.2 全基因組關(guān)聯(lián)分析

本研究共篩選出24個候選基因,其中部分可能是影響肉牛生長發(fā)育的關(guān)鍵候選基因。由胴體重與出欄重性狀共同定位的基因為ENTPD8、NSMF、NOXA、PNPLA7、EXD3。其中ENTPD8在雞的相關(guān)研究中被認為是參與調(diào)控雞肉風(fēng)味物質(zhì)5′肌苷單磷酸(inosine 5′-monophosphate,IMP)沉積的關(guān)鍵基因[32];NSMF基因表達水平的增高能增強成肌細胞的增殖,在基因敲除小鼠試驗中,NSMF的敲除不影響肌肉細胞的分化,但細胞增殖速率降低[33],本研究在11號染色體上定位到該基因,或與肉牛肌肉發(fā)育有關(guān);NOXA作為下游轉(zhuǎn)錄凋亡靶基因影響p53家族蛋白的表達,參與氧化應(yīng)激誘導(dǎo)的細胞凋亡[34]。Beyfuss和Hood[35]的研究表明,p53家族蛋白參與禁食等應(yīng)激反應(yīng)中骨骼肌的氧化應(yīng)激過程,維持骨骼肌細胞的健康狀態(tài)及功能。NOXA的表達或與肉牛應(yīng)激后體重變化有關(guān);PNPLA7編碼的蛋白質(zhì)能夠促進肝臟極低密度脂蛋白的分泌,從而參與調(diào)控脂質(zhì)代謝[36]。由胴體重與后腱子共同定位的基因為ADAP2、RNF135、TEFM、ATAD5、RHOT1。其中RNF135基因微缺失的病例被發(fā)現(xiàn)存在過度生長的現(xiàn)象[37],但還未在動物生長發(fā)育相關(guān)領(lǐng)域研究中有所報道;TEFM編碼的蛋白質(zhì)為線粒體轉(zhuǎn)錄延長因子,通過調(diào)控線粒體DNA的復(fù)制與表達來影響細胞能量轉(zhuǎn)換過程中的氧化磷酸化[38],該基因在美利奴綿羊中被認為與生長發(fā)育有關(guān)[39]。由屠宰率性狀3個顯著SNPs位點共同鑒定到的候選基因SLC22A11及SLC22A12均屬于有機陰離子轉(zhuǎn)運蛋白,參與尿酸代謝[40],而另一項研究表明這與印度人的特異型肥胖有關(guān)[41]。

由軀干部分分割肉塊性狀定位到的候選基因有APP、RF00001、LAPTM5、MATN1、MKX、SERBP1、FOXO1。背最長肌性狀2號染色體上的一個顯著位點在區(qū)間內(nèi)鑒定到兩個基因,分別是LAPTM5與MATN1,其中MATN1編碼軟骨基質(zhì)蛋白,該蛋白分布于細胞外基質(zhì),特別是軟骨細胞附近,調(diào)控軟骨細胞以及骨骼的正常發(fā)育,該基因或與牛的成年體型有關(guān),骨架大小也將最終影響育肥效果[42];MKX基因編碼的蛋白質(zhì)參與肌肉分化與損傷修復(fù)肌肉再生過程,當(dāng)MKX的表達被抑制時會造成跟腱骨化以及早期的肌肉發(fā)育不良[43],推測該基因與肉牛育成期肌肉及骨骼發(fā)育有關(guān);SERBP1編碼產(chǎn)物為mRNA結(jié)合蛋白,已有研究表明其可能與豬肌內(nèi)脂肪沉積有關(guān)[44];FOXO1編碼的蛋白叉頭框轉(zhuǎn)錄因子(forkhead box protein O1)為胰島素信號傳導(dǎo)的主要靶標(biāo),相關(guān)研究表明該轉(zhuǎn)錄因子的下調(diào)將促進脂肪沉積[45]、細胞增殖[46]以及骨骼肌再生[47],該基因位于牛12號染色體,由腹肉性狀鑒定得到,是潛在的與肉牛生長發(fā)育有關(guān)的候選基因。

4 結(jié) 論

本研究對華西牛胴體和原始分割肉塊重量性狀進行了遺傳參數(shù)估計及全基因組關(guān)聯(lián)分析,結(jié)果表明出欄重、胴體重、屠宰率、凈肉率、頸肉、上腦、眼肉、外脊、里脊、背最長肌、腹肉屬于中等遺傳力性狀(0.19～0.39),前腱子、米龍、黃瓜條、牛霖、臀肉及后腱子屬于高遺傳力性狀(0.41～0.57)。各分割肉塊之間,及其與出欄重、胴體重性狀間均為高遺傳相關(guān),而屠宰率、凈肉率與上述性狀之間的遺傳相關(guān)較之稍低。對這些性狀的全基因組關(guān)聯(lián)分析共檢測到26個顯著性SNPs位點,定位到24個候選基因,其中NSMF、NOXA、TEFM、MATN1、MKX,FOXO1可能是影響肉牛肌肉生長的候選基因。以上結(jié)果為后續(xù)華西牛育種策略的細化提供參考。