程 玲，張心菲，張鑫鑫，張衛(wèi)華，林元震 (1華南農業(yè)大學 林學與風景園林學院，廣東 廣州 510642：2 廣東省森林植物種質創(chuàng)新與利用重點實驗室，廣東 廣州510642： 廣東省林業(yè)科學研究院，廣東 廣州 510640)

與農作物相似，林木選育的良種也存在地域適用性，因此林木遺傳試驗應進行多地點試驗(Multi-environment trial,MET)。多地點試驗分為3個階段：初試階段、中試階段和大規(guī)模試驗階段，其中參試材料數(shù)一般隨著試驗規(guī)模的增大而減少，但試驗地點數(shù)恰好相反[1]。在多地點試驗中，一個突出的問題是林木基因型與環(huán)境之間往往存在顯著的交互作用(Genotype by environment interaction, GEI)，因此如何準確評估這種交互作用對于后續(xù)林木良種的選育和推廣至關重要。目前，農業(yè)上大多使用聯(lián)合回歸法[2]、主效可加互作可乘模型(Additive main effects multiplicative interaction,AMMI)[3]和基因型主效加基因型-環(huán)境互作效應(Genotype main effect plus genotype-by-environment interaction,GGE)雙標圖[4]分析多點試驗結果，并據(jù)此來進行品種評價、試驗點評價和品種生態(tài)區(qū)劃分，其中GGE雙標圖越來越受到科研人員的關注[5]。迄今，GGE雙標圖在林木良種選育上的應用仍然很少，僅在少數(shù)樹種如輻射松[6]、楊樹[7]和樂昌含笑[8]上有報道。

雖然AMMI和GGE雙標圖在農作物上應用廣泛，但這些模型均直接通過表型數(shù)據(jù)進行分析，而且存在一些限制[5]：1)上述方法均只限定于固定效應模型；2)假定各地點誤差同質；3)要求數(shù)據(jù)是平衡的。但對于林木試驗地而言，首先試驗地點誤差同質幾乎不可能，現(xiàn)有研究證實林地存在各種空間變異[1]；其次，試驗結果往往會存在缺株或缺區(qū)數(shù)據(jù)等現(xiàn)象；再者，作為重要遺傳參數(shù)之一的育種值,需要隨機效應模型才可估算。因此，這些問題制約著AMMI和GGE雙標圖在林木研究中的應用。針對這些問題，本研究提出采用空間變異結合因子分析法獲取林木基因型在各試驗地的無偏預測值(Best linear unbiased prediction,BLUP)，然后再利用GGE雙標圖進行林木基因型評估、試驗地評估和生態(tài)區(qū)劃分，以避免固定效應模型、試驗地點同質和缺失數(shù)據(jù)的約束，為GGE雙標圖在林木良種選育上的應用提供借鑒。

1 方 法

1.1 BLUP的統(tǒng)計模型

多地點試驗的混合模型為：

yijk=μ+Si+SGi(j)+eijk。

(1)

式中:yijk為個體性狀表型值；μ為總體均值；Si為第i個地點的固定效應；SGi(j)為第i個地點與第j個基因型的隨機交互效應；eijk為剩余殘差。

使用因子分析法[9]擬合SGi(j)效應時，其方差協(xié)方差矩陣可寫為：

G=(Γ×?！?Ψ)?I。

(2)

式中：G為SGi(j)效應的方差矩陣；Γ為因子載荷矩陣；?！錇橐蜃虞d荷矩陣的轉置矩陣；Ψ為特殊方差矩陣；I為單位矩陣；?為矩陣的Kronecker乘法。

(3)

∑ρ自相關矩陣為：

(4)

上述BLUP過程采用ASReml軟件[11]進行分析，具體方法參照文獻[1]。

1.2 GGE雙標圖的統(tǒng)計模型

參照文獻[6]，GGE雙標圖的統(tǒng)計模型為：

yij-μ-βj=λ1γi1δj1+λ2γi2δj2+εij。

(5)

式中：yij為第j個地點第i個基因型的均值；μ為總體均值，為試驗地點的主效應；βj為第j個地點所有基因型的均值；λ1、λ2為第1個和第2個主成分的特征值；γi1、γi2為第i個基因型在第1個和第2個主成分的特征向量；δj1、δj2為第j個地點在第1個和第2個主成分的特征向量；εij為剩余殘差。

