# Stability Analysis to Select the Stable and High Yielding of Black Soybean (Glycine max (L) Merril) in Indonesia

Authors: Gatut Wahyu Anggoro Susanto; Haris Maulana; Pratanti Haksiwi Putri; Runik Dyah Purwaningrahayu; Acep Atma Wijaya; Binastya Anggara Sekti; Agung Karuniawan
Journal: International Journal of Agronomy
 (2023)

---

## Abstract

Selecting high-yielding and stable genotypes is the main goal of black soybean breeding. Testing new lines in multiple environments is one of the development processes in black soybean breeding. The aims of the research were (i) to identify the effects of genotype, environment, and their interactions (GEIs) on the grain yield of soybean lines in Java Island; (ii) to select stable and high-yielding soybean lines; (iii) to determine the discriminative environments; and (iv) to determine the concept of stability measurements on black soybean grain yields. Field trials of 10 newF8 promising lines and three check varieties were conducted under eight different environments during four years (2016–2019). The measurement results showed that the grain yield was influenced by genotype (8.35%), environment (59.49%), and GEIs (32.16%). Grain yield stability measurements showed that the four newly identified lines had high yields and were stable in eight environments; they were A-5A-PSJ (S2), DB-96-CTY (S5), UP 161 (S6), and UP 162 (S7). The Ngawi (2017), followed by Bogor (2019) and Banyuwangi (2016), has the strongest interactive capabilities and is suitable for use as a trial environment. Grain yield (Y) was identified as having a positive and significant correlation (p<0.05) with S(3), S(6), NP(2), NP(3), NP(4), KR, and YSI stability measurements, which indicated that they were included in the concept of dynamic stability measurement. The selected genotypes can be proposed for the release of new high-yielding black soybean varieties in Indonesia.


---

## Body

## 1. Introduction

Black soybeans contain many nutrients that are beneficial for humans. The current research of black soybean leads to the results of secondary metabolites that are beneficial to health. There are many secondary metabolites of black soybeans that are beneficial to humans, one of which is isoflavones. Isoflavones are a group of flavonoid compounds that produce natural antioxidants [1]. The benefits of isoflavones for humans include preventing cancer, heart attacks, and climate syndrome [2, 3]. Isoflavones in black soybeans were found in all plant organs, especially soybean seeds. Sumardi et al. [4] reported that the two largest types of isoflavones contained in black soybean seeds were daidzein and genistein. The important role of black soybeans for food and human health makes this commodity increasingly needed by the community, so that its demand is predicted to increase.In Indonesia, black soybeans were still needed as raw materials for production in soy sauce industries. The advantage as a raw material for soy sauce was to improve the quality of the color to that of dark chocolate. Soy sauce made from black soybeans “Merapi” has a protein content of 3.20%, higher than the soy sauce from yellow soybean “Argomulyo” (2.37%) [5]. The large role of black soybeans in functional foods makes the demand for black soybeans increase. Triyanti [6] reported that soybean production in 2020 had increased to 632,333 tons from 2019 which produced 424,190 tons. However, this amount was still not enough to meet the higher demand for black soybeans, so it was predicted that the import value of s black soybeans to meet domestic needs was still high.One of the efforts that can be made is by assembling new superior black soybeans that have high grain yields. The stages that must be passed in the assembly of promising lines are testing in a variety of environments (multilocation trial). Multilocation trials were carried out to obtain information on stable and adaptive lines in all environments as well as specific adapted lines with high yields. The more environments used for testing (mega-environment), the higher the accuracy of stable lines [7, 8]. According to Srivastava et al. [9], mega-environment analysis using the GGE biplot provides information on suitable sites for extensive adaptation and suitable sites for specific adaptation testing. Silva et al. [10] reported that mega-environment trials could lead to complex GEIs and obtained information about the relationship between trial environments. Information obtained from mega-environment trials was very helpful for plant breeders in managing GEIs and extracting them into the same agro-climatic concept [11]. According to Zdziarski et al. [12], mega-environment trial provides information on specific adapted lines based on regional targets and regional climate. In addition, based on these environment groupings, lines that show high yield stability and high yielding in suitable locations for recommended lines. Based on this description, the objectives of this study were (i) to identify the effects of genotype, environment, and GEIs on the grain yield of soybean lines in Java Island; (ii) to select stable and high-yielding soybean lines; (iii) to determine the discriminative environments; and (iv) to determine the concept of stability measurements on black soybean grain yields. The result information was useful in the development of soybean lines in the future.


## 2. Materials and Methods

### 2.1. Plant Material and Field Experiments

Planting materials used include 10 new promising F8 black soybean lines, namely, A-4A-PSJ (S1), A-5A-PSJ (S2), A-7A-PSJ (S3), A-8A-PSJ (S4), DB- 96-CTY (S5), UP 161 (S6), UP 162 (S7), UP 163 (S8), UP 164 (S9), and UP 165 (S10), as well as three check varieties as a control, namely, Detam 1 (S11), Detam 3 (S12), and GH Detam 5 (S13). The experiment was carried out in eight environments on Java Island, Indonesia, during four years (2016–2019) including Ngawi (2016), Banyuwangi (2016), Banyuwangi (2017), Ngawi (2017), Probolinggo (2018), Malang (2018), Bantul (2019), and Bogor (2019). The climate conditions of eight environments are presented in Table 1. The field experiment used a randomized completed block design which was repeated four times. Each line was planted on a plot measuring 3 × 5 m2 with a spacing of 20 cm × 30 cm.Table 1
Climate conditions in the eight environments.
LocationsNgawi (2016)Banyuwangi (2016)Banyuwangi (2017)Ngawi (2017)Probolinggo (2018)Malang (2018)Bantul (2019)Bogor (2019)T avg (°C)27.9528.2827.5326.1026.5823.2426.2621.38RH (%)67.5092.0079.0085.0077.6071.4278.7590.13DR (mm)18.1831.9633.0415.3117.504.652.4414.67ss (hour/day)6.377.917.547.114.717.897.632.19Altitude (masl)55.00168.00168.0055.0085.00335.0045.00120.00T avg, temperature average (°C); RH, humidity (%); DR, daily rainfall (mm); ss, solar radiation (hour/day); masl, meters above sea level.


### 2.2. Grain Yield Measures

A total of 13 black soybean lines were observed for their grain yield following the standard descriptor for soyabean [13]. The grain yield was harvested when the pods are brown. The grain yield data per plot was converted to ton/ha.


### 2.3. Statistical Analysis

