中國荷斯坦牛繁殖性狀的基因組預測效果比較

師 睿，蘇國生，陳紫薇，李 想，羅漢鵬，劉 林，郭 剛，張 毅， 王雅春*，張勝利，張 勤，5

(1.中國農業(yè)大學動物科技學院，北京 100193； 2.丹麥奧胡斯大學，切勒 DK-8830； 3.北京奶牛中心，北京 100192； 4.北京首農畜牧發(fā)展有限公司，北京 100176； 5.山東農業(yè)大學動物科技學院，泰安 271018)

實際生產(chǎn)中，牛群的繁殖效率至關重要，對牛場的利潤會產(chǎn)生直接影響，體現(xiàn)在獸醫(yī)、繁殖管理以及牛群結構(淘汰率)等方面。奶牛繁殖性狀具有重要的經(jīng)濟價值，而其與生產(chǎn)性狀有一定的拮抗關系。因此，繁殖性狀在實際育種指數(shù)中占有一定權重。繁殖性狀為低遺傳力性狀，由微效多基因調控，在使用常規(guī)遺傳分析手段時難以獲得較好的選擇效果。

自基因組選擇(genomic selection，GS)的概念提出以來，其逐漸被廣泛應用于我國各種動物實際育種過程中，也成為奶牛育種策略的重要組成部分。GS是利用覆蓋于個體全基因組的單核苷酸多態(tài)性(single nucleotide polymorphism，SNP)標記代表個體的基因組信息，結合常規(guī)遺傳評估信息來計算個體的基因組估計育種值(genomic estimated breeding value，GEBV)。隨著GS理論的發(fā)展，人們使用高密度標記構建個體間的親緣關系矩陣，用于替代傳統(tǒng)育種中由系譜構建的親緣關系矩陣。其中，VanRaden基于SNP信息構建了矩陣，結合最佳線性無偏估計(best linear unbiased prediction，BLUP)方法用于估計個體GEBV，稱為GBLUP。Aguilar等提出，將系譜信息及基因組信息同時利用，構建了混合親緣關系矩陣陣用于遺傳評估，該方法稱為一步法GBLUP(single-step GBLUP，ssGBLUP)。在ssGBLUP方法中，可選擇不同的參數(shù)進行陣的構建，不同參數(shù)組合(α、β、τ和ω；定義詳見材料方法部分)可能會影響最終的預測效果。Misztal等對這些參數(shù)進行梯度賦值，評估了參數(shù)變化對荷斯坦牛體型性狀基因組預測效果的影響。研究表明，當τ的取值線性升高時，基因組預測準確性提高了0.14左右，且各性狀無偏性(unbiasedness，b)取值也會有不同程度的增加；當ω的取值線性升高時，基因組預測準確性上升了0.17左右，但各性狀b值會有較大程度的下降。因此，尋找ssGBLUP的最佳參數(shù)組合能夠保證基因組預測時無偏性最優(yōu)的同時達到較高的準確性。

GS較傳統(tǒng)奶牛育種的優(yōu)勢在于節(jié)省后裔測定成本、進行早期選擇以減少世代間隔，提高育種值估計的準確性，并加速遺傳進展。很多研究表明，將基因組信息納入育種過程后，能夠顯著提高估計育種值的可靠性，對低遺傳力性狀(例如繁殖性狀)更加有效。但是，目前國內尚未有成熟的奶牛繁殖性狀GS體系，且各類GS方法對中國荷斯坦群體繁殖性狀的基因組預測效果，以及預測模型不同參數(shù)選擇對GS效果的影響都是未知的。

本研究基于BLUP、GBLUP和ssGBLUP方法對北京地區(qū)中國荷斯坦牛群的部分繁殖性狀進行了基因組預測，并比較了系譜信息和基因組的權重對ssGBLUP預測的影響， 旨在填補我國奶牛繁殖性狀GS體系的空白，并為后續(xù)選種及育種策略改良提供理論依據(jù)。

