Crop variety development is a fundamental pillar of modern agriculture, with advancements over the past century playing a vital role in enhancing global food security, improving farmer livelihoods, and promoting sustainable farming practices (Qaim, 2020; Begna and Begna, 2021). By creating high-performing, adaptable cultivars that are resilient to biotic and abiotic stresses, researchers have enabled stable and abundant crop yields, bolstered food supplies, and supported the economic wellbeing of farming communities, while also facilitating the adoption of other sustainable agricultural innovations (Zsögön et al., 2022; Swarup et al., 2021). Continued investment and innovation in this field will be crucial as the world navigates the complex challenges of ensuring long-term food security and environmental sustainability (Blakeney, 2022).

Multi-environment trials (METs) are a crucial element of the crop variety development pipeline. In this process, newly bred crop genotypes are evaluated across a range of agro-ecological environments (Smith et al., 2001b, 2021; Brown et al., 2020). The primary purpose of METs is to capture the influence of diverse environmental factors on the expression of genotypic performance. This allows for the assessment of G × E effects- a critical consideration in identifying high-performing and stable crop varieties that can perform well under a range of conditions (Verbyla, 2023; Lee et al., 2023; Lisle et al., 2021).

Understanding G × E is a fundamental aspect of plant breeding, as it provides insights into how different genotypes respond to diverse environmental conditions. This interaction is critical because it directly affects important agronomic outcomes, such as yield and the ability of plants to adapt to varying climates and soil types (Van Eeuwijk et al., 2016; Malosetti et al., 2013, 2016; Bustos-Korts et al., 2009). This is also demonstrated in studies by Cooper et al. (2020); Mark et al. (2021) and Cooper et al. (2022) that focus on implementing management strategies to maximize crop productivity through combinations of Genotype–Management (G–M) technologies.

By testing new genotypes across multiple environments, breeders can gain valuable insights into how a variety's traits manifest and interact with the local environmental context. This information is essential for selecting cultivars that exhibit both high productivity and reliability, making them suitable for deployment across a wide geographical area (Smith et al., 2001b, 2019, 2021).

The conventional approach for analyzing MET data has relied on ANOVA-based methods, such as AMMI (additive main effects and multiplicative interaction) and GGE (genotype and genotype-by-environment interaction) analysis (Beeck et al., 2010; Zhang et al., 2020). These techniques have enabled researchers to obtain ANOVA results and gain insights into interactions. Additionally, they have facilitated data visualization through the use of biplot techniques.

However, these methods exhibit several inherent limitations that have become increasingly apparent as the complexity of MET experiments has grown. One key issue is that ANOVA-based methods may struggle to handle unbalanced and incomplete data structures, where some observations are missing or not all genotypes are present in each environment (Piepho, 1997; Smith et al., 2001a; Kelly et al., 2007). This is a common challenge in MET studies. Additionally, these methods may fail to adequately account for non-genetic sources of variance associated withenvironmental factors or experimental design effects, which can significantly influence the observed phenotypic performance (Smith et al., 2005). Furthermore, the ANOVA-based analysis commonly follows a two-stage modeling approach. This involves first analyzing the data within each environment, followed by a second stage of combining the results across environments for the G × E interaction analysis (Piepho et al., 2012; Smith et al., 2005). This approach can lead to a loss of information and less effective extraction of insights from the complex MET datasets.

Linear mixed model-based approaches have emerged as a more efficient methodology for the analysis of MET data. Recently, these mixed model approaches have become predominant, as they provide a flexible framework that can easily handle incomplete data and appropriately model the non-genetic variances between and within environments. This includes accounting for spatial variation within environments and error variance heterogeneity across environments (Smith et al., 2001a). The error variation within environments can be modeled using the approach of Gilmour et al. (1997), which appropriately models the three patterns of spatial trends associated with the field of trials: local, extraneous, and global trends.

Smith et al. (2001a, 2005) extended the G × E analysis by employing Factor Analytic Multiplicative Mixed (FAMM) models. These models assume random genotype effects and fixed environment effects, and use a one-stage analysis approach, where the models for residual effects are estimated simultaneously with the models for G × E effects. A key feature of the Factor Analytic (FA) model for MET data is its importance for the estimation of the associated variance structure for G × E effects. The FA model provides a good and parsimonious approximation to the unstructured form, and is generally more computationally robust (Kelly et al., 2007). Additionally, the best linear unbiased predictions (BLUPs) of the G × E effects and estimations of loadings and scores can be obtained, enabling bi-plot analysis to better understand the G × E interaction patterns.