The combined ANOVA statistical model to estimate GEIs follows the following equation:(1)Yopqr=μ+Go+Ep+GEop+Rqp+Brq+εopqr,where Yopqris the value of line o in plot r, and the value in location p of each replication q; μ is the grand mean; Go is the effect of line o; Ep is the effect of the environment; GEop is the effect of genotype by environment interactions online o and environment p; Rqp is the effect of replicate q on location p; Brq is the effect of replication q on plot r; and εopqr is the error effects from line o in plot r and repeat q of location p, respectively. To calculate the combined ANOVA, Genstat 12th is used.The grain yield stability was estimated using parametric and nonparametric measurements, including linear regression [14], mean variance component (θi) [15], GE variance component (θ(i)) [16], ecovalence value (Wi2) [17], Shukla’s stability variance σ2i [18], coefficient of variance (CVi) [19], yield Stability Indexed (YSI) [20], Nassar and Huehn [21, 22], Thennarasu [23], and Kang rank [24] (Table 2). To calculate grain yield stability based on parametric and nonparametric measurements, STABILITYSOFT is used [25].Table 2
Parametric and nonparametric stability measurements.
Parametric stability measurementsLinear regressionbi−1=∑ixij−X¯i.−X¯.j+X¯…X¯.j−X¯…/∑jX¯.j−X¯…2[14]Sdi2=1/N−2∑iX¯ij−X¯i.−X¯.j+X¯…−bi−22∑jX¯.j+X¯…2Mean variance component (θi)θi=p/2p−1q−1∑j−1qxij−X¯i.+X¯.j2+SSGE/2p−2q−1[15]GE variance component (θ(i))θi=−p/p−1p−2q−1∑j−1qxij−X¯i.−X¯.j+X¯..2+SSGE/p−2q−1[16]Eco valence value (Wi2)Wi2=∑Xij−X¯i.−X¯.j+X¯..2[17]Shukla’s stability varianceσ2iσi2=p/p−2q−1Wi2−∑Wi2/p−1p−2q−1[18]Coefficient of variance (CVi)CVi=SDl/X¯x100[19]Yield stability indexed (YSI)YSI=RASV+RGY[20]Nonparametric stability measurementsNassar and HuehnSi1=2∑jn−1∑j′=j+1nrij−rij′/Nn−1[21, 22]Si2=∑j=1nrij−r¯i.2/N−1Si3=∑j=1nrij−r¯i.2/r¯iSi6=∑j=1nrij−r¯i./r¯i.ThennarasuNP1=∑j=1nrij∗−Mdi∗/N[23]NP2=∑j=1nrij∗−Mdi∗/Mdi/NNP3=∑rij∗−ri.∗2/N/r¯i.NP6=2x∑j=1n−1∑j′=j+1nrij∗−ri.∗/r¯i./NN−1Kang rankKR = RGY + Rσ2i[24]Xij, grain yield total in of the ith line in jth environment; X¯i., average of the grain yield total from ith line at all (three) environments; X¯.j, mean of the grain yield in the jth environment; X¯…, average of the grain yield total; p and q, numbers of environments and lines; SDl, standard deviation of GEIs. rij, stability rank of the ith line in the jth environment; r¯i., average rank if ith line in all environments; and N, number of environment. rij∗, stability rank of the ith line in the jth environment (adjusted data); Mdi∗, adjusted data (median rank); Mdi, unadjusted data (median ranks of the same parameters). RGY, rank of grain yield; Rσ2i = rank of Shukla stability measurement.Identification of stable lines among the selected genotypes was conducted using AMMI following the study of [26]:(2)Yijk=µ+Gi+Ej+∑k=1mλkαikγjk+ρijr,where Yijk is the yield in location j from line i of replication k, μ is the average of grand yield, Gi is the influence of line i, Ej is the influence of the location j, λk is the value of primer component k, αik and γjk were the vector score for the line i and location j to component k, and ρijr is a mistake from line i and location j.While ASV was estimated following the study of [27](3)ASV=ssIPCA1ssIPCA2IPCA12+IPCA22,ss IPCA1 and ss IPCA2 were the sum of squares in IPCA 1 and 2, which shows the score of the main component because of the high contribution in genotype by location interactions. IPCA1 and IPCA2 were the first and second from IPCA scores for each genotype from the AMMI analysis.The model for GGE biplot following [28] with the formula is(4)Ῡmn−μm=βn+∑k=1tλoαmoγno+εmn,where Ῡmn; μm; βn; k; λo; αmo and γno; and εmn are the performance in location “n” from line “m”; overall average grain yield; the influence of location “n”; number of primer components; the singular value from primer component “o”; value of line “m” and location “n” for primer component “o”; and the error of the line “m” in location “n”, respectively. Visualization of AMMI and GGE biplot using the R program.


## 2.1. Plant Material and Field Experiments

Planting materials used include 10 new promising F8 black soybean lines, namely, A-4A-PSJ (S1), A-5A-PSJ (S2), A-7A-PSJ (S3), A-8A-PSJ (S4), DB- 96-CTY (S5), UP 161 (S6), UP 162 (S7), UP 163 (S8), UP 164 (S9), and UP 165 (S10), as well as three check varieties as a control, namely, Detam 1 (S11), Detam 3 (S12), and GH Detam 5 (S13). The experiment was carried out in eight environments on Java Island, Indonesia, during four years (2016–2019) including Ngawi (2016), Banyuwangi (2016), Banyuwangi (2017), Ngawi (2017), Probolinggo (2018), Malang (2018), Bantul (2019), and Bogor (2019). The climate conditions of eight environments are presented in Table 1. The field experiment used a randomized completed block design which was repeated four times. Each line was planted on a plot measuring 3 × 5 m2 with a spacing of 20 cm × 30 cm.Table 1
Climate conditions in the eight environments.
LocationsNgawi (2016)Banyuwangi (2016)Banyuwangi (2017)Ngawi (2017)Probolinggo (2018)Malang (2018)Bantul (2019)Bogor (2019)T avg (°C)27.9528.2827.5326.1026.5823.2426.2621.38RH (%)67.5092.0079.0085.0077.6071.4278.7590.13DR (mm)18.1831.9633.0415.3117.504.652.4414.67ss (hour/day)6.377.917.547.114.717.897.632.19Altitude (masl)55.00168.00168.0055.0085.00335.0045.00120.00T avg, temperature average (°C); RH, humidity (%); DR, daily rainfall (mm); ss, solar radiation (hour/day); masl, meters above sea level.


## 2.2. Grain Yield Measures

A total of 13 black soybean lines were observed for their grain yield following the standard descriptor for soyabean [13]. The grain yield was harvested when the pods are brown. The grain yield data per plot was converted to ton/ha.


## 2.3. Statistical Analysis

