INTRODUCTION

Soybean (Glycine max [L.] Merr.) is grown in an extensive cultivation area in Brazil, and is the main commodity of the national agricultural sector. World estimated production for the 2016-2017 crop season was approximately 336 million tons, and Brazil was responsible for over 107 million tons (^{Conab, 2017}).

The search for genotypes with high grain yield and stability is the key aspect of plant breeding programs, which aim to obtain genotypes with high performance for different cultivation regions. The Genotype × Environment (GE) interaction strongly influences genotype performance, making new cultivar selection difficult. One way to take advantage of this interaction is to identify genotypes with high production capacity and phenotypic stability in different environments (^{Bornhofen et al., 2017}).

Identifying genotypes with high grain yield and stability is the main goal of plant breeding programs. The literature describes numerous methods to study the stability and productive performance of genotypes. These methods differ in the evaluated parameters and statistics. Studies comparing different methods have already been described for different crops, such as wheat (^{Bornhofen et al., 2017}), maize (^{Cargnelutti Filho et al., 2009}), common beans (^{Pereira et al., 2009}), and soybean (^{Polizel et al., 2013}).

Among the main methods used to study phenotypic stability of genotypes are univariate parametric models (^{Wricke, 1965}; ^{Eberhart and Russell, 1966}), multivariate parametric models (^{Zobel et al., 1988}; ^{Yan, 2001}), non-parametric models (^{Lin and Binns, 1988}), and mixed models (^{Resende, 2006}). The ability to explain the sum of squares of the GE interaction is the main difference among methods (^{Bornhofen et al., 2017}).

Phenotypic stability can be divided into two principal types, that is, stability in the biological and agronomic sense. Stability in the biological (static) sense refers to the ability of genotypes to maintain constant production in different environments, with low variation between them, that is, genotypes exhibit "homeostasis". Stability in the agronomic (dynamic) sense indicates that the genotype positively responds to improvements in edaphoclimatic conditions of the environment and can perform above the mean in different locations (^{Sabaghnia et al., 2015}). This behavior is of interest to plant breeders and farmers. It is important for the plant breeder to adopt methods in which genotype stability is associated with high grain yield average. Among the different currently available methods, the genotype main effect + genotype x environment (GGE) biplot enables simultaneous analysis of genotypes by considering high grain yield average and stability. The combination of these two concepts is known as the "ideal genotype" (^{Yan et al., 2007}; ^{Yan, 2016}).

To compare the results obtained between different methods, many authors have used Spearman's correlation. This comparison methodology is based on ranks. However, it is not known if the methods provide a consistent result regarding the correlation between the different years under study. The new contribution of the present study is the evaluation of the consistency of Spearman's rank correlation between different years using multivariate methods, such as genotype main effect + genotype × environment (GGE) interaction effect, additive main effects and multiplicative interaction (AMMI1), and mixed models. Therefore, the objective of this study was to evaluate the consistency of association between grain yield and stability methods in multi-environment trials and identify which of the evaluated methods is more suitable for selecting superior genotypes.

MATERIALS AND METHODS

Soybean grain yield (GY) data were obtained from trials of Value for Cultivation and Use (VCU) at 10 locations in the 2013-2014, 2014-2015, and 2015-2016 crop seasons. The locations are in the soybean macro-region 2, situated in the states of Paraná and Mato Grosso do Sul, Brazil. In Paraná, trials were conducted at Bela Vista do Paraíso (23°00' S, 51°11' W), Cafelándia (24°37' S, 53°19' W), Floresta (23°36' S, 52°05' W), Palotina (24°17' S, 53°50' W), Rolándia (23°19' S, 51°22' W), Santa Terezinha de Itaipu (25°27' S, 54°24' W), and Sertanópolis (23°04' S, 51°02' W). In Mato Grosso do Sul, trials were conducted at Dourados (22°13' S, 54°49' W), Maracaju (21°37' S, 55°10' W), and Naviraí (23°04' S, 54°11' W).