GGE雙標圖的分析采用R軟件程序包GGEBiplotGUI[12]實現(xiàn)。參數(shù)設置時Scaled選擇0(非標準化)，Centerd選擇G+GE，SVP選擇Symmetrical。

1.3 測試數(shù)據(jù)

測試數(shù)據(jù)來自文獻[1]?；鹁嫠傻幕蛐陀?6種(1～36)，試驗地點6個(S1～S6)，試驗設計為拉丁方設計，每個試驗地設3次重復，每個重復設6個區(qū)組，測量性狀為種子產量。其中，除了第3個地點為9行12列外，其余均為6行18列。

1.4 原始數(shù)據(jù)和BLUP數(shù)據(jù)的比較

原始數(shù)據(jù)和BLUP數(shù)據(jù)的熱圖采用R軟件程序包AAfun[13]生成，方差分析采用R軟件aov函數(shù)[1]完成。

2 結果與分析

2.1 原始數(shù)據(jù)的BLUP分析結果

為克服試驗地誤差異質對GGE雙標圖分析的影響，本研究采用空間分析模型擬合每個試驗地的R殘差，目的是：(1)求證每個試驗地誤差是否異質；(2)獲取每個試驗地測量性狀的BLUP值，以用于后續(xù)的GGE雙標圖分析?？臻g變異擬合R殘差的參數(shù)BLUP結果如表1所示。

表1 火炬松36個基因型的種子產量原始數(shù)據(jù)的BLUP分析結果Table 1 BLUP result of parameters for original data of seed yield of 36 Pinus teada genotypes

注：行、列自相關值的上標值是t檢驗統(tǒng)計量，為相關估計值與其標準誤的比值。

Note:The superscript iststatistic autocorrelation,the ratio of estimated values to standard error.

由表1可知，各地點的環(huán)境誤差差異比較大，對于隨機誤差，地點S1最大，而地點S2和S5的估計值為0；對于空間誤差，所有地點均存在，其中地點S2和S4的比較大。此外，除地點S5以外，其余地點都存在顯著的行或列自相關性(t>1.5時，代表相關值在α=0.05水平上顯著)，說明試驗地存在顯著的空間變異模式。上述結果表明，各試驗地點環(huán)境誤差不同，空間變異的擬合結果有效，獲取的測量性狀BLUP值可靠。

2.2 原始數(shù)據(jù)與BLUP數(shù)據(jù)的比較

原始數(shù)據(jù)和BLUP數(shù)據(jù)的熱圖結果如圖1所示。

S1～S6.表示6個試驗地；R1～R3.分別為重復1～3；圖中的小方框代表區(qū)組S1-S6.Six trail sites;R1-R3.repeats 1-3.Small squares of the pictures represent blocks圖1 火炬松種子產量原始數(shù)據(jù)(左)和BLUP數(shù)據(jù)(右)的熱圖Fig.1 Heatmap of original data (left) and BLUP data (right) of Pinus teada seed yield

圖1顯示，試驗地間的種子產量差異明顯，而且原始數(shù)據(jù)和BLUP數(shù)據(jù)之間的差異也比較明顯。經(jīng)過空間變異的校正后，BLUP數(shù)據(jù)中數(shù)據(jù)點之間的變異幅度縮小，而且BLUP數(shù)據(jù)的空間變異模式更為清晰，例如地點S2，顏色接近的幾乎連成片(存在空間變異的一種典型特征)，但原始數(shù)據(jù)并未呈現(xiàn)出這種趨勢。上述結果進一步表明，各地點的確存在空間變異，與表1的分析結果一致。

原始數(shù)據(jù)和BLUP數(shù)據(jù)的方差分析結果(表2)顯示，對于原始數(shù)據(jù)，基因型和地點的效應均達極顯著水平(P<0.001)，基因型與地點的互作效應達顯著水平(P<0.05)，地點與重復的互作效應達極顯著水平(P<0.001)，說明上述因子對種子產量影響顯著；對于BLUP數(shù)據(jù)，基因型、地點及其互作效應的F值均大于原始數(shù)據(jù)的F值，而且均達極顯著水平(P<0.001)，但地點與重復的互作效應不顯著(P>0.05)，這是由于空間變異的擬合削弱了重復因子的效應，此外，誤差的均方顯著降低?；蛐?、地點及其互作效應對平方和的解釋百分比之和，原始數(shù)據(jù)為67.81%，BLUP數(shù)據(jù)為88.52%，說明基因型、地點及其互作效應引起的變異是種子產量變異的主要來源，即可以使用GGE雙標圖進行后續(xù)分析。上述結果表明，對于產量變異的解釋能力，BLUP數(shù)據(jù)比原始數(shù)據(jù)更高，表明BLUP數(shù)據(jù)用于GGE雙標圖的分析更為可靠。