1 材料與方法

1.1 試驗數(shù)據(jù)

1.1.1 表型數(shù)據(jù) 原始繁殖記錄來源于北京首農畜牧有限公司33個規(guī)?；膛霭⒎平鸸芾硐到y(tǒng)。數(shù)據(jù)收集時間為1990年1月至2020年12月，共保留出生記錄381 640條，產(chǎn)犢記錄509 054條，配種記錄1 789 650條和妊檢記錄663 759條。系譜數(shù)據(jù)來源于中國奶牛數(shù)據(jù)中心(https://www.holstein.org.cn)，通過對牛只出生記錄和產(chǎn)犢記錄進行整理追溯，去除信息不明和錯誤的個體，最終包含211 131頭個體，其中包括6 171頭公牛。根據(jù)繁殖性狀的研究報告及生產(chǎn)實際，本研究主要對3個重要繁殖性狀進行分析：1)產(chǎn)犢至首次配種間隔(interval from calving to first service，ICF，d)，即產(chǎn)犢日期與該胎次首次配種日期的間隔天數(shù)；2)青年牛配種次數(shù)(number of services for heifers，NSH)，即青年牛第一次產(chǎn)犢前的配種次數(shù)；3)成母牛配種次數(shù)(number of services for cows，NSC)，即成母牛每個胎次的配種次數(shù)。其中，成母牛的胎次范圍為1～3胎。

各繁殖性狀表型質量控制范圍參考北歐遺傳評估中心(Nordic cattle genetic evaluation，NAV)制定的標準。各性狀經(jīng)過最終質控后的描述性統(tǒng)計結果如表1所示。

表1 各性狀描述性統(tǒng)計Table 1 Descriptive statistics for each trait

1.1.2 基因組數(shù)據(jù) 芯片測定數(shù)據(jù)來自于課題組數(shù)據(jù)庫，共有12 195頭，包括8 718頭母牛和3 477頭公牛。所有芯片測定牛中，有3 119頭母牛和81頭公牛采用GeneSeek公司牛150K芯片(Illumina Bovine 150K SNP)，剩余牛只采用50K芯片(Illumina Bovine 50K SNP)，后續(xù)通過Beagle v5.0軟件填充至150K用于本研究?；蚪M數(shù)據(jù)通過PLINK1.9軟件質控，標準如下：1)個體檢出率大于90%；2)剔除最小等位基因頻率(minor allele frequency，MAF)小于5%的位點；3)設定哈迪-溫伯格平衡檢驗(Hardy-Weinberg equilibrium test)閾值為10。最終保留了111 181個SNPs位點。

1.2 逆回歸育種值計算

為了計算繁殖性狀的逆回歸育種值(De-regressed proof，DRP)，首先使用動物(重復力)模型BLUP估計各性狀的個體EBVs，然后基于VanRaden等和Guarini等提出的方法計算DRPs：

其中，為逆回歸育種值；為父母均值(parent average，PA)的有效記錄貢獻值(effective record contribution，ERC)；為個體的ERC；為PA的可靠性；為EBV的可靠性；=(1-)。個體DRP的可靠性()則由該公式估計：(+(1-))。

1.3 統(tǒng)計模型

本研究所使用的BLUP、GBLUP和ssGBLUP方法的數(shù)學模型形式相同：

=b+a+(c)+

式中，矩陣的構造為：α+β，由VanRaden提出的方法構建，本研究將參數(shù)α和β取默認值0.95和0.05；是有基因型個體的系譜親緣關系矩陣；τ和ω為和矩陣的權重(weight)。文獻表明τ和ω的取值會影響基因組預測的效果，且具有性狀特異性。因此，為了尋找最佳的權重，本試驗設定τ從0.5至1.5、ω從0至1的梯度變化。

