基于GBLUP等模型對(duì)梅花鹿(Cervus Nippon)生長(zhǎng)相關(guān)性狀基因組選擇的預(yù)測(cè)準(zhǔn)確性比較

李浩東，閔祥玉，周 雅，張禾垟，鄭軍軍，劉琳玲，王 平，王艷梅，楊福合*，王桂武*

(1.中國(guó)農(nóng)業(yè)科學(xué)院特產(chǎn)研究所，長(zhǎng)春 130112； 2.吉林省梅花鹿產(chǎn)業(yè)研究院，長(zhǎng)春 130600)

我國(guó)不僅擁有悠久的梅花鹿飼養(yǎng)歷史，還擁有最大的梅花鹿飼養(yǎng)群體，是梅花鹿遺傳育種的主戰(zhàn)場(chǎng)[1]。特別是在2020年，梅花鹿被列入《國(guó)家畜禽遺傳資源品種目錄》[2]，使得養(yǎng)鹿迎來了大發(fā)展。其主產(chǎn)品鹿茸具有極高的醫(yī)學(xué)和藥用價(jià)值，明代李時(shí)珍在《本草綱目》上就有關(guān)于鹿茸的記載，稱其“善于補(bǔ)腎壯陽，生精益血，補(bǔ)髓健骨”。而且鹿肉、鹿鞭等也能很好的契合健康的主題。但長(zhǎng)期以來，梅花鹿的選育都以鹿茸為主要經(jīng)濟(jì)目的，重茸輕肉，而忽略了對(duì)肉用鹿品種的培育。鹿肉性溫和，有補(bǔ)脾益氣、溫腎壯陽的功效，瘦肉多、結(jié)締組織少，蛋白質(zhì)含量豐富、脂肪少，具有極高的營(yíng)養(yǎng)價(jià)值。目前，鹿肉在國(guó)內(nèi)擁有廣闊的市場(chǎng)，但肉用鹿尚無專門化的品種，在市場(chǎng)及飯店出售的鹿肉多為老弱病殘的淘汰鹿，品質(zhì)不佳，安全性難以保證[3]。

在傳統(tǒng)育種中，畜禽的選配主要是利用表型以及經(jīng)驗(yàn)來進(jìn)行的，后過渡到以最佳線性無偏預(yù)測(cè)(best linear unbiased prediction，BLUP)為主通過評(píng)估畜禽育種值來進(jìn)行選擇的方法[4]。隨著高通量測(cè)序技術(shù)的快速發(fā)展和測(cè)序成本降低使得畜禽育種已經(jīng)進(jìn)入分子育種階段，分子標(biāo)記輔助選擇(marker assisted selection，MAS)的應(yīng)用使畜禽遺傳改良有了長(zhǎng)足的進(jìn)步。但由于經(jīng)過功能驗(yàn)證用于畜禽遺傳改良的功能基因與分子標(biāo)記較少，且較多分子標(biāo)記對(duì)遺傳變異比例解釋太低，限制了分子輔助標(biāo)記的應(yīng)用[5-6]。2001年，Meuwissen等[7]提出了基因組選擇的概念，全基因組選擇(genomic selection，GS)被作為是全基因組范圍內(nèi)的標(biāo)記輔助選擇，是當(dāng)下畜禽遺傳育種的主流策略[8]，相較于基于系譜的估計(jì)育種值(estimated breeding value, EBV)，利用覆蓋全基因組的標(biāo)記位點(diǎn)估計(jì)的基因組估計(jì)育種值(genomic estimated breeding value, GEBV)更加準(zhǔn)確[9]。統(tǒng)計(jì)模型是基因組選擇的核心，極大地影響了基因組預(yù)測(cè)的準(zhǔn)確性和效率。根據(jù)統(tǒng)計(jì)模型的不同，可以分為以GBLUP為代表的直接法[10]和以貝葉斯為代表的間接法[11]。牛[12-13]、豬[14-15]、雞[16-17]等動(dòng)物育種相繼進(jìn)入基因組選擇育種時(shí)代。尤其是在奶牛育種中，基因組選擇已在美國(guó)[18-19]、澳大利亞[20-21]、英國(guó)[22-23]、新西蘭[24-25]等國(guó)家廣泛應(yīng)用，并取得了明顯的遺傳進(jìn)展。

