The analyses were divided into three stages. In the first stage, individual analyses were conducted for each experiment to assess the existence of genetic variability and estimate selective accuracy and means for each experiment. In the second stage, different FAMM models were tested in multi-environment analysis, considering only the experiments that showed genetic variability (significant genotypic variance component of clones). In the third stage, the scores from the best-fitting FAMM model were used in the selection process.

Estimates of variance components were obtained using the REML (residual maximum likelihood) method (Patterson and Thompson, 1971); and the EBLUEs (empirical best linear unbiased estimates) of fixed effects and EBLUPs (empirical best linear unbiased predictions) of random effects were calculated using the equations of Henderson's mixed models (Henderson et al., 1959). The analyses using the mixed models methodology were performed using the Echidna Mixed Models software (Gilmour, 2021). Images were obtained using the 'matplotlib.pyplot' package in Python software.

## 2.2.1 Individual analyses

Individual analyses were performed based on the model presented in Eq. (1). For experiments designed as randomized complete blocks (RCB), the effect of incomplete blocks (ub) was removed from the model.

y = 1µ + Xrτr + Zcuc + Zbub +e (1)

where y (*Nj* × 1) is the vector of phenotypic observations and *N**j* is the number of plots in experiment j; 1 (*Nj* × 1) is a vector of ones; µ (1 × 1) is the intercept; τ*r* (*r* × 1) is the vector of fixed effects of repetitions associated with the X*r* (*Nj* × *r*) matrix (assuming a full-rank matrix) and r is the number of repetitions; u*c* (*c* × 1) is the vector of random genotypic effects of clones associated with the Z*c* (*Nj* × *c*) matrix and c is the number of clones evaluated in experiment j; u*b* (*b* × 1) is the vector of random effects of blocks within each repetition in the experiment associated with the Z*b* (*Nj* × *b*) matrix and b is the number of blocks within each repetition in experiment j; and e (*Nj* × 1) is the vector of random errors for experiment j.

It was assumed that the random effect vectors u*c*, u*b* and e are mutually independent and follow a multivariate Gaussian distribution with zero mean and (co)variance var(u*c*) = \({\sigma }_{c}^{2}\)I*c*, var(u*b*) = \({\sigma }_{b}^{2}\)I*b*, and var(e) = \({\sigma }^{2}\)I*Nj*, where \({\sigma }_{c}^{2}\), \({\sigma }_{b}^{2}\), and \({\sigma }^{2}\) are the variance components associated with the genotypic, block, and residual effects, respectively, and I*c* (*c* × *c*), I*b* (*b* × *b*), and I*Nj* (*Nj* × *Nj*) are identity matrices.

The significance of the variance component \({\sigma }_{c}^{2}\) was tested using the likelihood ratio test (LRT) presented in Eq. (2).

where ln() is the natural logarithm, L1 is the maximum point of the residual likelihood function for the reduced model, and L2 is the maximum point of the residual likelihood function for the full model.

accuracy of the genotypic effects of clones was estimated for each experiment using the expression in Eq. 3.

$${r}_{\widehat{c}c}=\sqrt{1-\frac{{\upsilon }_{c}}{{\sigma }_{c}^{2}}}$$

3

where \({\upsilon }_{c}\) is the average prediction error variance (PEV) associated with the vector of genotypic effects of clones (u*c*).

The mean of each experiment was also estimated using the linear combination presented in Eq. (4).

$$\stackrel{-}{\text{y}}=\mu +{\stackrel{-}{{\tau }}}_{r}$$

4

where \({\stackrel{-}{{\tau }}}_{r}\) is the average effect of repetitions or complete blocks.

## 2.2.2 Multi-environment analyses

Considering the set of experiments selected in section "2.2.1. Individual analyses", a multi-environment analysis was conducted using the model presented in Eq. (5). Comparison of the model of Eq. (1) with the model of Eq. (5) shows that Eq. (5) includes the vector of fixed effects of the experiment (τ*e*) and the vector of random effects of the clone-by-experiment interaction (u*ce*); it also replaces the vector of fixed effects of repetition (τ*r*) with the vector of fixed effects of repetition within experiments (τ*re*).

y = 1µ + X*e*τ*e* + X*re*τ*re* + Z*c*u*c* + Z*ce*u*ce* + Z*b*u*b* + e (5)