表2 火炬松種子產量原始數(shù)據(jù)和BLUP數(shù)據(jù)的方差分析Table 2 Analysis of variance for original data and BLUP data of Pinus teada seed yield

注：*.效應達顯著水平(P<0.05)，**.效應達極顯著水平(P<0.001)。

Note: *. significant (P<0.05), **. very significant (P<0.001).

2.3 原始數(shù)據(jù)和BLUP數(shù)據(jù)與GGE雙標圖結合分析結果的比較

2.3.1 試驗地分組結果 GGE雙標圖可以展示試驗地點間的關系，試驗地點向量之間的夾角余弦值為試驗點的遺傳相關系數(shù)，夾角越小，則相關值越大，當夾角小于90°時為正相關，而大于90°時為負相關。地點間關系的分析結果(圖2)顯示，對于原始數(shù)據(jù)，地點S1、S2和S5之間高度相關，其中地點S1、S2的相關系數(shù)接近1，地點S4和S6間也高度相關；對于BLUP數(shù)據(jù)，地點S1、S2和S5之間仍高度相關，但地點S1和S2間的相關系數(shù)明顯小于1，地點S4和S6間的相關程度也比原始數(shù)據(jù)弱。由此可見，BLUP數(shù)據(jù)可以校正試驗地點間的相關關系。此外，兩個主成分的方差解釋百分比和，原始數(shù)據(jù)為70.51%，BLUP數(shù)據(jù)為87.50%，與表2結果一致，進一步說明BLUP數(shù)據(jù)GGE雙標圖的分析結果更為可靠。

GGE雙標圖的Which Won Where/What功能將最外圍的基因型連成一個多邊形，通過原點對多邊形每條邊做垂線，據(jù)此將試驗地點分組，并得到各分組內的優(yōu)秀基因型。分組結果(圖3)顯示，不論是原始數(shù)據(jù)還是BLUP數(shù)據(jù)，6個試驗地點均分為2組，地點S4和S6為第1組，其余4個地點為第2組，分組結果與因子分析的多點結果[1]一致，說明GGE雙標圖的試驗地分組結果是可靠的。優(yōu)秀基因型篩選結果顯示，2種數(shù)據(jù)在第1，2組試驗地的種子最高產基因型分別是10和28。上述結果表明，雖然2種數(shù)據(jù)的試驗地點分組一致，但分組內的優(yōu)秀基因型存在差異。

2.3.2 試驗地的區(qū)分力和代表性 試驗地的選擇與良種選育的可靠性直接相關，一個理想的試驗地應具備較強的區(qū)分力和代表性。GGE雙標圖的Discrimitiveness VS.representativness功能可以展示試驗地點的區(qū)分力和代表性。試驗地評估結果見圖4，圖中各試驗點與原點的虛線向量長度代表試驗點的區(qū)分力，試驗點向量與平均環(huán)境軸的夾度反映的是試驗點的代表性，夾角越小表示試驗點的代表性越強，如果夾角為鈍角，則該試驗地不適合作為試驗點。由圖4可知，對于原始數(shù)據(jù)，區(qū)分力最好的是地點S4和S3，代表性最好的是地點S5、S1和S2；對于BLUP數(shù)據(jù)，區(qū)分力最好的是地點S3、S4和S1，代表性最好的是地點S5、S1和S2。綜合起來，原始數(shù)據(jù)最好的試驗點是地點S5，BLUP數(shù)據(jù)最好的是地點S1。