The objective of this research is to explore the potential of linear mixed model-based approaches, particularly FAMM models, as a high-performing alternative to the traditional randomized complete block (RCB) designs analysis for effectively analyzing MET data and extracting meaningful insights from the complex G × E interactions. By leveraging linear mixed model-based methodology, this research aims to demonstrate how the non-genetic variances within and between environments can be appropriately modeled, including the use of the Gilmour et al. (1997) approach for capturing spatial trends and error variation. Furthermore, the application of FAMM models is explored, as these approaches provide advantages in estimating the variance structure of G × E effects and enabling more informative visualizations. The ultimate goal is to provide researchers and breeders with a robust and efficient analytical tool for extracting meaningful insights from MET data, ultimately supporting the development of high-performingand adaptable crop varieties.

The spatial analysis the individual trial data revealed several key insights. The presence of significant local, extraneous and global spatial variation, as detected by the Wald and REMLR tests, highlights the critical importance of accounting for these spatial effects when analyzing data from multi-environment trials. The positive estimated AR1 correlations, indicating that nearby plots tend to be more similar than those located further apart, underscores the need to incorporate appropriate spatial covariance structures in the linear mixed models used for these analyses (Cullis et al., 1998; Stefanova et al., 2009).

The systematic spatial trends observed in the residual plots (Figure 1A) for the 2020MZNVT-AS trial clearly demonstrate the value of including random row and column effects, as well as fixed linear global trends, in the modeling approach. This aligns with the recommendations from previous studies on the analysis of multi-environment trials, which have emphasized the importance of capturing both local and global spatial patterns to improve the accuracy and reliability of genotypic evaluations (Möhring and Piepho, 2009; Piepho et al., 2012; Smith et al., 2005).

The findings from the current study are consistent with the work of Qiao et al. (2000), who highlighted the importance of evaluating experimental designs and spatial analyses in wheat breeding trials. Similarly, the application of enhanced statistical models for the analysis of multi-environment trials in finger millet (Tesfaye et al., 2023) and common bean (Argaw et al., 2024) underscores the broader relevance and applicability of the spatial modeling approaches employed in this study across various crop varietal evaluation trials.

Furthermore, the use of spatial mixed models for varietal selection field trials, as discussed by Woldemeskel and Fenta (2022), reinforces the practical utility of the methods demonstrated here. The insights gained from the spatial analysis, such as the identification of significant local and global spatial variation, can directly inform the selection of high-performing genotypes and the optimization of experimental designs in future trials. This can lead to more accurate assessments of genotypic performance and ultimately contribute to the development of improved crop varieties (Kelly et al., 2007).

The results of the spatial + G × E analysis across the multiple datasets provide valuable insights into the underlying patterns of the G × E effects. One of the key findings is the consistent superiority of the factor analytic (FA) modeling approach over the conventional RCB analysis in capturing the G × E variance. As the order of the FA models increased, meaning more factors were incorporated, the total percentage of G × E variance explained also increased (Table 3). This observation aligns with the theoretical foundations and empirical evidence presented in previous studies (Smith et al., 2001a, 2005; Burgueño et al., 2008; Meyer, 2009). This supports the notion that higher-order FA models are better able to effectively represent the complex and multidimensional nature of G × E effects.

However, the optimal FA model order was found to be dataset-dependent, as demonstrated by the ACB14-16 dataset, where the intermediate FA-3 and FA-4 models did not reach statistical significance, while the lower-order FA-2 and higher-order FA-5 models were significant (Table 1). This observation suggests that the appropriate model complexity may vary based on the specific characteristics of the data, such as the magnitude and structure of the G × E effects. This scenario could be associated with the implications of trial locations. Even with the best analytical methods, poorly distributed trials can complicated the MET data analysis. This leads us to the concept of the target population of environments (TPE). The TPE framework emphasizes the need for trials to be representative of the environments where the genotypes will be deployed. Misalignment between trial locations and the intended target environments can result in inaccurate assessments of genotypic performance (Cooper et al., 2020, 2021; Hajjarpoor et al., 2022). However, the identification of the optimal FA model order is a crucial step, as emphasized by De Resende and Thompson (2004) and Smith et al. (2015) as it ensures the most parsimonious yet informative representation of the G × E structure.