本研究所有計算均基于BLUPF90軟件。先使用AIREML90模塊，平均信息約束最大似然法(average information restricted maximum likelihood，AIREML)估計各性狀的方差組分，再使用BLUPF90模塊構建和矩陣進行基因組預測。

1.4 基因組預測驗證群

對于每個性狀，本研究將使用3套數(shù)據(jù)集進行基因組預測及驗證：1)總數(shù)據(jù)集①：經(jīng)過NAV標準質控后的數(shù)據(jù)，用于估計EBVs和計算DRPs，后者作為基因組育種值預測反應變量和預測模型驗證的及基準值；2)母牛驗證數(shù)據(jù)集②：基于數(shù)據(jù)集①和母牛出生時間劃分參考群和驗證群，驗證群為母牛；3)公牛驗證數(shù)據(jù)集③：在數(shù)據(jù)集①的基礎上剔除最近4年的母牛(女兒)記錄，再按個體出生時間劃分參考群和驗證群，驗證群為公牛。

在數(shù)據(jù)集②中，根據(jù)Su等的方法剔除了驗證群母牛在參考群中的近親，避免信息二次利用。在數(shù)據(jù)集③中，驗證公牛需滿足在參考群內無女兒和近親，但在數(shù)據(jù)集①中至少有10個女兒。各數(shù)據(jù)集的參考群都包括了公牛和母牛，不同性狀參考驗證群劃分情況具體見表2和表3。

表2 GBLUP方法參考驗證群劃分情況Table 2 Structure of training and validation data sets for GBLUP method

表3 BLUP和ssGBLUP方法參考驗證群劃分情況Table 3 Structure of training and validation data sets for BLUP and ssGBLUP

驗證群體通過各方法所預測的(基因組)估計育種值((G)EBVs)和基準值(DRPs)的皮爾遜相關系數(shù)(Pearson correlation coefficient)作為準確性，(G)EBVs對DRPs的回歸系數(shù)作為無偏性。以這兩個參數(shù)作為評價基因組預測效果的標準。準確性越高，說明預測值與基準值越接近；無偏性越接近于1，說明(G)EBVs的尺度與DRPs越接近。

2 結 果

2.1 中國荷斯坦牛繁殖性狀遺傳參數(shù)

通過AIREML估計得3個繁殖性狀的遺傳力分別為0.03、0.05和0.08。其中NSC的遺傳力最高，NSH次之，ICF的遺傳力最低。所有性狀皆為低遺傳力性狀，所估計的方差組分及遺傳力估計標準誤見表4。

表4 各性狀遺傳力和方差組分估計結果Table 4 Heritabilities and variance components for each trait

2.2 信息矩陣權重取值對ssGBLUP預測效果的影響

針對研究所選3個繁殖性狀，在不同的驗證群下篩選ssGBLUP方法的最優(yōu)權重組合(τ和ω)，預測準確性結果如圖1所示。對于所有性狀，不同驗證群下的預測結果有一定差異。當τ為0.5時，所有性狀使用任意驗證群皆可獲得較高的預測準確性，且隨著τ的增大，所有性狀的預測準確性皆有不同程度的下降。對于ω最優(yōu)取值，除了基于公牛驗證群的NSH外(圖1e，ω=1準確性最高)，其余預測皆在ω=0時達到最高預測準確性。

a.母牛驗證群下ssGBLUP對ICF的預測準確性；b.母牛驗證群下ssGBLUP對NSH的預測準確性；c.母牛驗證群下ssGBLUP對NSC的預測準確性；d.公牛驗證群下ssGBLUP對ICF的預測準確性；e.公牛驗證群下ssGBLUP對NSH的預測準確性；f.公牛驗證群下ssGBLUP對NSC的預測準確性a.Prediction accuracies of ssGBLUP for ICF of validation cows; b. Prediction accuracies of ssGBLUP for NSH of validation cows; c. Prediction accuracies of ssGBLUP for NSC of validation cows; d. Prediction accuracies of ssGBLUP for ICF of validation bulls; e. Prediction accuracies of ssGBLUP for NSH of validation bulls; f. Prediction accuracies of ssGBLUP for NSC of validation bulls圖1 不同權重組合及驗證群下各繁殖性狀ssGBLUP預測準確性Fig.1 Prediction accuracies of ssGBLUP for each fertility trait using various combinations of weights and validation populations