目前，梅花鹿的選育工作還停留于傳統(tǒng)的育種方法。由于梅花鹿尚未完全馴化，仍具有一定野性，導(dǎo)致生產(chǎn)性能測(cè)定難度極大，并且大多數(shù)養(yǎng)殖場(chǎng)系譜記錄混亂甚至無系譜記錄，而基因組選擇的優(yōu)勢(shì)在于能夠降低測(cè)定難度以及利用全基因組的標(biāo)記位點(diǎn)構(gòu)建親緣關(guān)系替代系譜，提高育種值估計(jì)的準(zhǔn)確性，縮短時(shí)代間隔。隨著梅花鹿由小戶散養(yǎng)向集約化養(yǎng)殖的轉(zhuǎn)變，基因組選擇是提高整體生產(chǎn)水平的最佳選擇。因此，本試驗(yàn)基于梅花鹿的生長(zhǎng)相關(guān)性狀來評(píng)估各種模型預(yù)測(cè)的準(zhǔn)確性，以探究在梅花鹿群體中實(shí)施基因組預(yù)測(cè)的可行性。

1 材料與方法

1.1 試驗(yàn)動(dòng)物

本研究動(dòng)物群體來自吉林省長(zhǎng)春市雙陽某梅花鹿繁育場(chǎng)2014—2019年所產(chǎn)純種梅花鹿，共計(jì)261只。其中公鹿96只，母鹿165只。測(cè)定的生長(zhǎng)性狀包括體重(body weight)、體高(body height)、體斜長(zhǎng)(body length)、管圍(shin circumference)、胸圍(chest girth)、臀端高(high hips)。利用Excel對(duì)表型數(shù)據(jù)進(jìn)行預(yù)處理，刪除表型的缺失值和異常值(平均值±3倍標(biāo)準(zhǔn)差)。通過R語言中的GLM函數(shù)校正性別和出生年份。梅花鹿表型描述性統(tǒng)計(jì)見表1。

表1 梅花鹿體重與體尺指標(biāo)的表型統(tǒng)計(jì)量

1.2 基因型數(shù)據(jù)的獲取與質(zhì)量控制

采取頸靜脈采血的方式抽取血液用于基因組DNA的提取，將樣品送至博瑞迪生物技術(shù)有限公司，利用本團(tuán)隊(duì)開發(fā)的梅花鹿液相40K育種芯片[26]進(jìn)行基因分型。芯片分型共得到了92 246個(gè)SNPs位點(diǎn)，使用PLINK(V1.90)軟件[27]對(duì)分型結(jié)果進(jìn)行質(zhì)量控制，質(zhì)量控制標(biāo)準(zhǔn)：1)個(gè)體SNP檢出率>90%；2)SNP檢出率>90%；3)最小等位基因頻率≥5%；4)哈代溫伯格平衡P<1×10-6。最終獲得261個(gè)個(gè)體的84 059個(gè)SNPs位點(diǎn)用于后續(xù)分析。

1.3 基因組選擇模型

使用ASReml4.1.0軟件[28]通過約束最大似然法進(jìn)行方差組分的估計(jì)與遺傳力的計(jì)算?；谌蚪MSNP信息構(gòu)建親緣關(guān)系G矩陣，通過加性遺傳方差除以總方差計(jì)算性狀遺傳力。計(jì)算公式如下：

1.3.1 基因組最佳線性無偏預(yù)測(cè)(GBLUP) GBLUP通過SNP信息構(gòu)建分子親緣關(guān)系矩陣(G矩陣)替代基于系譜信息構(gòu)建的關(guān)系矩陣(A矩陣)，該方法由 Vanraden[29]提出。該方法的統(tǒng)計(jì)模型如下：