The combined ANOVA statistical model to estimate GEIs follows the following equation:(1)Yopqr=μ+Go+Ep+GEop+Rqp+Brq+εopqr,where Yopqris the value of line o in plot r, and the value in location p of each replication q; μ is the grand mean; Go is the effect of line o; Ep is the effect of the environment; GEop is the effect of genotype by environment interactions online o and environment p; Rqp is the effect of replicate q on location p; Brq is the effect of replication q on plot r; and εopqr is the error effects from line o in plot r and repeat q of location p, respectively. To calculate the combined ANOVA, Genstat 12th is used.The grain yield stability was estimated using parametric and nonparametric measurements, including linear regression [14], mean variance component (θi) [15], GE variance component (θ(i)) [16], ecovalence value (Wi2) [17], Shukla’s stability variance σ2i [18], coefficient of variance (CVi) [19], yield Stability Indexed (YSI) [20], Nassar and Huehn [21, 22], Thennarasu [23], and Kang rank [24] (Table 2). To calculate grain yield stability based on parametric and nonparametric measurements, STABILITYSOFT is used [25].Table 2
Parametric and nonparametric stability measurements.
Parametric stability measurementsLinear regressionbi−1=∑ixij−X¯i.−X¯.j+X¯…X¯.j−X¯…/∑jX¯.j−X¯…2[14]Sdi2=1/N−2∑iX¯ij−X¯i.−X¯.j+X¯…−bi−22∑jX¯.j+X¯…2Mean variance component (θi)θi=p/2p−1q−1∑j−1qxij−X¯i.+X¯.j2+SSGE/2p−2q−1[15]GE variance component (θ(i))θi=−p/p−1p−2q−1∑j−1qxij−X¯i.−X¯.j+X¯..2+SSGE/p−2q−1[16]Eco valence value (Wi2)Wi2=∑Xij−X¯i.−X¯.j+X¯..2[17]Shukla’s stability varianceσ2iσi2=p/p−2q−1Wi2−∑Wi2/p−1p−2q−1[18]Coefficient of variance (CVi)CVi=SDl/X¯x100[19]Yield stability indexed (YSI)YSI=RASV+RGY[20]Nonparametric stability measurementsNassar and HuehnSi1=2∑jn−1∑j′=j+1nrij−rij′/Nn−1[21, 22]Si2=∑j=1nrij−r¯i.2/N−1Si3=∑j=1nrij−r¯i.2/r¯iSi6=∑j=1nrij−r¯i./r¯i.ThennarasuNP1=∑j=1nrij∗−Mdi∗/N[23]NP2=∑j=1nrij∗−Mdi∗/Mdi/NNP3=∑rij∗−ri.∗2/N/r¯i.NP6=2x∑j=1n−1∑j′=j+1nrij∗−ri.∗/r¯i./NN−1Kang rankKR = RGY + Rσ2i[24]Xij, grain yield total in of the ith line in jth environment; X¯i., average of the grain yield total from ith line at all (three) environments; X¯.j, mean of the grain yield in the jth environment; X¯…, average of the grain yield total; p and q, numbers of environments and lines; SDl, standard deviation of GEIs. rij, stability rank of the ith line in the jth environment; r¯i., average rank if ith line in all environments; and N, number of environment. rij∗, stability rank of the ith line in the jth environment (adjusted data); Mdi∗, adjusted data (median rank); Mdi, unadjusted data (median ranks of the same parameters). RGY, rank of grain yield; Rσ2i = rank of Shukla stability measurement.Identification of stable lines among the selected genotypes was conducted using AMMI following the study of [26]:(2)Yijk=µ+Gi+Ej+∑k=1mλkαikγjk+ρijr,where Yijk is the yield in location j from line i of replication k, μ is the average of grand yield, Gi is the influence of line i, Ej is the influence of the location j, λk is the value of primer component k, αik and γjk were the vector score for the line i and location j to component k, and ρijr is a mistake from line i and location j.While ASV was estimated following the study of [27](3)ASV=ssIPCA1ssIPCA2IPCA12+IPCA22,ss IPCA1 and ss IPCA2 were the sum of squares in IPCA 1 and 2, which shows the score of the main component because of the high contribution in genotype by location interactions. IPCA1 and IPCA2 were the first and second from IPCA scores for each genotype from the AMMI analysis.The model for GGE biplot following [28] with the formula is(4)Ῡmn−μm=βn+∑k=1tλoαmoγno+εmn,where Ῡmn; μm; βn; k; λo; αmo and γno; and εmn are the performance in location “n” from line “m”; overall average grain yield; the influence of location “n”; number of primer components; the singular value from primer component “o”; value of line “m” and location “n” for primer component “o”; and the error of the line “m” in location “n”, respectively. Visualization of AMMI and GGE biplot using the R program.


## 3. Results

### 3.1. Combined Analysis of Variance (ANOVA) to Identify GEIs

The AMMI ANOVA revealed that genotype, environment, and their interactions (GEIs) had a significant effect (P<0.05) on soybean grain yield variations (Table 3). The main effect of genotype accounted for 8.35% of the total variation, while environment and GEIs accounted for 59.49% and 32.16%, respectively. These results indicated that grain yields were significantly affected by environmental changes, followed by the effects of GEIs and genotypes (Table 3). Highly significant environmental effects and high variance components can be attributed to large differences between test environments such as elevation and differences in the amount and distribution of annual rainfall (Table 1). The GEIs were further divided into seven IPCAs, of which the first two IPCAs showed significant contributions. IPCA1 contributed 42.66% to the number of squared GEIs. The first two IPCAs accounted for 69.30% of the total variation in GEIs, indicating the accuracy of the data. The magnitude of the GEIs effects was also seen in the box plot (Figure 1) and the GEIs heat map (Figure 2). In Figure 1 the distribution of the average grain yield in each environment showed high variation. The GEIs heatmap (Figure 2) showed the relationship between each environment and the average yield. In Figure 2, the environment was divided into three main clusters, namely, environmental groups that have high average grain yields (E4 and E5), low grain yield (E2, E6, and E7), and medium grain yield (E1, E3, and E8). This data showed that the influence of the environment was very strong on soybean grain yields.Table 3
Combined ANOVA and AMMI analysis of new promising soybean lines in eight environments.
SourcedfSSMSF% treatment% GEIsTotal415108.390.26∗∗Treatments10350.660.492.97∗∗Genotypes124.230.352.13∗8.35Environments730.134.3110.21∗∗59.49Block2410.120.422.55∗∗Interactions8416.290.191.17∗32.16IPCA1186.950.392.33∗∗42.66IPCA2164.340.271.64∗26.64IPCA3142.930.211.2717.99IPCA4121.320.110.668.10IPCA5100.460.050.282.82IPCA680.230.030.171.41IPCA760.060.010.060.37Error28847.620.17Minimum (ton/ha)0.72Mean (ton/ha)2.11Maximum (ton/ha)3.68SD0.51CV (%)19.25df, degree freedom; SS, sum of square; MS, mean of square; GEIs, genotype by environment interactions;∗p<0.05; ∗∗p<0.01; SD, standard of deviation; CV, coefficient of variation.Figure 1
GEIs box plot from eight environments on 13 soybean lines in Indonesia.E1 = Ngawi (2016), E2 = Banyuwangi (2016), E3 = Banyuwangi (2017), E4 = Ngawi (2017), E5 = Probolinggo (2018), E6 = Malang (2018), E7 = Bantul (2019), and E8 = Bogor (2019).Figure 2
GEIs heatmap from eight environments on 13 soybean lines in Indonesia.E1 = Ngawi (2016), E2 = Banyuwangi (2016), E3 = Banyuwangi (2017), E4 = Ngawi (2017), E5 = Probolinggo (2018), E6 = Malang (2018), E7 = Bantul (2019), and E8 = Bogor (2019).


### 3.2. Parametric and Nonparametric Stability Measurements on Soybean Grain Yields