各性狀在不同驗證群及權重組合下，ssGBLUP的無偏性(b)結果如圖2所示。對于ICF和NSC兩個經(jīng)產(chǎn)牛性狀，母牛驗證的b值相比公牛驗證更加接近1，而NSH在兩種驗證群下b值均在1附近波動。當τ減小時，所有性狀母牛驗證下的b值呈增加趨勢。其中，NSH的增加幅度最大，ICF的增加幅度最小，且b值最大可能超過1(圖2a和2b)。對于公牛驗證群，各性狀b值隨τ增加的變化趨勢不一致(圖2d-2f)。當τ取值最小時，NSH和NSC的b值最接近1，而ICF則完全相反。對于ω，無論驗證群和性狀，b值隨增大皆呈下降趨勢?；赽值，ICF、NSH和NSC在母牛驗證群下最優(yōu)權重組合分別為：τ=1.3和ω=0，τ=0.5和ω=0.4以及τ=0.5和ω=0；在公牛驗證群下最優(yōu)權重組合分別為：τ=1.5和ω=0，τ=1.3和ω=0.8以及τ=0.5和ω=0。

a.母牛驗證群下ssGBLUP對ICF的預測無偏性；b.母牛驗證群下ssGBLUP對NSH的預測無偏性；c.母牛驗證群下ssGBLUP對NSC的預測無偏性；d.公牛驗證群下ssGBLUP對ICF的預測無偏性；e.公牛驗證群下ssGBLUP對NSH的預測無偏性；f.公牛驗證群下ssGBLUP對NSC的預測無偏性a.Unbiasedness of ssGBLUP for ICF of validation cows; b. Unbiasedness of ssGBLUP for NSH of validation cows; c. Unbiasedness of ssGBLUP for NSC of validation cows; d. Unbiasedness of ssGBLUP for ICF of validation bulls; e. Unbiasedness of ssGBLUP for NSH of validation bulls; f. Unbiasedness of ssGBLUP for NSC of validation bulls圖2 不同權重組合及驗證群下各繁殖性狀ssGBLUP預測無偏性Fig.2 Unbiasedness of ssGBLUP for each fertility trait using various combinations of weights and validation populations

2.3 不同驗證群及方法下基因組預測效果比較

通過不同驗證群比較傳統(tǒng)BLUP、GBLUP和ssGBLUP 三種方法對不同繁殖性狀的基因組預測準確性和無偏性的結果見表5，其中ssGBLUP選擇了使無偏性最接近1的權重組合。母牛驗證群下，BLUP、GBLUP和ssGBLUP對3個繁殖性狀的平均估計準確性分別為0.25、0.24和0.36，平均無偏性離差分別為0.263、0.267和0.017。對于各性狀，ssGBLUP較BLUP和GBLUP方法準確性分別提升了35.7%～65.2%和27.6%～65.2%。ssGBLUP的平均b值與1相差0.017，優(yōu)于BLUP(0.263)和GBLUP(0.267)。對于ICF，GBLUP的準確性低于BLUP，但b值更接近1；對于NSH，BLUP與GBLUP的準確性相同，且BLUP的b值更優(yōu)；對于NSC，GBLUP的預測準確性及無偏性都優(yōu)于BLUP。公牛驗證群下的結果也類似，GBLUP和ssGBLUP的預測準確性和無偏性均優(yōu)于BLUP，且ssGBLUP的準確性較GBLUP提升24.4%～70.6%。除了ICF外，ssGBLUP的無偏性離差均較GBLUP有所下降。