y=Xb+Zg+e

1.3.3 嶺回歸最佳線性無偏預(yù)測(cè)法(RRBLUP) RRBLUP是間接法的代表。Whittaker等[32]為了解決標(biāo)記效應(yīng)誤差的問題，因而提出的一種改良的最小二乘法，類似于GBLUP法。其將標(biāo)記效應(yīng)作為隨機(jī)效應(yīng)，假定標(biāo)記效應(yīng)服從正態(tài)分布，利用混合線性模型求解標(biāo)記效應(yīng)的總和即為個(gè)體育種值。RRBLUP模型的計(jì)算使用R語言環(huán)境下的RRBLUP包進(jìn)行[33]。

1.4 基因組預(yù)測(cè)準(zhǔn)確性的評(píng)估

對(duì)以上6種方法通過5-fold交叉驗(yàn)證的方法來評(píng)價(jià)GEBV的估計(jì)準(zhǔn)確性，先將梅花鹿群體隨機(jī)分成5個(gè)獨(dú)立的子集，然后依次將每個(gè)子集作為驗(yàn)證群體，并假定其表型未知，同時(shí)利用其余4個(gè)子集做參考群體來預(yù)測(cè)驗(yàn)證群體的GEBV。預(yù)測(cè)能力為表型與GEBV間的相關(guān)性。為避免遺傳力對(duì)預(yù)測(cè)準(zhǔn)確性的影響，將表型與GEBV間的皮爾遜相關(guān)系數(shù)除以相應(yīng)性狀遺傳力算術(shù)平方根作為 GEBV 準(zhǔn)確性的評(píng)價(jià)指標(biāo)[34]。計(jì)算公式為：

其中,y表示校正后的表型值,h2表示性狀遺傳力。對(duì)每個(gè)性狀的驗(yàn)證進(jìn)行10次重復(fù)，最后用相關(guān)系數(shù)的平均值作為模型準(zhǔn)確性的評(píng)估指標(biāo)。

2 結(jié) 果

2.1 梅花鹿生長(zhǎng)相關(guān)性狀遺傳參數(shù)的估計(jì)

利用全基因組SNP信息構(gòu)建親緣關(guān)系矩陣，建立混合線性模型，通過約束最大似然法對(duì)梅花鹿生長(zhǎng)相關(guān)性狀的加性遺傳方差、環(huán)境方差、遺傳力以及標(biāo)準(zhǔn)誤進(jìn)行估計(jì)。結(jié)果表明，梅花鹿群體重與體尺性狀的遺傳力范圍為0.15(胸圍)～0.50(臀端高)。其中，管圍與臀端高的遺傳力分別為0.43、0.50，屬于高遺傳力；體重、體高與體斜長(zhǎng)的遺傳力分別為0.22、0.30、0.27，屬于中等遺傳力，而胸圍的遺傳力為0.15，屬于低遺傳力。各性狀方差組分和遺傳力及標(biāo)準(zhǔn)誤估計(jì)見表2。

表2 梅花鹿生長(zhǎng)相關(guān)性狀的方差組分和遺傳力估計(jì)

2.2 梅花鹿生長(zhǎng)相關(guān)性狀基因組選擇不同模型的預(yù)測(cè)準(zhǔn)確性比較

通過5-fold交叉驗(yàn)證，利用GBLUP、Bayes A、Bayes B、Bayes C、Bayes Lasso、RRBLUP六種基因組選擇模型對(duì)生長(zhǎng)相關(guān)性狀的預(yù)測(cè)準(zhǔn)確性進(jìn)行評(píng)估，將表型與GEBV間的皮爾遜相關(guān)系數(shù)除以相應(yīng)性狀遺傳力算術(shù)平方根作為 GEBV 準(zhǔn)確性的評(píng)價(jià)指標(biāo)。