where y (*N* × 1) is the vector of phenotypic observations from all experiments and *N* is the total number of plots; 1 (*N* × 1) is a vector of ones; µ (1 × 1) is the intercept; τ*e* (*e* × 1) is the vector of fixed effects of experiments associated with the X*e* (*N* × *e*) matrix (assuming a full-rank matrix), where e is the number of experiments; τ*re* (*r*× 1) is the vector of fixed effects of repetitions within experiments associated with the X*re* (*N* × *r*) matrix (assuming a full-rank matrix), where *r* is the number of repetitions; u*c* (*c* × 1) is the vector of random genotypic effects of clones (predicted based on the mean of e experiments) associated with the Z*c* (*N* × *c*) matrix, where *c* is the total number of clones evaluated; u*ce* (*ce* × 1) is the vector of random effects of the clone-by-experiment interaction associated with the Z*ce* (*N* × *ce*) matrix; u*b* (*b* × 1) is the vector of random effects of blocks within repetitions in each experiment associated with the Z*b* (*N* × *b*) matrix, where b is the total number of blocks; and e (*N* × 1) is the vector of random errors.

It was assumed that the random effect vectors u*c*, u*ce*, u*b*, and e are mutually independent and follow a multivariate Gaussian distribution with zero mean and (co)variance matrices var(u*c*) = \({\sigma }_{c}^{2}\)I*c*, var(u*ce*) = \({\sigma }_{ce}^{2}\)I*c*⊗I*e*, var(u*b*) = \({\oplus }_{j=1}^{e}{\sigma }_{bj}^{2}\)I*bj*, and var(e) = \({\oplus }_{j=1}^{e}{\sigma }_{j}^{2}\)I*Nj*, where \({\sigma }_{c}^{2}\), \({\sigma }_{ce}^{2}\), \({\sigma }_{bj}^{2}\), and \({\sigma }_{j}^{2}\) are the variance components associated with the average genotypic effects of clones, clone-by-experiment interaction, blocks within repetitions in experiment *j*, and residuals in experiment *j*. I*c* (*c* × *c*), I*e* (*e* × *e*), I*bj* (*bj* × *bj*), and I*Nj* (*Nj* × *Nj*) are identity matrices; and ⊗ and ⊕ are the Kronecker product and direct sum operators, respectively. The block and residual structures were considered heterogeneous to accommodate different patterns of experimental precision due to variation in the number of repetitions and dimension of the experiments (Table 1) and to the presence of different experimental designs.

From the model presented in Eq. (5), the vector of genotypic effects of clones capitalizing on the G×E interaction (u*g*) can be predicted through the combination of the u*c* and u*ce* vectors (Eq. 6). Thus, the variance of the u*g* vector (Eq. 7) results in a more flexible (co)variance structure, as the (co)variance matrix between experiments (G*e*) can take on different forms, ranging from a simpler form like compound symmetry, which has two parameters, to the more general (unstructured) form with e(e + 1)/2 parameters (Smith et al., 2015).

u*g* = (1*e*⊗I*c*)u*c* + u*ce* (6)

var(u*g*) = G*e*⊗I*c* (7)

To ensure greater parsimony without losing the ability to accommodate the heterogeneities of variances and covariances typically observed in a G×E interaction scenario, the FAMM model was used (Smith et al., 2001). Thus, the (co)variance matrix between experiments can be written as G*e* = ΛΛ*T* + ψ, where Λ (*e* × *k*) is the factor loading matrix, with *k* being the number of factors used (1, 2, ..., *k*); and ψ (*e* × *e*) is the diagonal matrix of specific variances for each experiment.

From the FAMM model, it is also possible to predict the performance of clones capitalized by the G×E interaction on the original scale of the trait. To do this, it is necessary to re-parameterize Eq. (6) as presented in Eq. (8).

u*g* = (Λ⊗I*c*)f + δ (8)

where f (*ck* × 1) is the vector of factor scores from the factor analysis, and δ (*ce* × 1) is the vector of residual effects from the factor analysis.

A compound symmetry (CS) model (Smith et al., 2015) and a sequence of FAMM models, ranging from FA1 to FA*k*, were used. The goodness-of-fit of the models was assessed using the Akaike information criterion (AIC) (Akaike, 1974), presented in Eq. (9), as well as the LRT (Eq. 2) test, with the latter used to compare the FA*k* order FAMM models.

