一種適用于單胎動物多品種選擇的基因組關系矩陣

張洪志，任端陽，安麗霞，喬利英，劉文忠

(山西農業大學動物科學學院，太谷030801)

最佳線性無偏預測(best linear unbiased prediction, BLUP)使用系譜信息構建分子血緣相關矩陣(A陣)，結合表型信息計算常規估計育種值(estimated breeding value, EBV)[1]。對于易測量、中高遺傳力的性狀，BLUP取得了較好的應用效果[2-3]；但是對于低遺傳力性狀、限性性狀、難以測定及生命晚期才能測定的性狀等，其應用效果較差[4-5]。

單核苷酸多態性(single nucleotide polymorphism, SNP)芯片可提供生物個體全基因組范圍的基因型信息。隨著SNP芯片技術的發展及其制作成本的降低[6]，SNP芯片正被廣泛應用于基因組選擇(genomic selection, GS)[7-8]。GS是一種全基因組范圍的標記輔助選擇[9]，用來計算畜禽個體的基因組估計育種值(genomic estimated breeding value, GEBV)。最初實施GS使用“多步”法[10]，包括基因組最佳線性無偏預測(genomic BLUP, GBLUP)和貝葉斯及其改進方法[11-12]。多步法使用EBV[13]或逆回歸育種值[14-15]作為偽表型值，通過偽表型值和基因組信息預測GEBV，從而用于選種[16]。由于使用偽表型值及一些參數估值有偏，多步法的可靠性較低[17]。一步法(single-step GBLUP, ssGBLUP)將基因組信息添加到由系譜信息構建的A陣中構建基因組關系矩陣(H陣)，結合表型信息直接預測未基因分型個體的EBV及基因分型個體的GEBV。ssGBLUP的計算步驟較簡單且準確性不比多步法低[17-18]，已在奶牛育種中廣泛應用[19-20]，這得益于在奶牛中可以組建足夠規模的參考群體[21]。

對于肉牛、綿羊等家畜，通常群體規模較小，實施單品種GS(within-breed GS, wbGS)效果較差[22]。針對這一科學問題，提出并逐漸開始研究多品種GS(across-breeds GS, abGS)[23-24]。但是，目前缺少針對abGS的統計方法，使得abGS的應用效果較差[25]。

本研究旨在提出一種針對abGS的基因組關系矩陣構建方法，并通過模擬研究對該方法的使用效果進行評估，為abGS的實際應用提供新思路。

1 材料與方法

1.1 數據模擬

1.1.1 群體結構 使用QMsim軟件[26]模擬數據。模擬4個數量性狀基因座(quantitative trait loci, QTL)數水平(50、100、500、1 000)及3個遺傳力(heritability,h2)水平(0.1、0.3、0.5)，共產生12個 性狀。模擬牛的29對常染色體的50K SNP標記信息。對整個模擬過程重復10次，以減少抽樣誤差。模擬的群體結構和基因組參數見表1和表2。

表1 群體結構模擬參數

表2 基因組模擬參數

第一步，模擬歷史群體。起始群體規模為3 000頭，至500代時，群體規模遞增到40 000頭，至800代時，群體規模遞減到30 000頭，以得到與實際肉牛群體相似的連鎖不平衡(linkage disequilibrium, LD)程度和突變-漂變平衡[27]。整個歷史群體創建階段，公母數目保持一致，采用隨機交配。

第二步，模擬5個當代群體，每個群體模擬13個世代，視為5個品種。5個當代群體相互獨立，分別從歷史群體最后一代隨機選擇10頭公牛，100頭母牛作為當代群體的初始群體。模擬群體的前9代作為擴繁群體，公牛的淘汰率為60%，公、母牛的增長率為35%；10～13代公、母牛淘汰率分別為60%和20%，公母牛數量不再增長。公、母牛隨機交配，每頭母牛每年生產1頭后裔，公、母后裔各半。選擇估計育種值(estimated breeding value, EBV)高的個體留種。

使用Henderson動物模型[28]BLUP估計EBV。模擬性狀表型方差標準化為1，個體的真實育種值為所有QTL加性遺傳效應之和，個體的表型值為真實育種值與隨機殘差之和。模擬數據不考慮固定效應。