2.2.1 梅花鹿體重基因組選擇預(yù)測(cè)準(zhǔn)確度 如表3所示:體重的預(yù)測(cè)準(zhǔn)確度分別是0.231、0.143、0.121、0.123、0.086、0.077。在體重中，GBLUP預(yù)測(cè)準(zhǔn)確性最高，Bayes A、Bayes C、Bayes B和Bayes Lasso次之，RRBLUP最低。

2.2.2 梅花鹿體高基因組選擇預(yù)測(cè)準(zhǔn)確度 如表3所示：體高的預(yù)測(cè)準(zhǔn)確度分別是0.205、0.324、0.339、0.329、0.310、0.336。在體高中，Bayes B預(yù)測(cè)準(zhǔn)確性最高，RRBLUP、Bayes C、Bayes A和Bayes Lasso次之，GBLUP最低。

2.2.3 梅花鹿體斜長(zhǎng)基因組選擇預(yù)測(cè)準(zhǔn)確度 結(jié)果如表3所示：體斜長(zhǎng)的預(yù)測(cè)準(zhǔn)確度分別是0.256、0.336、0.333、0.338、0.314、0.314。在體斜長(zhǎng)中，Bayes C預(yù)測(cè)準(zhǔn)確性最高，Bayes A、Bayes B、RRBLUP和Bayes Lasso次之， GBLUP最低。

2.2.4 梅花鹿管圍基因組選擇預(yù)測(cè)準(zhǔn)確度 如表3所示：管圍的預(yù)測(cè)準(zhǔn)確度分別是0.197、0.222、0.230、0.222、0.223、0.217。在管圍中，Bayes B預(yù)測(cè)準(zhǔn)確性最高，Bayes Lasso、Bayes A、Bayes C和RRBLUP次之， GBLUP最低。

2.2.5 梅花鹿胸圍基因組選擇預(yù)測(cè)準(zhǔn)確度 如表3所示：胸圍的預(yù)測(cè)準(zhǔn)確度分別是0.190、0.601、0.593、0.579、0.529、0.583。在胸圍中，Bayes A預(yù)測(cè)準(zhǔn)確性最高，Bayes B、RRBLUP、Bayes C和Bayes Lasso次之， GBLUP最低。

2.2.6 梅花鹿臀端高基因組選擇預(yù)測(cè)準(zhǔn)確度 如表3所示：臀端高的預(yù)測(cè)準(zhǔn)確度分別是0.192、0.318、0.319、0.318、0.311、0.325。在臀端高中，RRBLUP預(yù)測(cè)準(zhǔn)確性最高，Bayes B、Bayes A、Bayes C和Bayes Lasso次之，GBLUP最低。

表3 基于6種基因組選擇模型對(duì)梅花鹿生長(zhǎng)相關(guān)性狀的預(yù)測(cè)能力與準(zhǔn)確性

2.3 梅花鹿生長(zhǎng)相關(guān)性狀遺傳力對(duì)基因組預(yù)測(cè)準(zhǔn)確性的影響

如表3所示，預(yù)測(cè)能力的范圍為從0.036(RRBLUP法中體重)到0.233(Bayes A法中胸圍)。在GBLUP中，各性狀的預(yù)測(cè)能力由高到低依次為臀端高、體斜長(zhǎng)、管圍、體高、體重、胸圍；在Bayes A與Bayes B中，各性狀的預(yù)測(cè)能力由高到低依次皆為胸圍、臀端高、體高、體斜長(zhǎng)、管圍、體重；在Bayes C、Bayes LASSO與RRBLUP中，各性狀的預(yù)測(cè)能力由高到低依次皆為臀端高、胸圍、體高、體斜長(zhǎng)、管圍、體重。有大量研究表明，基因組選擇的預(yù)測(cè)能力與性狀的遺傳力估計(jì)值呈正相關(guān)關(guān)系[35-37]。為了探索遺傳力與預(yù)測(cè)能力之間的關(guān)系，對(duì)遺傳力與預(yù)測(cè)能力進(jìn)行回歸分析。如圖1所示，6種基因組選擇模型的相應(yīng)回歸系數(shù)為0.144(GBLUP)、0.074(Bayes A)、0.101(Bayes B)、0.1(Bayes C)、0.153(Bayes Lasso)、0.132(RRBLUP)。