Eleven soybean cultivars were evaluated (coded as C1 to C11): 5958RSF IPRO (C1), 6563RSF IPRO (C2), NA 5909 RG (C3), BMX Potencia RR (C4), DM 5.9i (C5), NK 7059 RR (C6), 6160RSF IPRO (C7), 6458RSF IPRO (C8), 68I70RSF IPRO (C9), 61I59RSF IPRO (C10), and 63I64RSF IPRO (C11). Trials were carried out in a randomized complete block design with three replicates. Each experimental plot consisted of four rows, 5 m long, and 0.5 m row spacing. The two central rows were mechanically harvested from each plot (5 m²). Sowing density was standardized at 350 000 seeds ha^-1 for all cultivars. Crop management followed the technical recommendations for soybean cultivation. Grain moisture was adjusted to 13% and expressed in kg ha^-1.

Data were tested for normality according to the Lilliefors test. Joint and individual ANOVA were performed each year. The significance of genotype (G), environment (E), and G×E (GE) interaction effects were determined by the F test. Selective accuracy was calculated according to the method described by ^{Resende and Duarte (2007)}. Genotype performance for the group of environments was evaluated each year by the Scott-Knott method at 5% error probability. The method described by ^{Annicchiarico (1992)} was used to divide the environments as favorable and unfavorable in 2013, 2014, and 2015. After this division, ANOVA was performed for the environments classified as favorable and unfavorable. A new group of means was then performed by the Scott-Knott method. This procedure aimed to identify differences in cultivar grouping for the group of environments and favorable and unfavorable environments.

In order to evaluate stability, the following methods were tested: ^{Wricke (1965)} (coi), ^{Eberhart and Russel (1966)} (E-R), Lin and Binns modified by ^{Carneiro (1998)} (L-B/C), additive main effects and multiplicative interaction (AMMI1) (^{Zobel et al., 1988}), and the harmonic mean of the genotypic values (HMGV) across environments via mixed models by restricted maximum likelihood (REML)/best linear unbiased prediction (BLUP), model 54 (^{Resende, 2006}). In addition, for multivariate analysis via GGE, the biplot graph and output (Log) values of the analyses were used to identify the most stable genotypes. The GGE biplot (^{Yan, 2001}), Selegen (^{Resende, 2006}), and Genes (^{Cruz, 2013}) software were used to perform the analyses.

Spearman's rank correlation coefficient was used to compare the described methodologies. Cultivars were ranked by considering stability and productive performance. Thus, the genotypes were first ordered by GY based on mean performance for the complete set of environments; the best genotype performance was allotted number 1 and the worst number 11.

For the ci parameter, genotype ranking was carried out by the ecovalence index. In this method, when the modular ecovalence of the genotype is lower, stability is higher. The genotype that exhibited the lowest modular ecovalence value was allotted number 1 and the genotype with the highest value was number 11 in the ranking.

In the methodology described by E-R, genotypes were arranged according to the deviation of the regression. Genotypes that showed no significant deviation of the regression (S²d = 0) were classified as stable. In this group, the genotype with the lowest deviation of the regression in the module was ranked as 1 and up to the g^th deviation of the regression. Then, genotypes with significant deviation of the regression at 5% and 1% were ranked following the same procedure.

The L-B/C method is based on the P_i parameter where the most stable genotype has the lowest P_i value. The modification proposed by ^{Carneiro (1998)} enables the identification of stable genotypes in favorable, unfavorable, and general environments. Ranking genotypes by P_i (g), P_i (f), and P_i (u) parameters was carried out for general, favorable, and unfavorable environments, respectively. The genotype with the lowest P_i value was classified as the most stable and ranked as 1. The process was repeated for the 11 genotypes and three parameters.

The HMGV parameter was used for classification by the mixed model method because it estimates stability. In the multivariate methods, AMMI1 and GGE for stability (GGE STA), genotype stability was evaluated by the graphs and log values of the GGE biplot software. In these methods, vector projection in relation to the AEC-axis, regardless of the direction, represents stability, and when the vector is smaller, stability is higher. In the ideal genotype analysis (GGE IG), a value of 1.5 was used as weight for mean GY in relation to stability.