S1～S6.表示6個試驗地；第一主成分、第二主成分得分括號內的數(shù)據(jù)分別表示第一主成分和第二主成分解釋的方差百分比。下圖同S1-S6.Six trail sites;Data in brackets indicate the proportions explained by the first and the second principal components.The same below.圖2 基于火炬松種子產量原始數(shù)據(jù)(左)和BLUP數(shù)據(jù)(右)的試驗地點間關系Fig.2 Site relationship of original data (left) and BLUP data (right) of Pinus teada seed yield

1～36.表示36個基因型。下圖同1-36. 36 genotypes. The same below圖3 基于火炬松種子產量原始數(shù)據(jù)(左)和BLUP數(shù)據(jù)(右)的試驗地分組Fig.3 Site grouping result of original data (left) and BLUP data (right) of Pinus teada seed yield

2.3.3 供試基因型的高產性和穩(wěn)產性 GGE雙標圖的Mean VS.Stability功能可以展示供試基因型的高產性和穩(wěn)產性。不同基因型火炬松的高產性和穩(wěn)產性結果如圖5所示，圖中各基因型到平均環(huán)境軸的垂直虛線段代表各基因型在所有試驗地的平均產量和穩(wěn)產性，虛線段越長代表產量越不穩(wěn)定；與平均環(huán)境軸垂直的實線為產量總體均值，在其左側的基因型產量低于總體均值，距其越遠產量越低，在其右側的基因型產量高于總體均值，距其越遠產量越高。由圖5可知，對于原始數(shù)據(jù)，基因型28種子產量最高，其次是20，33，21和25；基因型9種子產量最低，其次是基因型32，2和31；基因型26和12種子產量接近總體均值；種子產量最不穩(wěn)定的是基因型10，3，11和28，而基因型21，30和16產量比較穩(wěn)定，其中16屬于穩(wěn)定但低產基因型。對于BLUP數(shù)據(jù)，基因型28種子產量最高，其次是33，20，21和19；基因型9種子產量最低，其次是32，2和31；基因型23種子產量接近總體均值；種子產量最不穩(wěn)定的是基因型10，其次是3和28，而基因型21，30，25，22和16產量比較穩(wěn)定，其中22和16屬于穩(wěn)定但低產基因型。綜合來看，原始數(shù)據(jù)中最理想的基因型是21，BLUP數(shù)據(jù)的理想基因型也是21，該基因型種子產量既高又穩(wěn)定。

圖4 基于火炬松種子產量原始數(shù)據(jù)(左)和BLUP數(shù)據(jù)(右)的試驗地的區(qū)分力和代表性Fig.4 Discrimination and representativeness of trial sites for original data (left) and BLUP data (right) of Pinus teada seed yield

圖5 基于火炬松種子產量原始數(shù)據(jù)(左)和BLUP數(shù)據(jù)(右)的基因型高產性和穩(wěn)定性Fig.5 Genotypic mean and stability of original data (left) and BLUP data (right) of Pinus teada seed yield

3 討 論

