Application of BLUP-GGE Biplot in Mega-Environment Analysis and Test Location Evaluation of Wheat Regional Trials in the Huanghuai Winter Wheat Region in China

1. Introduction

Wheat (Triticum aestivum L.) is the third-largest food crop in China, with a planting area second only to rice and maize, and it serves as a cornerstone of national food security [1]. Its production stability is critical for safeguarding China’s agricultural economy and food self-sufficiency [2]. The Huanghuai winter wheat region (HWWR) is the most strategically important wheat-producing area in China, contributing over 60% of the nation’s total wheat cultivation area and yield, and its ecological and agronomic characteristics directly determine national wheat production capacity [1,3]. Stretching approximately 1200 km east-west from the coastal plains of Sheyang (Jiangsu Province) to the Guanzhong Plain of Baoji (Shaanxi Province), the HWWR spans Shaanxi, Henan, Anhui, and Jiangsu provinces, with an average inter-test location distance of 385 km, forming a naturally extensive and heterogeneous experimental network.

The HWWR exhibits profound ecological heterogeneity, driven by complex topographical variations and a climatic gradient transitioning from subtropical in the south to warm-temperate in the north. This gradient results in significant spatial variations in annual precipitation, temperature regimes, and soil typologies [4], which inevitably induce strong genotype-by-environment (G×E) interactions. Such interactions cause cultivar performance—including yield formation, growth duration, and stress resistance—to vary substantially across locations, posing major challenges for the precise evaluation of new varieties and their rational spatial deployment. Accurate ME delineation based on robust G×E interaction analysis is therefore a practical necessity for crop breeding and regional trial research, as it enables the maximization of cultivar yield potential by matching superior genotypes to their optimal growing conditions [5], thereby enhancing overall production efficiency and sustainability.

Multi-environment variety trials (METs) are the fundamental basis for assessing the yield potential, stability, adaptability, and stress resistance of new wheat cultivars. The scientific rigor and statistical reliability of METs directly determine the quality of variety registration decisions and subsequent promotion strategies [6,7]. A critical component of successful MET analysis is the systematic evaluation of test locations and scientifically defensible delineation of MEs. Historically, ecological subregions in the HWWR were defined based on administrative boundaries or broad climatic classifications [1], which reflect geographical proximity rather than the specific G×E interaction patterns that govern cultivar performance. This limitation frequently results in suboptimal cultivar recommendations, as micro-environmental variations and soil-climate interactions can alter cultivar performance within traditionally defined zones [2]. Advanced statistical models that analyze MET data to uncover biological patterns of cultivar response to environmental cues provide a more objective, data-driven alternative for ME delineation and test location evaluation.

Statistical methodologies for MET analysis have advanced significantly in recent decades [8,9,10,11]. Traditional methods such as analysis of variance (ANOVA) and the GGE biplot have provided valuable insights, with the GGE biplot being widely adopted for its ability to visually integrate genotype main effects (G) and genotype-by-environment interaction effects (GE). This tool facilitates cultivar comparison, test location evaluation, and ME identification via the intuitive “which-won-where” pattern [12,13]. However, conventional GGE biplot analysis relies on raw phenotypic means and is grounded in fixed-effect model assumptions, including homogeneous error variances across test locations and balanced datasets—conditions that are rarely met in real-world regional trials due to practical constraints [4,14].

Real-world regional trial data often deviates from these statistical assumptions in two key ways: first, the “rolling” participation of cultivars (where poorly performing varieties are replaced with new candidates annually) creates unbalanced data structures with missing values across years and locations; second, the ecological diversity of test locations inherently leads to heterogeneous error variances, as environmental variability differentially impacts trial precision at each site. These complexities render traditional fixed-model analyses less accurate and potentially biased for evaluating test location representativeness, discriminating ability, and precise ME delineation.

To address these limitations, the integration of BLUP with the visualization strengths of the GGE biplot has emerged as a superior analytical framework. Based on linear mixed models, the BLUP method efficiently handles unbalanced data without missing value deletion or imputation, accounts for heterogeneous variances across environments, and provides accurate genotypic value estimates by treating genotypes as random effects—enabling shrinkage of estimates toward the mean and improved predictive accuracy [15,16]. Fitting multi-year data with mixed models allows BLUP to calculate adjusted means for each genotype-by-location combination, overcoming the challenges of missing data and producing values closer to true breeding values by correcting for environmental noise and trial error [17,18]. When these robust BLUP values are used to construct a cultivar × location two-way table for GGE biplot analysis, the results enable more reliable test location evaluation and scientifically sound ME delineation. A growing body of research on wheat[19], peanut[18], and forest tree[14] has confirmed that the BLUP-GGE biplot explains a greater proportion of total phenotypic variation and yields more reliable results for resource allocation and cultivar recommendation than analyses based on raw phenotypic means [20,21,22].