表5 各性狀基因組預測準確性和無偏性Table 5 Accuracies and unbiasedness of genomic prediction for each trait

對于BLUP方法的結果，使用母牛驗證群時準確性高于使用公牛驗證群的準確性，且無偏性更加接近于1?；诠ｒ炞C群的GBLUP和ssGBLUP所得平均準確性與基于母牛驗證群的平均準確性都相差0.07。除了ICF外，NSH和NSC的準確性都是在使用公牛驗證群時較高。就無偏性而言，母牛驗證下GBLUP和ssGBLUP方法的平均b值較公牛驗證群更接近于1。

3 討 論

3.1 中國荷斯坦牛繁殖性狀遺傳評估分析

本研究系統(tǒng)地分析了不同驗證群及預測方法下，荷斯坦牛3個繁殖性狀的基因組預測效力。首先，使用動物模型和重復力模型BLUP對各性狀的遺傳參數(shù)進行了估計，結果表明遺傳力估計值在0.03～0.08之間，符合該類性狀的正常范圍。對于配種次數(shù)性狀，本研究中青年牛群體(NSH)的遺傳力估計值要低于經(jīng)產(chǎn)牛群體(NSC)。Guarini等對加拿大荷斯坦牛群繁殖性狀的研究也表明NSC的遺傳力(0.08)要高于NSH的遺傳力(0.02)，與本研究結果一致。

3.2 權重取值對基因組預測結果的影響

使用ssGBLUP方法構建陣時，使用不同參數(shù)組合改變陣和陣的組合，能夠影響基因組預測的準確性，目前主要有3種構造陣的方法。其一為本研究所評估的方法，參數(shù)包括α、β、τ和ω(本研究中參數(shù)α和β取值為0.95和0.05)；方法二為方法一的特殊形式，即將τ和ω的取值都定為1，僅考慮參數(shù)α和β的取值；方法三則進一步對陣進行調整。構造陣時，不同的參數(shù)取值會影響基因組預測結果的準確性和無偏性，甚至影響算法的收斂容易度。針對τ和ω的取值，本研究結果表明，盡管各性狀最佳權重取值不同，但整體而言τ越小，所有性狀的基因組預測準確性和b值都越高，且大部分性狀在ω=0時能獲得更高的準確性和b值(圖1、圖2)。其中，ω對準確性和b值的變化規(guī)律與前人研究相同，但τ的變化趨勢則與之相反。這可能是由于該先前研究的性狀為產(chǎn)奶性狀，與本研究的繁殖性狀不同，同時其系譜和表型記錄數(shù)量(超過百萬)遠高于本研究。另外，其使用的SNP芯片密度為50K，而本研究中的150K芯片大部分是由50K芯片填充所得，部分SNP標記填充錯誤可能會影響矩陣的構造，進而影響最佳權重系數(shù)的取值。Guarini等使用公牛驗證群對加拿大荷斯坦群體繁殖性狀的GS研究表明，當τ=1且ω=0.6或0.7時，ICF、NSH和NSC的基因組預測效果最好。導致結果不一致的部分原因可能是其研究所使用的芯片密度為50K，且參考群體中擁有系譜個體與基因型數(shù)據(jù)個體的比例遠高于本研究。由于陣是由基因組和系譜構造，所以研究中參考群體的數(shù)量、基因測定個體的數(shù)量、系譜深度和SNP芯片密度皆會影響權重的取值。綜上，使用ssGBLUP方法之前，確定τ和ω的取值能夠在一定程度上提升基因組預測的效力，但該權重的取值會因研究群體和所用SNP芯片的不同而產(chǎn)生差異。

3.3 不同驗證群及方法對基因組預測結果的影響

