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 12, reproductive fitness, harvest time, and biotic and abiotic resistance 13. The emergence of GEIs effects in multi-environment experiments makes it difficult to select the superior line 12,14. 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 1). The box plot (Fig. 1) and heatmap GEIs (Fig. 2) also showed that the tested environments and lines had high grain yield variations. Since combined ANOVA cannot explain the effect of GEIs, other statistical models such as AMMI are more useful. Therefore, parametric and non-parametric 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 15. 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 2). According to Mandic et al. 16, rainfall was related to soybean grain yields. In other studies, rainfall also affects in sweet potato yields 17. In addition, the humidity during the experiment also showed suitable conditions for growth and development of soybean.
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, non-parametric, AMMI, and GGE biplot. The use of various measurements aims to increase accuracy in selecting stable and high yielding lines tested 12,18,19. The parametric and non-parametric measurements presented in Table 4, showed 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, σ²ᵢ, and θ₍ᵢ₎. In other studies, it was also reported that there were parametric measurements that had similar results, including Vaezi et al. 12 on barley in iran. According to Karuniawan et al.17, if there were stability measurements that have similar outputs, one of them can be chosen to estimate the stability of the yield. Ruswandi et al. 20 showed different things, all parametric and non-parametric 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 non-parametric measurements. The approach commonly used in cases like this was Hierarchical Cluster Analysis (HCA) 12,17. 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 (Fig. 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 21, a stable and high yielding lines is one of the targets for soybean breeding. In addition, lines with high and stable yields are expected to increase farmers' income 19. Thus, this group was strategically developed for a further sustainable soybean breeding program.
To classify parametric and non-parametric stability measurements, the principal component analysis (PCA) was used. The results of PCA (Fig. 4) and correlation analysis (Table 4) showed that the grain yield (Y) in the same group with NP(2), NP(3), NP(4), S(3), S(6), KR, and YSI. This group belongs to the concept of dynamic stability 12,18. Another group that correlates with each other were S(1) with S(2); bi with NP(1); S2di with ASV, θ₍ᵢ₎, σ²ᵢ, and Wi2; While CVi and θᵢ were not grouped with other measurements. These groups were included in the concept of static stability 12. Thus, dynamic stability groups can be used to select lines in favorable environments 18. While the static stability group can be used to select lines in a unfavorable environment.
AMMI biplots showed that the first two principal components (IPCA1 and IPCA2), explained 69.2% of the total variation of GEIs (Fig. 5). In Fig. 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 22,23. 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 24, 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 (Fig. 6), 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 19,23. Ngawi (2017) environment was classified as had the highest average yield compared to other environments, and the lines that 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) has a small negative IPCA value that was close to zero with an average grain yield of less than the overall average grain yield. These two environments were not suitable for testing and development of soybeans. Lines that were close to the horizontal line have stable grain yields 24. 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 can provide economic income 19,23.
The mega-environment trials uses the GGE biplot “which won where” (Fig. 7). 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 12,25. 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 (Fig. 7). Based on 4 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 12,25. 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 non-parametric) and graphical (AMMI and GGE biplot) measurements resulted the same pattern for the identification of stable soybean lines. For example, parametric and non-parametric 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. 17 used a numerical measurements, AMMI, and GGE biplot to select stable sweet potato in West Java, Indonesia. Several other studies also reported the relative contribution of these various measurements to identify the ideal line, including Vaezi et al. 26 on barley, Oyekunle et al. 27 on maize, Jamshidmoghaddam & Pourdad 28 on safflower. 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 were can be further proposed as candidates for superior black soybean lines in Indonesia.