圖1 6種基因組選擇模型的預(yù)測(cè)能力與遺傳力回歸分析

3 討 論

本研究利用基因組數(shù)據(jù)對(duì)梅花鹿體重體尺性狀進(jìn)行了遺傳力估計(jì)，并系統(tǒng)的探討了以GBLUP為代表的直接法模型和以Bayes類方法為代表的間接法對(duì)生長(zhǎng)相關(guān)性狀基因組選擇的預(yù)測(cè)能力與預(yù)測(cè)準(zhǔn)確性。利用基因組信息基于GBLUP通過約束最大似然法對(duì)方差組分進(jìn)行估計(jì)，進(jìn)而計(jì)算遺傳力。遺傳力的估計(jì)受動(dòng)物群體種類、群體大小以及標(biāo)記密度等因素的影響。但在本研究中，動(dòng)物群體以及大小等因素是固定存在的，因此均已無法探討，但已將性別與年齡作為固定效應(yīng)加入遺傳力的估計(jì)中。在本研究中，體尺性狀的遺傳力均為首次估計(jì)。其中，管圍與臀端高屬于高遺傳力；體重、體高與體斜長(zhǎng)屬于中等遺傳力，而胸圍屬于低遺傳力?；谌蚪MSNP信息構(gòu)建的親緣關(guān)系估計(jì)的梅花鹿體重遺傳力遠(yuǎn)低于周世朗等[38]基于系譜構(gòu)建的親緣關(guān)系估計(jì)的遺傳力，這可能有兩方面原因：1)模型高估了顯性遺傳方差，因此低估了加性遺傳方差，導(dǎo)致遺傳力的估計(jì)偏低[39]；2)全基因組SNP信息構(gòu)建的親緣關(guān)系矩陣相較于系譜構(gòu)建的親緣關(guān)系矩陣能夠更加準(zhǔn)確的反映個(gè)體間的遺傳關(guān)系[40]。在基于系譜構(gòu)建的親緣關(guān)系中，共同的環(huán)境組分會(huì)導(dǎo)致高估加性遺傳方差，使得遺傳力的估計(jì)偏高[41]。因此，在對(duì)梅花鹿生長(zhǎng)相關(guān)性狀進(jìn)行遺傳評(píng)估時(shí)，利用基因組信息構(gòu)建親緣關(guān)系矩陣會(huì)更加準(zhǔn)確。