The stability of soybean grain yields in multienvironmental trials based on parametric and nonparametric measurements is presented in Table4. Linear regression measurements on environmental indexed showed that grain yields of the two stability parameters bi and S2di were inconsistent in assessing the line response to various environmental conditions. Almost all lines showed a regression coefficient value (bi) which was not significantly different from one (1); on the contrary, some lines also showed a significant deviation from the regression value (S2di) greater than zero (Table 4). Based on the regression coefficients, all lines had an average response in all test environments. Thus, the lines UP161 (S6), UP 162 (S7), and UP 163 (S8) were declared the most stable according to these measurements. The other lines, Detam 1 (S11) and Detam 3 (S12), had a regression deviation (S2di) which was significantly greater than zero and a bi value that was not different from one, indicating that these two lines were better adapted to environments that produce high yields (favorable). The Wi2, σ2i, dan θ₍i₎ stability measurements showed the same results in determining the stability of soybean grain yields. Based on these three stability measurements, the A-5A-PSJ (S2) was declared the most stable line, followed by the UP 161(S6), DB-96-CTY (S5), and UP163 (S8) lines. CVi measurement indicated A-8A-PSJ (S4) as the most stable, followed by A-4A-PSJ (S1) and UP 162 (S7). θi stability measurements indicated the Detam 3 (S12) line as the most stable, followed by A-7A-PSJ (S3) and GH Detam 5 (S13). These three lines have a high average yield. ASV measurement revealed UP 161 (S6), UP 162 (S7), and A-5A-PSJ (S2) lines as the most stable. YSI measurement indicated UP 162 (S7) as the most stable line, followed by DB-96-CTY (S5) and A-5A-PSJ (S2). All three lines had high average grain yields across the eight test environments. In nonparametric stability measurements, DB-96-CTY (S5) was declared the most stable by S(3) and KR measurements. The UP 161 (S6) line was declared the most stable by NP(1) measurement. The UP 162 (S7) line was declared the most stable by S⁽1⁾, S⁽⁶⁾, NP⁽2⁾, NP⁽³⁾, and NP⁽⁴⁾ measurements. The GH Detam 5 (S13) line was declared the most stable by S(1) and S(2) measurements.Table 4
The result of stability analysis on soybean lines based on parametric and nonparametric measurements.
GenotypesYs2dibiWi2σ2iCViθ₍i₎θiASVYSIS⁽1⁾S⁽2⁾S⁽³⁾S⁽⁶⁾NP⁽1⁾NP⁽2⁾NP⁽³⁾NP⁽⁴⁾KRA-4A-PSJ2.010.030.940.230.0340.630.040.040.2715.004.6717.0329.906.733.561.190.811.0217.00A-5A-PSJ2.150.011.050.120.0141.830.050.030.259.004.3312.8613.043.442.560.300.380.557.00A-7A-PSJ2.280.081.090.620.0942.590.040.070.8314.004.2213.3612.033.034.110.340.510.4813.00A-8A-PSJ1.980.050.890.380.0539.640.040.050.4120.003.729.7514.634.253.330.780.790.7021.00DB-96-CTY2.200.031.070.200.0341.850.050.040.308.003.569.008.002.223.110.310.420.406.00UP 1612.080.031.000.180.0241.390.050.040.1710.003.9410.7512.903.502.110.380.430.5911.00UP 1622.190.031.030.220.0340.730.040.040.226.003.3910.508.692.072.560.270.320.359.00UP 1632.080.030.990.200.0341.370.050.040.3514.004.1112.7815.594.073.110.570.530.6312.00UP 1642.220.041.060.330.0441.710.040.050.5911.004.3913.7512.223.113.330.360.430.4910.00UP 1652.050.030.960.240.0340.770.040.040.4017.004.7817.8624.265.554.001.000.660.8117.00Detam 12.110.061.000.420.0642.250.040.050.6518.004.2813.1117.164.072.440.540.550.7017.00Detam 32.180.091.010.620.0941.810.040.070.7517.004.1111.8613.343.502.330.350.440.5818.00GH Detam 51.920.060.910.420.0642.370.040.050.6223.003.398.2516.505.503.781.001.100.8524.00RankYs2dibiWi2σ2iCViθ₍i₎θiASVYSIS⁽1⁾S⁽2⁾S⁽³⁾S⁽⁶⁾NP⁽1⁾NP⁽2⁾NP⁽³⁾NP⁽⁴⁾KRSRARSDA-4A-PSJ11696626848121213131013121381729.053.32A-5A-PSJ617119113331086542252894.683.39A-7A-PSJ11212121213122136810331347371538.054.24A-8A-PSJ129139919581243810810119121628.533.15DB-96-CTY3210331031152321263321753.953.00UP 16193222621214555617575894.682.77UP 1624555535921142141113623.262.10UP 163843445410676798698861226.421.98UP 1642888878695111144864441256.582.46UP 16510767747779131312121211101181739.112.59Detam 171111010111041111991193891081628.532.82Detam 3513413138131121066772566111487.793.74GH Detam 5131011111112113101311101111111312131889.903.71Y, yield; SR, sum of rank; AR, average of rank; SD, standard deviation.The grouping of soybean lines tested in this study is presented in Figure3. The grouping was the output of hierarchical clustering analysis (HCA) on the stability ranks of the tested lines based on parametric and nonparametric measurements. In Figure 3, soybean lines were divided into three main groups, namely, the unstable low yield cluster containing the A-4A-PSJ (S1), UP 165 (S10), A-8A-PSJ (S4), and GH Detam 5 (S13) lines; the unstable high yield cluster containing Detam 1 (S11), Detam 3 (S12), A-7A-PSJ (S3), and UP 164 (S9) lines; and a stable high yield cluster containing the DB-96-CTY (S5), UP 162 (S7), UP 161 (S6), UP 163 (S8), and A-5A-PSJ (S2). The third group was preferred because it had more stable and high grain yields than the overall average.Figure 3
Clustering of soybean lines based on parametric and nonparametric measurements in Indonesia.


### 3.3. The Relationship between Parametric and Nonparametric Stability Measurements with Soybean Line Grain Yields

The relationship between stability measurements was estimated using Spearman’s correlation to the stability ranks of each measurement. The results of the correlation analysis are presented in Table5. The grain yield (Y) was identified as having a positive and significant correlation (p<0.05) with the S(3), S(6), NP(2), NP(3), NP(4), KR, and YSI stability measurements. Another positive and significant correlation was shown by S(1) with S(2) (p<0.01); NP(1) with bi (p<0.01); KR with Wi2, σ2i, S2di, θ₍i₎, ASV, YSI (p<0.01), and θi (p<0.05); and CVi with ASV (p<0.05). To classify the relationship between each measurement, principal component analysis (PCA) was used on the grain yield stability ranks. The results of the PCA measurement showed that the first four PCs identified produced an eigenvalue >1, with an effect of 93.14% of the total variation (Table 6). PCA biplots regarding the relationship between stability measurements and grain yield were drawn based on the PC1 and PC2 values (Figure 4). Based on Figure 4, the grain yield (Y) was in one group with NP(2), NP(3), NP(4), S(3), S(6), YSI, and KR measurements. S(1) is grouped with S(2). NP(1) is grouped with bi. ASV is grouped with S2di, θ₍i₎, σ2i, and Wi2. The CVi measurement was weakly correlated with the last group, while the θi measurement was not correlated with all other measurements.Table 5
Correlation of parametric and nonparametric stability measurements to soybean grain yield.
YS⁽¹⁾S⁽²⁾S⁽³⁾S⁽⁶⁾NP⁽¹⁾NP⁽²⁾NP⁽³⁾NP⁽⁴⁾KRWi²σ²is²dibiCViθ₍i₎θiASVYSIY1.00S⁽¹⁾−0.061.00S⁽²⁾−0.160.951.00S⁽³⁾0.760.450.381.00S⁽⁶⁾0.860.350.260.961.00NP⁽¹⁾0.150.270.320.260.331.00NP⁽²⁾0.810.290.240.860.940.441.00NP⁽³⁾0.760.150.120.820.900.550.921.00NP⁽⁴⁾0.870.310.220.960.980.290.940.901.00KR0.63−0.09−0.100.620.720.310.670.830.691.00Wi²−0.01−0.020.070.230.250.350.240.490.210.751.00σ²i−0.01−0.020.070.230.250.350.240.490.210.751.001.00s²di−0.010.040.130.250.260.280.250.480.230.740.990.991.00bi0.05−0.13−0.18−0.150.020.710.100.27−0.030.220.250.250.151.00CVi−0.40−0.11−0.13−0.25−0.300.05−0.33−0.10−0.21−0.050.310.310.280.051.00θ₍i₎−0.01−0.020.070.230.250.350.240.490.210.751.001.000.990.250.311.00θi0.010.02−0.07−0.23−0.25−0.35−0.24−0.49−0.21−0.75−1.00−1.00−0.99−0.25−0.31−1.001.00ASV−0.170.110.140.160.150.370.150.400.130.590.910.910.900.240.520.91−0.911.00YSI0.640.060.000.740.800.310.740.880.780.940.700.700.700.140.030.70−0.700.631.00Bold numbers significantly correlate at the 5% level.Table 6
Eigenvalue, percent, and cumulative of stability measurements.
PC1234Eigenvalue9.134.792.211.56Percent (%)48.0725.2311.628.22Cumulative (%)48.0773.3084.9293.14PC, principal component.Figure 4
Principal component analysis (PCA) of parametric and nonparametric stability measurements on soybean lines in Indonesia.