Spearman's rank correlation coefficient was used to compare the methods to measure stability. Consistency of the methods between years was evaluated by the significance of the associations and the coexistence of similar correlations in the 3 years. When results were similar in the correlations for the 3 years, the correlated methods were accepted as consistent, which were significant or nonsignificant by the t test for n = 11. The magnitude of the coefficients of Spearman's correlation was according to the description by ^{Monteiro et al. (2015)} as follows: null magnitude (0.0 to 0.20), weak magnitude (0.21 to 0.40), average magnitude (0.41 to 0.70), strong magnitude (0.71 to 0.90), and extremely strong magnitude (0.91 to 1.00).

RESULTS AND DISCUSSION

The joint ANOVA for locations in each year indicated that the sources of variation of G, E, and GEI were significant (p < 0.01) for all the evaluated years (Table 1). This indicates that genotypes exhibited different behavior in the different environments, and that there was a change in the magnitude of the responses for GY due to environmental variation. Similar results were observed by ^{Cargnelutti Filho et al. (2007)}. The GY mean of the trials was 3529.34 kg ha^-1 in 2013, 3333.03 kg ha^-1 in 2014, and 3486.42 kg ha^-1 in 2015. The coefficient of variation (CV) was 7.96%, 9.74%, and 10.23%, and selective accuracy was 0.70, 0.86, and 0.87 in 2013, 2014, and 2015, respectively. These results indicate high experimental precision, providing security and reliability for selecting superior genotypes.

The GY means for the cultivars in the 3 years and 10 locations are shown in Table 2. The 11 cultivars were present in the 10 locations for the 3 years, and the inferences on stability were made for the same cultivars and locations in the different years. The difference between the cultivars was nonsignificant for the 10 locations as a group and for each of the three evaluated years. This was probably due to the complex interaction between cultivars and the experimental locations that changed the genotype ranking between locations evaluated within the same year (^{Costa et al., 2015}). The genetic correlation between performances in the evaluated environments (rgloc) was 0.10, 0.26, and 0.32 in 2013, 2014, and 2015, respectively (data not shown), indicating extremely complex GE interaction. Values of rgloc that are close to 0.50 indicate complex GE interaction with alterations in the ranking of cultivars between environments (^{Carvalho et al., 2016}). This shows the need for statistical methods that enable the identification of superior genotypes with high GY and stability.

Table 1 Joint analysis of variance (ANOVA) for soybean grain yield (GY) of 11 cultivars at 10 locations in 2013, 2014, and 2015. 

The absence of significant difference between genotypes for the group of 10 locations indicates the need to discriminate environments. The method by ^{Annicchiarico (1992)} allows classifying environments as favorable and unfavorable based on grain yield at each location (Table 3). After designating environments as favorable and unfavorable, a significant difference was observed between years (Table 2). Among the environments classified as favorable, differentiation between cultivars was found only in 2013 and cultivars C2, C7, C8, and C10 showed superior performance. Regarding unfavorable environments, both the 2013 and 2014 seasons exhibited significant differences between cultivars. In 2013, cultivars C2, C4, C6, C8, C9, C10, and C11 were superior, while cultivar C9 had superior performance in 2014. In 2015, classifying environments as favorable and unfavorable was ineffective to identify superior cultivars for GY.

The comparison of statistical methods is shown in Table 4. Spearman's correlation coefficient varied from -0.89 to 1.00. This result indicates different consistency categories between the parameters of cultivar classification. Significant positive associations were observed between GY and stability by the L-B/C methodology (Table 4). The magnitude of the associations of the L-B/C method with GY for general environments (P_i) and unfavorable environments (P_i(u)) were 0.65 and 0.83, 0.91 and 0.77, and 0.99 and 0.92 in 2013, 2014, and 2015, respectively. For the favorable environments (P_i(f)), the correlation with GY was significant in 2014 (0.85) and 2015 (0.97). In 2013, the correlation between GY and P_i(f) was nonsignificant (0.27). Therefore, the P_i(f) parameter was not consistent between years. The evaluation of productive stability by non-parametric methods, such as the ^{Lin and Binns (1988)} method, is advantageous since it combines GY and stability (^{Zali et al., 2011}). However, the lack of consistency between years can generate unreliable information and is associated with GY only in some years.

