The soybean (Glycine max) is an important oil seed crop that fulfills the demands of oil and proteins of millions of people around the world. Its cultivation area covered 130.94 million hectares during 2021/2022, with a net production of 355.59 million metric tons in the world (https://www.fas.usda.gov/data/world-agricultural-production). However, in Pakistan, its cultivation area is negligible, with a production of 1,000 metric tons. To meet the increasing need for food, it is essential for breeders to develop cultivars that have high yield and yield stability and are also resistant to biotic and abiotic stresses (Eltaher et al., 2021). In any crop species, grain yield is the most important factor. It is a complex quantitative trait that is influenced by multiple genes and environmental factors (Yongchun et al., 2008). Hence, it is necessary to dissect the underlying genetics of grain yield and other related traits for manipulating alleles at relevant loci to get maximum benefits (Yongchun et al., 2008). The selection of genotypes carried out in a single environment on the basis of their performance is not suitable for the development of varieties (Shrestha et al., 2012). So, the selection of the genotypes on the basis of yield stability evaluation is more important than their mean performance in multiple environmental conditions (Tariku et al., 2013; Islam et al., 2016). For crops like soybean, which grows in a wide range of ecological conditions, it is very important to select the genotypes for adaptability and stability before recommending any environmental condition. The photoperiod is the main climatic factor in soybean that determines its adaptability to different ecological conditions. Because of photoperiod sensitivity, each soybean cultivar is restricted to cultivation in a narrow range of latitudes to get maximum yield (Debebe et al., 2014). Although soybean grows in a wide range of latitudes (50°N–35°S) across the world, identification of traits that help to determine the performance of the most stable cultivars at different latitudes is very important (Li et al., 2020).

Genotype × environment interaction (GEI) has limitations in the study of important agronomic traits like yield and its components, as it complicates the understanding of genetic experimentations and restricts the selection of varieties adaptive to specific conditions (Farshadfar and Sutka, 2003). Normally, in plant breeding programs, the selection of genotypes for a specific environment is conducted by multi-environmental trials (METs) for the evaluation of genotypes based on their performance across environments (Li et al., 2020). Numerous research studies have been conducted using several statistical modeling approaches for checking the effect of GEI on yield and other agronomic traits (Grüneberg et al., 2005). These approaches mainly utilize a generalized linear model (GLM) to measure the variation caused by genotype, environment, and GEI for each variable by linear regression and joint analysis of variance (ANOVA) (Arif et al., 2021). GLMs lower the supposition of dependent variables (Olsson, 2002).

The additive main effect and multiplicative interaction (AMMI) model is mostly used in crop breeding programs for evaluating the genotypes for variety approval. First, the AMMI model uses ANOVA to divide variations into the main effect of genotype (G), main effect of environment (E), and effect caused by genotype-by-environment interaction (GEI). Second, it performs principal component analysis (PCA) by singular value decomposition for genotype and environment (Gauch Jr et al., 2008).

The most important method that visually helps to examine the relationship among genotypes, environments, and genotype-by-environment interaction and plays a significant role in the selection of the most stable and high-performing genotype for a specific environment in mega-environmental trials is the GGE biplot (Tiwari, 2019). GGE biplots play a major role in the selection of the most stable genotypes and discard those genotypes that are unstable across environments and/or have less yield (Li et al., 2020). In the past, many studies have been conducted to check the stability of soybean across environments. GGE biplots have been used to check the stability and adaptability of soybean genotypes that were cultivated in multiple environments and select the varieties that were highly stable and performed better across the environments (Mukuze et al., 2020; Carvalho et al., 2021). Other than soybean, GGE biplots were used in oat, sugarcane, rice, wheat, and maize for screening the stable genotypes in mega-environmental trials. Hence, it has been established that in agricultural research programs, the GGE biplot is the most effective method for the selection of suitable cultivars for specific environments in mega-environmental trials (Donoso-Ñanculao et al., 2016; Tena et al., 2019).

In general, genetic make-up (G), environment (E), and their interaction G × E influence the expression of any physiological and morphological trait. Due to their polygenic nature, yield and other quantitative traits are continuously controlled and affected by quantitative trait loci, genomic regions with associated genes, and environment (Said et al., 2022). As a result, genes that affect the yield and its components are highly sensitive to the environment and show QTL–environment interaction. This interaction between QTL and the environment can facilitate or constrain the responses toward artificial selection (Falconer, 1952). Therefore, breeding programs need to take these effects seriously and address them properly (El-Soda and Sarhan, 2021). Traditional QTL mapping and genome-wide association mapping are two methods that can be used for identifying the genes with an underlying natural variation that affects the genotypes.

In traditional breeding programs, the selection of genotypes is mostly carried out on the basis of phenotypic performance. Breeders mostly select the genotypes that perform well in a specific environment, costing time and resources (Baenziger, 2016). Association mapping (AM) is an alternative to conventional breeding and is considered an effective approach for dissecting the genomic location of genes or quantitative trait loci (Verdeprado et al., 2018). Based on the association between markers and traits, it performs rapid and fine mapping of the target locus (Mackay and Powell, 2007). In previous studies, association mapping conducted in soybean by using genome-wide SSR markers showed a successful marker–trait association (Ghione et al., 2021). The use of functional molecular markers, especially those derived from expressed sequence tags (ESTs), facilitates the association between phenotype and genotype by providing direct access to the population variation in genes of agronomically important traits (Mulato et al., 2010). In the past, linkage mapping was frequently used to check the effect of genotype, environment, and GEI (Ma et al., 2009). However, association mapping is rarely used for dissecting GEI (Lü et al., 2011). So, the aim of this study was 1) to check the effect of G × E interaction on the performance of soybean genotypes in multiple-environmental trials and select the genotypes that are most stable and have high adaptability across the environments, 2) identify the genotypes that give high yields in different environments, and 3) find the association between markers and traits for important agronomic traits.