### 3.4. AMMI Biplot Analysis of Soybean Grain Yield

The first two principal components (IPCA1 and IPCA2), accounted for 69.2% of the total variation in GEIs (Figure5). The length of the environment vector from the center of the biplot was proportional to the number of GEIs. An environment with a longer vector has a strong interaction, while an environment with a shorter vector has a weak interaction. The Ngawi (2017) environment has the strongest interactive power, followed by Bogor (2019) and Banyuwangi (2016). Malang (2018) environment has the weakest interactive. The lines that clustered close to the biplot axis had stable yields across the environment. Thus, the A-4A-PSJ (S1), A-5A-PSJ (S2), A-7A-PSJ (S3), DB-96-CTY (S5), and UP 162 (S7) lines had almost the same grain yield appearance (stable). Lines and environments with the same arrow pointing from the center of the biplot had positive GEIs, and vice versa. Thus, the lines UP 163 (S8), UP 164 (S9), and UP 165 (S10) have a positive GEIs with the Banyuwangi (2016), Ngawi (2018), and Bantul (2019) environments, respectively.Figure 5
AMMI biplot “PC1 vs. PC2” on 13 soybean lines in Indonesia during 2016–2019.The first IPCA (IPCA1) was plotted against the line and environmental mean (Figure6). The lines were distributed below and above the mean IPCA1 values of −40 to +60. The environments and lines on the left hand side of the vertical line have grain yields below the overall average, while those on the right side of the origin have grain yields above the overall average. The environments of Ngawi (2016), Ngawi (2017) and Probolinggo (2018) have large positive IPCA1 scores and have above-average grain yields. The Ngawi (2017) environment was classified as having the highest average yield compared to other environments, and the lines that had the same IPCA sign and were close to Ngawi (2017) were S5 and S7. The Malang (2018) environment was identified as having the lowest average grain yield and a small positive IPCA value close to zero, while Banyuwangi (2017) had a small negative IPCA value close to zero with an average grain yield less than the overall average. Lines that are close to the horizontal line have stable grain yields. Overall, of the 13 lines tested, six lines (46%) produced grain yields above the overall average grain yield, with three lines being close to the horizontal line, and seven lines (54%) produced grain yields below the overall average grain yields.Figure 6
AMMI biplot “PC1 vs. yield” on 13 soybean lines in Indonesia during 2016–2019.


### 3.5. Mega-Environments Analysis of 13 Soybean Lines in Eight Environments

The genetic correlations between the eight environments are presented in Table5 and reflected in the GGE biplot (Figure 7). In the GGE biplot analysis in 2016–2019, PC1 explained 46.28% of the total variation and PC2 explained 23.56% of the total variation, accounting for 67.84% of the total variation for grain yields (Figure 7). Biplots based on average data collected over four years showed that the eight test environments had five sectors with different winning lines (vertex). The vertex lines were A-7A-PSJ (S3), UP 164 (S9), GH Detam 5 (S13), A-8A-PSJ (S4), and Detam 3 (S12). Figure 7 represents the four mega-environments. The first mega-environment (ME-I) consisted of Ngawi (2016), Bogor (2019), Ngawi (2017), and Malang (2018) with the A-7A-PSJ (S3) line as the winner. The second mega-environment (ME-II) consisted of Banyuwangi (2016) and Probolinggo (2018), with the Detam 3 (S12) line as the winner. Bantul (2019) and Banyuwangi (2017), respectively, were the two individual mega-environments (ME-III and ME-IV) with A-8A-PSJ (S4) and UP 164 (S9) as vertex lines. Vertex line GH Detam 5 (S13), without the environment in its sector, is not a line that has a high average grain yield in all environments.Figure 7
Mega-environments on 13 soybean lines in Indonesia during 2016–2019.


## 3.1. Combined Analysis of Variance (ANOVA) to Identify GEIs

The AMMI ANOVA revealed that genotype, environment, and their interactions (GEIs) had a significant effect (P<0.05) on soybean grain yield variations (Table 3). The main effect of genotype accounted for 8.35% of the total variation, while environment and GEIs accounted for 59.49% and 32.16%, respectively. These results indicated that grain yields were significantly affected by environmental changes, followed by the effects of GEIs and genotypes (Table 3). Highly significant environmental effects and high variance components can be attributed to large differences between test environments such as elevation and differences in the amount and distribution of annual rainfall (Table 1). The GEIs were further divided into seven IPCAs, of which the first two IPCAs showed significant contributions. IPCA1 contributed 42.66% to the number of squared GEIs. The first two IPCAs accounted for 69.30% of the total variation in GEIs, indicating the accuracy of the data. The magnitude of the GEIs effects was also seen in the box plot (Figure 1) and the GEIs heat map (Figure 2). In Figure 1 the distribution of the average grain yield in each environment showed high variation. The GEIs heatmap (Figure 2) showed the relationship between each environment and the average yield. In Figure 2, the environment was divided into three main clusters, namely, environmental groups that have high average grain yields (E4 and E5), low grain yield (E2, E6, and E7), and medium grain yield (E1, E3, and E8). This data showed that the influence of the environment was very strong on soybean grain yields.Table 3
Combined ANOVA and AMMI analysis of new promising soybean lines in eight environments.
SourcedfSSMSF% treatment% GEIsTotal415108.390.26∗∗Treatments10350.660.492.97∗∗Genotypes124.230.352.13∗8.35Environments730.134.3110.21∗∗59.49Block2410.120.422.55∗∗Interactions8416.290.191.17∗32.16IPCA1186.950.392.33∗∗42.66IPCA2164.340.271.64∗26.64IPCA3142.930.211.2717.99IPCA4121.320.110.668.10IPCA5100.460.050.282.82IPCA680.230.030.171.41IPCA760.060.010.060.37Error28847.620.17Minimum (ton/ha)0.72Mean (ton/ha)2.11Maximum (ton/ha)3.68SD0.51CV (%)19.25df, degree freedom; SS, sum of square; MS, mean of square; GEIs, genotype by environment interactions;∗p<0.05; ∗∗p<0.01; SD, standard of deviation; CV, coefficient of variation.Figure 1
GEIs box plot from eight environments on 13 soybean lines in Indonesia.E1 = Ngawi (2016), E2 = Banyuwangi (2016), E3 = Banyuwangi (2017), E4 = Ngawi (2017), E5 = Probolinggo (2018), E6 = Malang (2018), E7 = Bantul (2019), and E8 = Bogor (2019).Figure 2
GEIs heatmap from eight environments on 13 soybean lines in Indonesia.E1 = Ngawi (2016), E2 = Banyuwangi (2016), E3 = Banyuwangi (2017), E4 = Ngawi (2017), E5 = Probolinggo (2018), E6 = Malang (2018), E7 = Bantul (2019), and E8 = Bogor (2019).


## 3.2. Parametric and Nonparametric Stability Measurements on Soybean Grain Yields

