華西牛胴體及原始分割肉塊重量性狀遺傳參數估計與全基因組關聯分析

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

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

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

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

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

1 材料與方法

1.1 試驗群體

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

1.2 基因型數據獲得與質控

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

1.3 影響因素統計與顯著性檢驗

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

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

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

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

1.4 遺傳參數估計模型

遺傳參數估計模型如下:

y=Xβ+Wu+e

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

1.5 全基因組關聯分析

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

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

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

隨機效應模型如下:

yi=ui+ei

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

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

1.6 基因功能注釋

使用Ensembl數據庫對顯著位點所在位置上、下游50 kb范圍內的基因進行檢索,收集該范圍內的基因信息,使用NCBI網站對基因進行查詢以及功能的注釋,對于能夠翻譯蛋白質的基因,利用基因名對應的蛋白名稱在uniprot數據庫中進行蛋白質的功能注釋與生物功能相關信息的檢索。

2 結 果

2.1 研究性狀的基本統計

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

2.2 固定效應及協變量顯著性檢驗

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

表2 固定效應不同水平的觀察個數Table 2 The number of observations at different levels of fixed effets

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

2.3 遺傳力估計

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

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

2.4 表型相關與遺傳相關

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

上三角為表型相關數值,數值與字體越大,相關性越強,下三角為性狀兩兩之間相關性擬合分布,中間為各性狀表型值分布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 華西牛胴體性狀與原始分割肉塊重量性狀表型相關熱圖Fig.1 Heatmap of phenotypic correlation for carcass traits and primal cuts weight traits in Huaxi cattle

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

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

2.5 全基因組關聯分析

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

A.出欄重、胴體重、屠宰率、凈肉率QQ圖與曼哈頓圖;B.軀體部分分割肉塊上腦、眼肉、外脊、背最長肌、腹肉QQ圖與曼哈頓圖;C.四肢及臀部區域分割肉塊米龍、牛霖、臀肉、后腱子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 華西牛胴體性狀與原始分割肉塊性狀關聯分析結果的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 討 論

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

3.1 遺傳參數估計

本研究利用770K芯片數據與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]利用大數據計算14個品種的胴體重遺傳力均值分別為0.40以及0.13～0.45。對上腦、眼肉、外脊、背最長肌的遺傳力計算結果為0.22～0.37。Naserkheil等[10,28]對韓牛的上腦、眼肉、外脊及其之和的遺傳力估計結果為0.21～0.42,其中韓牛的屠宰分割對于上腦性狀的分割起始位置較國內的更靠前,這部分分割肉塊與國內的頸肉部位有重合,因此影響了計算的遺傳力結果。Zhu等[29]對中國肉用西門塔爾牛經濟性狀的一項研究表明,上腦、眼肉、外脊性狀的遺傳力為0.24～0.27。以上報道與本研究中的相關結果0.22～0.37較為一致。Berry等[30]對腹肉與里脊部位重量性狀的遺傳力計算結果分別為0.14及0.37,Pabiou等[5]的研究結果為0.25與0.29。本研究對腹肉及里脊部位重量的遺傳力計算結果為0.25與0.39。此外,本研究對前腱子、后腱子、米龍、臀肉、黃瓜條、牛霖等四肢及臀部區域分割部位重量性狀進行計算得到遺傳力為0.41～0.57,結果表明這些性狀屬于高遺傳力性狀,這與此前不同肉牛品種研究中Zhu等[29]的結果0.40～0.62,Berry等[30]的結果0.45～0.75及Naserkheil等[10]的結果0.50～0.52較為一致。

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

3.2 全基因組關聯分析

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

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

4 結 論

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