The main objective of any breeding program is to develop genotypes that are resistant to biotic and abiotic stresses with high production. Environmental interaction has a significant impact on complex quantitative traits with several contributing factors such as yield. The existence of a substantial genotype main effect and G × E interactions revealed that genotypes respond differently in different environmental conditions. Studies on stability and GEI are crucial for effective breeding and adaptability in a wide range of environmental conditions (Liang et al., 2015). For the allocation of best resources in breeding or cultivar evaluation programs, the mega-environment concept is helpful (Gauch Jr and Zobel, 1997). A mega-environment is a collection of places that regularly use the same top cultivars (Yan and Rajcan, 2002). Because of the significant impact of genotype by mega-environment interaction, the evaluation of cultivars should be carried out separately for each mega-environment before the cultivar recommendation (Yan et al., 2011). Genotypes selected from the ideal test environment were mostly those with exceptional mean performance and greater adaptability (Yan et al., 2000). For identification of lines with high homeostasis in multilocation trials and coordinated variety testing programs, stability analysis models such as YSi statistics, AMMI, and GGE biplots were used. The main issue for plant breeders is to get the relevant knowledge concealed in multi-environment data and then to understand it for successful utilization. For mega-environment and cultivar evaluation, and assessment of varietal stability, GGE biplots have mostly been used (Rakshit et al., 2012; Zimmer et al., 2016). The GGE biplot was more beneficial when the mega-environment was used to evaluate a large set of genotypes, as the pattern of GEI could make the genotype evaluation more challenging (Krishnamurthy et al., 2017). In other words, environmental variation was inconsistent with the superiority of genotype, which restricted the selection of cultivars. The quality of selection can be improved by the selection of superior genotypes, with little stability variance produced through simultaneous selection for high mean and stability. In many crops, this technique has been effectively used, particularly for determining grain yield.