實際生產(chǎn)中，GS一般用于提前選擇青年公牛以減少育種的世代間隔，進而加快遺傳進展。以此為前提，一般會對具有后裔測定信息的公牛和青年公牛進行基因組測定，進而預測青年公牛的基因組育種值，達到提前選種的目的。但近年來，隨著芯片測定技術的發(fā)展，許多國家也對母牛進行基因組測定，并納入GS的參考群體來擴充參考群，同時利用基因組育種值選擇優(yōu)秀母?；蛘吖Ｄ赣H。本研究采用驗證群，公牛驗證群體個體數(shù)較少，但DRP可靠性較高；而驗證母牛群體個體數(shù)遠遠大于公牛驗證群體，但DRP可靠性遠低于公牛群體(表2)。驗證結果顯示,預測準確性和無偏性在兩種驗證群體有差別，最主要原因是DRP的準確性。本研究用(G)EBV和DRP的相關系數(shù)表示預測準確性，但該相關系數(shù)又與DRPs本身的估計準確性相關，而公牛驗證群體和母牛驗證群體DRPs準確性差別很大。實際上，本研究所使用的預測準確性僅解釋了對表型(DRPs)的預測效果，若要計算對估計育種值的預測效果，還需進一步除以表型值(例如DRPs)的準確性。然而，對于低遺傳力性狀(如繁殖性狀)，DRPs準確性較小(表2)，作為被除數(shù)時可能會因為DRPs準確性估計誤差導致異常結果(例如超過1)。因此，對低遺傳力性狀基因組預測的相關研究中，經(jīng)常使用表型值與預測值的皮爾遜相關系數(shù)作為準確性參數(shù)。因此，本研究兩種驗證群結果不能直接比較，但可以在同一驗證群體比較不同模型的預測能力，以及不同模型預測能力在兩種驗證群體的排序是否一致。

在奶牛實際育種中，受到成本限制，一般僅有少部分個體會進行基因組測定。使用GBLUP方法時，需要被評估個體有基因組信息，因此大量的表型數(shù)據(jù)及系譜信息無法利用。而ssGBLUP則可有效地將基因組、系譜和表型信息同時利用，無基因組信息個體也可以借助基因型信息進行評估。本研究結果表明，ssGBLUP方法的基因組預測效力最好，表現(xiàn)為較高的預測準確性及更好的無偏性。在部分性狀中，GBLUP方法的準確性和無偏性可能略遜于BLUP方法(ICF和NSH)，其原因可能是GBLUP僅利用少部分基因組測定個體，預測值誤差較大(特別是針對低遺傳力性狀)。該結果反映出ssGBLUP方法在GS體系中具有一定的優(yōu)勢。在不同群體及不同物種的研究中，ssGBLUP較BLUP和GBLUP方法都得到了較高的預測準確性和無偏性。因此，在本研究比較的3種方法中，ssGBLUP是對中國荷斯坦牛繁殖性狀進行基因組選擇的最佳方法。

4 結 論

本研究以中國荷斯坦牛為對象，利用BLUP、GBLUP和ssGBLUP方法對部分重要繁殖性狀進行了基因組預測效果分析。研究發(fā)現(xiàn)構建陣時，權重組合對ssGBLUP的預測效果存在影響。因此，正式進行預測前嘗試不同權重參數(shù)取值，并選擇最優(yōu)的權重組合能夠在一定程度上提升ssGBLUP方法的預測性能。另外，研究結果表明所使用的3種方法里，ssGBLUP的預測結果具有較高準確性和無偏性，因此是中國荷斯坦牛繁殖性狀基因組選擇的最佳方法之一。

感謝中國農業(yè)大學動物科技學院奶牛分子數(shù)量遺傳學團隊提供部分芯片數(shù)據(jù)，以及團隊老師在性狀數(shù)據(jù)整理分析及構建基因組預測模型等方面提供的幫助。感謝中國農業(yè)大學動物科技學院459實驗室同學胡麗蓉、張海亮、婁文琦、?，?、竹磊和安濤等人以及新疆農業(yè)大學動物科技學院魏趁和張夢華博士參與血樣采樣，并在文章撰寫過程中提出寶貴的意見。