The stability of soybean grain yields in multienvironmental trials based on parametric and nonparametric measurements is presented in Table4. Linear regression measurements on environmental indexed showed that grain yields of the two stability parameters bi and S2di were inconsistent in assessing the line response to various environmental conditions. Almost all lines showed a regression coefficient value (bi) which was not significantly different from one (1); on the contrary, some lines also showed a significant deviation from the regression value (S2di) greater than zero (Table 4). Based on the regression coefficients, all lines had an average response in all test environments. Thus, the lines UP161 (S6), UP 162 (S7), and UP 163 (S8) were declared the most stable according to these measurements. The other lines, Detam 1 (S11) and Detam 3 (S12), had a regression deviation (S2di) which was significantly greater than zero and a bi value that was not different from one, indicating that these two lines were better adapted to environments that produce high yields (favorable). The Wi2, σ2i, dan θ₍i₎ stability measurements showed the same results in determining the stability of soybean grain yields. Based on these three stability measurements, the A-5A-PSJ (S2) was declared the most stable line, followed by the UP 161(S6), DB-96-CTY (S5), and UP163 (S8) lines. CVi measurement indicated A-8A-PSJ (S4) as the most stable, followed by A-4A-PSJ (S1) and UP 162 (S7). θi stability measurements indicated the Detam 3 (S12) line as the most stable, followed by A-7A-PSJ (S3) and GH Detam 5 (S13). These three lines have a high average yield. ASV measurement revealed UP 161 (S6), UP 162 (S7), and A-5A-PSJ (S2) lines as the most stable. YSI measurement indicated UP 162 (S7) as the most stable line, followed by DB-96-CTY (S5) and A-5A-PSJ (S2). All three lines had high average grain yields across the eight test environments. In nonparametric stability measurements, DB-96-CTY (S5) was declared the most stable by S(3) and KR measurements. The UP 161 (S6) line was declared the most stable by NP(1) measurement. The UP 162 (S7) line was declared the most stable by S⁽1⁾, S⁽⁶⁾, NP⁽2⁾, NP⁽³⁾, and NP⁽⁴⁾ measurements. The GH Detam 5 (S13) line was declared the most stable by S(1) and S(2) measurements.Table 4
The result of stability analysis on soybean lines based on parametric and nonparametric measurements.
GenotypesYs2dibiWi2σ2iCViθ₍i₎θiASVYSIS⁽1⁾S⁽2⁾S⁽³⁾S⁽⁶⁾NP⁽1⁾NP⁽2⁾NP⁽³⁾NP⁽⁴⁾KRA-4A-PSJ2.010.030.940.230.0340.630.040.040.2715.004.6717.0329.906.733.561.190.811.0217.00A-5A-PSJ2.150.011.050.120.0141.830.050.030.259.004.3312.8613.043.442.560.300.380.557.00A-7A-PSJ2.280.081.090.620.0942.590.040.070.8314.004.2213.3612.033.034.110.340.510.4813.00A-8A-PSJ1.980.050.890.380.0539.640.040.050.4120.003.729.7514.634.253.330.780.790.7021.00DB-96-CTY2.200.031.070.200.0341.850.050.040.308.003.569.008.002.223.110.310.420.406.00UP 1612.080.031.000.180.0241.390.050.040.1710.003.9410.7512.903.502.110.380.430.5911.00UP 1622.190.031.030.220.0340.730.040.040.226.003.3910.508.692.072.560.270.320.359.00UP 1632.080.030.990.200.0341.370.050.040.3514.004.1112.7815.594.073.110.570.530.6312.00UP 1642.220.041.060.330.0441.710.040.050.5911.004.3913.7512.223.113.330.360.430.4910.00UP 1652.050.030.960.240.0340.770.040.040.4017.004.7817.8624.265.554.001.000.660.8117.00Detam 12.110.061.000.420.0642.250.040.050.6518.004.2813.1117.164.072.440.540.550.7017.00Detam 32.180.091.010.620.0941.810.040.070.7517.004.1111.8613.343.502.330.350.440.5818.00GH Detam 51.920.060.910.420.0642.370.040.050.6223.003.398.2516.505.503.781.001.100.8524.00RankYs2dibiWi2σ2iCViθ₍i₎θiASVYSIS⁽1⁾S⁽2⁾S⁽³⁾S⁽⁶⁾NP⁽1⁾NP⁽2⁾NP⁽³⁾NP⁽⁴⁾KRSRARSDA-4A-PSJ11696626848121213131013121381729.053.32A-5A-PSJ617119113331086542252894.683.39A-7A-PSJ11212121213122136810331347371538.054.24A-8A-PSJ129139919581243810810119121628.533.15DB-96-CTY3210331031152321263321753.953.00UP 16193222621214555617575894.682.77UP 1624555535921142141113623.262.10UP 163843445410676798698861226.421.98UP 1642888878695111144864441256.582.46UP 16510767747779131312121211101181739.112.59Detam 171111010111041111991193891081628.532.82Detam 3513413138131121066772566111487.793.74GH Detam 5131011111112113101311101111111312131889.903.71Y, yield; SR, sum of rank; AR, average of rank; SD, standard deviation.The grouping of soybean lines tested in this study is presented in Figure3. The grouping was the output of hierarchical clustering analysis (HCA) on the stability ranks of the tested lines based on parametric and nonparametric measurements. In Figure 3, soybean lines were divided into three main groups, namely, the unstable low yield cluster containing the A-4A-PSJ (S1), UP 165 (S10), A-8A-PSJ (S4), and GH Detam 5 (S13) lines; the unstable high yield cluster containing Detam 1 (S11), Detam 3 (S12), A-7A-PSJ (S3), and UP 164 (S9) lines; and a stable high yield cluster containing the DB-96-CTY (S5), UP 162 (S7), UP 161 (S6), UP 163 (S8), and A-5A-PSJ (S2). The third group was preferred because it had more stable and high grain yields than the overall average.Figure 3
Clustering of soybean lines based on parametric and nonparametric measurements in Indonesia.


## 3.3. The Relationship between Parametric and Nonparametric Stability Measurements with Soybean Line Grain Yields