Table 2 Ranking of 11 soybean cultivars evaluated at 10 locations during 3 years (2013, 2014, and 2015), as related to grain yield (GY) in all environments, in favorable (F) and unfavorable (U) environments, and parameters of stability methods. 

¹Means followed by the same letter in the same year do not differ according to the Scott-Knott test (p = 0.05).

C: Cultivar; GY: grain yield; Parameters to evaluate adaptability and stability: mi: ecovalence (^{Wricke, 1965}); S²d: ^{Eberhart and Russell (1966)}; P_i, P_i(_f) and P_i(_u): ^{Lin and Binns (1988}) method modified by ^{Carneiro (1998)}; HMGV: stability by mixed models (REML/BLUP); AMMI1: additive main effects and multiplicative interaction, first principal component of the AMMI1 analysis; GGE: genotype main effect + genotype x environment interaction, stability by GGE biplot analysis; i, (f), and (u): performance in general, favorable environments, and unfavorable environments, respectively.

Cultivars: C1: 5958RSF IPRO, C2: 6563RSF IPRO, C3: NA 5909 RG, C4: BMX Potencia RR, C5: DM 5.9i, C6: NK 7059 RR, C7: 6160RSF IPRO, C8: 6458RSF IPRO, C9: 68I70RSF IPRO, C10: 61I59RSF IPRO, and C11: 63I64RSF IPRO.

Table 3 Classification (Class) of locations as favorable (F) or unfavorable (U) according to method by ^{Annicchiarico (1992}) for soybean grain yield in 2013, 2014, and 2015. 

Table 4 Repeatability of Spearman's correlation coefficients between evaluation methods of grain yield and stability concepts using rankings of 11 soybean cultivars evaluated in 2013, 2014, and 2015. 

*, **Significant at the 0.05 and 0.01 probability levels according to Student's t-test, respectively; ^ns nonsignificant according to Student's t-test at 0.05 probability level. GY: grain yield; Stability parameters: mi (^{Wricke, 1965}); S²d (^{Eberhart and Russell, 1966}); P_i, P_i(_f), and P_i(_u) for stability of environments in general, favorable environments, or unfavorable environments according to the method by Lin and Binns modified by ^{Carneiro (1998)}; additive main effects and multiplicative interaction - AMMI1 (^{Zobel et al., 1988}); additive main effect of genotype plus multiplicative interaction of genotype main effect + genotype x environment interaction (GGE) for stability (GGE STA) and ideal genotype (GGE IG) (^{Yan, 2001}); mixed models - REML/BLUP, model 54 (^{Resende, 2006}); HMGV: stability obtained by harmonic mean of the genotypic values across environments.

Correlation coefficients between GY and the HMGV parameter were 0.88, 0.97, and 1.00 in 2013, 2014, and 2015, respectively. The observed correlations in all 3 years indicate that these methods enable the identification of stable genotypes with high GY, that is, they are associated with the dynamic concept of stability (^{Annicchiarico, 2002}). The HMGV method classifies genotypes by their estimated genetic values and also by stability. Therefore, this method is different because it evaluates genotypes in terms of estimated genetic values by best linear unbiased prediction (BLUP) and restricted maximum likelihood (REML) (^{Resende and Duarte, 2007}). On the other hand, GY was not positively associated with the E-R, mi, AMMI1, and GGE STA methodologies in any year. These results indicate that genotypes with high stability are not always the most productive; results concur with those reported by ^{Bornhofen et al. (2017}) when evaluating wheat genotypes.