While the “which-won-where” polygon view of the GGE biplot is widely used for ME delineation, it is heavily influenced by the performance of a few high-yielding “super cultivars”, which can disproportionately affect test location grouping and fail to reflect the overall similarity of environmental effects on all cultivars [23,24]. For clustering environments into biologically meaningful MEs, the ME biplot—based on the projection and clustering of test location coordinates— is more appropriate and stable, as it focuses on the similarity of environmental effects on the entire set of genotypes and is less sensitive to outlier cultivars [25].

This study utilized five consecutive years (2021–2025) of yield data from national wheat METs conducted at 21 strategically distributed test locations across the HWWR, covering core ecological subregions including the Guanzhong Plain, Northern/Central/Eastern Henan Plains, Huaibei Plain, and Northern Jiangsu Plain. The extensive spatial coverage captures a wide spectrum of environmental conditions, providing a robust basis for analyzing G×E interactions and delineating MEs. The primary objectives of this study were to: (1) quantify the relative contributions of genotype (G), environment (E), and G×E interaction to yield variation using a linear mixed model; (2) calculate BLUP values for each genotype-by-location combination to construct an adjusted two-way table corrected for data imbalance and environmental noise; (3) systematically evaluate the discriminating ability and representativeness of test locations using the BLUP-GGE biplot; (4) delineate scientifically based MEs in the HWWR via the ME biplot applied to BLUP-adjusted data, minimizing the undue effect of superior genotypes on environmental clustering; and (5) identify superior cultivars with specific adaptation to each delineated ME.

The findings of this study are expected to provide a robust scientific basis for optimizing the HWWR regional trial network by identifying core representative test locations, enabling targeted cultivar promotion strategies, and supporting breeding programs for ME-specific varieties. Ultimately, this work aims to contribute to sustainable yield improvement and production stability in this vital agricultural region, and the integrated BLUP-GGE biplot approach presented here serves as a reference paradigm for MET analysis in other crop production systems facing similar ecological and experimental challenges.

2. Materials and Methods

2.1. Test Locations and Dataset

This study was based on grain yield data from national multi-environment wheat trials conducted in the HWWR from 2021 to 2025. Over the five-year experimental period, 256 unique genotypes (including the check cultivar Zhoumai 36) were evaluated: 132 cultivars were tested for two consecutive years, and 124 cultivars were tested for one year before exclusion from further trials. Trials were conducted annually at the same 21 test locations (Table 1), establishing a consistent and stable multi-environment testing network that covers the core ecological and geographical range of the HWWR.

At each experimental site, trials were established in representative fields using a randomized complete block design (RCBD) with three replications. Plot area was no less than 13.3 m² (actual sown area: 13.3–17 m²), and grain yield was determined by whole-plot harvesting to ensure accurate evaluation of actual yield potential. All test sites implemented a unified technical protocol in accordance with the National Wheat Variety Regional Trial Implementation Scheme. Sowing was conducted with a plot seeder between 10 and 20 October at a seeding density of 240–270 plants m⁻². Field management practices (tillage, irrigation, fertilization, pest and disease control) followed local high-yield agronomic standards to ensure trial results reflected the genetic potential of each cultivar, free from confounding effects of inappropriate management.

The test locations span four provinces in the HWWR and exhibit distinct ecological gradients: altitude ranges from 2.5 m (Sheyang) to 612.9 m (Baoji); topographic types include coastal flats, alluvial plains, and low hills; soil types comprise fluvo-aquic soil, Shajiang black soil, Lou soil, saline soil, and paddy soil; mean annual temperatures range from 13.1 to 15.3 °C; and annual precipitation ranges from 570 to 1075 mm. This diversity ensures the trial network captures the full range of environmental conditions governing wheat growth and development in the HWWR.

2.2. Statistical Analysis

