Vol.:(0123456789) Theoretical and Applied Genetics (2025) 138:241 https://doi.org/10.1007/s00122-025-05023-6 ORIGINAL ARTICLE Prediction of Australian wheat genotype by environment interactions and mega‑environments Nick S. Fradgley1   · Guillermo S. Gerard2   · Velu Govindan2   · Julie M. Nicol3   · Amit Singh3   · Wuletaw Tadesse4 · Alexander B. Zwart1 · Richard Trethowan3   · Ben Trevaskis1   · Alex Whan1   · Jessica Hyles1  Received: 18 February 2025 / Accepted: 6 August 2025 / Published online: 4 September 2025 © The Author(s) 2025 Abstract Key message  Latent environmental effects of genotype by environment interactions could be predicted from observed environmental covariates. Predictions into the wider target population of environments revealed greater insights. Abstract  Wheat is grown across a diverse range of environments in Australia with contrasting environmental constraints. Targeted breeding to optimise genotypes in target environments is hindered by large and ubiquitous genotype by environment interactions (GEI). Common GEI in multi-environment trial experiments, which sample the target population of environ- ments, can be efficiently modelled using latent environmental effects from factor analytic mixed models. However, gener- alised prediction into the full target population of environments is difficult without a clear link to observed environmental covariates (ECs) that are defined from high-resolution weather and soil data. Here, we used a large wheat multi-environment trial dataset and demonstrated that latent environmental effects can be associated with and predicted from observed ECs. We found GEI-based environment classes could be defined by combinations of key ECs. Prediction of main and latent effects in a wider set of environments covering the full TPE across the Australian grain belt over 13 years revealed the complex trends of environmental effects and GEI over regional scales demonstrating high year-to-year variability. Regional environment types often shifted year-to-year. Cross-validation of forward genomic prediction into untested year environments demonstrated that increased accuracy is possible if estimated genetic effects are also accurate and ECs of new environments are known. These findings may guide Australian wheat breeders to better target specifically adapted material to mega-environments defined by static GEI while also considering broad adaptability and non-static GEI resulting from year-to-year variability. Introduction Australia is a major wheat producing nation. Between 2010 and 2022, Australia contributed between 5.1 and 15.4% of the 145.7 to 200 million tonnes of global exports (FAO 2024). Along with improved agronomy and crop manage- ment, significant genetic gains achieved by plant breeders have increased wheat yields in key Australian growing regions (Kirkegaard 2019; Sadras Lawson 2011), but pro- duction environments are highly variable across locations and seasons. Abiotic factors such as water limitation and high temperature stress (Flohr et al. 2017) combined with biotic factors, including soil-borne pathogens (Thompson et al. 1999) limit yields, and extensive genotype by envi- ronment interactions (GEI) limit breeders' selections across multiple environments (Basford and Cooper 1998; Cullis et al. 2000). Increased temperatures due to changing cli- mates are likely to reduce the potential of widely adapted varieties and reduce selection efficiency, hence requiring greater diversity of adapted breeding materials (Xiong et al. 2024). Wheat breeders in Australia work with a relatively refined gene pool and a few major varieties often dominate the acre- age and market share (Brennan and Fox 1998; Cockram Communicated by Philomin Juliana. * Nick S. Fradgley nick.fradgley@csiro.au 1 CSIRO Agriculture and Food, GPO Box 1700, Canberra, ACT​ 2601, Australia 2 International Maize and Wheat Improvement Center (CIMMYT), KM 45 Carretera Mexico‑Veracruz, 56237 Texcoco, Mexico 3 The Plant Breeding Institute, University of Sydney, 107 Cobbity Road, Cobbity, NSW 2750, Australia 4 International Center for Agricultural Research in the Dry Areas (ICARDA), PO Box 6299, Rabat, Morocco http://crossmark.crossref.org/dialog/?doi=10.1007/s00122-025-05023-6&domain=pdf http://orcid.org/0000-0002-7868-8795 http://orcid.org/0000-0002-9112-3588 http://orcid.org/0000-0001-9502-4352 http://orcid.org/0009-0009-4578-4105 http://orcid.org/0000-0002-5984-3622 http://orcid.org/0000-0003-0105-875X http://orcid.org/0000-0001-7277-4994 http://orcid.org/0000-0002-6839-2915 http://orcid.org/0000-0002-1836-4334 Theoretical and Applied Genetics (2025) 138:241241  Page 2 of 24 2024), which is in contrast to other major wheat producing nations such as the US (Chai et al. 2022). The International Maize and Wheat Improvement Centre (CIMMYT) and the International Centre for Agricultural Research in the Dry Areas (ICARDA) are large international wheat breeding pro- grammes that collaborate with Australian wheat breeders and researchers to introduce and evaluate new sources of yield potential and adaptation to biotic and abiotic stresses as part of the CIMMYT Australia ICARDA Germplasm Evalu- ation (CAIGE) programme (Trethowan et al. 2024). Unpredictability of GEI is a long-standing problem that has limited targeted breeding for specific environments as well as resilience to year-to-year variability. Many statistical modelling approaches have been developed to understand and predict GEI in multi-environment trial (MET) experi- ments (Crossa et al. 2024; Malosetti et al. 2016). Whereas early methods simply regressed genotype performance against an environmental index defined by the environment mean yields (Finlay and Wilkinson 1963; Yates and Cochran 1938), more recent approaches have sought to replace this environmental index with biologically meaningful observed environmental covariates (ECs) so that prediction into untested environments is possible. As an analogy to genotyping, methods for ‘envirotyping’ derive large numbers of ECs from high-resolution weather data to characterise crop growing environments (Costa-Neto et al. 2021; Resende et al. 2021; Xu 2016), and are often informed by mechanistic crop growth models (Cooper et al. 2020; Zhang et al. 2012). Factorial regression models can incorporate genotype-specific responses to individual ECs (Boer et al. 2007; Piepho et al. 1998), but are problematic to fit when more than a few ECs need to be used. Rogers and Holland (2022) used principal component dimension reduc- tion of large EC matrices, while multiple ECs can be linearly combined to define multiple synthetic covariates for GEI terms in similar mixed model frameworks (Piepho 2022; Piepho and Blancon 2023; Tadese et al. 2024). In a con- trasting approach, Jarquín et al. (2014) developed reaction norm models that used environmental relationship matrix (ERM) kernels, defined by many ECs, to efficiently model similarity among environments so that GEI was modelled as a single term by multiplication of genomic and environ- mental relationship matrices. These reaction norm models generally assumed equal weighting among all ECs, but using partial least squared regression as a preliminary step to learn environment-phenotype association was shown to increase predictive ability (Costa-Neto et al. 2023), while Rincent et al. (2019) developed methods to determine the optimum subset of ECs needed to calculate the ERM for main effects and GEI terms separately. Factor analytic models approximate common environ- ment-specific genetic effects to small numbers of latent environmental factors and allow efficient analysis of highly unbalanced or sparse METs (Burgueño et al. 2008; Piepho and Williams 2024; Smith et al. 2001; Smith and Cullis 2018). Latent factor loadings from factor analytic models can then be used to define interaction classes of environments (iClasses) with minimal genotype cross-over interactions (Smith et al. 2021). In comparison to traditional methods of visualising biplots of single vector decompositions of two- way tables of genotype by environment means (Kempton 1984), and assigning ‘which one where’ mega-environment groups of environments with the same single best genotype (Yan et al. 2000), the iClass method of environment clus- tering is advantages in that it is flexible to highly unbal- anced MET datasets, the factor loadings are identified in a single-stage model rather than genotype means, combines more than two principal component loadings and separates out common and specific GE effects. The genotype score response slopes against each latent factor loading can then be used to identify adapted and stable genotypes (Smith et al. 2021). Smith et al. (2023) further extended this approach to separately model additive and non-additive genetic variance for application in predictive breeding for untested chickpea genotypes. Tolhurst et al. (2022) integrated ECs into fac- tor analytic models using cotton breeding trial data, while other authors have recently investigated correlations between ECs and factor analytic derived latent factors to characterise environmental drivers of GEI and define parameters of target environment types in several crop and animal breeding spe- cies (Bakare et al. 2022; Callister et al. 2024; Fairlie et al. 2024; Oliveira et al. 2020; Rogers and Holland 2022; Sae- Lim et al. 2014). Two-stage multivariate prediction of latent factor loadings using ECs in partial least squares regression models was shown to improve predictive ability of tested rice and soybean genotypes in new environments (Araújo et al. 2024). Methods and models to predict GEI and characterise envi- ronment types using ECs have proven useful for predicting into new environments. However, this assumes that ECs for future environments are already known, so application of these approaches is only possible for expected static GEI effects, which are consistent over years at particular loca- tions, rather than due to non-static year-to-year GEI variabil- ity (Cullis et al. 2000). For breeders to develop genotypes that are specifically adapted to particular mega-environment regions that exhibit repeatable GEI (Yan et al. 2023), and that are broadly adaptable to year-to-year variability within mega-environments, large numbers of environments rep- resenting the target population of environments across space and time should be considered. Here, we use a large MET dataset of CAIGE bread wheat genotypes tested over 13 years in Australian growing environments to: (i) define latent environmental factors and iClass environment types; (ii) demonstrate that latent environmental effects can be related to, and predicted from, observed ECs; (iii) predict out Theoretical and Applied Genetics (2025) 138:241 Page 3 of 24  241 to a large number of extended untested environments cover- ing the Australian grain belt; (iv) evaluate if this framework enables more accurate genomic prediction of new genotypes in new year environments; and (v) evaluate the adaptation patterns of new breeding material introduced through the CAIGE programme. Methods Multi‑environment trial data A bread wheat MET dataset as part of the CAIGE project (Trethowan et al. 2024) included a total 34,293 trial plot observations of wheat grain yield for 3,355 genotypes (vari- eties, cultivars and breeding lines) tested across 114 trial environments (location by year combinations). This broadly covered the major Australian wheat growing regions. For the purposes of this study, we consider Basford and Cooper (1998) descriptions for trial environments in Western Aus- tralia (WA) as the western growing region, environments in South Australia (SA), Victoria (VIC) and southern New South Wales (NSW) as the southern growing region, and environments in northern NSW and Queensland (QLD) as the northern growing region (Fig. 1). Trial years spanned 2011 to 2023 with between three and 14 trials per year, and between 147 and 697 genotypes tested each year (Sup- plementary Table 1). Grain yield was the primary trait of interest but measures of days to heading and plant height were also obtained for a subset of key testing locations each year. Over 70% of trials were sown in May but sowing dates ranged from between the 17th of April to the 17th of August. Genetic material A total of 1,691 and 1,546 genotypes from the CIMMYT and ICARDA breeding programmes, respectively, were tested across all years. Different cohorts of breeding lines from each of CIMMYT and ICARDA were tested each year (Supplementary Table 1). Genotypes were partially repli- cated across trial locations each year and each trial environ- ment included between 114 and 475 entries. Very few breed- ing lines were tested in more than one trial year, but a set of 22 commercial varieties were tested across multiple years as check genotypes (Supplementary Table 2). On average, 188.5 (min = 12; max = 473) genotypes were in common between trials in the same year, but only two to 21 check entries were in common between trials in different years. The trial at Spring Ridge in 2020 was an exception to this where 227 entries from the 2019 trials cohort were included. An efficient sparse MET design with partial replication (Cul- lis et al. 2018, 2006) was used for each trial year. Across all trials, a median of 70.3% (min = 0%; max = 93.4%) of entries was single replicates while a median of 28.5% (min = 0%; max = 99.4%) of entries was replicated twice in each trial. The maximum number of replicates for any line in any single trial was between two (32.5% of trials) and eight (only at Narrabri in 2011). Genotypic analysis Genotypic data were available for 2,058 (61.3%) of the genotypes tested in the MET. These included data from two genotyping platforms. Firstly, 6,530 markers for 1,432 CIMMYT breeding lines were derived from genotype by sequencing data, using methods described by Poland et al. (2012). Secondly, 20,512 single-nucleotide polymorphism (SNP) markers for 1,042 CIMMYT and ICARDA breeding lines and control varieties were genotyped using the Infin- ium Wheat Barley 40 K SNP array (Keeble-Gagnère et al. 2021). There were 416 lines genotyped with both marker platforms, and missing marker data were imputed based on linkage disequilibrium for all genotypes using a custom algorithm (R code for which is provided in Fradgley 2025b). Briefly, for each marker a subset of up to 500 markers with the greatest correlation to the focal marker and with differ- ences in minor allele frequency < 0.1 were used to predict the categorical allele calls of the focal marker using random forest models via the ranger version 0.16.0 package (Wright and Ziegler 2017) in R version 4.4.0 (R Core Team 2023). This method did not require genetic or physical map marker positions. Only markers that could be predicted from a ran- dom forest model having out-of-bag prediction error < 0.2 were imputed. Accuracy of marker imputation was deter- mined by cross-validation within the subset of lines geno- typed on both platforms. 95.2% of marker calls were cor- rectly imputed when randomly masking 10% of marker calls. Cross platform imputation accuracy of randomly masking all marker data from one platform for 10% of lines demon- strated average per marker imputation accuracy of 85.4% and 81.1% for markers on the Infinium Wheat Barley 40 K SNP and genotype by sequencing platforms, respectively. After imputation, markers with > 5% missing data, < 0.025 minor allele frequency or heterozygous call frequency > 0.1 were removed. The combined markers were also pruned to remove redundancy and bias between closely linked markers using a custom algorithm (R code for which is provided in  Fradg- ley  (2025b). Each marker was sequentially evaluated, and all other linked markers (defined as having correlation r < − 0.8 or r > 0.8) were removed. This resulted in m = 7,381 imputed and pruned markers for n = 2,058 genotypes. The n × n genomic relationship matrix (GRM) between all pairs of genotypes was then calculated following Endelman and Jannink (2012) using the A.mat function using the rrBLUP package in R as GRM = MM′ c  , where Mik = Xik + 1 − 2ps , Theoretical and Applied Genetics (2025) 138:241241  Page 4 of 24 X is the n × m marker matrix, ps is the allele frequency at marker s, and c = 2 ∑ sps(1 − ps) . A principal coordinate analysis plot of the GRM is shown in Supplementary Fig. 1. Pedigree and selection history information associated with the lines and check varieties was also used to calculate the coefficient of parentage (COP) matrix among all pairs of 3,405 lines. Pedigrees formatted as text strings in Purdy format (Purdy et al. 1968) were converted to a three-column format with ID, P1 (parent 1) and P2 (parent 2) as column names using a custom script (R code for which is provided in Fradgley 2025b). Seven generations of inbreeding before crossing were assumed for each individual genotype except where a back-cross was specified by ‘*’ in the Purdy nota- tion, where it was assumed that the back-cross was made to an F1 generation plant without inbreeding. Where avail- able, selection histories for breeding lines were also used to determine partially shared inbreeding of lines within fami- lies from the same cross. The full pedigree in three-column format including intermediate generations of inbreeding was then used to calculate the COP matrix using the kin- ship function from the kinship2 version 1.9.6.1 package in R (Sinnwell et al. 2014). Environmental characterisation Environmental data were collected from multiple sources to characteris Fe each trial environment. Daily weather records were accessed from the ‘SILO Patched Point Dataset—Aus- tralian climate data from 1889 to yesterday’ public database as described by Jeffrey et al. (2001). Weather variables included total daily rainfall (mm), daily minimum tempera- ture ( Tmin ; °C), daily maximum temperature ( Tmax ; °C), total daily solar radiation (SR; MJ m−1) and daily vapour pressure deficit (VPD; hPa). Daily photothermal quotient ( PQ ; MJ m−2 day−1 °C−1) was also calculated as the ratio of photo- synthetically active radiation to mean temperature as described by Dreccer et al. (2018), including the conversion of solar radiation to photosynthetically active radiation as applied by Pinker and Laszlo (1992), using the formula PQ = SR×0.47 (Tmin+Tmax)∕2  . The day lengths (h) were derived using each trial location latitude and the daylength function as part of the ChillR version 0.75 package in R (Luedeling et al. 2023). As described by Heslot et al. (2014), daily weather data were processed to define biologically meaningful environ- mental stress covariates within estimated crop development stages. Similarly to Cane et al. (2013), key crop development stages were estimated based on cumulative crop thermal time ( TT  ) growing degree days calculated for each day from trial sowing date with base, optimum and maximum cardinal temperature thresholds of 0, 26 and 34 °C, respectively, as outlined in the Agricultural Production Systems sIMulator (APSIM) wheat crop growth model (McCown et al. 1996). Where dates of sowing were not available, the median value across all known trial environments was assumed. The daily mean temperature ( T  ) was calculated as T = Tmax+Tmin 2  , and thermal time ( TT  ) for each day was then calculated as: A modified version of crop growth stage estimates described in APSIM (McCown et al. 1996) was used. The crop emergence growth stage (Emer), corresponding to growth stage (GS) 10 (Zadoks et al. 1974) was estimated to occur 14 days after sowing if more than three mm of rain occurred in the seven days prior to sowing date or 14 days after the first day of rain after sowing (Sow). The end of juvenile growth stage (Juv; GS30) was estimated to occur when 500 cumulative TT after Emer was reached. Based on a subset of trial environments in which days to heading (He; GS55) was recorded, the cumulative TT after Sow that median ear emergence across all plots occurred was esti- mated as 1,374, 1,301 and 1,194 TT for environments in the North, South and West, respectively. Flowering (Flw; GS65) was estimated to occur when 250 cumulative TT after He was reached. Start of grain filling (Sgf; GS71) was estimated to occur when 250 cumulative TT after Flw was reached. End of grain filling (Egf; GS87) was estimated to occur when 250 cumulative TT after Sgf was reached. Finally, maturity (Mat; GS92) was estimated to occur when 400 cumulative TT after Egf was reached. Thus, dates of seven crop growth intervals including Sow2Emer, Emer2Juv, Juv2He, He2Flw, Flw2Sgf, Sgf2Egf and Egf2mat were defined for each trial environment and a summary of all ECs is provided in Sup- plementary Table 4. Growth interval specific environmental covariates (ECs) were then calculated from daily weather conditions between consecutive estimated crop growth stages. These included; the average values of T, Tmin and Tmax (Avtemp, Avmintemp and Avmaxtemp, respectively), day lengths (AveDL), solar radiation (AveSR) and vapour pressure deficit (AveVPD), in addition to the total number of days (Ndays), the total rain- fall (TotRain), the number of dry days with zero rain (Ndd), the number of frost days when Tmin < 0 °C (Ndays < 0), the number of warm days when Tmax > 26◦C (Ndays > 26) and the number of hot days when Tmax > 34◦C (Ndays > 34) for each interval between consecutive growth stages. For example, Avtemp_He2Flw denotes the average of daily mean temperature ( T  ) between the heading and flowering crop growth stages. Additionally, ECs over multiple intervals were calculated which included the total rain between Sow and Flw (TotRain_Sow2Flw), total rain between Flw and Egf TT = ⎧ ⎪⎨⎪⎩ T 0 < T ≤ 26 26 8 (34 − T) 26 < T ≤ 34 0 T ≤ 0orT > 34 Theoretical and Applied Genetics (2025) 138:241 Page 5 of 24  241 (TotRain_Flw2Egf), prior stored soil moisture estimated as the total rain between the 1st of January and the date of sow- ing (TotRain_priorSow), the number of days between Sow and Flw (Ndays_Sow2Flw) and the number of days between Flw and Egf (Ndays_Flw2Egf). The number of frost days when Tmin < 0 °C within seven days of the estimated flower- ing day (Mintemp < 0_Flw) was also calculated. ECs also included soil characteristics for each trial loca- tion derived from the Soil and Landscapes Grid of Australia dataset (Grundy et al. 2015) using the SLGACloud version 1.0.1 package in R. Point data were extracted for the 90 m2 grid closest to latitude and longitude coordinates of each trial location. Soil attributes were taken from the national soil attributes map at 0–5, 5–15, 15–30, 30–60, 60–100 and 100–200 cm soil depth layers. These included: soil texture percentages of clay, silt and sand; available water capacity (%); bulk density of whole earth including coarse fragments (g cm−1); effective cation exchange capacity (meq 100 g−1); organic carbon (%); mineral associated, particulate and pyrogenic organic carbon (%); coarse fragments (%); pH of water; fraction of both total nitrogen and phosphorus (%), available phosphorus (mg kg−1); as well as the volumetric drained upper and lower limits (%). Depth of soil profile A and B horizons (m) was also extracted. ECs with zero variance across trial environments were removed leaving a total of 96 weather stress ECs and 90 soil characteristic ECs. Finally, the total of w = 186 weather and soil ECs were combined into an p × w matrix (W). We also derived ECs in the same manner as described above for 369,161 extended environments in the Australian wheat belt, over 13 growing seasons from 2011 to 2023, using historical gridded weather data. The mean sowing date across all environments in the observed MET of the 24th of May was assumed for all extended environments. Locations were selected from a grid of 572,721 latitude and longi- tude coordinates covering Australia with 681 latitude values between -44 and -10, and 841 longitude values between 112 and 154. Of these, 28,397 locations were selected which had a fraction of wheat cropping area > 0.02 according to data from Monfreda et al. (2008) using the geodata version 0.6–2 package in R. Functions for deriving ECs are provided in the EC4MET R package (Fradgley and Whan 2025). Multi‑environment trial analyses Data from the multi-environment trials (MET) were ana- lysed with linear mixed models fitted using the ASREML-R version 4.2.0.276 software package (Butler et al. 2017) in R. One-stage MET analyses including all data from individual trial plots across all environments with specific random and residual spatial covariance structures for each trial environ- ment were considered superior in terms of statistical effi- ciency compared to two-stage models (Damesa et al. 2017; Gogel et al. 2018; Welham et al. 2010). To enable fitting of a one-stage model, we first progressively determined parsimo- nious model structures for each trial environment separately using single-trial and single-year models. Fitted parameters from these simpler models were used as average informa- tion algorithm parameter starting values for the one-stage MET model as detailed below. R code for MET analyses are provided in Fradgley (2025b). Single‑trial analyses Following Smith et  al. (2021), linear mixed mod- els with genetic and non-genetic effects were fitted as y = X� + Zgug + Zpup + � , where � is a vector of fixed effects with the associated design matrix X; ug is the vector of random genetic effects associated with the design matrix Zg ; up is the vector of random non-genetic peripheral effects associated with the design matrix Zp and � is the vector of residual errors. A total of 28 models were fitted for each trial, with genotype identity as a random factor effect ( ug ) and testing all combinations of trial row and column coordi- nates factor peripheral random effects, in addition to row and column independent (id), auto-regressive order 1 (ar1) and auto-regressive order 2 (ar2) residual structure non-genetic effects ( up ) (Supplementary Table S3). For each trial, the model with the lowest Akaike information criterion (AIC) was selected as the best fitting model and the estimated ran- dom and residual effects parameters for these were recorded and recycled for use in multi-environment models. Model residually were visually assessed for normality and to iden- tify outlier plot values. Prediction error variance (PEV) of each estimated genetic effect per line was calculated as PEV = SE2 , where SE denotes the standard error. Reliability (r) of each estimate was calculated as: r = 1 − PEV∕�2 G  , where �2 G is the genetic variance. Seven environments with low mean reliability of all genotype estimates per environment (r < 0.3) were removed from further multi-environment analyses. Single‑year MET analyses To determine suitable starting values for multi-year FA1 models, linear mixed models for all trials within each year were then fitted as detailed by Smith et al. (2023). Trial environment identity was included as a fixed main effect ( � ) and genetic random effects using either the GRM or COP to model genetic covariances among genotypes ( Gn ). Consid- ering p as the number of environments and n the number of individuals with genotype data that are present in the GRM or that have pedigree data and are included in the COP, the total genetic variance G is defined as G = Ge ⊗ Gn , where Ge is the p × p matrix of genetic covariance between environ- ments, Gn is the n × n GRM or COP numerator matrix and Theoretical and Applied Genetics (2025) 138:241241  Page 6 of 24 ⊗ indicates the Kronecker product. The Ge matrix was first structured for a baseline model having diagonal (DIAG) var- iance–covariance structure that assumes a specific genetic variance ( �2 ) at each of e environments but with independent genetic covariance (0) among all pairs of p environments: Then, reduced rank factor analytic order 1 (FA1) models that approximated heterogeneous genetic covariance among environments in Ge using a single factor (k = 1) were fitted for all trials within each year. Reduced rank factor analytic (FAk) models combine both reduced rank and diagonal com- ponents as Ge = �D� � +� , where � is a matrix of p × k factor loadings, D is a k × k matrix of factor score variances and � is the �2 1 , �2 2 ,… , �2 e diagonal (non-identity) matrix of environment-specific genetic variances. Therefore, �D� ′ accounts for common and repeatable genotype nested with environment effects (GE) occurring in more than one envi- ronment which can be approximated by the FA structure, while � accounts for environment-specific GE effects that are not repeatable across tested trial environments. Rather than fitting both additive and non-additive effects with sepa- rate FA structures as implemented by Smith et al. (2023) and Smith and Cullis (2018), we only fitted additive genetic effects considering the GRM or COP similarly to Tolhurst et al. (2022). To aid convergence of the mixed model average informa- tion (AI) algorithm, REML starting values for trial specific random and residual covariance structure row and column effects ( up ) for single-year DIAG models were recycled from the single-trial analyses best fitting models. Starting values for single-year FA1 models specific genetic variances ( � ), as well as non-genetic random effects variance parameters ( up ) were recycled from fitted parameters from converged single-year DIAG models. Auto-regressive row and column covariance effects for environments that prevented single- year models from converging were adjusted to independent ‘id(ROW):id(COL)’ for further analyses. Full multi‑year MET analyses Full one-stage models were finally fitted using all plot data across all environments with high mean reliability (r > 0.3) from all trial years including either the GRM or COP as Gn . Note that at this stage, the dimension size (n) of the Gn matri- ces were n = 2,058 and 3,405 for GRM and COP models, respectively. Like the single-year analyses described above, the one-stage DIAG model was first fitted using � fixed envi- ronmental main effects and up parameters starting values Ge = ⎡ ⎢⎢⎢⎣ �2 1 0 ⋯ 0 0 �2 2 ⋯ 0 ⋮ ⋮ ⋱ ⋮ 0 0 ⋯ �2 p ⎤ ⎥⎥⎥⎦ recycled from single-year DIAG models. Then reduced rank FA models with incrementally increased k factors were progressively fitted up to k = 3 when, as suggested by (Smith et al. 2015), more than 80% of the total variance was explained. For the one-stage FA1 model, � and up param- eter starting values were recycled from the one-stage DIAG model and ug parameter starting values including � and the first factor loadings were recycled from respective single- year FA1 models. All possible � , ug (including k > 1 factor loadings) and up starting values for subsequent FA2 and FA3 models were recycled from lower order (k-1) models. As described by Smith et al. (2021), the estimated � matrix of factor loadings and associated genotype scores were rotated to a principal component solution using the svd function in R. Ge covariance matrices estimated from rotated � could then be converted to genetic correlation matrices among all pairs of environments using the cov2cor func- tion in R. The percentages of variance explained by each r factor loading for each of p environments ( vre ) was calculated as vre = 100(Dkkdiag(�k� � k ))∕diag(Ge) , where Dkk is the value of factor score variance for factor k, �k is the vector of envi- ronment factor loadings for factor k and diag(Ge ) is the diag- onal vector component of the Ge matrix. The total variance explained by the full model including all factors ( v ) was calculated as v = 100tr(�D� � )∕tr(Ge) , where tr(Ge) com- putes the trace, i.e. the sum of the diagonal component of Ge. Each environment was assigned to one of several defined interaction classes (iClasses) environment types. These iClasses minimise cross-over interactions from common GE effects within classes (considering that changes in rank ordering of genotype scores for each factor only occur when each factor loading equals zero) (Burgueño et al. 2008; Smith et al. 2021, 2023). This method uses the vector of k rotated factor loadings for each of e environments ( �e ) where positive loadings are assigned as ‘p’ and negative loadings are assigned to ‘n’. These vectors of characters are then concatenated to iClass strings of length k with 2k pos- sible iClasses. As proposed by Smith et al. (2023), iClass groups containing only one environment were merged with the iClass with the highest mean genetic correlation. Envi- ronment iClasses were then compared across geographic regions and years to illustrate spatial and temporal drivers of GEI. Genotype scores ( ̃f  ) from FA models describe the regres- sion slopes responses of each genotype to each factor load- ing and are used to estimate environment-specific empirical best linear unbiased predictions (EBLUPs). As detailed by Smith and Cullis (2018), general overall performance (OP) across all iClasses for genotype i was defined as OPi = Λ1 f̃1i , where Λ1 is the average first factor loading across environ- ments, and f̃1i is the first factor genotype score for genotype Theoretical and Applied Genetics (2025) 138:241 Page 7 of 24  241 i. Note that the OP score will not apply a minority of envi- ronments with negative first factor loadings. Similarly, the iClass overall performance (iClassOP) for each genotype i in each iClass � was calculated as: where Λr� is the mean factor loadings for factor r of all environments in iClass � , and f̃ri is the estimated genotype scores for genotype i and factor r (Smith et al. 2021; Smith and Cullis 2018). A genotype stability index based on the root mean squared error (RMSD) of higher order factor scores was also then calculated as: where ∈̃ij is the deviations of the linear regression line of common GE effects against first factor loadings for genotype i in environment j. Environment-specific EBLUPs integrat- ing common GE effects could be calculated for any new environment where factor loadings ( � ) can be predicted. Predictability of environment main effects and factor loadings To increase interpretability and predictability of iClasses, we then investigated associations between both estimated environmental main effects and latent environmental fac- tors against observed crop stress ECs. Pearson’s correlation coefficients were calculated between all combinations of w ECs and each of �1∶k factor loadings vectors in addition to estimated environment main effects. One thousand per- mutations of randomly sampling each FA loadings vector were used to define 95% confidence interval thresholds for associations with ECs. We then tested several approaches for predicting �1∶k vec- tors of factor loadings and environmental main effects in a second-stage prediction model from w EC predictors. Note that this approach used fixed rather than random effect esti- mates of environmental main effects from the linear mixed model because they were used again to fit the second-stage model from which predictions were made to new environ- ments. With more predictor ECs than environments (w > p), we tested multivariate prediction models including: random forest (RF) regression using ensembles of decision trees (Breiman 2001); a regularisation approach including the least absolute shrinkage and selection operator (LASSO; Tibshirani 1996); and partial least squares regression (PLS; Helland 1990). RF models were fitted using the randomFor- est version 4.7–1.1 package in R (Liaw and Weiner 2022) OP�i = k∑ r=1 Λr� f̃ri RMSD = √√√√1 p p∑ j=1 ∈̃ 2 ij including default parameters of 500 trees, one-third of the variables randomly sampled at each split, and a maximum node size of five. Feature importance scores from full fitted models were also used to identify highly influential ECs. LASSO models were fitted using the glmnet version 4.1–8 package in R (Friedman et al. 2010) where optimum values of lambda for each model were selected based on tenfold cross-validation for 100 lambda values between 1 × 10–8 and 0.04. Predictive ability of these models was assessed for a leave-one-environment-out cross-validation scheme as the Pearson’s correlation coefficients and the root mean squared errors (RMSE), calculated as RMSE = �∑ (yobs−ypred) 2 p  , where yobs and ypred are the observed and predicted values of each FA loadings and environmental main fixed effects, respectively. To account for shrinkage in predictions of envi- ronment main effects, out-of-fold predictions were adjusted so that the mean and variance of predictions equalled those o f o b s e r v e d v a l u e s a s y��� = y−y√∑ (y2)∕e−1 × √∑ (y��� 2)∕e − 1 + yobs , where y is the vector of predictions, y is the mean of the predictions, yobs is the vector of observed values, yobs is the mean of the observed values and e is the number of environments. After predictive ability of models for environmental main effects and factor loadings influencing GE was validated, we then made predictions into the larger set of 369,161 extended environments as described above. To ensure realistic vari- ance for the yield predictions, bias and shrinkage in pre- dictions of environmental main effects were adjusted using intercept and slope linear regression coefficients from linear models fitted between observed and out-of-fold predicted main effect values as described above. Genomic prediction in untested environments To validate the approach of prediction of factor analytic mixed model environmental main effects and factor load- ings for untested environments, forward cross-validation of genotype performance prediction was applied. The cross-validation involved four rounds of forward predic- tion including trial environments between 2020 to 2023, 2021 to 2023, 2022 to 2023 and only 2023 as test sets of environments and using all environments from pre- ceding years as training datasets to fit models. First, a genetic main effects only model (MM) was fitted as a baseline model which included only a single vector of ug genetic effect (not with covariance structures nested within environments) so that genetic effects were the same across all environments. Effects for peripheral ( up ) effects were fitted with starting values from single-year models as described above. More complex models were Theoretical and Applied Genetics (2025) 138:241241  Page 8 of 24 then sequentially fitted from a DIAGGRM model to an FA3GRM model including lines present in the test years in the GRM, as detailed above. Genetic effects for each line from DIAG models were averaged across all environ- ments in the training set of environments, and these aver- aged effects were used for predictions in all test environ- ments. For all FA models, as also detailed above, random forest models were trained on the estimated factor load- ings and environmental main effects from models fitted on the training set of environments with observed ECs as predictors. Factor loadings and environmental main effects for test set year environments were then predicted, and EGBLUPs for all genotypes in all test environments were calculated. Therefore, MM and DIAG models made no predictions of GE in the test environments, whereas FA models did predict environment-specific genotype performance. Within-environment predictive ability of tested and untested genotypes for each test environment was cal- culated as the Pearson’s correlation and RMSE between Table 1   Summary statistics of all fitted one-stage models for multi-year MET analyses AIC is the Akaike Index Criterion, BIC is the Bayesian Index Criterion, LL is log likelihood and % var is the percentage variance explained Model Parameters Coefficients AIC BIC LL % var DIAGCOP 522 364,436  − 15,704.57  − 11,321.47 8,374.29 – FA1COP 742 729,515  − 20,244.21  − 14,937.47 10,754.11 53.53 FA2COP 852 732,804  − 21,285.98  − 15,062.99 11,383.49 71.68 FA3COP 962 736,093  − 21,641.58  − 14,512.74 11,669.79 76.54 DIAGGRM 522 227,926  − 12,609.19  − 8,226.09 6,826.59 – FA1GRM 742 455,254  − 14,805.28  − 9,498.53 8,034.64 65.39 FA2GRM 852 457,302  − 15,117.69  − 8,895.70 8,299.84 76.58 FA3GRM 962 459,350  − 15,531.96  − 8,403.12 8,614.98 86.94 Fig. 1   Map of CAIGE bread wheat multi-environment trial locations across three major wheat growing regions and the main Australian wheat growing areas in 2020 according to data from Monfreda et al. (2008) Theoretical and Applied Genetics (2025) 138:241 Page 9 of 24  241 predicted values and observed genetic effects estimated from single-trial analyses. Statistically significant differ- ences between mean values of within trial environment predictive ability for the DIAG and each FA models were compared to the baseline MM model using paired two- tailed t-tests with the t.test function in R. Results We analysed a large wheat MET dataset including 3,354 genotypes and 114 field trial environments to character- ise genotype by environment interactions and associated latent environmental effects. Data from weather and soil environmental covariates could then be related to latent environmental effects to characterise GEI-based environ- ment types (iClasses) for Australia. Model comparisons By sequentially fitting single-trial, then single-year analysis models, a full one-stage model incorporating 32,472 yield plot observations could be fitted. Linear mixed models included reduced rank factor analytic (FA) approximation of heterogeneous genetic covariance structure among envi- ronments up to FA order k = 3 with either pedigree (COP) or genomic (GRM)-based relationship matrices. The COP models estimated more coefficients had lower values for AIC and BIC and greater log likelihood values relative to the GRM models (Table 1). This reflects the greater number of lines with pedigree information compared to those with genetic data and therefore greater capacity to fit larger COP models. The lowest AIC and BIC values across FA model orders differed for both COP and GRM models so similarly to Smith et al. (2015) and Tolhurst et al. (2022), we took a pragmatic approach considering the percentage of vari- ance explained per environments. GRM models explained a greater proportion of the phenotypic variance than COP models. On average across all environments, 84.6% of the variance was accounted for per environment and more than two-thirds (69.0%) of environments had greater than 80% variance accounted for (Table 1). FA3GRM was the only model to explain more than 80% of the total variance and so was chosen as the best fitting model and used to further define clusters of environments with similar GE. The vari- ance in estimated coefficients of common GE effects that were approximated by the reduced rank FA structure was 0.148 and much larger than that of environment-specific genetic effects (0.007) which models non-repeatable GE. In fact, over half (51.8%) of the environments had near-zero (< 0.0001) specific genetic variance. Interaction classes in observed multi‑environment trials Estimates of environmental main effects in yield varied widely from 0.14 t ha−1 at Junee in 2019 to 6.97 t ha−1 at Balaklava in 2016. As expected, almost all (100 of 110) environments had positive loadings on the first factor which explained 50.0% of the total variance. This indicated that cross-over GEI was rare over the first factor while it mainly accounted for heterogeneous scaling of within-environment genetic variances across environments; higher yielding envi- ronments generally had greater genetic variance. Indeed, the correlation coefficient between the first factor loading and the estimated environmental main effects was 0.42. The sec- ond factor accounted for 26.3% of the GE variance with more balanced positive and negative loadings values (83 and 27 environments, respectively), representing the most important latent environmental effects influencing common cross-over GE. The second factor loadings were completely uncorrelated with environmental main effects (r = 0.01) sug- gesting that these environmental effects on yield are highly dependent on genotype with minimal overall yield impacts. Finally, the third factor accounted for 10.6% of the GE vari- ance had 57 and 53 environments with positive and negative factor loadings, respectively, and was more closely corre- lated with environment main effects (r = − 0.40). We used the iClass environment type clustering system to define groups of environments with the same sign (‘p’ = pos- itive, ‘n’ = negative) of factor loadings, and minimal cross- over interactions (similar genotype ranking) within classes. The npp iClass was rarely identified, including only three environments whereas ppp was the most frequent, including 39 environments. Genetic correlations among all environ- ments, estimated from the factor loadings and environment- specific genetic variances, ranged from − 0.90 to 0.99 with a mean of 0.40 across all pairs of environments (Fig. 2a). The two rarest iClasses (npp and npn) had small negative first factor loadings but still had relatively high genetic correla- tions with equivalent iClasses with positive first factor load- ings (ppp and ppn, respectively). Genetic correlations among environments within iClasses groups were greater, ranging from 0.59 to 0.97 with a mean of 0.70, whereas genetic cor- relations among environments in different iClasses ranged from − 0.45 to 0.77 with a mean of only 0.17 across all pairs of iClasses (Fig. 2b). These results highlight the large GEI effects among the broad range of Australian wheat growing environment types. Of note, the ppn iClass was the second most frequently characterised, in 34 environments. The ppn iClass exhibited high within-class genetic correlations (rang- ing from 0.35 to 0.999 with a mean of 0.67), and so is an example of a frequently occurring class of environments with relatively consistent genetic effects and low GEI. Theoretical and Applied Genetics (2025) 138:241241  Page 10 of 24 Comparison of iClasses in thirteen trial years across major growing regions revealed that iClasses were not specific to either regions, locations or years (Fig. 3). Nev- ertheless, some trends across regions within years were Fig. 2   Heatmaps of genetic correlations among all pairs of environments ordered by iClasses (a); and mean values of genetic correlations among pairs of environments within (diagonal) and between (off-diagonal) the six iClass groups (b) Fig. 3   Distributions of yield main effects across iClasses (a); interaction classes (iClasses) are defined from positive and negative factor loadings (b–d); patterns of estimated iClass identity for 110 trial environ- ments across regions, locations and growing years (e). The size of the points in plots b–d and within each cell in plot e indicates the percentage of the genetic variance explained by the FA3GRM model for each trial environment Theoretical and Applied Genetics (2025) 138:241 Page 11 of 24  241 apparent. For example, in 2019, all environments in the North and South regions were assigned the ppp iClass, except for Balaklava where only 6.7% of the variance was explained, and three of the four environments in the West were also assigned the ppp iClass. The pnn iClass was also never assigned in the West in any year. Single locations rarely had consistent iClasses even when tested over many years. Junee in the North, which was included in all thirteen years, was allocated to every iClass except npp in at least one year. Narrabri in the North was described as iClass ppn for six of the thirteen years. In contrast, Roseworthy in the South was assigned to iClass ppn only once in 2017 and instead was iClass ppp for eight of the thirteen years. These patterns are illustrated in (Fig. 3). Environmental covariates characterise environment interaction classes We found environmental main effects and factor loadings defined by FA models could be related to observed envi- ronmental covariates (ECs) to define generalised environ- ment types of each iClass. Key associated ECs were iden- tified by comparing correlations between ECs and both environmental main effects and factor loadings estimated from the FA3GRM model in addition to importance scores derived from fitted random forest multivariate predic- tion models. Environmental main effects on yield were positively correlated with rainfall throughout the season where total rainfall between sowing and flowering, juve- nile to heading, heading to flowering, flowering to end of grain fill, and start of grain fill to end of grain fill were all positively correlated (0.33 < r < 0.40) and were among the most important explanatory variables in random for- est models (Fig. 4a). Some soil characteristics were also positively associated with environmental main effects on yield, including organic carbon at the 0.6 to 1 m soil depth and available phosphorus at the 0.15 to 3 m soil depth. The six environments that received more than 200 mm of rainfall between sowing and flowering and had more than 0.5% organic carbon at the 0.6 to 1 m soil depth all had high yield main effects between 5.55 and 6.94 t ha−1, while the eight environments that received less than 100 mm of rainfall between sowing and flowering and had less than 0.4% organic carbon at the 0.6 to 1 m soil depth had much lower yield main effects between 0.61 and 4.07 t ha−1. The 29 ECs with the strongest positive or negative cor- relation with the first factor loading (0.25 < r < − 0.25) were soil property covariates from environments with strong posi- tive first factor loadings (Fig. 4b). For example, Kalkee and Breeza (2023) and Spring Ridge (2020) soil was described by high drained upper limit, available water capacity and cation exchange capacity, with heavier soil types and a high percentage of clay or silt and low percentage of sand. The second factor loading, which explained the largest proportion of cross-over GEI effects but was unrelated to environmental main effects, was more clearly associated with temperature and the duration of growth stage related ECs. Key ECs with both strong correlations and high impor- tance scores included the average daily mean temperature and average daily maximum temperature between the head- ing and flowering growth stages (Fig. 4c). Environments with strong positive second factor loadings, such as Rose- worthy in 2016 and 2019 as well as Balaklava in 2016, had generally cooler average temperatures (10.0 to 11.7 °C) and a long interval (21 to 25 days) between heading and flower- ing. In contrast, environments with strong negative second factor loadings, such as Breeza and Narrabri in 2023 and Bellata in 2022, were located in the Northern region, with greater average temperatures between heading and flowering Fig. 4   Correlation coefficients for each environmental covariate and factor analytic mixed model environmental main effects (a), factor loadings (b–d) Importance scores of each environmental covariate are derived from multivariate random forest prediction models. Grey areas indicate 95% confidence intervals based on one thousand per- mutations of random sampling Theoretical and Applied Genetics (2025) 138:241241  Page 12 of 24 (15.0 to 17.7 °C) and shorter duration between these two growth stages (14 to 15 days). Similarly to environmental main effects on overall yield, the third factor loading was strongly related to ECs quan- tifying rainfall throughout the growing season. Total rain- fall between flowering and the start of grain fill, between flowering and the end of grain fill and between the juvenile to heading stages all correlated negatively with the third factor loading (r < − 0.43) and had high importance scores (Fig. 4d). Environments with a strong negative third factor loading, such as Bellata, Breeza and Junee in 2022, had high rainfall between flowering and the end of grain fill (between 162.5 and 207.2 mm). Similarly, Junee and Horsham in 2016 had high rainfall between the juvenile and heading growth stages (between 174 and 289.7 mm). In contrast, environ- ments with strong positive third factor loadings, such as Roseworthy (2013, 2015, 2020 and 2021), Narrabri (2019) and Spring Ridge (2020), received lower rainfall (between 0.5 and 61 mm) between flowering and the end of grain fill as well as the period between the juvenile and heading growth stages (between 13.4 and 130.7 mm). These envi- ronments could therefore characterised as relatively drought stressed environment types. The environmental characteristics of all factor loadings could then be used to characterise environment types of each iClass group of environments. For example, the pnn iClass demonstrated the highest average yield main effects (5.42 t ha−1) and therefore more favourable environmental conditions. The pnn iClass associated factor loadings were defined by more fertile soils with better water holding capac- ity, which was indicated by the positive first factor loading. In addition, higher temperatures after heading associated with the negative second factor loadings and greater rainfall associated with the negative third factor loadings further described these favourable environments. In contrast, envi- ronments in the npp and ppp iClass had the lowest aver- age yields (2.55 and 2.89 t ha−1, respectively) due to less favourable environmental factors, particularly including lower rainfall and more severe drought stress (indicated by positive third factor loadings). Despite these trends, defini- tion of iClasses based on individual ECs was difficult with considerable overlap in ranges of ECs between iClasses. Prediction of environmental main effects and factor loadings We tested several multivariate models for the prediction of environment main effects on yield and estimated factor load- ings for test environments with known ECs in a leave-one- environment-out cross-validation strategy. Random forest models most accurately predicted environment main effects (r = 0.67) with the lowest error (RMSE = 1.23 t ha−1; Fig. 5). Noticeable shrinkage of predicted values was evident, particularly for random forest models, where the variance of the predictions (0.59) was much smaller than the variance of the observed environment main effects (2.31). Adjusted predictions with equal mean and variance as the observed environment main effects showed slight increase in error but similar differences among models (Fig. 5). For prediction of factor loadings influencing GE, random forest models performed best with the greatest predictive ability and the smallest RMSE. Interestingly, prediction accuracies were greater for higher order factor loadings in more complex FA3 models which explained only a small proportion of GE variation (Fig. 6). The sign (positive or negative) of the three factor loadings was correctly predicted for 90.9, 73.6 and 70.0% of environments, respectively, by a random forest model. Furthermore, iClass identities derived from the predicted factor loadings were correctly predicted for 44.6% of all environments. The two most frequent iClass environment types (ppn and ppp) had relatively high fre- quencies of 30.9 and 35.5%, respectively, while they were more frequently correctly predicted at 70.6 and 56.4% of the time, respectively. Conversely, models were unable to pre- dict the rarer npn and npp iClasses that were never correctly predicted, indicating significant generalisation of this stage in environment type predictions. Importantly, the multivari- ate predictive ability of both environmental main effects and factor loadings were considerably greater than the strongest correlations with any single EC, suggesting that multiple soil and weather ECs defined over multiple crop growth intervals act in combination to influence grain yield and are essential parameters to predict GEI. Definition of Australian wheat mega‑environments After demonstrating that environmental main effects and factor loadings for untested environments can be predicted with reasonable accuracy, we used the full fitted FA3GRM model to predict yield, GE factor loadings and iClass iden- tity for 28,397 locations covering the Australian wheat belt for 13 years from 2011 to 2023. In total, this represented prediction of genotypes into 369,161 new environments. Maps of predicted mean yield each year clearly illustrated that higher rainfall regions closer to the coast in the North- ern and Southern growing regions and the south-west of the Western growing region were associated with higher predicted yields (Fig. 7). Dramatic year-to-year variability was highlighted with the greatest predicted yields occurring in the years where actual rainfall was high (2016, 2022). The lowest prediction for yield was for 2019 and this year coincided with widespread drought in Australia (Wittwer and Waschik 2021). Compared to the inconsistent trends in iClasses across the tested trial environments in the MET dataset (Fig. 3), mapping the predicted iClass identity over the landscape Theoretical and Applied Genetics (2025) 138:241 Page 13 of 24  241 revealed more interpretable trends in iClass environment types within years, and shifts in distributions of iClasses due to year-to-year variability (Fig. 8). Patterns of iClass envi- ronment types were not separable at the scale and boundaries of major growing regions and often changed considerably at small scales within regions or between years. For example, when particularly low yields were predicted across most of the Western Australian region in the low rainfall year of 2019, a greater proportion of ppp iClass environments were predicted, whereas in the higher rainfall years of 2022 (all regions) and in the northern and southern regions (2016), the ppn iClass was more frequently predicted to occur. Despite shifts in predicted iClass identity across regions due to seasonal year-to-year weather variability, consid- ering the patterns of iClass frequency at each location across the period of 13 years revealed spatial distribu- tions of mega-environments where each iClass would be more likely to occur considering static average climate conditions and soil properties (Fig. 9). The ppn iClass was the most frequently predicted across locations and was particularly prevalent closer to the coast in all grow- ing regions that generally experienced higher rainfall. In contrast, the ppp iClass was more prevalent at drier inland locations (in particular, in the wheat growing areas of SA or the north-west of VIC) and was very rare in the northern region in QLD (Fig. 9). Overall, the pnn and pnp iClasses occurred rarely and only in the northern region in the north of NSW and QLD. Some locations, such as the south-west of the wheat Western growing region in WA, were highly likely to be associated with the ppn iClass, whereas locations such as the centre of the wheat belt in NSW or WA had much less certainty in the predicted iClass identity (Fig. 9). Fig. 5   Model comparisons of observed and predicted (a) and adjusted predicted (b) values of environment main effects on yield. RF = ran- dom forest; LASSO = least absolute shrinkage and selection operator; and PLS = partial least squares. RMSE = root mean squared error and r = predictive ability. Predicted values were adjusted to minimise bias and shrinkage of predicted values. Solid diagonal lines indicate x = y, and diagonal dashed lines indicate RMSE thresholds Theoretical and Applied Genetics (2025) 138:241241  Page 14 of 24 Genetic diversity adapted to Australian environments The different sets of genotypes used in this study to char- acterise GEI and iClasses included check varieties that are each well adapted to specific environment types that make up the wide range of Australian environments. Breeding lines from both the CIMMYT and ICARDA international wheat breeding programmes introduced through the CAIGE project were found to be more genetically diverse than the Australian check varieties (lower average genomic kinship estimates within groups) and were similarly distant between groups (Supplementary Fig. 2). The Australian check vari- eties Mace, Scepter, Rockstar, Vixen and Calibre are all highly connected through pedigrees and have closer kinship relationships than 99.89 and 99.78% of pairs of CIMMYT and ICARDA lines, respectively. These genetic resources may provide useful specifically or broadly adapted materials. Differences in genotype scores and predicted performance in each iClass environment illustrated where breeding has targeted specific environment types (Fig. 10). In general, the Australian check varieties had higher positive genotype scores for all three factors for the FA3GRM model (Fig. 10a). This indicates (as expected) high overall performance across iClasses, but particular adaptation to cooler, drier and lower yielding iClasses environment types such as npp and par- ticularly ppp (Fig. 10b), which both frequently occurred throughout all regions (Fig.  9). In contrast, ICARDA breeding lines demonstrated lower third factor loadings on average, and CIMMYT breeding lines demonstrated more negative third factor loadings (Fig. 10a), representing adap- tation to warmer and wetter environments (pnn iClass) that was predicted to occur more often in the wheat growing regions of northern NSW and QLD (Fig. 10b). Although the Australian check varieties on average had the highest overall performance in the ppn iClass, there was a wide range in overall performance of CIMMYT lines for the ppn iClass. Several of these showed the greatest overall performance, highlighting the potential of this material for breeding in Australia (Fig. 10b). Also of note, there was particularly large variance in third factor genotype scores among Australian check varieties (Fig. 10a) and these also often had high RMSD across all iClasses (Fig. 10c). This highlighted cultivars, such as Rockstar and Vixen, that had Fig. 6   Comparison of predictive ability and root mean squared error (RMSE) of three multivariate prediction models for predicting factor ana- lytic model factor loadings. RF = random forest; LASSO = least absolute shrinkage and selection operator; and PLS = partial least squares Theoretical and Applied Genetics (2025) 138:241 Page 15 of 24  241 high overall performing and high RMSD and so were par- ticularly adapted to specific environment types. These two checks were ranked first and second, respectively, for iClass overall performance for both the npp and ppp iClasses but displayed the third and second lowest iClass overall perfor- mance in the pnn iClass (Fig. 10b). Conversely, varieties such as Sunmaster and Suntop demonstrated similar patterns of adaptation relative to the CIMMYT material, performing relatively poorly in npp and ppp iClasses but particularly well in the pnn iClass compared to other checks (Fig. 10b). Of note, the genotypes with the greatest overall performance and low RMSD (Fig. 10c) were often CIMMYT derived CAIGE lines, underlining their strong potential for broad adaptability and yield stability. To further investigate mechanisms of genotype adapta- tion and GEI, we compared how genotype scores for each factor relate to phenology and plant height traits within each group of genotypes. Among all comparisons, only a weak but significant correlation was found between genetic effects for estimated days to heading and the third factor score for CIMMYT lines (r = 0.16, p < 0.001). This suggests that earlier maturity CIMMYT lines demonstrated lower third factor scores and greater overall performance in pnn and ppn iClass environments. Associations were also found for plant height where a small but significant negative correla- tion with the third factor scores was found for CIMMYT lines (r = − 0.15, p < 0.001), but this association was much stronger for the Australian checks (r = − 0.89, p < 0.001). This suggests that shorter varieties and genotypes within this height range from both groups had generally higher third factor scores and greater performance in the more drought stressed pnp and ppp iClass environments. This trend was also reflected by differences in reduced height gene status among lines where semi-dwarf Rht-1 (Rht-B1b and Rht- D1a) genotypes had a significantly (p < 0.001) lower mean third factor genotype score of 1.2 compared to a mean of Fig. 7   Maps illustrating predicted yield main effects across Australian wheat belt locations. Predictions were made from soil and weather envi- ronmental covariates from 2011 to 2023 Theoretical and Applied Genetics (2025) 138:241241  Page 16 of 24 3.0 for Rht-2 (Rht-B1a and Rht-D1b) genotypes and concurs with results from Eagles et al. (2014). Approximately 96 and 95% of lines tested from the CIMMYT and ICARDA breeding programmes, respectively, were semi-dwarf Rht-1. Genomic prediction into untested environments We further validated the above approach of prediction of factor analytic mixed model factor loadings with observed ECs by testing a forward genomic prediction cross-valida- tion scheme into untested environments with known ECs. Predictive ability for check varieties that were tested over several years was significantly increased for higher order FA models compared to the baseline genetic main effects (MM) model in all cross-validation folds. This increase was smallest in the first cross-validation fold, which included 9 years of trial data in the training set, where prediction accuracies increased from 0.32 for the MM model to 0.38 for the FA3 model. When more years of data were used in the training set for the cross-validation fold including only 2023 as the test year, a greater increase in predictive ability was observed from 0.29 for the MM model to 0.53 for the FA3 model. (Fig. 11). The increase in prediction accuracies was much lower for untested genotypes (CAIGE lines tested in single years). In most cases, there was no significant difference in predictive ability between DIAG or FA models compared to the MM baseline model. Nevertheless, prediction accuracies were slightly, but statistically significantly, increased from 0.14 and 0.13 for MM models to 0.16 and 0.15 for FA2 mod- els in the 2020 to 2023 and 2021 to 2023 cross-validation test folds, respectively (Fig. 11). Overall, these results con- firm that prediction of factor loadings with observed ECs facilitates increased prediction of GE effects when genotype scores are estimated from tested genotypes, or from accurate genomic predictions. Fig. 8   Maps illustrating predicted iClass identity across Australian wheat belt locations based on soil and weather environmental covariates from 2011 to 2023 Theoretical and Applied Genetics (2025) 138:241 Page 17 of 24  241 Discussion Wheat breeding in Australia is hindered by prevalent and unpredictable GEI. By making use of a large dataset of MET series and well-established factor analytic mixed model analysis methods, we demonstrated that environmental main effects on yield and latent factor loadings that drive GEI can be related to and are predictable from observed ECs. Taking a similar approach to Araújo et al. (2024) and Costa-Neto et al. (2020), second-stage models to predict these param- eters were used to predict out to a much larger expanded set of environments that cover Australian wheat growing environments at a much broader scale across locations and time. Importantly, this mapping approach revealed general- ised spatial and temporal trends in predicted GEI that were unclear when only considering the observed MET environ- ments. Yield main effects and environment types were both significantly impacted by year-to-year variability which emphasises the importance of multi-year, multi-environment testing to select broadly adaptable genotypes to these non- static environmental effects. Nevertheless, by predicting into a large number of environments that thoroughly represent the Australian wheat TPE over several years and considering the probability of predicted iClass environment types occur- ring at each location over several years, some trends in static GEI effects could be used to define mega-environments of locations for which breeding can target specifically adapted genotypes (Basford and Cooper 1998). Importantly, these results justify considering each year-location environment uniquely independent rather than modelling separate loca- tion and years random effects and demonstrates how con- sidering predictions for a long-term distribution of weather scenarios based on historical observed weather enables dif- ferentiation of static and non-static GEI effects that would Fig. 9   Maps of probability of predicted iClasses occurring between 2011 and 2023 across Australian wheat belt locations based on soil and weather environmental covariates. Upper plots show the pre- dicted distributions of probabilities of each iClass occurring over the multiple years while the lower plot shows the most frequently pre- dicted iClass at each location Theoretical and Applied Genetics (2025) 138:241241  Page 18 of 24 be likely to occur in future years. In practice, selection for high overall performance in one iClass would be the best approach for locations with a high probability of that iClass occurring. However, if multiple iClasses are predicted to occur with reasonable probability over several years, then a more complex approach of selection for high overall perfor- mance proportional to the multiple iClasses in combination with low RMSD stability values would be required. Heslot et al. (2014) considered the analogy of the Anna Karenina effect; that ‘Happy families are all alike; every unhappy family is unhappy in its own way’, for the unpre- dictability and difficulty in generalisation of GEI due the large number of independent environmental effects that may limit yield. We found that common GE effects variance was much larger than that for specific effects. This may be due to the large number of environments and seasons that the MET Fig. 10   Comparisons of genotype adaptation for check varieties and CAIGE lines from the CIMMYT and ICARDA breeding pro- grammes. Upper plots (a) show the genotype scores for each factor loading from the FA3GRM model. Middle plots (b) show overall per- formance scores for each genotype for the average environment of each iClass (iClass OP). The lower plot (c) compares the overall per- formance across all iClasses (OP) over the first factor and yield stabil- ity over higher order factors (RMSD) Theoretical and Applied Genetics (2025) 138:241 Page 19 of 24  241 effectively sampled from the target population of environ- ments (TPE) covering the diverse Australian wheat grain belt. Thus, contrary to the Anna Karenina effect, we found that generalisation of GEI and predictability into untested environments is possible to some extent. Nevertheless, we acknowledge that methods presented here inherently gen- eralise GEI, first through factor analytic approximation of common rather than environment-specific genetic effects considered random rather than fixed, and secondly through predictions in the second-stage model including ECs out to the wider set of environments with another degree of error and uncertainty. As a result, finer scale and rare GEI may well be overlooked. We observed significant shrink- age towards the mean of predictions of environmental main effects. Despite attempts to rescale the variance of predic- tions through cross-validation, predictions may be less accurate in extreme environments. We also note that gen- eralisation of GEI is inherent to the scale of the TPE; in that similar approaches using MET data for a more refined TPE would likely resolve finer scale GEI effects and mega- environments better than presented here. The genotypes used in this study also likely influenced the generalised environ- ment type effects. We found that the CAIGE wheat lines from both CIMMYT and ICARDA represent a more diverse and distant genepool than the Australian check varieties and demonstrated novel adaptation patterns across environment types. Different results, in terms of iClass mega-environment trends and boundaries, would likely be found using a dif- ferent panel of genotypes or lines within specific breeding programmes with different selection histories. We also anec- dotally observed that the Australian checks often displayed contrasting genetic correlations compared to CAIGE lines. Further work may investigate whether fitting different covar- iance structures to groups of related genotypes increases pre- dictive performance. The iClass environment types that more frequently occurred in northern NSW and QLD were identified as tar- get environments for CAIGE material (particularly from CIMMYT) adaptation. These findings agree with previous work that found high genetic correlations between environ- ments in the northern Australian wheat growing region and the CIMMYT selection environments in Obregon, Mexico Fig. 11   Boxplots comparing within trial environment predictive abil- ity of different models over four forward prediction cross-validation schemes. Each point indicates a predictive ability for a single envi- ronment, and large points indicate the mean predictive ability or each model. MM is the genetic main effects models, DIAG is diagonal genetic covariance structured models, and FA1 to FA3 are factor ana- lytic structured models. Asterisks above each boxplot indicate levels of significant differences between each model and the MM baseline model from paired two-tailed t-tests where ‘*’ indicates p < 0.05, ‘**’ indicates p < 0.01 and ‘***’ indicates p < 0.001 Theoretical and Applied Genetics (2025) 138:241241  Page 20 of 24 and that this is likely due to the similar temperature and period of the growing season as well as the deep soil profile (Cooper et al. 1993; Cooper and Woodruff 1993; Mathews et al. 2007). Historical use of CIMMYT material, particu- larly in Australian breeding programmes in NSW and QLD, also supports these trends (Brennan and Quade 2006). We also note that the combination of high overall performance and low RMSD of many CIMMYT lines in CAIGE trials compared to specifically adapted checks highlights the value of these broadly adaptable materials across diverse Austral- ian growing environments. We found that factor loadings from a first stage MET model could be associated with ECs and were predictable with moderate predictive ability from second-stage multi- variate models using many ECs. The two-stage modelling approach enabled a large and complex MET dataset to be fitted in the first stage. Fitting one-stage MET analysis models including ECs with latent or factorial regression structures, as applied by Tolhurst et al. (2022) or Piepho and Blancon (2023), may be a more direct approach and could enhance the predictive ability and generalisability of models, but increases the number of parameters to be esti- mated in a mixed model. In contrast, we found that the ran- dom forest machine learning approach in the second-stage model including ECs was highly computationally efficient and flexible to modelling complex interactions among ECs without the need for restrictive assumptions of variance structures. Standard methods to integrate ECs into mixed model frameworks would also only model linear additive effects. While multivariate models provided better predic- tive ability than single EC associations, inference of direct causal effects from any single EC is problematic due to the high collinearity and potentially confounding effects among ECs. For example, the strong north to south temperature gradient is confounded with day length, and soil constraints in South and Western Australia are also known to be more severe than in the northern growing regions of NSW and QLD, which have summer rather than winter dominated rainfall patterns. Nevertheless, combined modelling of ECs is pragmatic considering non-additive effects among ECs are expected. For example, multiplicatively more severe drought stress might be expected due to low rainfall at locations with lighter soils with poorer water holding capacity, and pat- terns of soil characteristics and the rainfall gradient across the Australian grain belt are largely independent. The better performance of random forest prediction models, which can model interaction effects among predictors, also suggests that non-additive effects among ECs are important. Further research to validate and refine the GEI effects may involve trial experiments to control particular environmental effects, such as heat or drought stress, with more directed popula- tions of genotypes for quantitative trait loci mapping, simi- larly to Bennett et al. (2012), or with near-isogenic lines that contrast for key known genetic effects, similarly to Chapman et al. (2007). An advantage of our approach is that it builds upon a well-established framework of MET analysis with fac- tor analytic linear mixed models that many crop breeders routinely apply in large commercial breeding programmes (Smith et al. 2021). While analysis of large METs is often limited by computational capacity and software for fitting complex linear mixed model covariance structures (Smith et al. 2023), further work could better integrate ECs into single-stage MET analyses or build assumptions of covari- ance structures, that are of particular relevance to plant breeding contexts, into computationally efficient machine learning approaches. Definition of ECs from weather and soil data makes many assumptions of environmental effects on crop growth, and these methods could be further devel- oped integrating insights from crop models and plant physi- ology. Standardised methods to process environmental data from multiple sources and define ECs for Australian growing environments will be facilitated by the EC4MET R package developed as part of this work. In the present study, we used ECs that defined abiotic environmental stresses using weather and soil characteris- tics. Further work may integrate additional covariates that quantify biotic stresses  to further explain components of GEI. Several soil-borne pathogens are of particular impor- tance across Australian wheat growing environments and can be quantified using DNA-based tests from soil samples (Ophel-Keller et al. 2008), and Thompson et al. (2021) found strong correlations between FA loadings and root-lesion nematode (Pratylenchus thornei) abundance in soil profiles in a wheat MET. Although many of the soil-borne pathogen species, such as crown rot (Fusarium pseudograminearum), are known to occur in defined regions (Alahmad et al. 2018), inference of such biotic ECs into many new and unobserved environments would be difficult compared to the well inter- polated weather and soil datasets used here. Foliar pathogen incidence, such as stripe rust (Puccinia striiformis f.sp. trit- ici), is known to vary greatly due to year-to-year due to weather (Park 1990) and with incursions of new pathotypes (Wellings 2007). The static and non-static effects of abiotic ECs as demonstrated here, as well as more difficult to predict and quantify biotic effects on GEI, which may be influenced by management practices should also be considered prag- matically by plant breeders looking to define target regions. Crop breeding programmes are subject to inertia; in that breeding cycles operate over long timespans and changes to breeding strategies or direction are slow to come to frui- tion. Breeders should therefore consider prevailing envi- ronment types and climates several decades into the future. Methods presented here to predict environment types into many new environments could well be applied using cli- mate change projections to anticipate near-term impacts of Theoretical and Applied Genetics (2025) 138:241 Page 21 of 24  241 climate change and ensure that genetic gain is prioritised and achieved in environment types that are likely to become more prevalent. Current climate projections indicate a high certainty of increased temperatures across Australia but pro- jected changes in rainfall patterns are less consistent across independent climate models (Grose et al. 2020). Therefore, the higher average temperature iClass environment types identified here that were predicted to occur most frequently in the northern growing regions (iClasses pnn and pnp) may be considered of high importance in future climate change scenarios and highlights the potential value of novel CAIGE genotypes with high overall performance in these environ- ments. These methods could also be applied to plant breed- ing for other crop species or for improved informatics for farmers’ variety selection using large variety testing trial networks where the farmer’s target environment is a season that has not occurred yet and at a location for which there have been no trial experiments. Conclusion Insights from this work demonstrate how high-resolution environmental data can be processed and used by plant breeders to gain a better understanding of GEI effects to improve selection decisions. The value of specific adaptation indicated by static GEI effects, compared to broad adapta- tion influenced by non-static year-to-year GEI effects, will be important to consider. Finally, the increasing importance of climate change will make these approaches particularly relevant to modern crop breeding. Supplementary Information  The online version contains supplemen- tary material available at https://​doi.​org/​10.​1007/​s00122-​025-​05023-6. Acknowledgements  We thank the many financial and in-kind contri- butions of the research and plant breeding industry members of the CAIGE project. Germplasm used in this work was developed and made available by breeding programme members of the CGIAR. Author contributions  NF was involved in conception and design of the work; data analysis and interpretation; writing manuscript origi- nal draft; and software development. GG, VG and WT took part in germplasm development. AS participated in data acquisition. AZ was responsible for data analysis and interpretation. AW participated in software development. RT and JH were responsible for funding acquisi- tion; conception and design of the work; and project supervision. All authors reviewed and approved the manuscript. Funding  Open access funding provided by CSIRO Library Services. This work was funded by the CSIRO Early Research Career programme in collaboration with the CAIGE project which is funded by the Grains Research and Development Corporation (GRDC) project UOS2203- 004RTX and the Consultative Group on International Agricultural Research (CGIAR). Data and code availability  MET yield data are available through CAIGE project website (www.​caige​proje​ct.​org.​au). Some of the genotype data used in this study are available at (https://​doi.​org/​10.​ 25919/​3x88-​4p26; Fradgley 2025a. Weather data are available through the SILO website (www.​longp​addock.​qld.​gov.​au/​silo/), and soil data are available through the SLGA website (https://​esoil.​io/​TERNL​ andsc​apes/​Public/​Pages/​SLGA/​index.​html). The EC4MET R package (https://​doi.​org/​10.​5281/​zenodo.​15400​245; Fradgley and Whan 2025) has been developed to facilitate workflows to define ECs from weather and soil data for integration into MET analyses. Example analysis R scripts for fitting MET models and other functions described in this manuscript are provided (https://​doi.​org/​10.​5281/​zenodo.​15400​263; Fradgley 2025b). Declarations  Conflict of interest  The authors have no relevant financial or non-fi- nancial interests to disclose. Open Access  This article is licensed under a Creative Commons Attri- bution 4.0 International License, which permits use, sharing, adapta- tion, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visithttp://​creat​iveco​mmons.​org/​licen​ses/​by/4.​0/. References Alahmad S, Simpfendorfer S, Bentley AR, Hickey LT (2018) Crown rot of wheat in Australia: Fusarium pseudograminearum taxonomy, population biology and disease management. Australas Plant Pathol 47(3):285–299. https://​doi.​org/​10.​1007/​s13313-​018-​0554-z Araújo MS, Chaves SFS, Dias LAS, Ferreira FM, Pereira GR, Bezerra ARG, Alves RS, Heinemann AB, Breseghello F, Carneiro PCS, Krause MD, Costa-Neto G, Dias KOG (2024) GIS-FA: an approach to integrating thematic maps, factor-analytic, and enviro- typing for cultivar targeting. Theor Appl Genet 137(4):80. https://​ doi.​org/​10.​1007/​s00122-​024-​04579-z Bakare MA, Kayondo SI, Aghogho CI, Wolfe MD, Parkes EY, Kula- kow P, Egesi C, Jannink J-L, Rabbi IY (2022) Parsimonious geno- type by environment interaction covariance models for cassava (Manihot esculenta). Front Plant Sci. https://​doi.​org/​10.​3389/​fpls.​ 2022.​978248 Basford KE, Cooper M (1998) Genotype×environment interactions and some considerations of their implications for wheat breeding in Australia this review is one of a series commissioned by the advi- sory committee of the journal. Aust J Agric Res 49(2):153–174. https://​doi.​org/​10.​1071/​a97035 Bennett D, Reynolds M, Mullan D, Izanloo A, Kuchel H, Langridge P, Schnurbusch T (2012) Detection of two major grain yield QTL in bread wheat (Triticum aestivum L.) under heat, drought and high yield potential environments. Theor Appl Genet 125(7):1473– 1485. https://​doi.​org/​10.​1007/​s00122-​012-​1927-2 Boer MP, Wright D, Feng L, Podlich DW, Luo L, Cooper M, van Eeuwijk FA (2007) A mixed-model quantitative trait loci (QTL) analysis for multiple-environment trial data using environmen- tal covariables for QTL-by-environment interactions, with an https://doi.org/10.1007/s00122-025-05023-6 http://www.caigeproject.org.au https://doi.org/10.25919/3x88-4p26 https://doi.org/10.25919/3x88-4p26 http://www.longpaddock.qld.gov.au/silo/ https://esoil.io/TERNLandscapes/Public/Pages/SLGA/index.html https://esoil.io/TERNLandscapes/Public/Pages/SLGA/index.html https://doi.org/10.5281/zenodo.15400245 https://doi.org/10.5281/zenodo.15400263 http://creativecommons.org/licenses/by/4.0/ https://doi.org/10.1007/s13313-018-0554-z https://doi.org/10.1007/s00122-024-04579-z https://doi.org/10.1007/s00122-024-04579-z https://doi.org/10.3389/fpls.2022.978248 https://doi.org/10.3389/fpls.2022.978248 https://doi.org/10.1071/a97035 https://doi.org/10.1007/s00122-012-1927-2 Theoretical and Applied Genetics (2025) 138:241241  Page 22 of 24 example in maize. Genetics 177(3):1801–1813. https://​doi.​org/​ 10.​1534/​genet​ics.​107.​071068 Breiman L (2001) Random forests. Mach Learn 45(1):5–32. https://​ doi.​org/​10.​1023/A:​10109​33404​324 Brennan JP, Fox PN (1998) Impact of CIMMYT varieties on the genetic diversity of wheat in Australia, 1973–1993. Aust J Agric Res 49(2):175–178. https://​doi.​org/​10.​1071/​a97065 Brennan JP, Quade KJ (2006) Evolving usage of materials from CIM- MYT in developing Australian wheat varieties. Aust J Agric Res 57(9):947–952. https://​doi.​org/​10.​1071/​AR054​00 Burgueño J, Crossa J, Cornelius PL, Yang R-C (2008) Using factor analytic models for joining environments and genotypes with- out crossover genotype × environment interaction. Crop Sci 48(4):1291–1305. https://​doi.​org/​10.​2135/​crops​ci2007.​11.​0632 Butler DG, Cullis BR, Gilmour AR, Gogel BJ, Thompson R (2017) ASreml-R reference manual version 4. VSN International Ltd, Hemel Hempstead, HP1 1ES, UK Callister AN, Costa-Neto G, Bradshaw BP, Elms S, Crossa J, Brawner JT (2024) Enviromic prediction enables the characterization and mapping of Eucalyptus globulus Labill breeding zones. Tree Genet Genomes 20(1):3. https://​doi.​org/​10.​1007/​s11295-​023-​01636-4 Cane K, Eagles HA, Laurie DA, Trevaskis B, Vallance N, Eastwood RF, Gororo NN, Kuchel H, Martin PJ (2013) Ppd-B1 and Ppd-D1 and their effects in southern Australian wheat. Crop Pasture Sci 64(2):100–114. https://​doi.​org/​10.​1071/​CP130​86 Chai Y, Pardey PG, Silverstein KAT (2022) Scientific selection: a cen- tury of increasing crop varietal diversity in US wheat. Proc Natl Acad Sci U S A 119(51):e2210773119. https://​doi.​org/​10.​1073/​ pnas.​22107​73119 Chapman SC, Mathews KL, Trethowan RM, Singh RP (2007) Relation- ships between height and yield in near-isogenic spring wheats that contrast for major reduced height genes. Euphytica 157(3):391– 397. https://​doi.​org/​10.​1007/​s10681-​006-​9304-3 Cockram J (2024) Genome-wide resources for genetic locus discovery and gene functional analysis in wheat. In: Appels R, Eversole K, Feuillet C, Gallagher D (eds) The wheat genome. Springer Inter- national Publishing, pp 289–320 Cooper M, Woodruff DR (1993) Predicting grain yield in Australian environments using data from CIMMYT international wheat per- formance trials 3. Testing predicted correlated response to selec- tion. Field Crops Res 35(3):191–204. https://​doi.​org/​10.​1016/​ 0378-​4290(93)​90153-E Cooper M, DeLacy IH, Byth DE, Woodruff DR (1993) Predicting grain yield in Australian environments using data from CIMMYT international wheat performance trials. 2. The application of clas- sification to identify environmental relationships which exploit correlated response to selection. Field Crops Res 32(3):323–342. https://​doi.​org/​10.​1016/​0378-​4290(93)​90040-T Cooper M, Tang T, Gho C, Hart T, Hammer G, Messina C (2020) Integrating genetic gain and gap analysis to predict improvements in crop productivity. Crop Sci 60(2):582–604. https://​doi.​org/​10.​ 1002/​csc2.​20109 Costa-Neto GMF, Morais Júnior OP, Heinemann AB, de Castro AP, Duarte JB (2020) A novel GIS-based tool to reveal spatial trends in reaction norm: upland rice case study. Euphytica 216(3):37. https://​doi.​org/​10.​1007/​s10681-​020-​2573-4 Costa-Neto G, Galli G, Carvalho HF, Crossa J, Fritsche-Neto R (2021) EnvRtype: a software to interplay enviromics and quan- titative genomics in agriculture. G3 Genes|genomes|genetics 11(4):jkab040. https://​doi.​org/​10.​1093/​g3jou​rnal/​jkab0​40 Costa-Neto G, Crespo-Herrera L, Fradgley N, Gardner K, Bentley AR, Dreisigacker S, Fritsche-Neto R, Montesinos-López OA, Crossa J (2023) Envirome-wide associations enhance multi-year genome-based prediction of historical wheat breeding data. G3 13(2):jkac313. https://​doi.​org/​10.​1093/​g3jou​rnal/​jkac3​13 Crossa J, Montesinos-Lopez OA, Costa-Neto G, Vitale P, Martini JWR, Runcie D, Fritsche-Neto R, Montesinos-Lopez A, Pérez-Rod- ríguez P, Gerard G, Dreisigacker S, Crespo-Herrera L, Pierre CS, Lillemo M, Cuevas J, Bentley A, Ortiz R (2024) Machine learning algorithms translate big data into predictive breeding accuracy. Trends Plant Sci. https://​doi.​org/​10.​1016/j.​tplan​ts.​2024.​09.​011 Cullis BR, Smith A, Hunt C, Gilmour A (2000) An examination of the efficiency of Australian crop variety evaluation programmes. J Agric Sci 135(3):213–222. https://​doi.​org/​10.​1017/​S0021​85969​ 90081​63 Cullis BR, Smith AB, Coombes NE (2006) On the design of early generation variety trials with correlated data. J Agric Biol Environ Stat 11(4):381–393. https://​doi.​org/​10.​1198/​10857​1106X​154443 Cullis B, Cocks N, Smith A, Butler D (2018) Sparse multi-environment trial designs for early stage selection experiments in plant breed- ing programmes. National Institute for Applied Statistics Research of Australia Working Paper Series, University of Wollongong. https://​Niasra.​Uow.​Edu.​Au/​Worki​ngpap​ers/​Index.​Html. Damesa TM, Möhring J, Worku M, Piepho H-P (2017) One step at a time: stage-wise analysis of a series of experiments. Agron J 109(3):845–857. https://​doi.​org/​10.​2134/​agron​j2016.​07.​0395 Dreccer MF, Fainges J, Whish J, Ogbonnaya FC, Sadras VO (2018) Comparison of sensitive stages of wheat, barley, canola, chickpea and field pea to temperature and water stress across Australia. Agric for Meteorol 248:275–294. https://​doi.​org/​10.​1016/j.​agrfo​ rmet.​2017.​10.​006 Eagles HA, Cane K, Trevaskis B, Vallance N, Eastwood RF, Gororo NN, Kuchel H, Martin PJ (2014) Ppd1, Vrn1, ALMT1 and Rht genes and their effects on grain yield in lower rainfall environ- ments in southern Australia. Crop Pasture Sci 65(2):159–170. https://​doi.​org/​10.​1071/​CP133​74 Endelman JB, Jannink J-L (2012) Shrinkage estimation of the realized relationship matrix. G3 Genes|genomes|genetics 2(11):1405– 1413. https://​doi.​org/​10.​1534/​g3.​112.​004259 Fairlie W, Hughes D, Cullis B, Edwards J, Kuchel H (2024) Genotype- by-environment interaction for wheat falling number performance due to late maturity α-amylase. Crop Sci. https://​doi.​org/​10.​1002/​ csc2.​21348 FAO (2024) FAOSTAT. Crops and Livestock Products. https://​www.​ fao.​org/​faost​at/​en/#​data/​TCL Finlay KW, Wilkinson GN (1963) The analysis of adaptation in a plant- breeding programme. Aust J Agric Res 14(6):742–754. https://​doi.​ org/​10.​1071/​ar963​0742 Flohr BM, Hunt JR, Kirkegaard JA, Evans JR (2017) Water and tem- perature stress define the optimal flowering period for wheat in south-eastern Australia. Field Crops Res 209:108–119. https://​doi.​ org/​10.​1016/j.​fcr.​2017.​04.​012 Fradgley N, Whan A (2025) NickFlagleaf/EC4MET: V1.1 (Version v1.1) [Computer software]. Zenodo. https://​doi.​org/​10.​5281/​ zenodo.​15400​245 Fradgley N (2025a) 40K SNP array genotype data for CAIGE wheat lines (Version 1). CSIRO. https://​doi.​org/​10.​25919/​3x88-​4p26 Fradgley N (2025b) NickFlagleaf/GxE-FA-EC-Prediction: V1.0 (Ver- sion v1,0) [Computer software]. Zenodo. https://​doi.​org/​10.​5281/​ zenodo.​15400​263 Friedman JH, Hastie T, Tibshirani R (2010) Regularization paths for generalized linear models via coordinate descent. J Stat Softw 33:1–22. https://​doi.​org/​10.​18637/​jss.​v033.​i01 Gogel B, Smith A, Cullis B (2018) Comparison of a one- and two- stage mixed model analysis of Australia’s National Variety Trial Southern Region wheat data. Euphytica 214(2):44. https://​doi.​org/​ 10.​1007/​s10681-​018-​2116-4 Grose MR, Narsey S, Delage FP, Dowdy AJ, Bador M, Boschat G, Chung C, Kajtar JB, Rauniyar S, Freund MB, Lyu K, Rashid H, Zhang X, Wales S, Trenham C, Holbrook NJ, Cowan T, Alexan- der L, Arblaster JM, Power S (2020) Insights from CMIP6 for https://doi.org/10.1534/genetics.107.071068 https://doi.org/10.1534/genetics.107.071068 https://doi.org/10.1023/A:1010933404324 https://doi.org/10.1023/A:1010933404324 https://doi.org/10.1071/a97065 https://doi.org/10.1071/AR05400 https://doi.org/10.2135/cropsci2007.11.0632 https://doi.org/10.1007/s11295-023-01636-4 https://doi.org/10.1071/CP13086 https://doi.org/10.1073/pnas.2210773119 https://doi.org/10.1073/pnas.2210773119 https://doi.org/10.1007/s10681-006-9304-3 https://doi.org/10.1016/0378-4290(93)90153-E https://doi.org/10.1016/0378-4290(93)90153-E https://doi.org/10.1016/0378-4290(93)90040-T https://doi.org/10.1002/csc2.20109 https://doi.org/10.1002/csc2.20109 https://doi.org/10.1007/s10681-020-2573-4 https://doi.org/10.1093/g3journal/jkab040 https://doi.org/10.1093/g3journal/jkac313 https://doi.org/10.1016/j.tplants.2024.09.011 https://doi.org/10.1017/S0021859699008163 https://doi.org/10.1017/S0021859699008163 https://doi.org/10.1198/108571106X154443 https://Niasra.Uow.Edu.Au/Workingpapers/Index.Html https://doi.org/10.2134/agronj2016.07.0395 https://doi.org/10.1016/j.agrformet.2017.10.006 https://doi.org/10.1016/j.agrformet.2017.10.006 https://doi.org/10.1071/CP13374 https://doi.org/10.1534/g3.112.004259 https://doi.org/10.1002/csc2.21348 https://doi.org/10.1002/csc2.21348 https://www.fao.org/faostat/en/#data/TCL https://www.fao.org/faostat/en/#data/TCL https://doi.org/10.1071/ar9630742 https://doi.org/10.1071/ar9630742 https://doi.org/10.1016/j.fcr.2017.04.012 https://doi.org/10.1016/j.fcr.2017.04.012 https://doi.org/10.5281/zenodo.15400245 https://doi.org/10.5281/zenodo.15400245 https://doi.org/10.25919/3x88-4p26 https://doi.org/10.5281/zenodo.15400263 https://doi.org/10.5281/zenodo.15400263 https://doi.org/10.18637/jss.v033.i01 https://doi.org/10.1007/s10681-018-2116-4 https://doi.org/10.1007/s10681-018-2116-4 Theoretical and Applied Genetics (2025) 138:241 Page 23 of 24  241 Australia’s future climate. Earths Future 8(5):e2019EF001469. https://​doi.​org/​10.​1029/​2019E​F0014​69 Grundy MJ, Rossel RAV, Searle RD, Wilson PL, Chen C, Gregory LJ, Grundy MJ, Rossel RAV, Searle RD, Wilson PL, Chen C, Gregory LJ (2015) Soil and landscape grid of Australia. Soil Res 53(8):835–844. https://​doi.​org/​10.​1071/​SR151​91 Helland IS (1990) Partial least squares regression and statistical mod- els. Scand J Stat 17(2):97–114 Heslot N, Akdemir D, Sorrells ME, Jannink J-L (2014) Integrating environmental covariates and crop modeling into the genomic selection framework to predict genotype by environment interac- tions. Theor Appl Genet 127(2):463–480. https://​doi.​org/​10.​1007/​ s00122-​013-​2231-5 Jarquín D, Crossa J, Lacaze X, Du Cheyron P, Daucourt J, Lorgeou J, Piraux F, Guerreiro L, Pérez P, Calus M, Burgueño J, de los Campos G (2014) A reaction norm model for genomic selec- tion using high-dimensional genomic and environmental data. Theor Appl Genet 127(3):595–607. https://​doi.​org/​10.​1007/​ s00122-​013-​2243-1 Jeffrey SJ, Carter JO, Moodie KB, Beswick AR (2001) Using spatial interpolation to construct a comprehensive archive of Australian climate data. Environ Model Softw 16(4):309–330. https://​doi.​ org/​10.​1016/​S1364-​8152(01)​00008-1 Keeble-Gagnère G, Pasam R, Forrest KL, Wong D, Robinson H, Godoy J, Rattey A, Moody D, Mullan D, Walmsley T, Daetwyler HD, Tibbits J, Hayden MJ (2021) Novel design of imputation-enabled SNP arrays for breeding and research applications supporting multi-species hybridization. Front Plant Sci. https://​doi.​org/​10.​ 3389/​fpls.​2021.​756877 Kempton RA (1984) The use of biplots in interpreting variety by envi- ronment interactions. J Agric Sci 103(1):123–135. htt