Significant positive associations were found between the mi and AMMI1 (0.69, 0.67, and 0.68) and E-R (0.95, 0.95, and 0.99) methodologies in 2013, 2014, and 2015, respectively. Other authors have also found positive associations between these methods (^{Paula et al., 2014}; ^{Bornhofen et al., 2017}). This occurs because these methods are associated with the concept of static stability, which is obtained independently of genotype mean yield (^{Annicchiarico, 2002}). The selection exclusively made by these methods can result in selecting stable but not necessarily productive genotypes.

The E-R (S²d) methodology showed a significant association with the P_i(f) parameter (0.66) from the L-B/C method only in 2013. However, ^{Dehghani et al. (2008)} report a negative association (-0.55) between the L-B/C and E-R methods. In addition, an association was observed between L-B/C and the mi and AMMI1 methodologies only in 2013, and the association with the GGE STA parameter was absent. There is therefore no consistency between years, indicating differences between methods because they are based on different concepts of stability. Representatives of both groups should be used simultaneously in practice to consider the different concepts of stability and enable more reliable selection and recommendation of genotypes.

Significant positive associations between different methods indicate that genotypes are similarly classified as to stability. Methods may be supplying redundant information, and only one would be sufficient to select the best genotypes (^{Cargnelutti Filho et al., 2009}). However, even with strong association between methods, genotype ranking can be different between methodologies. This is the case, for example, of the association between GY and HMGV in 2013, that is, the correlation is high (0.88) but the methods differ as to the best cultivar; cultivar C2 was superior in HMGV, while cultivar C10 was superior for GY. In the HMGV method, cultivar C10 was only the fourth best genotype. Moreover, the HMGV method shows genotypic variation, whereas GY shows phenotypic variation. It is therefore necessary to use two or more methodologies to ensure more efficient selection.

The methods based on multivariate statistics make it possible to visually evaluate genotype behavior and more clearly identify genotypes with characteristics of interest. The GGE method consists in analyzing biplot graphs to evaluate genotype stability and GY (Figures 1 A, B, and C). The most stable cultivar for 2013 was C2, followed by C11, C6, C8, C3, and C1. In 2014, cultivars C5 and C3 were the most stable. The most stable cultivars in 2015 were C8, C3, and C1. However, the stability analysis of this method did not show a significant association with GY in all 3 years. Considering only the stability parameter, the absence of association of the GGE STA method in 2013, 2014, and 2015, with associations of 0.27, -0.49, and -0.01, respectively, indicates that this method is consistent, but not correlated, with GY. Only using stability results is insufficient for selecting genotypes with superior GY performance. However, the analysis of mean versus stability on the GGE biplot allows selecting superior genotypes in terms of both yield and stability (^{Yan et al., 2007}). When the genotype is more to the right of the biplot, it is more productive, whereas it is more stable when it is closer to the single arrow x-axis.

The GGE biplot software (^{Yan, 2001}) allows identifying the best genotypes with "ideal genotype" analysis (GGE IG) by GY and stability statistics (Figures 1 D, E, and F). In this analysis, it was possible to attribute different weights to genotype GY in relation to stability, and this value varied from 1 to 4. Lower values indicated higher weight for GY, that is, productive performance has to be high for the position to be at the extreme right of the biplot. On the other hand, higher values (e.g., 4) allow easily identifying an ideal genotype because there is a greater demand for genotype stability. In this study, the genotype identified as ideal was the same as the one that had the most productive performance for the evaluated period because 1.5 GY/stability weight was used. Weights close to 4 can alter the identification of the ideal genotype, especially when highly productive genotypes are very unstable.

The GGE IG method was efficient for selecting stable and productive genotypes, and showed a significant association with GY in 2013, 2014, and 2015 (0.73, 0.90, and 0.97, respectively). These results validate the efficiency of the method for selecting superior genotypes and confirm its consistency between years.

