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

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

(山西農(nóng)業(yè)大學(xué)動物科學(xué)學(xué)院，太谷030801)

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

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

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

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

1 材料與方法

1.1 數(shù)據(jù)模擬

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

表1 群體結(jié)構(gòu)模擬參數(shù)

表2 基因組模擬參數(shù)

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

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

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

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

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

1.2 統(tǒng)計模型

使用ssGBLUP模型：

y=Xb+Zu+e

其中，y為表型值向量；b為固定效應(yīng)向量；u為育種值向量；e為隨機殘差效應(yīng)向量。X和Z分別為固定效應(yīng)和育種值向量的關(guān)聯(lián)矩陣。

混合模型方程組為：

H陣為：

H=

H陣的逆表示為：

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

(1)常規(guī)/標(biāo)準(zhǔn)G陣：

(2)新型G陣(newG)：

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

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

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

1.3 GEBV的預(yù)測準(zhǔn)確性

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

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

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

2 結(jié) 果

2.1 參考群體和候選群體

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

表3 模擬數(shù)據(jù)的當(dāng)代群體結(jié)構(gòu)

2.2 基因分型個體系譜無缺失時各模型的預(yù)測準(zhǔn)確性

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

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

2.3 基因分型個體系譜部分缺失時各模型的預(yù)測準(zhǔn)確性

基因分型個體系譜不同比例缺失(20%、40%、60%)時，QTL數(shù)和h2變化時各模型的預(yù)測準(zhǔn)確性見表5?？梢姡鱍TL數(shù)和h2情況下，newG陣的準(zhǔn)確性最高，newGw陣與Gw陣準(zhǔn)確性有略微差異(-0.4%～0.3%)，G陣準(zhǔn)確性最低。隨著基因分型個體系譜缺失比例增加，在各QTL數(shù)及h2情況下，除G陣外，其他3種G陣準(zhǔn)確性逐漸下降，相對于newG陣，newGw陣和Gw陣準(zhǔn)確性降低幅度較大。系譜缺失對G陣無明顯影響，隨著缺失比例增加，G陣準(zhǔn)確性甚至?xí)黾印?/p>

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

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

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

2.4 基因分型個體系譜全部缺失時各模型的預(yù)測準(zhǔn)確性

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

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

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

3 討 論

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

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

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

4 結(jié) 論

本研究針對abGS的G陣構(gòu)建進行了探索，提出了一種G陣構(gòu)建方法。結(jié)果表明，新型G陣在不需要加權(quán)的情況下就能達到常規(guī)G陣加權(quán)時的GEBV預(yù)測準(zhǔn)確性，在系譜部分缺失時，應(yīng)用新型G陣較加權(quán)常規(guī)G陣的GEBV準(zhǔn)確性有一些提高。