In this study, soybean genotypes were analyzed in five different environmental conditions, and the features that substantially correlated with stability were discussed. With the help of stability analysis models, stable genotypes with high mean yields were identified. These models suggested that the most stable genotype for TY was G69, followed by G92, G85, and G40 (Figure 6A). A genotype with a high mean yield and great stability would be an ideal genotype. A genotype located closer to the mean environment’s direction and having a projection of zero onto the perpendicular AEC coordinate is considered an “ideal” genotype. For mean yield and stability across environments, lines G69 and G92 displayed high mean rankings and were determined to be the best-performing ones (Figure 7C). Results obtained from GGE biplots suggested that G92 and G69 are most suited to E1, G69 to E2, G18 to E3, G1 to E4, and G42 to E5 (Figure 7B). These results are in line with the findings of Arif et al. (2021), who reported that genotypes DCD, BRC-457, and D-14005 were ideal genotypes for E5 as they have high mean yield and stability.

A significant variation was observed in the yield of genotypes, which may be due to the presence of diversity across environments. Similar results were reported by many researchers (Kumar et al., 2014; Bhartiya et al., 2017). Carvalho et al. (2021) also reported the epistatic effect on yield. Results obtained by GGE biplot analysis recommended that E3 was the most suitable environment for the selection of high-yielding genotypes and general adaptability (Figure 7D). E4 was the most discriminating and least representative environment for genotype evaluation and would be helpful in choosing genotypes that are specifically adapted (Mulugeta et al., 2013). The environment with the least discrimination but the most representation for the majority of attributes was E3 (Figure 7D). The large environmental difference suggested that there were genotypic variations in adaptation (Krisnawati and Adie, 2018). Regarding the stability of yield attributes, genotypes also varied greatly. Genotype G69 stood out in all evaluated environments (Table 2). However, when average yield was considered, G92, G85, and G40 performed better. Results obtained from this study are aligned with those from the work of AbdulHamid et al. (2017) and Nagarajan et al. (2017), who showed the importance of cultivating soybean genotypes with yield-contributing traits.

With the advancement of phenotyping and genotyping technologies, it has become easy to analyze the genomic regions related to quantitative traits in larger populations. Considering the relatively important interactions between the environments and genotypes, association mapping was performed for yield-contributing traits with SSR and EST-SSR markers for each environment separately. For six agronomic traits, 26 marker–trait associations were found in five environments (Table 4). No common markers among environments were discovered, which may be due to absence or very weak significant relationships between environments for TY. Additionally, many studies have demonstrated that variations in the size and structure of the population might affect the outcomes of association mapping (Liu and Cheng, 2020). Given that a smaller population offers fewer allelic types, genetic drift may be the cause of this variation. Another possible source could be the type and number of SSRs, as a preference for single-locus SSR in the present study may miscalculate the value after lowering the complexity of genotyping (Vigouroux et al., 2002). A typical population for an association study should consist of multiple unrelated and independent individuals from the same origin (Porth et al., 2013). Therefore, in order to eliminate false positive associations, we still need to confirm the association results by targeting allelic variations in coding regions via molecular biology approaches such as knockout studies (Abdurakhmonov and Abdukarimov, 2008) that offer a high-precision estimation of allelic variation.

From this study, we concluded that significant genetic variation was present between the genotypes for yield in different environments. In total, 13 environment-specific genotypes showing their maximum grain yield in each environment were identified. Genotype G69 was an ideal genotype with higher grain yield and broad adaptation to environments E1 and E2. In the case of environments, E3 was a more ideal environment as it was plotted near the concentric circle and was very informative for the selection of genotypes with high yield. Furthermore, association mapping revealed a total of 26 marker–trait associations for six agronomic traits in five environments, with the highest significance for plant height and the lowest significance for hundred-grain weights. As the G × E interaction has a significant effect on yield, it is necessary to further evaluate the ideal location for introducing suitable genotypes with stable high yield. Plant breeders must concentrate on improving features with high heritability.