1.1.2 基因組 基于牛的真實基因組信息[29-30]，模擬29對常染色體，總長2 333 cM。為了模擬接近真實情況的SNP-QTL的LD，依據真實染色體長度及數目模擬基因組信息。模擬50 000個SNPs標記，假定SNP標記均勻分布在全基因組范圍內，標記無效應。對QTL數設(50、100、500、1 000)4個水平，QTL的最小等位基因頻率(minor allele frequency, MAF)為0.1。假定QTL在全基因組范圍內隨機分布，QTL效應值服從伽馬分布(形狀參數為0.4)。為了創建突變-漂變平衡，設置歷史群體中QTL和SNP標記的回復突變率為10-5。

1.1.3 參考群體和候選群體 選擇5個品種的第8～12代個體作為聯合參考群體，5個品種的第13代作為聯合候選群體。5個群體中挑選10～12代 有后裔的個體保留基因型信息，第13代隨機選擇250個個體(每個品種50個)保留基因型信息。對于基因分型個體，對其模擬不同比例(20%、40%、60%、100%)的系譜缺失，即個體的父號或母號未知。之前大部分GS的模擬研究只將無表型有基因分型的個體作為候選群體，Christensen用真實數據進行交叉驗證發現[31]，ssGBLUP模型對候選群體中基因分型和未進行分型的個體均能提高GEBV的準確性，在GS實施階段，可以對部分候選群體進行基因分型就能整體提高GEBV的準確性。本研究只對部分候選群體進行基因分型，更有實踐意義。

1.2 統計模型

使用ssGBLUP模型：

y=Xb+Zu+e

其中，y為表型值向量；b為固定效應向量；u為育種值向量；e為隨機殘差效應向量。X和Z分別為固定效應和育種值向量的關聯矩陣。

混合模型方程組為：

H陣為：

H=

H陣的逆表示為：

其中，Gw=(1-w)G+wA22，A陣為基于系譜構建的分子血緣相關矩陣，A22陣為基因分型個體對應的分子血緣相關矩陣，G為利用標記信息構建的基因組關系矩陣，w為對A22的加權系數，本研究中，w取值0和0.05(其中：newGw=0.95newG+0.05A22，Gw=0.95G+0.05A22)。使用兩種G陣構建H陣：

(1)常規/標準G陣：

(2)新型G陣(newG)：

其中，M陣的維度為m行n列，m為基因分型個體數，n為標記數。標記位點基因型用(0，1，2)表示，0和2代表兩種純合基因型，1代表雜合基因型。P陣中的元素為標記的第二等位基因頻率，Q陣的元素為標記的第二純合基因型頻率，pj為第j個標記的第二等位基因頻率。

對于聯合群體中符合哈代-溫伯格平衡(Hardy-Weinberg equilibrium, HWE)的位點，若Q0為該位點的第二純合基因型頻率，p0為第二等位基因頻率，則：

對于聯合群體的非HWE(Hardy-Weinberg disequilibrium, HWDE)位點，若Q1為第二純合基因型頻率，p1為第二等位基因頻率，則：

1.3 GEBV的預測準確性

使用Gmatrix軟件[32]構建常規G陣，使用R語言自編程序構建新型G陣，使用DMU軟件[33]計算GEBV。

在不同h2、QTL數及基因分型個體系譜不同缺失比例下，通過候選群體GEBV與真實育種值(true breeding value, TBV)的相關性(correlation coefficient,cor)來衡量應用不同G陣的預測準確性：

其中，Cov為協方差(covariance)，Var為方差(variance)。

2 結 果

2.1 參考群體和候選群體

本研究模擬了從相同的歷史群體得到的5個獨立的當代群體(品種)，群體結構模擬結果見表3。使用5個品種的第8～12代作為聯合參考群體，個體數為31 935，其中10～12代中共有5 865個個體有基因型信息。將5個群體的第13代作為聯合候選群體。候選群體共有7 445個個體，共有250個個體有基因型信息。

表3 模擬數據的當代群體結構

2.2 基因分型個體系譜無缺失時各模型的預測準確性

基因分型個體系譜無缺失時，QTL數和h2變化時各模型的預測準確性見表4。可見，QTL數為50，h2為0.3和0.5，及QTL數為100、500和1 000，各種h2時，newGw陣和Gw陣準確性相同，且較newG陣高0.0%～0.2%。QTL數為50，h2為0.1時，newG陣較newGw陣和Gw陣準確性高0.1%。在各QTL數及h2情況下，G陣準確性最低。

表4 基因分型個體系譜無缺失時不同QTL數和遺傳力下各模型的預測準確性(10次重復的均值)

