From 2014 to 2018, 51 fields were located in 16 different locations in the Hebei Province of China ( Table 1 , Figure 1 ). According to the annual trial results, new genotypes were added to the evaluation plan every year. A total of 14, 13, 14, 16, and 14 winter wheat genotypes were evaluated in 2014, 2015, 2016, 2017, and 2018, respectively. Information on the winter wheat genotypes evaluated each year is shown in Supplementary Tables S1–S5 . A randomized complete block design (RCBD) was adopted with three replications in a plot of size 13.33 m². The agronomic measures during the experiment were based on local field management, and the grain yield was measured in kg ha⁻¹, with a correction of 13%.

Firstly, 19 meteorological factors from 16 locations from 1989 to 2019 were collected using the R package EnvRtype (Costa-Neto et al., 2021). To better classify mega-environments (MEs), we filtered the data for 19 meteorological factors each year and only screened for meteorological data covering the winter wheat growing season (October to June). In addition, data for six local soil chemical factors were obtained using the SoilType package (Fritsche-Neto, 2023). A dataset containing 25 environmental covariates (ECs) is listed in Table 2 . The 25 ECs were used to construct an envirotype-covariable matrix (W) that further computed the environmental affinities W_matrix using functions in the EnvRtype package. With a nine-month period (October to June) to represent the temporal variation of the crop growth period, 6,750 (30 years × 25 ECs × 9 intervals) variables were obtained, and the envirotype covariable matrix (25 environmental rows × 6,750 climatic variables’ columns) was used to calculate the enviromic kernel (KE) using the following formula:

where K_E is the environmental similarity kernel based on environment “omics” and W is the envirotype matrix. To identify MEs, hierarchical clustering (average linkage method) was applied to the K_E .

Finally, to show the correlations between the 25 environmental variables, principal component analysis (PCA) was performed by creating a bidirectional table containing environmental variables using the ‘fviz_pca_biplot()’ function from the R package factoextra (Kassambara and Mundt, 2017).

All statistical analyses of the raw data involved in this study were performed using R software 4.3.1 (R Core Team, 2022) with the packages mentioned previously.

Based on 30 years (1989–2019) of climate and soil information, four mega-environments (MEs) were identified using similarity analysis of 25 environmental covariates (ECs), including 19 meteorological factors and six soil physicochemical factors. ME1 comprised locations—Shenzhou, Wuyi, Malan, Dacaozhuang, Xingtai, and Xinhe. ME2 includes locations such as Linxi, Yongnian, Feixiang, and Handan. Nanpi and Fucheng were included in ME3, whereas Gaocheng, Zhengding, Luquan, and Xinle were included in ME4 ( Figure 2 ). These MEs exhibited geographic proximity, with the first two principal components accounting for 68.8% of location variability, indicating significant differences in environmental variables across environments ( Supplementary Figure S1 ).

Principal component analysis based on 30-year climate data and soil composition information revealed that ME1 had higher TRANGE, PRECTOT, GDD, N, SAND, PHAQ, and SILT values. ME2 had higher RTA, WS2M, FRUE, SPV, TMAX, ASKLW, and PETP. ME3 had lower climate and soil variables, namely, VPD, ORGC, NITKJD, and CLAY, and ME4 was associated with higher ETP, ASKSW, and n ( Figure 3 ). RTA, T2M, TMIN, FRUE, and GDD were the climatic variables that contributed the most to the environmental scores ( Supplementary Figure S2 ). During the 2014–2018 METs, grain yield (GY) was found to be positively correlated with environmental variables such as ETP, VPD, TMAX, SPV, GDD, T2M, FRUE, and ORGC, and negatively correlated with environmental factors such as RH2M, PRECTOT, PETP, and T2MDEW ( Figure 4 ). This variation was predominantly attributed to TMIN, FRUE, T2M, SPV, and GDD during 2014–2018 ( Supplementary Figure S3 ).

The combined analysis of variance for each year showed that genotype (G), environment (E), and GE were all highly significant (P ≤0.001) for grain yield at all locations ( Table 3 ). Among the sources of variation, the environmental effect recorded the highest sum of squares, indicating that the highest degree of variation was 47.72%, 66.82%, 70.59%, 68.28%, and 69.26% in all five years, respectively. The proportions of GE interactions accounted for 15.54%, 9.51%, 7.76%, 9.70%, and 10.00%, respectively. From 2014 to 2018, genotypic main effects accounted for 12.64%, 5.46%, 3.24%, 3.96%, and 4.99% of the total variation. During 2014–2018, the proportion of total variation recorded by GE interactions was greater than that of the genotype, indicating its importance for variety performance across various locations.