AIC = -2ℓ + 2p (9)

where ℓ is the logarithm of the maximum point of the residual likelihood function, and *p* is the number of parameters.

## 2.2.2.1 Factor rotation

To interpret the factors, both the factor loading matrix (Eq. 10) and the factor scores vector (Eq. 11) of the selected FA*k* model were rotated using an orthogonal matrix T (Cullis et al., 2010; Smith et al., 2001, 2005, 2015, 2021a, 2021b; Smith & Cullis, 2018). This matrix was obtained in two ways: through singular value decomposition (SVD) (Smith & Cullis, 2018) and the varimax method (Kaiser, 1958), using the 'base' and 'stats' packages in the R software (R Core Team, 2016; R Core Team, 2022).

Λ* = ΛT (10)

f* = (T*T*⊗I*c*)f (11)

To enhance interpretation of the FAMM models, the percentage of genetic variation retained in the FA*k* order models was estimated for individual experiments (νj) and overall (ν) using Eqs. (12) and (13), respectively. Additionally, Pearson correlations between the genotypic effects of clones in experiment j and the rotated factor k scores (Eq. 14) were calculated. Furthermore, approximate genetic correlation matrices between experiments were also derived for the FA*k* models.

$$\nu j =\frac{100diag\left(\varLambda \varLambda T\right)j}{diag\left(\varLambda \varLambda T + \psi \right)j}$$

12

$$\nu =\frac{100tr\left(\varLambda \varLambda T\right)}{tr\left(\varLambda \varLambda T + \psi \right)}$$

13

$${\rho }_{jk}=\frac{cov\left({u}_{{g}_{j}}, {f}_{k}\right)}{\sqrt{var\left({u}_{{g}_{j}}\right)\times var\left({f}_{k}\right)}}$$

14

## 2.2.2.2 Selection of superior clones

For selection of clones based on adaptability and stability, two methodologies were used: OP-RMSD (overall performance and root mean square deviation) and genotype-ideotype distance (GI), based on the rotated scores from the selected FAMM model. For OP-RMSD, the SVD method was used to rotate the loadings and scores, where the first factor represents the overall behavior of the genotypes, while the remaining factors capture the G×E interaction (Smith & Cullis, 2018). For GI, in turn, the varimax method was used, which allows groups of environments to be associated with the factors and, consequently, with the performance of the clones. Thus, if the most important loadings for the k factors have predominantly the same sign, a selection direction can be defined based on the scores (positive or negative values).

For the GI distance, the scores of the *k* factors of the selected model (assuming that the vector f is ordered by clones) were combined using Eq. (16). adapted from Rocha et al. (2018), using a reference ideotype that reflects maximum performance (highest EBLUP) across all environments (Eq. 15). Geometrically, the elements of the PGI vector can be interpreted as the inverse of the normalized Euclidean distance within the range of 0 to 1. Based on the PGI vector, the clones were sorted in descending order, with higher PGI values corresponding to clones with greater adaptability and stability.

f(ideotype) = Λ+(ug(ideotype) - δ(ideotype)) (15)

PGI = w-1{(\({1}_{\varvec{k}}^{\varvec{T}}\)⊗I*c*)[f - (1c⊗Ik)f(ideotype)]2}-1/2 (16)

where **f**(ideotype) (*k* × 1) is the ideotype score vector, **Λ**+ (*k* × *e*) is the Moore-Penrose generalized inverse of the matrix **Λ** (*e* × *k*), **u***g*(ideotype) is the ideotype EBLUP vector (maximum performance across all environments), **δ**(ideotype) is the factor analysis residuals corresponding to the elements of the **u***g*(ideotype), vector, **P**GI is the spatial probability vector of the genotype-ideotype index, \({1}_{k}^{T}\) (1 × *k*) is a vector of ones, and *w* is a scalar representing the sum of the elements of the (\({1}_{k}^{T}\)⊗**I***c*)[**f** - (**1***c*⊗**I***k*)**f**(ideotype)]2}−1/2 vector. The Moore-Penrose generalized inverse was obtained using the 'MASS' package (Ripley et al., 2013) in the R software.