The relationship between stability measurements was estimated using Spearman’s correlation to the stability ranks of each measurement. The results of the correlation analysis are presented in Table5. The grain yield (Y) was identified as having a positive and significant correlation (p<0.05) with the S(3), S(6), NP(2), NP(3), NP(4), KR, and YSI stability measurements. Another positive and significant correlation was shown by S(1) with S(2) (p<0.01); NP(1) with bi (p<0.01); KR with Wi2, σ2i, S2di, θ₍i₎, ASV, YSI (p<0.01), and θi (p<0.05); and CVi with ASV (p<0.05). To classify the relationship between each measurement, principal component analysis (PCA) was used on the grain yield stability ranks. The results of the PCA measurement showed that the first four PCs identified produced an eigenvalue >1, with an effect of 93.14% of the total variation (Table 6). PCA biplots regarding the relationship between stability measurements and grain yield were drawn based on the PC1 and PC2 values (Figure 4). Based on Figure 4, the grain yield (Y) was in one group with NP(2), NP(3), NP(4), S(3), S(6), YSI, and KR measurements. S(1) is grouped with S(2). NP(1) is grouped with bi. ASV is grouped with S2di, θ₍i₎, σ2i, and Wi2. The CVi measurement was weakly correlated with the last group, while the θi measurement was not correlated with all other measurements.Table 5
Correlation of parametric and nonparametric stability measurements to soybean grain yield.
YS⁽¹⁾S⁽²⁾S⁽³⁾S⁽⁶⁾NP⁽¹⁾NP⁽²⁾NP⁽³⁾NP⁽⁴⁾KRWi²σ²is²dibiCViθ₍i₎θiASVYSIY1.00S⁽¹⁾−0.061.00S⁽²⁾−0.160.951.00S⁽³⁾0.760.450.381.00S⁽⁶⁾0.860.350.260.961.00NP⁽¹⁾0.150.270.320.260.331.00NP⁽²⁾0.810.290.240.860.940.441.00NP⁽³⁾0.760.150.120.820.900.550.921.00NP⁽⁴⁾0.870.310.220.960.980.290.940.901.00KR0.63−0.09−0.100.620.720.310.670.830.691.00Wi²−0.01−0.020.070.230.250.350.240.490.210.751.00σ²i−0.01−0.020.070.230.250.350.240.490.210.751.001.00s²di−0.010.040.130.250.260.280.250.480.230.740.990.991.00bi0.05−0.13−0.18−0.150.020.710.100.27−0.030.220.250.250.151.00CVi−0.40−0.11−0.13−0.25−0.300.05−0.33−0.10−0.21−0.050.310.310.280.051.00θ₍i₎−0.01−0.020.070.230.250.350.240.490.210.751.001.000.990.250.311.00θi0.010.02−0.07−0.23−0.25−0.35−0.24−0.49−0.21−0.75−1.00−1.00−0.99−0.25−0.31−1.001.00ASV−0.170.110.140.160.150.370.150.400.130.590.910.910.900.240.520.91−0.911.00YSI0.640.060.000.740.800.310.740.880.780.940.700.700.700.140.030.70−0.700.631.00Bold numbers significantly correlate at the 5% level.Table 6
Eigenvalue, percent, and cumulative of stability measurements.
PC1234Eigenvalue9.134.792.211.56Percent (%)48.0725.2311.628.22Cumulative (%)48.0773.3084.9293.14PC, principal component.Figure 4
Principal component analysis (PCA) of parametric and nonparametric stability measurements on soybean lines in Indonesia.


## 3.4. AMMI Biplot Analysis of Soybean Grain Yield

The first two principal components (IPCA1 and IPCA2), accounted for 69.2% of the total variation in GEIs (Figure5). The length of the environment vector from the center of the biplot was proportional to the number of GEIs. An environment with a longer vector has a strong interaction, while an environment with a shorter vector has a weak interaction. The Ngawi (2017) environment has the strongest interactive power, followed by Bogor (2019) and Banyuwangi (2016). Malang (2018) environment has the weakest interactive. The lines that clustered close to the biplot axis had stable yields across the environment. Thus, the A-4A-PSJ (S1), A-5A-PSJ (S2), A-7A-PSJ (S3), DB-96-CTY (S5), and UP 162 (S7) lines had almost the same grain yield appearance (stable). Lines and environments with the same arrow pointing from the center of the biplot had positive GEIs, and vice versa. Thus, the lines UP 163 (S8), UP 164 (S9), and UP 165 (S10) have a positive GEIs with the Banyuwangi (2016), Ngawi (2018), and Bantul (2019) environments, respectively.Figure 5
AMMI biplot “PC1 vs. PC2” on 13 soybean lines in Indonesia during 2016–2019.The first IPCA (IPCA1) was plotted against the line and environmental mean (Figure6). The lines were distributed below and above the mean IPCA1 values of −40 to +60. The environments and lines on the left hand side of the vertical line have grain yields below the overall average, while those on the right side of the origin have grain yields above the overall average. The environments of Ngawi (2016), Ngawi (2017) and Probolinggo (2018) have large positive IPCA1 scores and have above-average grain yields. The Ngawi (2017) environment was classified as having the highest average yield compared to other environments, and the lines that had the same IPCA sign and were close to Ngawi (2017) were S5 and S7. The Malang (2018) environment was identified as having the lowest average grain yield and a small positive IPCA value close to zero, while Banyuwangi (2017) had a small negative IPCA value close to zero with an average grain yield less than the overall average. Lines that are close to the horizontal line have stable grain yields. Overall, of the 13 lines tested, six lines (46%) produced grain yields above the overall average grain yield, with three lines being close to the horizontal line, and seven lines (54%) produced grain yields below the overall average grain yields.Figure 6
AMMI biplot “PC1 vs. yield” on 13 soybean lines in Indonesia during 2016–2019.


## 3.5. Mega-Environments Analysis of 13 Soybean Lines in Eight Environments

The genetic correlations between the eight environments are presented in Table5 and reflected in the GGE biplot (Figure 7). In the GGE biplot analysis in 2016–2019, PC1 explained 46.28% of the total variation and PC2 explained 23.56% of the total variation, accounting for 67.84% of the total variation for grain yields (Figure 7). Biplots based on average data collected over four years showed that the eight test environments had five sectors with different winning lines (vertex). The vertex lines were A-7A-PSJ (S3), UP 164 (S9), GH Detam 5 (S13), A-8A-PSJ (S4), and Detam 3 (S12). Figure 7 represents the four mega-environments. The first mega-environment (ME-I) consisted of Ngawi (2016), Bogor (2019), Ngawi (2017), and Malang (2018) with the A-7A-PSJ (S3) line as the winner. The second mega-environment (ME-II) consisted of Banyuwangi (2016) and Probolinggo (2018), with the Detam 3 (S12) line as the winner. Bantul (2019) and Banyuwangi (2017), respectively, were the two individual mega-environments (ME-III and ME-IV) with A-8A-PSJ (S4) and UP 164 (S9) as vertex lines. Vertex line GH Detam 5 (S13), without the environment in its sector, is not a line that has a high average grain yield in all environments.Figure 7
Mega-environments on 13 soybean lines in Indonesia during 2016–2019.


## 4. Discussion