The AMMI1 method (Figure 2) also simultaneously provides productive performance and stability in identifying superior cultivars. This method divides the main effects of genotype and environment and the effects of the GE interaction. It also combines ANOVA for the main effects of genotype and environment, represented on the x-axis, and the effect of the interaction of the principal components on the y-axis (^{Bose et al., 2014}). Cultivars C10, C2, C8, and C11 were the most productive in 2013 and cultivar C2 was able to combine GY and stability among the evaluated environments. In 2014, cultivar C9 had superior productive performance, although it was the most unstable among the evaluated genotypes. In the same year, cultivar C2 exhibited high GY and stability. Cultivar C11 was the most productive in 2015, and C2 was the most stable, in 2013 and 2014. The AMMI1 method had no significant association with GGE IG in any of the years under study, which confirms the differences between the methods and the possibility of using them together as complementary analyses.

The use of stability and GY estimation methods are fundamental to identify highly productive and stable genotypes. Several parameters were positively correlated one with the other for the different years, thus indicating consistency in the results and helping to suggest more accurate methods for GY and stability analysis of the soybean crop. Due to the different concepts of stability, the breeder should prioritize using a group of methods that are related or unrelated to GY. Furthermore, the methods can be differentiated by using phenotypic data, or by estimating genetic values. The group of methods to be used can therefore consist of the HMGV and GGE IG methods since they are consistently associated with GY in different evaluated years. In this case, the HMGV method is emphasized because it uses genetic value prediction and the concept of dynamic stability, while GGE IG is pointed out for evaluating stability and GY together. Moreover, genotype selection by one of the mi, AMMI1, or S²d parameters is recommended because they exhibit positive and consistent associations with each other in the different years under study. These methods use the concept of static stability (not associated with yield). Among the different studied methods, the GGE IG method is highlighted. Despite its stability parameter to use the static concept, it allows verifying GY and the concept of ideal genotype, as well as the ease provided by graphical analysis.

The evaluated cultivars were Cl: 5958RSF IPRO, C2: 6563RSF IPRO, C3: NA 5909 RG, C4: BMX Potencia RR, C5: DM 5.9i, C6: NK 7059 RR, C7: 6160RSF IPRO, C8: 6458RSF IPRO, C9: 68I70RSF IPRO, C10: 61I59RSF IPRO, and C11: 63I64RSF IPRO.

Figure 1 Mean, stability, and ideal genotype for the 11 soybean cultivars evaluated in 2013, 2014, and 2015 at 10 locations. A, B, and C correspond to analysis of mean and stability in 2013, 2014, and 2015, respectively, and D, E, and F correspond to analysis of ideal genotype in 2013, 2014, and 2015, respectively. 

The evaluated cultivars were C1: 5958RSF IPRO, C2: 6563RSF IPRO, C3: NA 5909 RG, C4: BMX Potencia RR, C5: DM 5.9i, C6: NK 7059 RR, C7: 6160RSF IPRO, C8: 6458RSF IPRO, C9: 68I70RSF IPRO, C10: 61I59RSF IPRO, and C11: 63I64RSF IPRO.

Figure 2 Additive main effects and multiplicative interaction 1 (AMMI1) analysis for the 11 soybean cultivars evaluated in 2013 (A), 2014 (B), and 2015 (C) at 10 locations. 

CONCLUSIONS

The evaluated methodologies show different results regarding the consistency of associations between grain yield and stability evaluation methods across years.

The Lin and Binns method modified by Carneiro (L-B/C) exhibits positive associations with the Eberhart and Russel (E-R), Wricke (mi), and additive main effects and multiplicative interaction (AMMI1) methodologies, as well as with grain yield (GY), but it is not consistent across years and generates unreliable information.

The harmonic mean of the genotypic values (HMGV) across environments parameter and the genotype main effect + genotype x environment interaction effect by ideal genotype (GGE IG) method has consistent results and strong or extremely strong associations with GY; it is the most adequate method to select productive and stable genotypes.

The mi methodology consistently and significantly correlates with AMMI1 and E-R, and it is recommended to prioritize one of these methodologies to avoid redundant results.

Genotype selection by the GGE biplot multivariate methodologies is valid when considering the concepts of mean, stability, and ideal genotype, not only the concept of stability in the biological sense without the interference of yield.