Accounting for spatial trends in multi-environment diallel analysis in maize breeding

Spatial trends represent an obstacle to genetic evaluation in maize breeding. Spatial analyses can correct spatial trends, which allow for an increase in selective accuracy. The objective of this study was to compare the spatial (SPA) and non-spatial (NSPA) models in diallel multi-environment trial analyses in maize breeding. The trials consisted of 78 inter-populational maize hybrids, tested in four environments (E1, E2, E3, and E4), with three replications, under a randomized complete block design. The SPA models accounted for autocorrelation among rows and columns by the inclusion of first-order autoregressive matrices (AR1 ⊗ AR1). Then, the rows and columns factors were included in the fixed and random parts of the model. Based on the Bayesian information criteria, the SPA models were used to analyze trials E3 and E4, while the NSPA model was used for analyzing trials E1 and E2. In the joint analysis, the compound symmetry structure for the genotypic effects presented the best fit. The likelihood ratio test showed that some effects changed regarding significance when the SPA and NSPA models were used. In addition, the heritability, selective accuracy, and selection gain were higher when the SPA models were used. This indicates the power of the SPA model in dealing with spatial trends. The SPA model exhibits higher reliability values and is recommended to be incorporated in the standard procedure of genetic evaluation in maize breeding. The analyses bring the parents 2, 10 and 12, as potential parents in this microregion.


Introduction
Maize (Zea mays L.) is the most cultivated crop worldwide [1]. Quantitative traits, such as grain yield and plant height, are controlled by several genes and are highly influenced by the environment [2]. In this scenario, the genotype-by-environment (G × E) interaction, also known as phenotypic plasticity [3,4], plays an essential role in phenotypic expression and can lead to difficulties in genetic selection [4].
Diallel mating designs are used for progeny tests [5] and are widely adopted in plant breeding [6][7][8][9]. These mating designs allow the evaluation of general and specific combining abilities, which are additive genetic effects based on general combining ability (GCA), and dominance genetic effect based on specific combining ability (SCA) [10,11]. The GCA is given by the mean of the performance of a particular individual in combination with many others, and the SCA is the genetic effect of a specific cross [11,12].
In diallel MET analyses, there are many sources of variation [38]. As main environmental sources, Burgueño et al. [39] cited soil structure, moisture, light interception, pathogen infections, and even crop management. To deal with spatial trends, one of the most adopted model structures is the separable first-order autoregressive (AR1�AR1), which considers a partitioned residual structure (e), composed of the correlated error (ξ) and, the independent random error (η) [40,41], where e = ξ+η. For MET analyses, Smith et al. [31,42] proposed the two-stage analysis, being stage 1 the modeling of AR1�AR1 spatial structure of residuals, to determine the best-fitted model. The output of stage 1 is imputed in stage 2, where the MET analysis is performed. This analysis has recently been applied in wheat [43] and soybean [44] breeding.
The two-stage analysis can present some issues, such as convergence failure and high computational demand, owing to the large amount of data [45]. Selective accuracy, heritability, and gain with selection can be used to compare different statistical models [46]. Selective accuracy indicates the precision of the predicted genotypic values [46], guides the strategies of breeding programs, and assists in decision-making [47].

Genetic material
A diallel mating design (13 × 13), with no reciprocal crossing, was implemented (Fig 1 and S1 Table in S1 File). The genetic material consisted of 78 interpopulational hybrids evaluated in four environments (E1 to E4). Six commercial hybrids (F 1 ), which were widely adopted in the region, were used as checks (S2 Table in S1 File).