The analytical strategy was designed to achieve two core objectives: (1) evaluate the discriminating ability and representativeness of test locations, and (2) delineate MEs based on G×E interaction patterns. The analysis employed a two-step approach integrating BLUP and GGE biplot analysis, with all statistical analyses performed using Genstat 24th edition (VSN International Ltd., Hemel Hempstead, UK) and GGEbiplot™ software (https://www.ggebiplot.com).

2.2.1. Estimation of Adjusted Means via BLUP

Given the unbalanced dataset arising from rolling cultivar participation across years, the BLUP method was used to derive adjusted means for each cultivar-by-location combination [14,15,19,26]. Based on a linear mixed model, this approach treats genotypic and environmental effects as random factors, effectively addressing unbalanced data and heterogeneous variances across environments [15]. BLUP values for each cultivar-location combination were extracted to construct an adjusted cultivar × location two-way table for subsequent GGE biplot analysis; this table contains predicted grain yields corrected for experimental design effects and unbalanced data structure.

The linear mixed model fitted via the restricted maximum likelihood (REML) algorithm is as follows [14]:

y_ijkl​ = μ + G_i​ + L_j + Y_k​ + GL_ij​ + GY_ik​ + LY_jk​ + GLY_ijk​ + R_l(jk)​ + ε_ijkl

(1)

Where, y_ijkl is the yield observation of the i-th cultivar in the j-th test location in the k-th year in the l-th replication, μ is the grand mean, G_i is the random effect of the i-th cultivar, the variety effect, L_j is the random effect of the j-th test location, Y_k is the random effect of the k-th year, GL_ij is the random interaction effect between the i-th cultivar and the j-th test location, GY_ik is the random interaction effect between the i-th cultivar and the k-th year, LY_jk is the random interaction effect between the j-th test location and the k-th year, GLY_ijk is the random interaction effect among the i-th cultivar, the j-th test location and the k-th year, R_l(jk) is the block effect within the location-year combination, and ε_ijkl is the residual error.

2.2.2. GGE Biplot Analysis

The adjusted two-way matrix of BLUP values for all cultivar-location combinations was subjected to GGE biplot analysis using GGEbiplot™ software, a specialized platform for visualizing MET data. The discriminating ability and representativeness of each test location were visualized and interpreted using the discriminating ability vs. representativeness view of the GGE biplot: discriminating ability was quantified by the length of the environmental vector extending from the biplot origin (a longer vector indicates a greater capacity to differentiate among cultivars) [23,25,27] , and representativeness was defined by the angle between the location vector and the average environment axis (AEA) (a smaller angle indicates higher representativeness of the average environmental conditions across the testing network) [28]. An “ideal” test location is characterized by the longest vector (highest discriminating ability) and a zero angle with the AEA (perfect representativeness). The desirability index of each location was depicted in the ideal tester view of the GGE biplot, where proximity to the ideal location marker correlates with higher desirability [28].

MEs were initially delineated using the traditional which-won-where polygon view of the GGE biplot [29], a method highly sensitive to the performance of top-yielding cultivars [24]. To alleviate this limitation, a test location clustering-based method was adopted for ME identification: test locations grouped into the same cluster were assigned to the same ME, as they induced consistent genotypic performance across all evaluated cultivars. This strategy provides a more stable and reliable foundation for ME delineation, essential for rational cultivar deployment. To clearly illustrate inter-ME relationships, an ME biplot was constructed [25]: the mean performance of each ME was represented by a vector from the biplot origin to the mean coordinate of all locations within the ME, and intra-ME variation was visualized by connecting the ME mean coordinate to the coordinate of each individual location

The GGE biplot model based on the cultivar-location matrix is expressed as[30]:

Y_ij=μ + β_j + λ₁ξ_i1η_j1 + λ₂ξ_i2η_j2 + ε_ij

(2)

Where, Y_ij is the trait mean value of genotype i in environment j; μ is the grand mean; β_j is the main effect of environment j; μ+β_j is the average yield of all genotypes in environment j; λ₁ and λ₂ are the singular values (SV) of the first and second principal components respectively; ξ_i1 and ξ_i2 are the eigenvectors of genotype i on the first and second principal components respectively; η_j1 and η_j2 are the eigenvectors of environment j on the first and second principal components respectively; ε_ij is the residual of genotype i in environment j. The GGE biplot takes the first principal component score (PC1) of varieties and test sites as the abscissa and the second principal component score (PC2) as the ordinate.

3. Results

3.1. Variance Component Analysis of Grain Yield via Linear Mixed Model

Variance components of grain yield in the multi-environment trials were decomposed using the REML method based on the linear mixed model (Table 2). Total phenotypic variation was primarily driven by environmental factors, with genotypic main effects and G×E interactions contributing comparatively little—a pattern consistent with previous crop regional trial studies [4,20]. Test location was the largest single source of variation, accounting for 40.5% of total variance with a significant effect (Z-ratio > 1.5), indicating that environmental heterogeneity (soil, climate, topography) among locations was the primary factor governing yield expression. The year effect accounted for 19.5% of total variation (marginally non-significant), while the year × location interaction was the second-largest source (30.7% of total variance, highly significant), implying that the combined year-location environment strongly modulated phenotypic performance. Collectively, the year effect and its interactions with location accounted for 55.1% of total variation, and the combined variance components of year, location, and their interactions exceeded 90% of total variation—confirming that temporal and spatial environmental effects were the dominant drivers of yield variation in the HWWR trials.

Notably, the Z-ratio for test location was 2.89 (significant), indicating substantial heterogeneous experimental error across locations, which was attributed to inherent ecological differences among sites (soil properties, climatic conditions, topography). This validates the rationality and necessity of using linear mixed models with heterogeneous variance structures for variance component analysis in this study.

For genotypic and interaction effects, the cultivar main effect accounted for 2.3% of total variation (highly significant), confirming genuine genetic divergence among the tested cultivars. The cultivar × location interaction contributed 2.1% of total variation (highly significant), demonstrating significant G×E interactions and specific environmental adaptability of cultivars across locations. The G×L×Y three-way interaction accounted for 4.8% of total variation (highly significant), further highlighting the complexity of environmental effects on yield. The block effect within location-year combinations constituted a small proportion of total variation, indicating that field management, plot layout, and experimental operations were standardized, random errors were effectively controlled, and experimental data quality was reliable.

3.2. Inter-Location Relationships: Raw Data vs. BLUP-GGE Biplots

The relation among testers view of GGE biplots was constructed for both raw grain yield data and BLUP values to compare their ability to resolve inter-location relationships (Figure 1).

The raw-data GGE biplot (Figure 1a) exhibited highly complex inter-location correlations, with both positive and negative vector relationships: vectors for LYG and SZ were nearly opposite to those for HY, XY, PY, LH, LY, and SOU; the BJ vector was almost diametrically opposed to those for YL, ZMD, SY, HA, and XZ. This indicates divergent cultivar performance between these opposing locations, with an apparent division into two tentative MEs (LYG, SZ, SQ, FY, ZMD; and all other locations), though relationships within each group remained intricate. SY exhibited an exceptionally long vector, suggesting strong influences from unique environmental conditions (e.g., coastal saline-alkali soil). PC1 and PC2 explained 36.0% and 26.4% of total variation, respectively (combined 62.4%), indicating a reasonable fit to the cultivar-location matrix but insufficient resolution for reliable ME delineation and test location evaluation. This suboptimal performance is attributed to the use of unadjusted phenotypic means, which fail to partition and eliminate year effects and other interaction effects, leading to environmental noise interference. For this reason, subsequent analyses exclusively used the BLUP-GGE biplot.

The BLUP-GGE biplot (Figure 1b) exhibited distinct characteristics relative to the raw-data biplot: all location vectors formed acute angles (< 90°), with the maximum angle (83.9°) between LYG and BJ, confirming positive correlations among all test locations. PC1 and PC2 explained 65.6% and 15.5% of total variation, respectively (combined 81.1%)—a substantial improvement in explanatory power, which enhances the reliability of cultivar selection, test location evaluation, and ME delineation. This improvement is due to the BLUP model’s ability to account for year, location, cultivar, and their interaction effects, producing phenotypic value estimates corrected for environmental noise and experimental error. While all locations showed positive correlations, weak correlations persisted between some sites, suggesting potential subregional clustering: relatively strong associations were observed among BJ, XY, YY, SY, and XZ; among LYG, SQ, FY, SZ, and HX; and among the remaining 12 locations. These clusters represent distinct potential MEs, providing a scientific basis for targeted cultivar selection in the HWWR.

3.3. Test Location Evaluation via BLUP-GGE Biplot

The discriminating ability vs. representativeness view (Figure 2a) and ideal tester view (Figure 2b) of the BLUP-GGE biplot were used to evaluate test location performance, with numerical values for representativeness, discriminating ability, and desirability index presented in Table 3.

In Figure 2a, ZMD lies exactly on the AEA, reflecting optimally high representativeness of the target ME; XX, HY, SOU, and LH also form small angles with the AEA, indicating excellent representativeness of HWWR ecological conditions. In contrast, BJ, XY, YY, SY, and XZ (below the AEA) and LYG, SQ, FY, SZ, and HX (above the AEA) form large angles with the AEA, suggesting poor representativeness. Discriminating ability was defined by vector length: SY, XY, YY, HX, BJ, and ZMD had long vectors (strong discriminating ability), while YL, HY, LY, and LH had short vectors (weak discriminating ability).

In the ideal tester view (Figure 2b), concentric circles serve as reference lines, with the central blue circle representing the ideal location marker. ZMD was the closest to the ideal marker, followed by SQU and PY, while LYG, BJ, SQ, and XY were the farthest, exhibiting moderately inferior comprehensive performance.

Table 3 presents the ranking of the 21 test locations for the three key indicators: (1) representativeness (best to worst): LY > ZMD > Y > XX > SQU > LH > PY > XMQ > HA > GY > YL > XZ > SZ > HX > SY > FY > YY > SQ > XY > LYG > BJ; (2) discriminating ability (strongest to weakest): SY > XY > YY > HX > BJ > ZMD > SZ > GY > SQU > PY > FY > XZ > SQ > HA > LYG > XX > XMQ > LH > LY > HY > YL; (3) comprehensive desirability index (best to worst): ZMD > SQU > PY > GY > XX > HA > LH > XZ > SY > XMQ > LY > HX > HY > SZ > YY > YL > FY > XY > SQ > BJ > LYG.

3.4. ME Delineation via the “Which-Won-Where” View of the BLUP-GGE Biplot

Despite positive inter-location correlations in the BLUP-GGE biplot (Figure 1b and Figure 2), weak correlations among some sites lead to inconsistent cultivar performance, justifying further ME delineation. The which-won-where polygon view (Figure 3a) was used for traditional ME delineation, with the corresponding ME biplot (Figure 3b) constructed to visualize inter-ME correlations and intra-ME variation.

The which-won-where view (Figure 3a) partitions the HWWR into three sectors (MEs), with the cultivar at each polygon vertex identified as the top-performing genotype for the corresponding sector (blue labels): (1) ME1 (14 locations: ZMD, SQU, GY, XX, HA, LH, XMQ, LY, HX, HY, SZ, YL, FY, SQ) – Guomai 169 as the superior cultivar; (2) ME2 (2 locations: PY, XZ) – Luomai 66 as the superior cultivar; (3) ME3 (4 locations: YY, SY, BJ, XY) – Lunxuan 148 as the superior cultivar. Remaining genotype markers are represented by “+” for clarity.

The ME biplot (Figure 3b) quantifies inter-ME correlations (vector angles) and intra-ME variation (distance from ME mean coordinate to individual locations): ME vectors are defined as rays from the origin to the mean coordinate of each ME’s locations, with smaller vector angles indicating stronger inter-ME correlations. The delineated MEs exhibited strong positive correlations: the ME2-ME3 angle was extremely small (16.5°), and the ME1-ME3 angle was 50.1° (maximum among all ME pairs). No significant differences in representativeness or desirability index were detected among the three MEs (Section 3.6), indicating insufficient inter-ME differentiation. This failure to meet the core ME delineation criterion—maximizing inter-ME variation and minimizing intra-ME variation—limits the practical value of the which-won-where view for targeted cultivar selection. Additionally, ME delineation via this method was highly sensitive to “super cultivars”: the inclusion or exclusion of top-performing genotypes significantly altered sector partitioning and ME classification, demonstrating poor robustness for multi-year multi-location data analysis.

Collectively, MEs delineated by the which-won-where view exhibit strong inter-ME correlations and insufficient differentiation, providing limited guidance for the targeted selection of cultivars with specific environmental adaptation— a core objective of wheat breeding and regional deployment in the HWWR.

3.5. ME Delineation via Test Location Clustering-Based Method

In contrast to the which-won-where view, test location coordinates in the BLUP-GGE biplot are determined by the aggregate performance of all evaluated cultivars, rendering them stable and less susceptible to distortion by extreme genotypes. A test location clustering-based method—using Euclidean distances among location coordinates in the BLUP-GGE biplot for hierarchical clustering—was therefore adopted for ME delineation, with the resulting ME biplot (Figure 4a) visualizing inter-ME relationships and intra-ME variation. An extended ME biplot (Figure 4b) was constructed to highlight cultivar-ME interaction effects, with locations depicted as small circles and selected representative cultivars labeled for clarity.

The clustering-based ME biplot (Figure 4a) stably partitions the HWWR into three distinct MEs, with ME3 consistent with that delineated by the which-won-where view: (1) ME1 (5 locations: HX, SZ, FY, SQ, LYG); (2) ME2 (12 locations: ZMD, SQU, GY, XX, HA, LH, XMQ, LY, HY, YL, PY, XZ); (3) ME3 (4 locations: SY, YY, XY, BJ). Inter-ME differentiation was markedly enhanced relative to the which-won-where view: the ME1-ME3 vector angle was 68.5° (near right angle, weak positive correlation), ME2 occupied a transitional position between ME1 and ME3, and the ME1-ME2 (31.4°) and ME2-ME3 (35.2°) angles were all larger than the corresponding angles in the which-won-where ME biplot (Figure 3b). These results confirm that the clustering-based method produces MEs with substantially improved inter-ME differentiation and greater robustness.

The extended ME biplot (Figure 4b) reveals clear cultivar-ME interaction effects, with each ME having an optimally adapted cultivar: Guomai 169, Zhengmai 33, Zhengmai 28, and Luomai 66 exhibited superior performance in ME2; Lunxuan 148 performed exceptionally well in ME3 but moderately in ME1; Tianyikemai 10 showed the inverse pattern (excellent in ME1, poor in ME3); Kejujin 6 performed suboptimally in ME1; and Sanfengmai 616 performed poorly in ME3. While most genotype markers are represented by “+” for clarity, numerous other cultivars exhibited analogous interaction patterns. These observations firmly confirm significant cultivar-ME interaction effects, which form the fundamental basis for ME delineation in wheat METs.

3.6. Comparative Analysis of ME Delineation Methods

To quantitatively compare the performance of the which-won-where view and the test location clustering-based method for ME delineation, we summarized the key quantitative indicators (discriminating ability, representativeness, desirability index) for the delineated MEs (Table 4) and calculated inter-ME correlation coefficients and vector angles (Table 5). Multiple comparisons of means were conducted using the least significant difference (LSD) method (p < 0.05 and p < 0.01).

For MEs delineated by the which-won-where view (Table 4), no significant differences were detected in representativeness or desirability index among the three MEs; only the discriminating ability of ME3 was significantly higher than ME2 (p < 0.05) and highly significantly higher than ME1 (p < 0.01). This indicates that representativeness and overall desirability are comparable across MEs, with only ME3 exhibiting relatively strong discriminating ability. Since discriminating ability is determined by vector length—a parameter weakly associated with sector partitioning in the which-won-where view—this method produces MEs with limited inter-ME differentiation, constraining its theoretical and practical utility for the HWWR.

In contrast, MEs delineated by the clustering-based method (Table 4) exhibited highly significant differences in all three key indicators: (1) discriminating ability: ME3 >> ME1 = ME2 (p < 0.01); (2) representativeness: ME2 >> ME1 = ME3 (p < 0.01); (3) desirability index: ME2 >> ME1 (p < 0.01) and ME2 > ME3 (p < 0.05). These results demonstrate clear functional characterization of the clustering-based MEs: ME2 (12 locations) has the highest comprehensive desirability and strongest representativeness (core HWWR production area); ME3 (4 locations) has the strongest discriminating ability (distinct environmental characteristics); ME1 (5 locations) has moderate overall performance (low representativeness, weak discriminating ability).

Table 5 further quantifies inter-ME differentiation via correlation coefficients and vector angles: the which-won-where view produces small vector angles (16.5° to 51.1°) and high correlation coefficients, confirming strong inter-ME correlations; the clustering-based method produces substantially larger vector angles (31.4° to 68.5°) and lower correlation coefficients, with the maximum ME1-ME3 angle approaching 70°. These results confirm that the clustering-based method achieves the core ME delineation criterion of maximizing inter-ME divergence, yielding more robust and reliable ME classification for wheat METs in the HWWR.

4. Discussion

4.1. Superiority of the BLUP-GGE Biplot in Analyzing Multi-Environment Variety Trial Data

Accurate dissection of G×E interactions is central to interpreting crop regional trial data, as the choice of statistical methodology directly dictates the reliability of test location evaluation and ME delineation efficacy [4,20]. Traditional GGE biplot analysis relies on raw phenotypic means and is constrained by fixed-effect model assumptions (homogeneous error variances, balanced datasets) that are rarely met in real-world METs [14,19,31]. In this study, the raw-data GGE biplot explained only 62.4% of total phenotypic variation and exhibited erratic inter-location correlations (positive and negative), which hinders reliable location evaluation and ME classification. This suboptimal performance arises because raw phenotypic means conflate year, location, cultivar, and their interaction effects, leading to severe environmental noise interference and obscuring true genetic differences among cultivars and ecological correlations among test locations—consistent with recent studies demonstrating that raw-data GGE biplots are prone to bias in complex multi-year multi-location trial systems [17,20].

The integration of BLUP with GGE biplot analysis effectively overcomes these limitations by leveraging the strengths of linear mixed models [15,16]. Treating genotypes and environments as random effects, the BLUP framework efficiently accommodates unbalanced data from rolling cultivar participation and fully decomposes/adjusts for year, location, and G×E interaction effects. In this study, the BLUP-GGE biplot explained 81.1% of total phenotypic variation (an 18.7% increase over the raw-data biplot) and all test location vectors formed acute angles (maximum = 83.9°), confirming consistent positive inter-location correlations and the ecological homogeneity of the HWWR. BLUP estimates attenuate environmental noise and experimental error, more accurately capturing inherent inter-location relationships and true genotypic performance—aligning with Liu et al.[21], who reported superior phenotypic variation explanatory ability for the BLUP-GGE biplot relative to conventional GGE biplots.

For test location evaluation, the BLUP-GGE biplot clearly discriminated the discriminating ability and representativeness of each 21 test locations, identifying ZMD as the most ideal location (highest desirability index), followed by SQU and PY, while LYG, BJ, SQ, and XY exhibited inferior comprehensive performance. This provides a clear quantitative basis for optimizing the HWWR regional trial network, including strengthening core location construction and adjusting/replacing low-desirability locations. Collectively, these results confirm that the BLUP-GGE biplot is more accurate and robust than the raw-data GGE biplot for multi-year multi-environment trial analysis, and it is the more appropriate statistical tool for test location evaluation in the ecologically complex HWWR.

4.2. Limitations of the “Which-Won-Where” View for Multi-Year Multi-Location ME Delineation

The “which-won-where” view is the most extensively applied method for ME delineation in crop METs, with successful applications in oat[32], wheat[33], cotton[25,34], maize[35], soybean[36,37], and other crops[38,39]. However, this study identifies clear limitations when the method is directly applied to multi-year multi-location trial data in the HWWR, stemming from its original design for single-year multi-location data analysis [23,24].

First, the which-won-where view cannot effectively eliminate year effects and other confounding sources of variation. The cultivar-by-location two-way table from multi-year trials integrates cumulative temporal effects, and the lack of year-effect decomposition obscures the genuine ecological characteristics of test locations. In this study, MEs delineated by this method exhibited strong inter-ME correlations (vector angles 16.5°–51.1°) and no significant differences in representativeness or desirability index, failing to meet the core ME delineation criterion of maximizing inter-ME variation. This renders the delineation outcomes of limited practical value for targeted cultivar selection and deployment in the HWWR.

Second, the method is highly susceptible to interference from “super cultivars” (top-yielding genotypes). Designating the polygon vertex cultivar as the top-performing genotype for each sector means the inclusion/exclusion of these cultivars can substantially alter test location sector partitioning and ME delineation results. In this study, ME classification via the which-won-where view varied with changes in core super cultivars, demonstrating poor robustness for multi-year multi-location data analysis.

Conventional improved methods for multi-year multi-location data analysis also have critical limitations: the “year-by-year analysis with multi-year induction” method [24] is extremely cumbersome in practical operation, and substantial annual environmental variability in the HWWR leads to poor repeatability and practicality of year-by-year ME delineation results; the location grouping method proposed by Yan(2019) [23] suffers from high annual variation and unclear regularity, with successful applications limited to a small number of studies [25,40], making it difficult to popularize in the ecologically complex HWWR.

The fundamental reason for these limitations is that the which-won-where view is predicated on a relatively homogeneous annual environmental background (single-year data), and it fails to accommodate the complex G×L×Y three-way interaction in multi-year datasets [23,24]. This leads to the confounding of multiple temporal and spatial effects, distorting ME delineation results. Collectively, direct application of the which-won-where view to multi-year multi-location wheat trial data in the HWWR does not yield reliable ME delineation outcomes, highlighting the need to optimize and refine ME delineation methods based on the unique characteristics of multi-year trial data.

4.3. Advantages of the Test Location Clustering-Based Method for ME Delineation

To address the limitations of the which-won-where view, this study proposed a test location clustering-based method for ME delineation in the HWWR: constructing an ME biplot via hierarchical clustering of test locations, using a Euclidean distance matrix derived from location coordinates in the BLUP-GGE biplot. This method is grounded in the BLUP-GGE biplot framework, inheriting the key advantage of BLUP values—effective attenuation of environmental noise and accurate capture of inherent ecological correlations among test locations. Most importantly, it abandons the sector-partitioning logic of the which-won-where view (susceptible to super cultivar interference) and adopts the overall similarity of test location responses to all evaluated cultivars as the clustering criterion. This fundamental methodological shift significantly enhances the stability and reliability of ME delineation outcomes for complex multi-year multi-location trial data.

The clustering-based ME biplot exhibited three distinct advantages in this study: (i) Enhanced inter-ME differentiation: The method produced substantially larger inter-ME vector angles (31.4°–68.5°) than the which-won-where view (16.5°–51.1°), with the ME1-ME3 angle approaching 70° (near right angle, weak positive correlation). Multiple comparisons confirmed significant differences in all key indicators (discriminating ability, representativeness, desirability index) among the three MEs, forming clearly functionally characterized MEs aligned with the core ME delineation criterion. (ii) High robustness to super cultivar interference: Clustering is based on the overall coordinate characteristics of test locations in the BLUP-GGE biplot, which integrate the comprehensive responses of all cultivars to the test location environment. This avoids the one-sidedness of the which-won-where view, which focuses solely on top-performing cultivars, and ensures ME delineation results are not distorted by extreme genotypes. (iii) Clear cultivar-ME interaction effects: The method explicitly reveals cultivar-ME interactions, with each ME having its own optimally adapted cultivars (Guomai 169/Zhengmai 33 in ME2, Lunxuan 148 in ME3, Tianyikemai 10 in ME1). This clear cultivar-ME matching relationship provides a direct scientific basis for the precise deployment and promotion of wheat cultivars in the HWWR, enabling the maximization of cultivar genetic potential and regional yield levels.

By utilizing BLUP values as the data foundation and hierarchical clustering as the core analytical approach, the test location clustering-based ME biplot effectively resolves the issue of confounding effects in multi-year multi-location data analysis and compensates for the deficiencies of the which-won-where view [29], Yan’s location grouping method [23], and the “year-by-year analysis with multi-year induction” method [24,41]. Beyond its applicability to wheat regional trials in the HWWR, this method offers a novel reference for ME delineation in other crops characterized by complex ecological conditions and long-term regional trial datasets [22,42], holding significant popularization and application value in crop breeding and regional trial research.

4.4. Practical Implications for Optimizing the HWWR Wheat Regional Trial System

The findings of this study hold important practical implications for optimizing the wheat regional trial system in the HWWR, a core wheat-producing area in China, with recommendations for test location layout, ME-based cultivar promotion, and statistical methodology application:

1. Test location layout optimization: The BLUP-GGE biplot identified core test locations with high discriminating ability and strong representativeness (ZMD, SQU, PY). We recommend strengthening the construction of these core locations, standardizing trial management protocols, and ensuring the stability and reliability of trial data to serve as benchmark sites for the HWWR trial network. For test locations with low comprehensive desirability (LYG, BJ, SQ, XY), we suggest three strategies: (i) adjust trial layout to improve ecological representativeness; (ii)optimize field management practices to enhance discriminating ability; or (iii)replace with more representative alternative locations. These measures will improve the efficiency and cost-effectiveness of the HWWR regional trial network.

2. ME-based precise cultivar promotion: The three MEs delineated by the clustering-based method encompass the main ecological types of the HWWR, enabling targeted cultivar deployment based on ME characteristics: (i) ME2 (12 core locations): the most favorable comprehensive ecological conditions, suitable for promoting high-yield, widely adaptable cultivars (Guomai 169, Zhengmai 33); (ii) ME3 (4 locations with distinct environmental characteristics e.g., coastal saline-alkali soil): strong discriminating ability, ideal for breeding and promoting stress-tolerant, specialized cultivars (Lunxuan 148); (iii) ME1 (5 locations with low representativeness): conduct targeted cultivar selection tailored to local ecological conditions (e.g., Tianyikemai 10 for ME1 adaptation). The precise matching of cultivars to their optimal MEs will effectively unlock cultivar genetic potential and elevate the overall wheat yield level in the HWWR.

3. Statistical methodology popularization: We recommend popularizing the BLUP-GGE biplot in HWWR wheat regional trials, replacing the conventional raw-data GGE biplot, to improve the accuracy of test location evaluation and cultivar selection. Additionally, the test location clustering-based ME biplot should be adopted as the preferred method for ME delineation of multi-year multi-location trial data in the HWWR, providing a more scientific basis for the rational spatial layout of wheat cultivars [4,20]. These methodological improvements will enhance the scientific rigor and practical value of the HWWR wheat regional trial system, supporting the development and promotion of ME-adapted wheat cultivars.

5. Conclusions

Based on national wheat regional trial yield data from 21 test locations in the Huanghuai winter wheat region (2021–2025), this study systematically compared the performance of raw-data GGE biplots and BLUP-GGE biplots for MET analysis, and explored the applicability of the traditional “which-won-where” view and a novel test location clustering-based method for mega-environment (ME) delineation. The BLUP-GGE biplot effectively eliminated environmental noise and experimental errors by decomposing year, location, and genotype-by-environment (G×E) interaction effects, explaining 81.1% of total phenotypic variation (a substantial improvement over the raw-data GGE biplot’s 62.4%) and providing more reliable test location evaluation results—confirming its superiority for analyzing unbalanced multi-year multi-environment trial data in the ecologically complex HWWR. ZMD was identified as the most desirable test location, followed by SQU and PY, while LYG, BJ, SQ, and XY exhibited relatively poor comprehensive performance, providing a quantitative basis for trial network optimization.

The traditional “which-won-where” view exhibited obvious limitations for multi-year multi-location data analysis, including failure to eliminate confounding temporal effects, high susceptibility to super cultivar interference, and insufficient inter-ME differentiation (vector angles 16.5°–51.1°), which constrains its practical utility for targeted cultivar selection. In contrast, the proposed test location clustering-based method—grounded in the BLUP-GGE biplot and Euclidean distance hierarchical clustering—effectively overcomes these limitations, stably delineating the HWWR into three distinct MEs (ME1: HX, SZ, FY, SQ, LYG; ME2: ZMD, SQU, GY, XX, HA, LH, XMQ, LY, HY, YL, PY, XZ; ME3: SY, YY, XY, BJ) with markedly enhanced inter-ME differentiation (maximum vector angle ≈70°). These MEs exhibited significant differences in key performance indicators (discriminating ability, representativeness, desirability index) and clear cultivar-ME interaction effects, with each ME having an optimally adapted cultivar (Guomai 169/Zhengmai 33 in ME2, Lunxuan 148 in ME3, Tianyikemai 10 in ME1).

This study establishes a comprehensive technical system for test location evaluation and ME delineation in the HWWR based on the BLUP-GGE biplot, optimizing the ME delineation method for multi-year multi-location trial data and offering a new path for the scientific optimization of wheat regional trial systems and precise cultivar promotion in the region. The proposed test location clustering-based ME biplot also has important reference significance for multi-environment trial analysis and ME delineation of other crops with complex ecological conditions, contributing to the improvement of crop breeding and regional trial research efficiency worldwide.