本研究基于5-fold交叉驗(yàn)證對(duì)GBLUP、Bayes A、Bayes B、Bayes C、Bayes Lasso、RRBLUP六種基因組選擇模型預(yù)測(cè)的準(zhǔn)確度進(jìn)行了比較，分析了不同模型對(duì)不同性狀的適用性。由于交叉驗(yàn)證分組存在隨機(jī)性，單次的交叉驗(yàn)證不具有代表性。驗(yàn)證過程中，重復(fù)的次數(shù)越多，預(yù)測(cè)的準(zhǔn)確性更高。本研究對(duì)各個(gè)性狀的驗(yàn)證進(jìn)行10次5-fold交叉驗(yàn)證重復(fù)，用以校正評(píng)估的準(zhǔn)確性。在體重中，GBLUP法的預(yù)測(cè)準(zhǔn)確度均優(yōu)于Bayes類與RRBLUP法。證明在樣本量較少的情況下，GBLUP的預(yù)測(cè)準(zhǔn)確度并非一定低于間接法[42]；Bayes B在各個(gè)性狀中的預(yù)測(cè)準(zhǔn)確性都較高，可能是由于模型假設(shè)更加符合數(shù)量性狀的遺傳結(jié)構(gòu)。它假設(shè)有效應(yīng)位點(diǎn)比例為1-π，通常為0.05左右；Bayes A在體重、體斜長(zhǎng)、胸圍等低遺傳力性狀中預(yù)測(cè)準(zhǔn)確度優(yōu)于Bayes B。因而在低遺傳力性狀中，首選Bayes A進(jìn)行預(yù)測(cè)是最佳選擇；體重、體高、體斜長(zhǎng)、管圍、胸圍、臀端高預(yù)測(cè)準(zhǔn)確度最高的分別為GBLUP、Bayes B、Bayes C、Bayes B、Bayes A、RRBLUP。因此，在基因組選擇中進(jìn)行模型選擇時(shí)并沒有統(tǒng)一的標(biāo)準(zhǔn)，要根據(jù)性狀的遺傳結(jié)構(gòu)選擇最適合的模型才能使得預(yù)測(cè)更加精準(zhǔn)。

本研究利用6種基因組選擇模型探討了對(duì)體重體尺性狀的預(yù)測(cè)能力，并對(duì)遺傳力與預(yù)測(cè)能力進(jìn)行了回歸分析。研究發(fā)現(xiàn)，性狀遺傳力的高低會(huì)影響基因組選擇預(yù)測(cè)的準(zhǔn)確性[25-27]。在GBLUP中，各個(gè)性狀遺傳力與預(yù)測(cè)準(zhǔn)確性大致呈正相關(guān)關(guān)系，這與相關(guān)研究結(jié)果相一致。但在Bayes類方法與RRBLUP中，低遺傳力的胸圍卻擁有較高的預(yù)測(cè)準(zhǔn)確度,原因可能是梅花鹿的特殊性導(dǎo)致，尚未完全馴化，具有較大的野性，在鹿群進(jìn)食過程中飼料營(yíng)養(yǎng)不均衡，使得鹿只胸圍的變異系數(shù)較大，導(dǎo)致在進(jìn)行遺傳力估計(jì)時(shí)高估了環(huán)境方差，導(dǎo)致遺傳力的估計(jì)偏低。因此，在梅花鹿的實(shí)際生產(chǎn)過程中，要更加注意飼養(yǎng)管理，以達(dá)到鹿只的體況均勻。既能保障生產(chǎn)水平的提高，同時(shí)也能提高基因組選擇預(yù)測(cè)的準(zhǔn)確性。

4 結(jié) 論

本研究以梅花鹿作為研究群體，對(duì)梅花鹿體重體尺性狀進(jìn)行了遺傳力估計(jì)，并基于5-fold交叉驗(yàn)證對(duì)GBLUP、Bayes A、Bayes B、Bayes C、Bayes Lasso、RRBLUP六種基因組選擇模型進(jìn)行了基因組選擇分析。研究發(fā)現(xiàn)，體重、體高、體斜長(zhǎng)、管圍、胸圍與臀端高的遺傳力分別為0.22、0.30、0.27、0.43、0.15、0.50，且性狀遺傳力的高低對(duì)基因組預(yù)測(cè)準(zhǔn)確性有重要影響。體重、體高、體斜長(zhǎng)、管圍、胸圍、臀端高預(yù)測(cè)準(zhǔn)確度最高的分別為GBLUP、Bayes B、Bayes C、Bayes B、Bayes A、RRBLUP。說明在實(shí)際生產(chǎn)中，沒有能夠完全適應(yīng)所有性狀的模型，必須根據(jù)預(yù)測(cè)的準(zhǔn)確性以及預(yù)測(cè)的時(shí)效性來特異的選擇最佳模型。本研究為梅花鹿生長(zhǎng)相關(guān)性狀的基因組選擇與實(shí)際應(yīng)用提供了重要的參考信息。