The AMMI analysis results for all five years are presented in Table 4 . The GE interaction matrix was partitioned to form a multiplicative component. The first three IPCAs (interaction principal component axes) obtained by singular value decomposition of GE interactions were statistically significant (p <0.01) using the F-test. IPCA1 accounted for 27.7%, 32.0%, 33.9%, 44.2%, and 55.5% of total GE interactions between 2014 and 2018, respectively.

The evaluation identified the optimal models for the period 2014–2018. Based on our multi-year datasets showing various genotype-environment interaction (GEI) patterns, our analysis indicated that the Best Linear Unbiased Prediction (BLUP) model provided the most accurate predictions. Furthermore, we observed that AMMI9, AMMI6, AMMI0, AMMI4, and AMMI2 demonstrated the highest accuracy among the AMMI models during the period 2014–2018 ( Figure 5 ).

A heat map ( Figure 6 ) was generated to compare various stability indices derived from the AMMI, BLUP, and WASSY methods based on their rankings ( Supplementary Tables S6 , S7 ). The rank correlation analysis of these indices demonstrated that the BLUP and AMMI indices formed two separate clusters with a strong correlation observed within each cluster, while correlations between clusters were weaker.

During the period from 2014 to 2018, the correlation between AMMI indices and grain yield varied between −0.53 to 0.42, −0.53 to −0.13, −0.33 to 0.15, −0.70 to −0.13, and −0.60 to −0.15, respectively, for each respective year, whereas the correlation of the BLUP index with grain yield was 1.00^*** to 1.00^***, 1.00^*** to 1.00^***, 0.99^*** to 1.00^***, 1.00^*** to 1.00^***, and 1.00^*** to 1.00^*** across the same time frame. Overall, this suggests that the BLUP indices exhibited a high correlation with grain yield, whereas the AMMI indices displayed a low or negative correlation.

The WAASY index was significantly and positively correlated with grain yield and BLUP indices. This suggests that the WAASY index considers both grain yield and GE interactions when identifying the best genotypes. The correlation between WAASY index and AMMI index in 2014–2018 was −0.48 to −0.86^***, −0.60^* to −0.85^***, −0.50 to 0.84^***, −0.68^** to −0.94^*** and −0.54^* to −0.91^***, respectively, while the correlation with BLUP index was 0.81^*** to 0.82^***, 0.84^*** to 0.86^***, 0.73^*** to 0.75^***, 0.78^*** to 0.79^*** and 0.83^*** to 0.84^***, respectively.

According to the GGE biplot analysis, the first two PCs (principal components) explained 66.36%, 62.60%, 64.43%, 65.75%, and 80.67% of the total genotype and GE variation from 2014 to 2018, respectively. The winning genotype at each location was identified using the which-won-where view of the GGE biplot. During 2014, genotypes JM196 exhibited high yield in locations Shenzhou, Zhengding, and Yongnian to become the champion genotype, while genotypes WM4176 and HN6119 were identified as superior genotypes in the remaining locations. During 2015, genotypes ZX4899 and H9966 became the general winners across six locations including Xinle, Feixiang, Shenzhou, Handan, Gaocheng, and Wuyi, and genotype LM22 became the universal winner across the remaining five locations. During 2016, genotypes BM7, KN8162, and KM3 were the champion genotypes on locations Fucheng and Xinhe, locations Nanpi and Gaocheng, locations Zhengding, Yongnian, Xingtai, and Shenzhou, respectively. During 2017, the champion genotype across locations Gaocheng, Linxi, Handan, Xingtai, Nanpi, and Luquan was HH14-4019, while the genotypes HM15-1 and HH1603 performed well in the locations Yongnian, Malan, Shenzhou, and Xinle as the universal winner. During 2018, the genotype S14-6111 had broad adaptability in the locations Malan, Xingtai, Dacaozhuang, Gaocheng, Wuyi, and Luquan as the universal genotype, while JM5172 performed best in the location Yongnian ( Figure 7 ).

Based on the discriminative and representative view of GGE biplot, during 2014, Yongnian, Malan, Gaocheng, Handan, Feixiang, and Luquan were the most discriminative locations. Moreover, the Handan and Feixiang locations were more representative. Thus, Handan and Feixiang were locations with both discrimination and representativeness. In 2015, Yongnian, Xinle, and Feixiang were the most discriminative locations, and Nanpi and Xingtai were more representative. In 2016, the highest discriminative locations were Fucheng, Yongnian, Xingtai, and Zhengding, whereas the representative location was Gaocheng. During 2017, locations Shenzhou and Xinle were higher discriminative, while Yongnian, Malan, Shenzhou, and Xingtai showed good representative ability among the test locations. Thus, Shenzhou was the most representative and best discriminating location for evaluating winter wheat. Likewise, during 2018, Shenzhou, Yongnian, Malan, and Xingtai were the higher discriminative locations, while Shenzhou, Luquan, Nanpi, and Handan showed good representative ability. Shenzhou was the location with higher discriminative and good representative ( Figure 8 ).