GGE雙標圖的開發(fā)者嚴威凱[5]指出，目前AMMI和GGE雙標圖在多點試驗分析中的主要限制有：僅限固定效應模型、各地點環(huán)境誤差同質和平衡數(shù)據(jù)，這3個因素可能正是GGE雙標圖在林木選育和推廣中應用較少的關鍵原因。但GGE雙標圖法在農業(yè)品種評價、試驗點評價和品種生態(tài)區(qū)劃分上已是先鋒方法，林木良種的選育和推廣也需要解決這3個問題，因此其在林業(yè)上的應用潛力很大?？紤]到林木試驗地環(huán)境差異往往比較大，GGE雙標圖要用于林木的遺傳分析，首先要解決試驗地環(huán)境的異質性問題?？臻g分析已成為解決環(huán)境異質性的主要方法[10]，因此本研究采用空間分析結合因子分析法獲取各基因型在每個地點中的BLUP數(shù)據(jù)，并對原始數(shù)據(jù)和BLUP數(shù)據(jù)的GGE雙標圖分析結果進行比較，結果表明空間分析可有效擬合試驗地的空間變異，BLUP數(shù)據(jù)明顯比原始數(shù)據(jù)有更強的空間變異趨勢，方差分析的結果進一步確認BLUP數(shù)據(jù)對產量變異的解釋能力要高于原始數(shù)據(jù)，這些結果說明BLUP數(shù)據(jù)已去除空間變異對目標性狀的影響，即BLUP數(shù)據(jù)基本解決了環(huán)境異質性問題。原始數(shù)據(jù)和BLUP數(shù)據(jù)的GGE雙標圖分析結果表明，2種數(shù)據(jù)對試驗地的分組結果一致，也與因子分析法的分組結果[1]一致，意味著GGE雙標圖對于試驗地的分組結果是可靠的，但是BLUP數(shù)據(jù)中地點間的遺傳相關被弱化了；同時，BLUP數(shù)據(jù)中兩個主成分的方差解釋百分比之和(87.50%)高于原始數(shù)據(jù)(70.51%)，由此可知，BLUP數(shù)據(jù)的GGE雙標圖結果更為可靠。對試驗地的評估結果表明，原始數(shù)據(jù)最好的試驗點是地點S2，而BLUP數(shù)據(jù)是地點S1，雖然地點S1和地點S2之間的相關性比較大，但在BLUP數(shù)據(jù)中，地點S2的區(qū)分力明顯偏弱，因此地點S1應為最理想的試驗地。試驗基因型的評估結果表明，雖然原始數(shù)據(jù)和BLUP數(shù)據(jù)中最理想的基因型均為21，但2種數(shù)據(jù)種子最高產的前5個基因型的一致性為80%，而產量穩(wěn)定性強的前5個基因型的一致性僅為40%，說明2種數(shù)據(jù)的基因型評估差異較大。綜上認為，GGE雙標圖可以用于林木多點試驗分析，但需要通過空間分析解決地點環(huán)境異質性，因此采用BLUP和GGE雙標圖相結合的模型，無論是在試驗地劃分、試驗地評估還是林木基因型評估上，均比原始數(shù)據(jù)的GGE雙標圖更為可靠。

[1] 林元震.R與ASReml-R統(tǒng)計學 [M].北京：中國林業(yè)出版社，2016:526-534.

Lin Y Z.R and ASReml-R statistics [M].Beijing:China Forestry Publishing House,2016:526-534.

[2] Finlay K W,Wilkinson G N.The analysis of adaptation in a plantbreeding programme [J].Aust J Agric Res,1963,14(6):742-754.

[3] Gauch H G,Zobel R W.Identifying mega-environments and targeting genotypes [J].Crop Sci,1997,37(2):311-326.

[4] Yan W.GGEbiplot:a Windows application for graphical analysis of multi-environment trial data and other types of two-way data [J].Agron J,2001,93(5):1111-1118.

[5] 嚴威凱.雙標圖分析在農作物品種多點試驗中的應用 [J].作物學報,2010,36(11):1805-1819.

Yan W K.Optimal use of biplots in analysis of multi-location variety test data [J].ActaAgron Sin,2010,36(11):1805-1819.

[6] Ding M,Tier B,Yan W,et al.Application of GGE biplot analysis to evaluate genotype (G), environment (E) and G×E interaction onPinusradiata:a case of study [J].New Zealand J For Sci,2008,38(1):132-142.

[7] Sixto H,Salvia J,Barrio M,et al.Genetic variation and genotype-environment interactions in short rotationPopulus,plantations in southern Europe [J].New Forests,2011,42(2):163-177.

[8] Wang R H,Hu D H,Zheng H Q,et al.Genotype×environmental interaction by AMMI and GGE biplot analysis for the provenances ofMicheliachapensisin South China [J].J For Res,2016,27(3):659-664.

[9] Smith A,Cullis B,Thompson R.Analyzing variety by environment data using multiplicative mixed models and adjustments for spatial field trend [J].Biometrics,2001,57(4):1138-1147.

[10] Dutkowski G W,Costae S J,Gilmour A R,et al.Spatial analysis methods for forest genetic trials [J].Can J For Res,2002,32(12):2201-2214.

[11] Gilmour A R,Gogel B J,Cullis B R,et al.ASReml user guide release 4.0 [R].London:Vsn International Ltd Hemel,2016.

[12] Frutos E,Galindo M P,Leiva V.An interactive biplot implementation in R for modeling genotype-by-environment interaction [J].Stoch Environ Res Risk Assess,2014,28(7):1629-1641.

[13] Lin Y Z.AAfun:ASReml-R Added Functions.R package version 2.6.1[CP/OL].(2016-06-18)[2016-12-20].https://github.com/ yzhlinscau/AAfun.