2.3 基因分型個體系譜部分缺失時各模型的預測準確性

基因分型個體系譜不同比例缺失(20%、40%、60%)時，QTL數和h2變化時各模型的預測準確性見表5。可見，各QTL數和h2情況下，newG陣的準確性最高，newGw陣與Gw陣準確性有略微差異(-0.4%～0.3%)，G陣準確性最低。隨著基因分型個體系譜缺失比例增加，在各QTL數及h2情況下，除G陣外，其他3種G陣準確性逐漸下降，相對于newG陣，newGw陣和Gw陣準確性降低幅度較大。系譜缺失對G陣無明顯影響，隨著缺失比例增加，G陣準確性甚至會增加。

表5 基因分型個體系譜不同缺失比例下不同QTL數和遺傳力下各模型的預測準確性(10次重復的均值)

不同基因分型個體系譜缺失比例下，newG陣較Gw陣的優勢(預測準確性差值)變化見圖1。可見，對于所有QTL數，隨著h2降低，newG陣優勢逐漸增加。這表明，在基因分型個體系譜部分缺失時，newG陣較Gw陣在低遺傳力時優勢更明顯。

圖1 基因分型個體系譜缺失(20%(A)，40%(B)，60%(C))時newG陣與Gw陣GEBV準確性差值

2.4 基因分型個體系譜全部缺失時各模型的預測準確性

基因分型個體系譜全部缺失，QTL數和h2變化時各模型的預測準確性見表6。可見，各QTL數及h2情況下，newG陣、newGw陣、Gw陣的準確性有輕微差異(±0.1%)，且均較G陣高。

表6 基因分型個體系譜完全缺失時不同QTL數和遺傳力下各模型的預測準確性(10次重復的均值)

結合表3～表6可見，在不同比例的基因分型個體系譜缺失情況下，對于所有G陣，隨著h2增加，預測準確性增加，而預測準確性隨QTL數并沒有明顯變化規律。

3 討 論

GS方法的模擬研究是其發展的重要過程，新方法實際應用之前，需要使用模擬數據對其電腦程序、算法及可靠性進行評估，為其實際應用提供參考。然而，數據模擬軟件的設計及運行軟件時的參數設置與真實情況往往有一定差異，使得模擬研究與實際應用的結果有偏差，如：在模擬研究中，非線性模型明顯優于GBLUP法[34]，然而使用真實數據驗證時，非線性模型較GBLUP只有輕微或沒有優勢[35]。本研究中，newG陣在模擬研究中顯示出一些優勢，而其還需要在實際數據中進一步驗證。本研究模擬肉牛的多品種聯合群體，對于其他單胎動物的多品種或同品種不同地區、不同國家的聯合群體理論上仍然適用。

試驗結果顯示，不同模型的準確性隨著h2的增大而提高，Luan等[36]使用實際數據應用GBLUP及Bayes模型時顯示出了相似的結果，由此推測，低遺傳力性狀需要更多的表型來提高GEBV預測準確性。本試驗中，不同模型的準確性隨QTL數的變化無明顯變化，在Daetwyler等[37]的模擬研究中，BayesB模型隨著QTL數減少而準確性增高，GBLUP在不同QTL數時，準確性基本無變化。由此可以推測，這可能是由于GBLUP及ssGBLUP假設所有標記解釋相同的效應方差，且無主效基因[16]，造成了本試驗中應用G陣和newG陣時對QTL數的變化不敏感。

阻礙肉牛、綿羊等畜種實施GS的因素有：1)參考群體規模少；2)群體中同品種個體較少，群體中多品種混雜[40]；3)信息記錄系統不完善，如：系譜記錄，表型記錄等[41]。因此，研究多品種聯合或同品種不同地區聯合實施GS是未來的趨勢[25]。對于abGS，各品種間等位基因頻率的差異會造成多品種聯合群體出現大量HWDE位點[42]，目前還沒有將HWDE考慮入GS的方法，本研究考慮HWDE因素構建基因組關系矩陣，在系譜部分缺失情況下具有優勢，然而，HWDE對abGS的影響還有待進一步研究。

4 結 論

本研究針對abGS的G陣構建進行了探索，提出了一種G陣構建方法。結果表明，新型G陣在不需要加權的情況下就能達到常規G陣加權時的GEBV預測準確性，在系譜部分缺失時，應用新型G陣較加權常規G陣的GEBV準確性有一些提高。