The results of this study confirm that GEI plays a significant role in determining winter wheat yield. The findings revealed that environmental factors, especially temperature, relative humidity, and vapor pressure deficit (VPD), significantly affected genotypic performance. These interactions are particularly pronounced in multi-year and multi-location trials, where the variability of the GEI exceeds that of genotype effects alone. This result aligns with those of previous studies (Özdoğan, 2011; Semenov et al., 2014; Bajwa et al., 2020), highlighting that environmental factors, especially those linked to climate change, cannot be overlooked in wheat production. The North China Plain (NCP), a region with dynamic agricultural conditions, serves as a prime example of how the GEI shapes genotype performance (Xiao et al., 2020; Zhang et al., 2020). Understanding the role of environmental changes in wheat growth is crucial for the selection of genotypes with superior adaptability. In this context, selecting genotypes with stable and high yield potential across various environments is paramount, as this can directly improve productivity in regions affected by climate variability (Nataraj et al., 2024).

The study revealed that warmer conditions, especially in locations such as Dacaozhuang, Xinhe, and Xingtai, led to elevated temperatures (T2M, TMAX, and TMIN), reduced relative humidity, and consequently, higher VPD. These conditions result in increased evapotranspiration (ETP), especially in Xinhe, suggesting that water loss through evapotranspiration is more pronounced in warmer environments. This positive correlation between environmental temperature and ETP supports the hypothesis that environmental warming increases the challenges associated with water loss in wheat production (Casagrande et al., 2024). These findings underscore the need for climate-resilient wheat genotypes that can withstand environmental stresses. Genotypes that exhibit stable yields under varying temperature and relative humidity are essential for maintaining productivity in the face of climate change. Focusing on these environmental factors when selecting drought-resistant or heat-tolerant traits could improve the resilience of future wheat crops (Wang et al., 2023).

Selecting appropriate statistical models is essential for improving the yield predictions in METs (Koundinya et al., 2021; Yue et al., 2022b). In this study, the Best Linear Unbiased Prediction (BLUP) model outperformed the Additive Main Effects and Multiplicative Interaction (AMMI) model for predicting yield. These findings support the use of BLUP as a preferable tool for large datasets, as it accounts for both genotype and environmental effects, and provides more reliable predictions for breeders (Piepho, 1994). While AMMI remains useful for modeling genotype–environment interactions and imputing missing data, our results suggest that BLUP is particularly valuable for predicting yield potential in trials with diverse environmental conditions. These insights can be directly applied to breeding programs that seek to optimize yield predictions and to select the best-performing genotypes for specific locations or climatic conditions (Piepho et al., 2008).

Wheat is a globally distributed and highly adaptable crop that not only has high nutritional value but also possesses unique gluten properties and excellent processing characteristics, making it suitable for the production of a wide variety of foods (Hyles et al., 2020; Yadav et al., 2022). Although China is the largest wheat producer in the world, with an annual production exceeding 137 million tons, accounting for approximately 17% of the global total, the mean yield in China is 5, 812 kg/ha, which is only approximately 80% of the yield levels seen in advanced agricultural countries like France. The analysis of numerous genotypes across diverse environments presents challenges in identifying consistent responses across environments, particularly without graphical representations of the data. The GGE biplot enabled visualization of genotype performance across multiple environments. By identifying mega-environments, the GGE biplot helps breeders focus on genotypes that are best suited to specific environmental conditions (Yan, 2001; Yan and Holland, 2010). Locations like Handan, Feixiang, Nanpi, and Shenzhou were consistently identified as representative environments over multiple years. These locations offer valuable insights for breeders looking to select genotypes that can thrive under specific climatic conditions, either for general or targeted adaptation. For example, genotypes such as JM196, WM4176, and ZX4899 performed consistently well in their respective mega-environments. These genotypes show promise for use in breeding programs that target specific regions. Moreover, environments with high discriminative power, such as Yongnian and Feixiang, are ideal for identifying stable genotypes with broad adaptability, contributing to the selection of high-yielding varieties suitable for various regions. Recognizing consistently similar environments aids in optimizing location selection for multi-environment trials (METs), thus reducing MET costs. Similar analyses have been applied to rice (Singh et al., 2023), maize (Yue et al., 2021), winter (Ferrante et al., 2021), and sugarcane (Mehareb et al., 2022), utilizing the GGE biplot to pinpoint the most discriminative test environments. This study also underscores the unpredictability of year-to-year variations at the same location, emphasizing the need for a stability analysis to better align specific genotypes with specific environments.