Compared to the standard RCB analysis, the spatial analysis improved estimates of genetic variance, while the spatial + G × E analysis led to even greater improvements in genetic variance estimates for most datasets. This potential can be attributed to the spatial + G × E analysis's ability to extract the genetic effect that is confounding with the residual effect and the effects stored in the G × E interactions. Notably, the spatial analysis reduced estimates of error variance, and the spatial + G × E analysis (Table 4) was able to substantially minimize residual variability, outperforming both the RCB and spatial-only methods.

The integration of spatial variability into the G × E analysis through the spatial + G × Emodeling approach has led to substantial improvements in heritability estimates across all datasets compared to the RCB and spatial analyses (Figure 2). This supports the findings of Smith et al. (2015), who demonstrated the ability of factor analytic mixed models to enhance the estimation of genotypic effects and heritability. Interestingly, the extent of improvement was more pronounced in datasets with relatively larger data volumes, such as DCP14-16 and LCB19-20, compared to smaller datasets like LCB14-16 and IHMZ20. This observation highlights the importance of data quality and quantity in effectively capturing the underlying patterns of G × E interactions and larger datasets provide more information about the complex G × E structure, allowing the factor analytic models to better estimate the genotypic effects and associated heritability. The advancement lies in demonstrating the superiority of the factor analytic model across different sizes of MET datasets, particularly in its ability to effectively investigate G × E interactions. This approach allows for a comprehensive understanding of how varying dataset sizes influence the model's performance to uncover underlying patterns in G × E relationships.

The visual representations, including the genetic correlation heat map and dendrogram, offer intuitive and informative insights into the trial relationships within each dataset. The observed patterns of strong positive, negative, and weak correlations in the heat maps (Figure 3), as well as the distinct trial clusters identified in the dendrograms (Figure 4), reflect the complex structure of the G × E interactions, as discussed in previous studies (Cullis et al., 2010; Argaw et al., 2024). These visual aids can assist breeders and researchers in better understanding the spatial and genotypic patterns, which can inform decision-making processes related to cultivar selection, target environment identification, and trial design optimization.

Our research constitute a noteworthy contribution to the field by demonstrating the effectiveness of mixed model-based approaches, particularly spatial + G × E analysis, in enhancing MET analysis. Unlike previous studies, we provide comparative insights across various MET datasets, revealing that the optimal order of the FA model for G × E effects varies depending on the specific dataset. Furthermore, our application of advanced linear mixed model methods to Ethiopia's crop MET data—characterized by notable G × E interactions and plot variability—offers critical insights that are highly relevant to crop breeding programs in the region.

The spatial analysis of the individual trial data revealed significant local, extraneous, and global spatial variation, underscoring the need to use mixed model spatial analysis for multi-environment trial data analyses. The positive estimated AR1 correlations indicate that nearby plots tend to be more similar, further emphasizing the importance of capturing local spatial patterns to improve the accuracy and reliability of genotypic evaluations. Employing these spatial modeling techniques in future trials is recommended to account for the inherent spatial variability, which can enhance the genetic gain achieved through variety selection.

The results of the spatial + G × E analysis showcased the value of factor analytic (FA) modeling in effectively representing the complex and multidimensional nature of G × E effects. The consistent superiority of the FA approach over the traditional RCB analysis supports the use of FA models, though the optimal model order may vary across datasets. The approach is significant because it applies the FA model to MET datasets of varying sizes, offering insights into GxE effects, and improves our understanding of how dataset size impacts model performance, enabling more accurate identification of underlying patterns in G × E relationships.

Researchers are encouraged to explore the application of these spatial + G × E modeling techniques across a broader range of crop breeding evaluation programs to further validate their benefits in increasing the efficiency of variety selection and improving genetic gain.

The integration of spatial variability into the G × E analysis led to substantial improvements in heritability estimates across all datasets. This highlights the ability of factor analytic mixed models to enhance the estimation of G × E effects and heritability, with more pronounced improvements observed in larger datasets. By accurately partitioning the sources of variation, these advanced statistical models can lead to more reliable predictions of genotypic performance, ultimately contributing to greater genetic gain in crop improvement programs.

The insights gained from this study can contribute to the development of improved crop varieties and enhanced experimental designs in terms of establishing good trailing system that optimize statistical precision associated with the size future trials. Continued research and adoption of advanced spatial and G × E modeling techniques in crop breeding and evaluation programs will be crucial for driving progress in the creation of high-performing agricultural cultivars. By accounting for the complex spatial and G × E structures inherent in multi-environment trials, researchers and breeders can improve the estimation of genetic gain and make more informed decisions to deliver high-performing crop varieties to farmers and consumers.