Experimental data
The trials were carried out during the 2018 winter season at four locations in the southwestern Goiás State, Brazil: the environments, E1 and E2, in the municipality of Jataí, and the environments, E3 and E4, in the municipalities of Caiapônia and Mineiros, respectively. The planting and harvesting dates, geographic coordinates, altitudes, average temperatures, and precipitation data for each environment are presented in S3 Table and S1 Fig in S1 File [58].
According to Alvares et al. [59], the weather in the southwestern region Goiás State is wet and temperate, with dry winters and hot summers (Cwa). The average annual temperature is around 21.5˚C and the average rainfall is between 1.400 and 2.000 mm year -1 . Agricultural practices for maize crops in Brazil [60] were adopted in this study and irrigation was not applied in the field (S3 Table in S1 File). The experiments were conducted in a randomized complete block design with three replications and 44 plants per plot. Plots were 4m long, with four rows spaced 0.45 m apart and a plot size of 7.2 m 2 .
To evaluate the grain yield (GY) trait, the ears at physiological maturity were harvested and shelled. Then, the grain weight and grain moisture percentage were recorded, and the GY trait was calculated at 13% moisture. All phenotypic measurements were taken from the middle rows, leaving the two border rows.

Statistical analyses
The estimation of variance components and the prediction of genotypic values for the GY trait were made using the REML/BLUP procedure [16,17], according to Gilmour et al. [61].
The individual analyses, considering the randomized complete block design with one observation per plot, was given by the following equation: where y is the vector of phenotypic data, τ is the vector of replications and checks (assumed to be fixed) added to the general mean, u g is the vector of the GCA effect (assumed as random), u s is the vector of the SCA effect (assumed as random), and e is the vector of residuals (random). Uppercase letters refer to the incidence matrices for these effects. In this model, u g Nð0; s 2 g Þ; u s Nð0; s 2 s Þ; and e N(0,R); where s 2 g is the GCA variance (related to additive genetic variance), s 2 s is the SCA variance (related to dominance genetic variance), and R refers to the residual variance matrix. Models accounting for different residual variance structures considering the correlations among observations (rows and columns) were tested.
The joint analysis, considering the randomized complete block design with one observation per plot and four environments, was performed using the following equation: where y is the vector of phenotypic data, τ is the vector of block-locations-checks combinations (assumed to be fixed), which comprises the effects of environment, replication within the environment, and checks added to the overall mean; u g is the vector of additive genetic effects (assumed as random); u s is the vector of dominant genetic effects (assumed as random); u ge is the vector of the interaction between additive genetic effects and environments (random), u se is the vector of the interaction between dominance genetic effects and environments (random), and e is the vector of residuals (random). Uppercase letters refer to the incidence matrices for these effects. In this model: u g Nð0; Is 2 g Þ; u s Nð0; Is 2 s Þ; u ge Nð0; Is 2 ge Þ; u se Nð0; Is 2 se Þ, and eN(0,R); where s 2 g is the GCA variance (related to additive genetic variance), s 2 s is the SCA variance (related to Interpopulational hybrids development process. The process occurred as following: (I) the most recommended commercial genetic materials (F 1 ) were selected as parents 1 to 13; (II) these parents (F 1 ) were cultivated in isolated fields; (III) the isolated fields were harvested and generated 12 F 2 hybrids and one rustic genetic material; (IV) the F 2 genetic materials were crossed artificially (controlled pollination) among themselves, in a diallel mating design; (V) as a result it was generated 78 interpopulational hybrids; and, (VI) the 78 interpopulational hybrids were cultivated in 4 environments (E1 to E4). https://doi.org/10.1371/journal.pone.0258473.g001 dominance genetic variance), s 2 ge is the GCA by environment interaction variance, s 2 se is the SCA by environment interaction variance, and R refers to the residual variance matrix. Models accounting for different residual variance structures considering the correlations among observations (rows and columns) were tested.
SPA and NSPA analyses were performed for comparison purposes.
Modelling non-genetic effects. In the individual analyses, the residual (co)variance structure could account for heterogeneity within the trial. The best fitted models were tested by the inclusion or the absence of the spatial information (given by rows and columns) in the fixed and random parts of the model. The joint analysis accounted for heterogeneity across environments in residual and genetic (co)variance matrices in two different ways: (i) by assuming spatially independent observations and (ii) by allowing spatial autocorrelations among observations, indexed by rows and columns.
The presence of autocorrelation among rows and columns was considered by the inclusion of first-order autoregressive matrices (AR1) for rows (S r (ρ r )) and columns (S c (ρ c )). There are three possible models beyond the NSPA model (R ¼ s 2 e I r � I c ), where s 2 e is the residual variance, I r is the identity matrix of rows (r × r), and I c is the identity matrix of columns (c × c). The baseline model did not consider any correlations among the rows and columns. The following two residual structures account for the correlation among rows (M2) and columns (M3), respectively, by the inclusion of the AR1 structure. M2 and M3 contemplate the residual structure as R ¼ s 2 x S r ðr r Þ � I c and R ¼ s 2 x I r � S c ðr c Þ, respectively. The last model accounts for correlations among rows and columns (M4), being R ¼ s 2 x S r ðr r Þ � S c ðr c Þ, where s 2 x is the SPA variance among columns and among rows, � is the Kronecker product, and, AR1x, or, S x , represent first-order autoregressive correlation matrix, for x (row or column).
After selecting the best residual structure for each trial, the SPA information (row and column specifications) was included in the fixed and random parts of the model. First, the selected residual structure (MRx) of each trial was adopted as the first model, named Mx.1. Second, we tested the presence of fixed linear effects for rows and columns (Mx.2 and Mx. 3). Third, we tested the presence of random linear effects for rows and columns (Mx.4 and Mx.5). Finally, the measurement error was added to the random part of the model (Mx.6).
The fixed part of the model (rows and columns) deals with rows and columns as covariates, and was developed as a regression of the phenotypic responses as a function of the row and column positions. In the random part of the model, the rows and columns effects were included to quantify the variability among rows and columns. Finally, we also tested the inclusion of an independent error regarding the measurement error, which is an independent residual variance, s 2 Z . Fig 2 presents the flowchart of all steps to encounter the best SPA models for each trial by modeling the fixed and random effects and to select the best genetic structure when analyzing all trials jointly in a MET analysis.
Model selection, significance test, and gain with selection. To select the best-fit model, the conditional Bayesian information criterion (BIC c ) [62] was used. The comparison among models was based on the parameter vector (b k ), which stands for the covariance matrix to calculate the corrected LogL (LogL c ), as follows: where LogL c is the corrected LogL byb kbk ; p is the number of random parameters, q is the number of fixed parameters, and n is the model's degrees of freedom. The variances of the random effects were assessed using the likelihood ratio test (LRT) as follows [63]: where LogL R is the logarithm of the maximum residual likelihood function of the reduced model. For the LRT test, chi-square statistics with 1 degree of freedom and 5% probability of Type I error were considered. The additive genetic variance (ŝ 2 A ), dominance genetic variance (ŝ 2 D ), and broad-sense (H 2 ) heritability (given by the additive genetic variance among hybrids) were estimated for each trial using the following equations [64][65][66]: Flowchart of the spatial analysis in each trial and jointly. All steps of the spatial analysis by considering the spatial information (rows and columns) in the fixed and random parts of the model, in each trial. Following, it is considered all individual models to encounter the best joint model, and the most adequate genetic covariance structure in the multienvironment analysis. https://doi.org/10.1371/journal.pone.0258473.g002

PLOS ONE
whereŝ 2 g andŝ 2 s are the estimates of the general and specific combining ability variances, respectively.
In addition to the genetic parameters by each trial, the additive genetic variance (ŝ 2 A ), dominance genetic variance (ŝ 2 D ), additive by environment interaction variance (ŝ 2 AE ), dominance by environment interaction variance (ŝ 2 DE ), and broad-sense (H 2 ) heritability, were estimated considering the four environments in a joint analysis using the following equations [64][65][66]: In addition, the genotypic correlation across environments (r g ) and the overall mean (μ g ) for the GY trait was also calculated as follows: Selective accuracies (rÂ A and rD D ) were estimated for individual and joint (rĜ G ) analyses using the following expression [67]: ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffiffi ; and; ½11� rĜ G ¼ ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffiffi ffi ffi ffiffi where PEV is the prediction error variance, according to the respective effect obtained by the diagonal elements of the inverse of the coefficient matrix of the mixed model equations.
The correlations among the NSPA and SPA models were evaluated using the Pearson correlation coefficient. Cohen's Kappa coefficient (K) [68] was applied to estimate the agreement between the NSPA and SPA models. The agreement among the selected interpopulational hybrids was calculated as follows: where A is the number of matching selected interpopulational hybrids between the NSPA and SPA model rankings, C is the number of selected interpopulational hybrids due to chance (C = bD, where b is the selection intensity = 0.13), and D is the number of selected interpopulational hybrids (10). Contrasting SPA and NSPA analyses. After selecting the individual and joint SPA models, further analyses were conducted using the NSPA model to expose the relative importance of the SPA model. Therefore, all analyses performed under the spatial approach were also remade under a non-spatial approach, considering the independence among rows and columns.

Modelling non-genetic effects
In trials E1 and E2, the NSPA model (MR1), which did not consider any correlation among the rows, and among the columns, was the best fitted model (Table 1). In trial E3, the best fitted model (MR4) considered the existence of correlation among rows and columns in the residual effects, while in trial E4, the best fitted model (MR3) considered the correlation among columns in the residual effects (Table 1).

Modelling fixed and random effects
The E1, E2, and E3 trials did not consider any other effect beyond the selected residual in the previous step; the E4 trial was the only trial that considered the inclusion of a linear effect of rows (M3.2), which improved the model's goodness-of-fit (Table 2).

Modelling the genetic effects in multi-environment trials
After selecting the best-fitted model for each trial, a joint analysis was performed, considering the best-fitted models, indicated by BIC c in the previous steps. In the joint analysis, the genetic The variance model includes just the Natural, which means the residual effects. Number of fixed parameters (p), number of variance parameters (q), the full loglikelihoods (Full-LogL), and Conditional Bayesian information criteria (BIC c ), are given for the spatial models fitted for each trial. # : best-fitted model; s 2 R and s 2 x : variance components of units and correlated residuals; I r and I c : identity matrices of rows and columns; S r (ρ r ) and S c (ρ c ): correlation matrices for the row model (order r, autocorrelation parameter ρ r ) and column model (order c, autocorrelation parameter ρ c ).

Estimates of variance components, heritability, and coefficient of determination of the dominance effects
The LRT results for both the SPA and NSPA analyses are presented in Fig 3. We observed differences in the significance of the random effects, where the SPA and NSPA models were considered, beyond the additive by environment interaction and dominance genetic effects, the additive genetic effects (GCA) were significant. When analyzing each trial (Fig 3), the same pattern was observed for trials E3 and E4, where, after SPA, genetic effects were considered significant and were not considered by NSPA.
The SPA and NSPA models used for the joint analyses fit the compound symmetry structure ( Table 4). The NSPA model fitted a homogeneous residual structure, while SPA Number of fixed parameters (p), number of variance parameters (q), the full log-likelihoods (Full-LogL) and the conditional Bayesian Information Criteria (BIC c ), are given for the spatial models fitted for each trial. # : best-fitted model; s 2 r ; s 2 c ; s 2 R , and s 2 x : variance components of rows, columns, units, and correlated residuals, respectively; I r , I c , and I n : identity matrices of rows, columns, and units, respectively; β r and β c : liner regression coefficient for rows (R) and columns (C), respectively; S r (ρ r ) and S c (ρ c ): correlation matrices for the row model (order r, autocorrelation parameter ρ r ) and column model (order c, autocorrelation parameter ρ c ).
https://doi.org/10.1371/journal.pone.0258473.t002 considered all required spatial information for each trial in the residual structure, as well as fixed and random effects.
Differences between SPA and NSPA analyses were observed for the estimates of additive and dominance genetic variances and residual variance in each trial (Table 4). In the joint analyses, where the additive by environment interaction (ŝ 2 AE ) and dominance by environment interaction (ŝ 2 DE ) variances were estimated, the NSPA and SPA also conducted different estimates (Table 4). Trials E1, E3, and E4 presented higher dominance genetic variances than additive genetic variances (Table 4). Conversely, trial E2 presented a higher additive genetic variance ( Table 4). The joint analyses conducted to similar genetic variance estimates among NSPA and SPA, except for dominance by environment interaction variances that were 50% lower in NSPA ( Table 4). The residual structure, which was homogeneous in NSPA, became heterogeneous in SPA, and it was possible to obtain residual variances per environment. Regarding the residual variance, E3 and E4 trials presented reduced values in SPA, achieving a  reduction of 20% in the E4 trial when adopting the indicated SPA models over the NSPA model (Table 4). In individual analyses, heritability estimates varied across trials. Considering the broadsense heritabilities, trials E1 and E2 presented estimates superior to trials E3 and E4. In trial E3, SPA showed an increase of approximately 1.2% over the NSPA estimate. Broad-sense heritability presented estimates that were superior to 15% in all trials. In addition, trials E1 and E2 had the same estimates in the SPA and NSPA analyses, since they did not require any spatial information. In trials E3 and E4, the broad-sense heritability increased by 1.1% and 3%, respectively, due to the adoption of SPA instead of NSPA models.
The selective accuracies (rĜ G ; rÂ A , and, rD D ) varied from NSPA to SPA analyses in each trial. E1 and E2 presented the same values in the NSPA and SPA analyses, which were superior to 0.80 the selective accuracy of genotypic effects and selective accuracy for genetic additive effects. The selective accuracies for dominance genetic effect were 0.60 and 0.52 in the trials E1 and E2, respectively. The trials that required spatial information (E3 and E4) had increased accuracy values from the NSPA to SPA analyses. The highest difference was observed in the Table 4 trial E3, where the selective accuracy for genetic additive effects increased from 0.62 to 0.73 when considered the spatial information. The same behaviors were observed in the accuracy estimates in trials E3 and E4. The joint analyses showed improvement of the genetic parameter estimates (heritability, accuracy, and genotypic correlation) from NSPA to SPA analyses. First, the NSPA analyses presented broad-sense heritability estimates of approximately 5%. For SPA analysis, the broadsense heritability estimates were split per environment, and their values varied according to each environment, from 4% (E1 and E3) to over 7% (E4). Compared to the individual analysis, the broad-sense heritability decreased from 13% in trial E1 to 4.8% in joint analysis, but in trial E4, it increased from 6.6% to 7.5%. The genotypic correlation across environments increased from 0.71, in NSPA analysis, to 0.75, in SPA analysis.

Contrasting additive and dominance genetic effects rankings
We contrasted the additive and dominance genetic effects by NSPA and SPA analyses, in each trial, to check the differences when the spatial information was used. Considering the three best parents selected, the same selected parents were observed in trials E1 and E2 (P1, P10, and P12) (S4 Table in S1 File). In trial E3, only one parent (P6) was identified among the top three highest additive genetic values in both NSPA and SPA. In trial E4, two parents were the same, P6 and P12, in the NSPA and SPA analyses. Based on the predicted additive and dominance genetic effects, it was also possible to predict the genotypic effects of the interpopulational hybrids in each trial (S5 Table in S1 File) and the gains with selection (Fig 4) by NSPA and SPA analyses.
The predicted dominance genetic effects of the parents' combinations (S6 Table in S1 File), represented by the interpopulational hybrids (Hi,j, where i and j the parents crossed). The top ten genetic materials were compared and detached. The E3 and E4 trials presented different crossings into the top 10 dominance genetic values. The Pearson correlation coefficients were equal to 1.0 in the E1 and E2 trials, while in the E3 and E4 trials, they were higher than 0.8 (S6 Table in S1 File). The Cohen's Kappa coefficients presented the same results in the E1 and E2 trials (1.0), while in E3 and E4 trials they were around 0.60 (S6 Table in S1 File).

Discussion
The use of statistical methods in breeding programs enables the breeder to access all information available from simpler, as phenotypic data, to complex inputs, such as genomic relationships and spatial correlations [71][72][73]. The accuracy tends to increase by modeling more sources of information, leading to more efficient use of resources [74].
According to Cooper et al. [75], there are some variations among plots that the breeder cannot predict a priori, even when the soil fertility is well controlled and the experimental design is appropriate. These spatial trends contribute to a large number of non-target signals [76]. Indifference or proper accounting for exogenous sources of variation decreases the accuracy [73,77] and selection mistakes [41,78].
The occurrence of diseases and insects within and across blocks, missing plots, irrigation misdistribution, and micro soil fertility spots can affect the residual independence among plots [79]. In this manner, trials E3 and E4, which required first-order autoregressive (AR1) modeling for rows and/or columns, also presented the lowest values for broad-sense heritability, selective accuracy, and phenotypic means for the GY trait.
So and Edwards [29] found, in a simulated balanced dataset study, similar results compared compound symmetry and unstructured models. This statement is relevant because of the genotypic correlation estimate across environments (0.75), by SPA analysis, where the compound symmetry structure, even being a more conservative structure [36], is capable of dealing with the level of G×E interaction effect without loss of predictive capability.
When comparing the NSPA and SPA models, the results demonstrated the importance of spatial analysis. Relevant studies on potato [80] and soybean [44] also presented satisfactory results when spatial trends were considered. Higher selective accuracies were encountered in the SPA models, reinforcing the importance and potential of the SPA analyses in crop breeding. Another simulation study of soybean breeding considering the accounting of spatial information in the model also presented gains in selective accuracy when autoregressive structures were used for spatial information [76].
The significance of the random effects tested by the LRT differed when spatial information was considered in both individual and joint analyses. The selective accuracies and genetic gains demonstrate the capacity of data correction, when necessary, by spatial analyses. Similar results were reported by Resende [11], who affirmed that the presence of spatial dependence is required when there is any kind of heterogeneity in the experiment, or even inside the block. Furthermore, by the adoption of SPA analyses, when applicable, it was possible to estimate higher values of genetic variances and parameters, which proves the importance of considering the best-fitted model. For instance, in trial E3, after modeling the residuals and including spatial information in the model, it was possible to obtain a higher portion of the additive and dominance genetic effects, turning them into significant effects. They were non-significant in the NSPA analyses and became significant in the SPA analyses. In addition, the genetic variances, heritabilities, and selective accuracies were maximized using SPA analyses. The genotype ranking also varied between the NSPA and SPA analyses. This information can be used in breeding programs for genetic selection purposes.
Selective accuracy demonstrates the correlation between the predicted and true genetic values. As mentioned previously, there was an increase from NSPA to SPA analyses. These results reinforce the superiority of SPA analyses, which maximize the selective accuracy when compared to NSPA analyses. Residual modeling considering the dependence among rows and/or columns decreased the residual variance within the trials. The consideration of two errors, dependent and independent, allows modeling of the dependence part and mitigates the residue in explained portions. This error investigation maximized the selective accuracies in the environments, which indicated the spatial models as the best-fitted model because of the possibility of modeling more environmental effects that were not considered in the NSPA analyses.
One of the main goals of plant breeding programs is to select the best genotypes to increase genetic gains [81]. In addition, to maximize the frequency of favorable alleles for a specific trait, plant breeders use the best-ranked genotypes in the selection process to develop hybrids [81]. Following this idea, the selection gains based on the predicted genotypic effects by the BLUP method confirm the corrective capacity of SPA analyses over the NSPA analyses. The predicted genetic values can be used to predict the breeding value of individuals [82]. The joint analysis allows the calculation of the additive genetic effect and its interaction with the environment when considering all four environments. This fact increases the accuracy of the selection gain, as it is based on the additive genetic effect free of the GxE interaction effect. Herein, we confirmed the corrective capacity of spatial analyses when considering the selection gains, as well the genotype ranking, which led to the selection of the right genotypes in the SPA analyses ranking.
Furthermore, the differences between SPA and NSPA analyses rankings, in trials E3 and E4 (which spatial information was used) reinforce the importance of this analysis. It was observed that the spatial information was useful for selecting different parents and interpopulational hybrids using SPA over NSPA analyses. In addition, using the Pearson correlation coefficient demonstrates the existence of different genetic materials when contrasting NSPA and SPA analyses, and the Cohen's Kappa coefficient, which focuses on the selected genetic materials, confirms the difference among the selected.
By observing the differences between the additive and dominance genetic effects, some interesting conclusions can be drawn for use in maize breeding programs. For instance, in trial E1, the H10,12, from the best parents (by additive effects) also presented a high dominance genetic effect; however, in trial E2, this effect was negative, demonstrating the presence of dominance by environment interaction effect, and its importance to calculate it. P6, which showed one of the lowest additive genetic effects, when combined with P12 (one of the highest additive genetic effects) presented one of the highest dominance genetic effects. These characteristics were observed in all the trials. In trial E3, a significant difference was noted between variance estimates and genetic parameter values from NSPA and SPA analyses, which demonstrates the importance of the SPA analyses in this study. Focusing on P3 and P12 as an example, the P3 was into three better parents in NSPA analyses, and decreased significantly in SPA analyses, while P12 arose from a median rank (in NSPA) to the third better in SPA analyses. Following, when looking into the dominance genetic effects just of the same two parents cross (H3,12), it was in the top 10 in the NSPA analyses, and dropped to the 16 th position in SPA analyses. As demonstrated, the spatial trend corrections affected the additive and dominance genetic effects and must be considered by the breeder in breeding program strategies. The selection will be made based on these genotypic, additive genetic, and dominance genetic effects values, and their accuracy will directly affect the subsequent cycles of the breeding process.
Identifying the best genetic material is the main goal of any breeding program. The main goal of this study was to identify the best potential parents and their combinations to form heterotic groups. The interpretation individually and jointly brings the parents 2, 10 and 12, as potential parents in this microregion composted by the four environments. Overall, they presented a high additive genetic effect (meaning they have a high proportion of favorable alleles to GY trait) and a high dominance genetic effect (seem to be contrasting and complementary), which would be very interesting in a hybrid composition, as the heterosis effect is expected in maize.
In this way, this study highlighted the importance of considering the dependence among rows and columns when required for specific trials and applied all additional information into the joint analysis, combined with the best genetic covariance structure. The adoption of SPA analyses led to better results, as observed in BICc, selective accuracy, heritability, and selection gains. This study confirmed the importance of SPA diallel MET analyses in maize breeding, which is a pioneer in the literature.

Conclusion
The adoption of MET permits a better understanding of the genetic effect, as it enables the breeder to account for the GxE interaction effect. This makes it possible to isolate the genetic effect or further explore the GxE interaction effect, depending on the goal. In addition, the diallel matching design enriches the study as it accesses the effects within the genotypic effect, both additive and dominance genetic effects, by crossing all parents among themselves. We found that the SPA analyses maximize the selective accuracy when compared to NSPA analyses, confirming the superiority of the SPA models. In addition, the LRT confirmed that the mismodeling process can lead to the misinterpretation of model effects, where some significant effects can be considered non-significant due to inappropriate analyses.
Based on the fitted information, it was possible to identify the genotypes with the highest potential for future heterotic groups that can be used as potential commercial hybrids. The parents 2, 10 and 12, were indicated as potential parents in this microregion. Owing to their high additive genetic effects, bringing favorable genes to the future parental lines, their positive combination, given by the dominance genetic effect, will be expressed as a strong heterosis in future hybrids.