The development of stable and high-yielding black soybean lines is one of the goals of the current breeding program. Effects of GEIs are frequent for many quantitative traits including grain yield [29], reproductive fitness, harvest time, and biotic and abiotic resistance [30]. The emergence of GEIs effects in multienvironment experiments makes it difficult to select the superior line [29, 31, 32]. In this study, the combined ANOVA showed that the effects of environment (E) and GEIs explained 59.49, and 32.16% of the total variability of the grain yield, respectively (Table 3). The box plot (Figure 1) and heatmap GEIs (Figure 2) also showed that the tested environments and lines had high grain yield variations. Since a combined ANOVA cannot explain the effect of GEIs, other statistical models such as AMMI are more useful. Therefore, parametric and nonparametric measurements were used to estimate the adaptability of the lines to a wide range of environments. The magnitude of the environmental and GEIs effects showed that the test environment had a significant effect on grain yields. Differences in rainfall, humidity, and light in each environment during the experiment had a significant effect on soybean grain yields [33]. In this study, the Ngawi (2017) had the highest average grain yield (2.54 t/ha) compared to other environments. The average daily rainfall during the experiment in Ngawi (2017) was 15.31 mm (Table 1). According to Mandić et al. [34], rainfall was related to soybean grain yields. In other studies, rainfall also affects sweet potato yields [35]. In addition, the humidity during the experiment also showed suitable conditions for the growth and development of soybeans.Stability and adaptability analyzes were used when the results of the combined ANOVA indicated the presence of GEIs. In this case, we used various stability measurements, namely parametric, nonparametric, AMMI, and GGE biplot. The use of various measurements aims to increase accuracy in selecting stable and high-yielding lines tested [29, 36, 37]. The parametric and nonparametric measurements presented in Table 4 show that each measurement selects a different stable line. This is because each measurement has different assumptions in terms of calculating stability and adaptability. However, there were several measurements showing uniform output, namely Wi2, σ2i, and θ(i). In other studies, it was also reported that there were parametric measurements that had similar results on barley in Iran [29]. According to Karuniawan et al. [31], if there are stability measurements that have similar outputs, one of them can be chosen to estimate the stability of the yield. Ruswandi et al. [8] showed different results; all parametric and nonparametric measurements showed differences in selecting stable maize. This difference in output indicates that another approach was needed to classify the stable lines based on parametric and nonparametric measurements. The approach commonly used in cases like this was hierarchical cluster analysis (HCA) [29, 31]. Based on the HCA results, the tested lines were divided into three main groups, namely, unstable low yield, unstable high yield, and stable high yield groups (Figure 3). Lines belonging to the unstable low-yield group include A-4A-PSJ (S1), UP 165 (S10), A-8A-PSJ (S4), and GH Detam 5 (S13). This first group was less favored because of low grain yields, which will result in low economic income. The second group, namely unstable high yield, consists of the Detam 1 (S11), Detam 3 (S12), A-7A-PSJ (S3), and UP 164 (S9) lines. The second group can be used for specific adapt. Meanwhile, the third group, namely stable high yields, consists of DB-96-CTY (S5), UP 162 (S7), UP 161 (S6), UP 163 (S8), and A-5A-PSJ (S2). The third group was the expected ideal group. According to Goksoy et al. [38], stable and high-yielding lines are one of the targets for soybean breeding. In addition, lines with high and stable yields are expected to increase farmers’ income [37]. Thus, this group was strategically developed for a further sustainable soybean breeding program.To classify parametric and nonparametric stability measurements, the principal component analysis (PCA) was used. The results of PCA (Figure4) and correlation analysis (Table 5) showed that the grain yield (Y) in the same group as NP(2), NP(3), NP(4), S(3), S(6), KR, and YSI. This group belongs to the concept of dynamic stability [29, 36]. Another group that correlates with each other were S(1) with S(2); bi with NP(1); S2di with ASV, θ₍i₎, σ2i, and Wi2. While CVi and θi were not grouped with other measurements. These groups were included in the concept of static stability [29]. Thus, dynamic stability groups can be used to select lines in favorable environments [36]. While the static stability group can be used to select lines in an unfavorable environment.AMMI biplots showed that the first two principal components (IPCA1 and IPCA2) explained 69.2% of the total variation of GEIs (Figure5). In Figure 5, the lines A-4A-PSJ (S1), A-5A-PSJ (S2), A-7A-PSJ (S3), DB-96-CTY (S5), and UP 162 (S7) were closest to the biplot axis point. According to several researchers, the lines that were close to the biplot axis point were the most stable [39, 40]. The Ngawi (2017) environment has the longest vector length, followed by Bogor (2019) and Banyuwangi (2016), while Malang (2018) has the shortest vector. According to Gauch [41], the environmental vector length from the center of the biplot was proportional to the number of GEIs. This shows that the Ngawi (2017), Bogor (2019), and Banyuwangi (2016) have a strong interaction, while the Malang (2018) environment has a weak interaction.In the AMMI biplot PC1 vs. Yield (Figure6), the first IPCA biplot (IPCA1) was plotted against the line and environmental mean. The environments and lines tested on the left side of the vertical line had grain yields below the overall average grain yield, while those on the right had above the overall average grain yield [37, 40]. Ngawi (2017) environment was classified as having the highest average yield compared to other environments, and the lines that were close to this environment were DB-96-CTY (S5) and UP 162 (S7). The Malang (2018) environment was identified as had the lowest average yield and a small positive IPCA value close to zero while Banyuwangi (2017) had a small negative IPCA value that was close to zero with an average grain yield that was less than the overall average grain yield. These two environments were not suitable for the testing and development of soybeans. Lines that were close to the horizontal line had stable grain yields [41]. Overall, of the 13 lines tested, six lines (46%) produced grain yields above the overall average grain yield, namely A-5A-PSJ (S2), A-7A-PSJ (S3), DB-96-CTY (S5), UP 162 (S7), UP 164 (S9), and Detam 3 (S12), with three lines that were close to the horizontal line (stable), namely A-5A-PSJ (S2), DB-96-CTY (S5), and UP 162 (S7). While the other seven lines (54%) had grain yields below the overall average grain yield. According to several researchers, the lines that were on the right side and close to the horizontal line were potential lines and could provide economic income [37, 40].The mega-environment trials use the GGE biplot “which won where” (Figure7). The results of the GGE biplot analysis provided an opportunity to detect differences between environments with different characteristics related to the performance of the tested lines in these environments [29, 42]. The line at the top of the sector (vertex) has the highest yield for the environment in that sector. One of the important features of this biplot was the clustering of environments, which suggests the possibility of different mega-environments [7]. Our results showed that during 2016–2019, the environments fell into different clusters with varied environmental clustering patterns. The first two PCs accounted for 67.84% of the total variability attributable to the effects of genotype (G), environment (E), and their interactions (GEIs) over 4 years (Figure 7). Based on four years of average data, we found four mega-environments with different vertex lines. This suggests a specific adaptation of the line to the mega-environment and positive exploitation of GEIs [29, 42]. Our findings revealed that some of the new lines were better adapted to different environments on the Java island (stable) than other lines and check varieties. Based on our results, both numerical (parametric and nonparametric) and graphical (AMMI and GGE biplot) measurements revealed the same pattern for the identification of stable soybean lines. For example, parametric and nonparametric stability measurements identified the DB-96-CTY (S5), UP 162 (S7), UP 161 (S6), UP 163 (S8), and A-5A-PSJ (S2) lines as the most stable and high yielding than others. The AMMI biplot identified A-5A-PSJ (S2), DB-96-CTY (S5), UP 161 (S6), and UP 162 (S7) lines, as the most stable and high yielding. The GGE biplot identified A-5A-PSJ (S2), DB-96-CTY (S5), UP 161 (S6), UP 162 (S7), and Detam 1(S11) as the most ideal lines. Similarly, Karuniawan et al., [31] used numerical measurements, AMMI, and GGE biplot to select stable sweet potatoes in West Java, Indonesia. Several other studies also reported the relative contribution of these various measurements to identify the ideal line, including on barley [43], on maize [44, 45], on safflower [46], on stevia [47], and soybean [48, 49]. From our test results, there were four new soybean lines that were identified as stable and high-yielding by various stability measurements, namely A-5A-PSJ (S2), DB-96-CTY (S5), UP 161 (S6), and UP 162 (S7). They can be further proposed as candidates for superior black soybean lines in Indonesia.


## 5. Conclusions

The effects of genotype (G), environment (E), and their interactions (GEIs) have an influence on black soybean grain yield variations, with contributions of 8.35%, 59.49%, and 32.16% of the total variation, respectively. Four new black soybean lines were identified as having high and stable grain yields in diverse environments on the Java island, namely A-5A-PSJ (S2), DB-96-CTY (S5), UP 161 (S6), and UP 162 (S7). The Ngawi (2017), followed by Bogor (2019) and Banyuwangi (2016) could be considered as environments with the strongest interaction. The stability measurements NP(2), NP(3), NP(4), S(3), S(6), KR, and YSI were included in the concept of dynamic stability that can select the best lines in a favorable environment.



---

*Source: 7255444-2023-01-20.xml*