Multi-trait multi-environment models in the genetic selection of segregating soybean progeny

At present, single-trait best linear unbiased prediction (BLUP) is the standard method for genetic selection in soybean. However, when genetic selection is performed based on two or more genetically correlated traits and these are analyzed individually, selection bias may arise. Under these conditions, considering the correlation structure between the evaluated traits may provide more-accurate genetic estimates for the evaluated parameters, even under environmental influences. The present study was thus developed to examine the efficiency and applicability of multi-trait multi-environment (MTME) models by the residual maximum likelihood (REML/BLUP) and Bayesian approaches in the genetic selection of segregating soybean progeny. The study involved data pertaining to 203 soybean F2:4 progeny assessed in two environments for the following traits: number of days to maturity (DM), 100-seed weight (SW), and average seed yield per plot (SY). Variance components and genetic and non-genetic parameters were estimated via the REML/BLUP and Bayesian methods. The variance components estimated and the breeding values and genetic gains predicted with selection through the Bayesian procedure were similar to those obtained by REML/BLUP. The frequentist and Bayesian MTME models provided higher estimates of broad-sense heritability per plot (or heritability of total effects of progeny; hprog2) and mean accuracy of progeny than their respective single-trait versions. Bayesian analysis provided the credibility intervals for the estimates of hprog2. Therefore, MTME led to greater predicted gains from selection. On this basis, this procedure can be efficiently applied in the genetic selection of segregating soybean progeny.

soybean crop management classification [27], aiming to exploit genetic variability for the selection of productive progeny. Seventy-two, 71, and 60 F 2 plants of populations 1, 2, and 3, respectively, were separately bulk-harvested and threshed. One sample was collected from each F 2:3 progeny to compose the 203 F 2:4 progeny that were used in this study by the withinprogeny bulk method [28,29].
The experiments were set up as a randomized complete block design with three replicates per site. Plots consisted of two 3.0-m rows spaced 0.5 m apart, with a plant stand density of 13 seeds per meter, totaling a density of 256,000 plants ha -1 . All plant management operations were undertaken in accordance with the requirements of the crop in the region [30].
Three target agronomic traits were evaluated, namely: number of days to maturity (DM, days), 100-seed weight (SW, g), and average seed yield per plot (SY, g). The first variable, DM, corresponds to the number of days before 95% of the pods were mature, as indicated by their color [31]. To evaluate the SW and SY traits, the grains were dried to 13% moisture. All data and code to run the main models used in this study are available in S1 Table and S2 Table, respectively.

Frequentist statistical analyses
The Restricted Maximum Likelihood/Best Linear Unbiased Prediction (REML/BLUP) procedure was adopted for statistical analyses under a frequentist approach (Patterson and Thompson [32] and Henderson [33]). The frequentist single-trait multi-environment (FSTME) statistical model associated with the evaluation of segregating progeny in a randomized block design in two environments, with one observation per plot, is given by: where y is the vector of phenotypes; b is the vector of block effects added to the overall mean (assumed fixed); g is the vector of progeny effects (assumed random), in which g � Nð0; s 2 g Þ; i is the vector of the G×E interaction effects (random), in which i � Nð0; s 2 int Þ; and e is the vector of residuals (random), in which e � Nð0; s 2 e Þ. The capital letters (X, Z and W) represent the incidence matrices for the effects of b, g, and i, respectively. The b vector encompasses all replicates of all locations. In this case, this vector comprises the effects of locations and of replicates within locations. The goodness of fit was obtained using Akaike information criteria (AIC [34]), defined by AIC = −2LogL+2p, in which LogL is the log-likelihood function and p is the number of estimated parameters, according to Little et al. [35]; and the Likelihood Ratio Test (LRT), following Wilks [36], using Chi-square Statistics with one degree of freedom, which is calculated by the following equation: λ = 2[Log e L p+1 −Log e L p ], in which L p+1 and L p are the maximum likelihood associated with the full and the reduced models, respectively.
In the frequentist multi-trait multi-environment (FMTME) approach, this model was transformed into g~N(0,∑ g �I), i~N(0,∑ int �I) and e~N(0,∑ e �I), where ∑ g is the progeny covariance matrix; ∑ int is the G×E interaction covariance matrix; ∑ e is the residual covariance matrix; I is an identity matrix with order appropriate to the respective random effect; and � denotes the Kronecker product. ∑ g , ∑ int , and ∑ e are also unstructured covariance structures (US) [18,19]. The following Y 1 , Y 2 , and Y 3 are the vectors of observed responses for each trait To perform the statistical analyses of the FSTME and FMTME models and obtain the variance components and breeding values, we applied the ASReml 4.1 package [21] of integrated R software (Development Core Team- [37]).

Bayesian statistical analyses
The single-trait multi-environment and multi-trait multi-environment models as well as the covariance matrix structures (US) of the frequentist analyses were used according to the Bayesian approach (BSTME and BMTME, respectively) via Markov Chain Monte Carlo (MCMC) to estimate the variance components and genetic parameter. To gather further information about representation of model, matrix, and vectors structures as well as priori probability distributions, we recommend reading Mrode [18] and Gilmour et al. [21] and recent publications Junqueira et al. [23], Torres et al. [25] and Mora and Serra [38].
We assumed that the variance-covariance matrices follow an Inverse-Wishart (IW) distribution and independent Inverse-Gamma (IG) distributions, which were used as a priori to model the variance-covariance matrix [39][40][41]. The ensuing covariance matrix distribution is such that all standard-deviation parameters have Half-t distributions and the correlation parameters have uniform distributions on (-1,1) for a particular choice of the IW shape parameter [13]. The advantage of this approach is that it allows us to choose the shape and scale parameters that achieve high non-arbitrary information of all standard deviations and correlation parameters [40].
The full (considering the genotype and G × E interaction effects) models were compared with the reduced (disregarding the genotype or G × E interaction effects) models by the deviance information criterion (DIC) proposed by Spiegelhalter et al. [42]: yÞ is a point estimate of the deviance obtained by replacing the parameters with their posterior mean estimates in the likelihood function and p D is the effective number of parameters in the model. Models with lower DIC should be preferred over models with higher DIC.
The Bayesian models (BSTME and BMTME) were implemented in the "MCMCglmm" package [43,44] of R software (R Development Core Team- [37]). A total of 1,000,000 samples were generated, and assuming a burn-in period and sampling interval of 500,000 and 5 iterations, respectively, this resulted in a final total of 100,000 samples. The convergence of MCMC was checked by the criterion of Geweke [45], which was performed using the "boa" [46] and "CODA" (Convergence Diagnosis and Output Analysis) [47] R packages.

Genetic and non-genetic components
Broad-sense heritability per plot (or heritability of total effects from progeny) for the frequentist and Bayesian models were computed based on the approximated estimators, as discussed in Piepho et al. [48], using the following expression: where s 2 g : variance of progeny; s 2 int : variance of the progeny × environment interaction; s 2 e : error variance; n: number of locations; and r: number of replicates. For the Bayesian models, the posterior estimates were calculated from the posterior samples of the variance components obtained by the model.
The accuracy of progeny selection for the frequentist models (FSTME and FMTME) was estimated based on the following expression [49]: 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 where s 2 g is the genotypic variance and PEV is the prediction error variance extracted from the diagonal of the generalized inverse of the coefficient matrix of the mixed-model equations.
For the Bayesian approach (BSTME and BMTME models), the accuracy of progeny selection was estimated according to Resende et al. [19] from the posterior distribution, given by:rĝ where sðg Þ is the standard deviation of the predicted breeding value. According to Resende et al. [19], in Bayesian inference, the variance of the very parameter that is assumed as a random variable is computed.
The "boa" [46] and "bayesplot" [50,51] packages of R software were used to calculate and plot the highest posterior density (HPD) intervals for all parameters, respectively. Estimates of the coefficient of experimental variation (CVe) and selective accuracy (rĝ g ) were used to evaluate the experimental quality of the models [52].

Genetic correlation
To determine the genetic covariance by the frequentists and Bayesian single-trait models (FSTME and BSTME, respectively), a pairwise analysis of the sum of phenotypic values of the traits was performed. Thus, the covariances were obtained, as proposed by Resende et al. [53], using the following expression: where s 2 g ðtrait i þtrait j Þ is the variance of the sum of phenotypic values of traits i and j; s 2 gðtrait i Þ is the genotypic variance of trait i; and s 2 gðtrait j Þ is the genotypic variance of trait j. Genetic covariances by the multi-trait multi-environment models (FMTME and BMTME models) were obtained directly by the mixed-model output from each applied methodology.
The genetic correlation coefficients between the DM, SW, and SY traits were obtained, as suggested by Piepho et al. [54], using the expression below for all models: r ðtrait i ;trait j Þ ¼ s gðtrait i ;trait j Þ 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 s 2 gðtrait i Þ s 2

Progeny selection
The Spearman rank correlation coefficient was calculated between the BLUP (breeding value) of the progeny from the analyses of FSTME, FMTME, BSTME, and BMTME, and its significance was verified using nonparametric bootstrap in the R package 'boot' [55,56]. The agreement between the selected progeny was also checked by the coincidence index (CI) proposed by Hamblin and Zimmermann [57], as shown below: where A is the number of coincident progeny in the two methods, M is the number of selected progeny, and C is the number of progeny coincident due to chance (C = bM, where b is the selection intensity = 0.15; i.e., 15%). Selection gain was predicted for each trait (DM, SW, and SY) based on the expression below: where GV i is the predicted genotypic value of progeny i and n is the number of selected progeny (30).
In order to perform a simultaneous selection and infer about the efficiency of selection gain for each evaluated trait between the frequentist (FSTME and FMTME) and Bayesian (BSTME and BMTME) models, we applied the additive genetic index using Selegen REML/BLUP software (AGI- [58]). The predicted breeding values of the selected progeny were thus compared by the frequentist and Bayesian approaches. The weights for each trait were defined based on the coefficients of genetic variation [59]. For all traits, the progeny were selected to increase the phenotypic expression, or provide the highest BLUP. After the direction of selection was defined, the genotypic values for each progeny [weighted by the pre-established weights (CV g ) for each trait] were summed, generating the AGI value. Subsequently, they were organized in descending order.

Analysis of deviance and model fitting
The significance of progeny (G) and G×E interaction effects of the FSTME model were evaluated. Significant G and G×E interaction effects (P � 0.01) were detected in the LRT for the DM, SW, and SY traits (Table 1). According to the AIC from the results obtained with the FSTME model, the model including the G and G×E interaction effects (full model) showed the best fit (lowest AIC value) for all traits (Table 2). Thus, according to the two methodologies, the full FSTME model was the most suitable to estimate the genetic parameters and predict the genotypic values. Likewise, the FMTME model was also appropriate, as it showed the lowest AIC value. In the Bayesian models, all chains achieved convergence by the criterion of Geweke [45]. Overall, the DIC values were smaller when using the full Model (considering genotype and G × E interaction effects), being the difference in relation to full Model higher than 2 ( Table 2), which according to Spiegelhalter et al. [42] it's enough to suggest that the use of full Model can lead to higher accuracy in estimating the parameters (Table 2). Therefore, since this model component is important the "best" genotypes measured in different environments couldn't the same.

Genetic and non-genetic components
The estimates of genotypic variance ðs 2 g Þ and G×E interaction variance ðs 2 int Þ were higher and lower, respectively, when obtained via FMTME and BMTME, for all traits. By contrast, the residual variance (s 2 res ) and individual phenotypic variance (s 2 phen ) estimates were similar regardless of the procedure used for the frequentist models, whereas in the Bayesian models these estimates were higher for all traits, even for the genotypic component ðs 2 g Þ. However, the s 2 res and s 2 phen estimates were similar between the BSTME and BMTME models. As a consequence, both the frequentist and Bayesian MTME models led to higher estimates of h 2 prog , CV g , and Ac prog compared with the frequentists and Bayesian single-trait models. The exception was the SW trait, for which no significant differences were observed (Table 3) between single-and multi-trait multi-environment models. For the other components, both procedures led to results proportional to the estimates of s 2 g and s 2 int . Compared with FSTME, the FMTME model provided increases of the orders of 1.81, 0.22, and 21.72% in the estimate of s 2 g and reductions of 6.63, 0.26, and 17.61% in the estimate of s 2 int for the respective traits DM, SW, and SY. As a result, the h 2 prog estimates increased by 1.03, 0.1, and 15.94%; CV g , by 0.88, 0.14 and 10.28%; and Ac prog , by 4.38, 0.25, and 42.29% for the respective traits.
Compared with BSTME, the BMTME model led to 1.56, 0.02, and 31.81% higher estimates of s 2 g and 8.56, 0.4, and 32.57% lower estimates of s 2 int for the respective traits DM, SW, and SY. The reduced ranging of h 2 prog (lower and upper difference) with Bayesian MTME credible intervals (probability of 95%) were 9.66, 0.22, and 42.44% for the DM, SW, and SY traits, Table 2 respectively (Table 4). Thus, the mean increased by 1.08, 0.00, and 24.01%; CV g increased by 0.88, 0.14, and 10.82%; and Ac prog increased by 1.04, 0.28, and 21.37% for the respective traits.

Akaike information criteria (AIC)
The posterior density intervals of heritability genetic parameters (Table 4, and Figs 1-3) were accessed to assist in the selection of genotypes in the Bayesian models. Thus, the breeding values and their HPD intervals obtained from the Bayesian STME and MTME models for each trait can be useful tools in progeny selection. Posterior density intervals of estimates of genotypic variance for BSTME and BMTME were showed in S1-S4 Figs.

Genetic correlation
Genetic correlations between the DM, SW, and SY traits obtained by the frequentist and Bayesian models are given in Table 5. For the DM-SW and SW-SY pairs, low correlations were detected in every comparison, indicating absence of linear associations. However, a high and positive correlation was found between DM and SY in both methodologies. Additionally, the correlation between this pair of traits estimated via FSTME and BSTME exceeded the parameter space. The same was not true when it was estimated via FMTME and BMTME, which showed similar values within the parameter space and were thus more realistic (not biased). Table 3. Estimates of variance components and genetic and non-genetic parameters for number of days to maturity (DM), 100-seed weight (SW) (grams), and average seed yield per plot (SY) (grams) evaluated in 203 soybean F 2:4 progeny via frequentist single-trait multi-environment (FSTME) and multi-trait multi-environment (FMTME) and Bayesian single-trait multi-environment (BSTME) and multi-trait multi-environment (BMTME) models.  Table 4. Posterior inferences for mode, mean, median, and higher posterior density (HPD) interval of the broad-sense heritability per plot, considering the proposed Bayesian single-trait multi-environment (BSTME) and multi-trait multi-environment (BMTME) models for number of days to maturity (DM), 100-seed weight (SW) (grams), and average seed yield per plot (SY) (grams). Genetic selection of segregating soybean progeny

Progeny selection
The Spearman rank correlations between breeding values via single-and multi-trait multienvironment models (FSTME and FMTME; BSTME and BMTME) were significant in all comparisons. This correlation was medium for SY (79.16 and 70.92) in the frequentist and Bayesian models, respectively, and high for the other comparisons. These results were confirmed by the agreement (Table 6).
In the frequentist and Bayesian analyses, the predicted selection gains were equivalent for all traits. However, the MTME models showed greater gains. For the SW trait, both the frequentist and the Bayesian procedures generated very similar results (Table 6). For DM and SY, however, there was less agreement between the selected progeny, especially for the SY trait, which culminated in greater discrepancy between the gains predicted from selection.
For the DM trait, the MTME models led to increased gains predicted from selection: 4.01 and 5.13% for the frequentist and Bayesian methodologies, respectively. Considering the SY trait, for which the agreement between the selected progeny was lower than 50%, the MTME models showed increases of 20.26% and 28.43% (frequentist and Bayesian models, respectively) in the gain predicted from selection, as compared with their respective single-trait models (Table 6).
When the AGI was used for the simultaneous selection of the 30 best soybean progeny, higher index gains were found for the FMTME and BMTME models (Tables 7 and 8) compared with the gain of the overall mean of AGI. The gains predicted from selection were similar for both models (FMTME and BMTME), for all traits. Moreover, greater gains were observed for the SY variable (15.03 and 15.71% for FMTME and BMTME, respectively).

Analysis of deviance
The LRT for the FSTME model revealed that the progeny and G×E interaction effects are significant (P < 0.01) for the DM, SW, and SY traits. Consequently, the respective variance Genetic selection of segregating soybean progeny components are significantly different from zero and so are the respective coefficients of determination ( Table 1). The fit of the frequentist models was checked by AIC. This criterion indicated the full model as the most suitable to estimate the variance components and predict the genotypic values (Table 2). For the Bayesian approach, the full (considering genotype and G × E interaction effects) and reduced (only genotype or G x E interactions effects) models were compared through DIC (Deviance Information Criterion), which suggests that models

Table 5. Genetic correlations between number of days to maturity (DM), 100-seed weight (SW) (grams), and average seed yield per plot (SY) (grams) evaluated in 203 soybean F 2:4 progeny via frequentist single-trait multi-environment (FSTME) and multi-trait multi-environment (FMTME) and Bayesian single-trait multienvironment (BSTME) and multi-trait multi-environment (BMTME) models.
Correlation FSTME FMTME BSTME BMTME with smaller DIC are better supported by the data. According to Spiegelhalter et al. [42], models with differences in DIC values lower than 2 need to be considered as equally well. Therefore, since DIC values obtained were higher than 2, it is possible to indicate the superiority of full model over the restrict models. The generalization of AIC is the most common method of assessing the fit of a statistical model estimated via Bayesian inference (DIC). The effects of Bayesian models can be used as an inference to the test of hypothesis [19]. In Bayesian statistics, the lowest expected deviance has the HPD [60], and this was observed in the present study for the MTME model (Tables 2 and 3).
In both criteria (AIC and DIC) for the choice of statistical models, the obtained results revealed that the MTME models showed the best fit, explaining the genetic variability of the experiment and selection considering the genetic (progeny) and environmental interaction effects ( Table 2).

Variance components
Variance components are the variances associated with the random effects of a model. Knowing them is of great importance in genetics and breeding, since the population and the breeding method to be used depend on information that can be obtained from these components. The solution of mixed-model equations depends on knowledge of the variance and covariance matrix, whose structure is known, but its components often are not. At present, the standard method for the estimation of variance components is REML, developed by Patterson and Thompson [32].
The BLUP method [33] maximizes the correlation between the predicted and the true genotypic value; i.e., it minimizes the prediction error variance (PEV). Additionally, it is not biased, as we expect the predicted genotypic value to be equal to the true genotypic value [61]. Further, BLUP allows for the simultaneous use of several sources of information as well as information originating from experiments carried out in one or various locations and evaluated in one or various harvests [62].
Although the mixed-model methodology by the frequentist approach has several desirable characteristics [49], the adoption of Bayesian statistical inference for genetic evaluation in the breeding of crop species has shown to be advantageous. Bayesian models have been used since 1986 [63] and further exploited in recent years [23][24][25]64,65] due to the great computational advancements and new methodological applications and elucidations.
Bayesian analysis is based on the knowledge of the posterior distribution of the parameters to be estimated. This allows for the construction of exact credibility intervals for the estimates of random variables, variance components, and fixed effects [66]. Higher values for the interval with 95% credibility of distribution for the broad-sense heritability parameter found in this study (Table 4) were also presented in the study of Torres et al. [25] to estimate genetic Table 6

Trait
Predicted selection gain Agreement (%) FSTME FMTME BSTME BMTME FSTME × FMTME BSTME × BMTME parameters for N-uptake efficiency and N-utilization efficiency under contrasting N levels in the soil via BMTME models. The difference between mean, mode, and median of broad-sense heritability estimates (Table 4) reflects some lack of symmetry in the posterior distribution estimates [38]. However, for the SY trait, differences between the BSTME and BMTME are clear when we analyze the posterior densities, mainly because the posterior MTME resulted in a more narrow and symmetric distribution, confirming the increase in precision (Fig 3).
When the prior distribution is informative, the credibility interval tends to be narrower than the confidence intervals. When the mixed-model parameters are assigned non-informative distributions, Bayesian and frequentist inferences should be equivalent [67]. Mathew et al. [68] showed that Bayesian inference is superior to frequentist inference when the posterior distribution of a variance component is bimodal. Mathew et al. [68], Waldmann et al. [69], Schenkel et al. [70], and Harville et al. [71] did not find relevant differences between the breeding values predicted by frequentist or Bayesian approach. Schenkel et al. [70] also observed that the breeding values presented the same bias and accuracy. Silva et al. [72] found results from noninformative analyses and results from REML/BLUP analyses (frequentist) for some components of variance and heritability and for breeding values. The specific results obtained by the frequentist and Bayesian approaches were similar (Table 3). This was expected, Table 8 since non-informative prior distributions were used in Bayesian analysis. The modes of the marginal posterior distributions of the genetic parameters were similar to the corresponding REML estimates. From the Bayesian point-of-view, the estimates obtained via REML correspond to the modes of the combined posterior distributions of the variance components, obtained by Bayesian approach, given the use of uniform priors for the fixed effects and variance components [66]. The frequentist and Bayesian MTME models provided higher s 2 g estimates and lower s 2 int estimates, which resulted in higher h 2 prog for all evaluated traits. The highest h 2 prog estimates were found for DM and SW and the lowest for SY, which confirms that DM and SW are less complex traits and are thus less influenced by the environment than SY [73][74][75].

. Order, progeny (Prog), breeding value (u+g), and Additive Genetic Index (AGI) of the 30 progeny selected simultaneously via Bayesian single-trait multienvironment (BSTME) and multi-trait multi-environment (BMTME) models for number of days to maturity (DM), 100-seed weight (SW; grams), and average seed yield per plot (SY; grams) evaluated in 203 soybean
The genotypic coefficient of variation (CV g ) quantifies the magnitude of genetic variation available for selection, and thus high values are desirable [76]. In this way, the increase seen in this parameter with the use of the MTME models is important for breeding programs. The residual coefficient of variation (CV e ) is a measurement of experimental precision of statistical and non-genetic nature. According to Resende and Duarte [52], CV e is of moderate magnitude for the SY trait and low magnitude for the DM and SW traits, indicating good experimental precision. Moreover, as expected, there were no alterations in the CV e estimates when the FSTME and FMTME models were used (Table 3).

Genetic correlations
Studies on genotypic, phenotypic, and environmental correlations in soybean involve traits that are evaluated from flowering to maturity; notably, yield and its components [77][78][79][80]. Our results corroborate those reported by Cober et al. [81], who also obtained a high correlation between DM and SY and no linear associations between the DM-SW and SW-SY pairs. The authors argued that the genes controlling maturity in soybean have pleiotropic effects with grain yield. Ablett et al. [82] and Lee et al. [83] reported that late maturity was associated with high yield. Liu et al. [84] investigated the genetic architecture of three growth period traits and confirmed that the soybean growth stages are highly correlated with grain yield. Li et al. [85] and Zhang et al. [75] observed that the association between the BARC-016957-02165 marker and seed yield was located in the same region as a QTL controlling pod maturity on chromosome 6, which explains the high correlation observed between these traits, as the QTL were closer to each other.
According to Pollak et al. [14], selection biases may occur when traits are analyzed individually. This bias was observed in the present study, where the genetic correlation value between the DM and SY traits exceeded the parameter space (value higher than 1) ( Table 5). Viana et al. [15] evaluated two traits in popcorn and also found that the genetic correlation obtained with the single-trait model exceeded the parameter space.
According to Thompson and Meyer [86], the increase in accuracy obtained with the use of multi-trait BLUP analysis compared with single-trait analysis is proportional to the difference between the genetic and environmental correlations of the analyzed traits. In the context of whole-genome prediction, Jia et al. [87], Guo et al. [88], and Jiang et al. [22] found that joint prediction of multiple traits benefits from genetic correlations between traits and significantly improves prediction accuracy compared to single-trait methods, specifically for low-heritability traits that are genetically correlated with a high-heritability trait. This fact was observed in our study, in which the DM variable showed high heritability and high correlation with SY, consequently generating significant increases in selection accuracy for SY. However, for both methodologies-frequentist and Bayesian-there was no significant increase in selection accuracy for the SW trait, as verified by its low correlation with the other evaluated traits.

Progeny selection
The observed differences in the genotypic values predicted by the Bayesian and BLUP/REML procedures were small, leading to a slight alteration in the ranking of the progeny selected by both procedures. This finding was confirmed by Rank Spearman correlation (Table 6). However, compared with the FSTME and BSTME models, the respective MTME models showed higher genetic gains for DM and, especially, for SY. According to Resende et al. [66], these are the conditions for there to be correspondence between the frequentist and Bayesian methodologies for the fixed and random effect parameters: attribution of non-informative priori for the fixed effects, normal priori for the random effects, and normal likelihood for the observations vector. These promises were used in the present study, which explains the obtained results.
Despite the high agreement between the progeny selected for the DM and SW traits by both procedures, there was little agreement for the SY trait, which resulted in greater gains predicted from selection via MTME. Piepho et al. [3] and Piepho et al. [62] recommended the use of multiple-trait models to predict breeding values in annual crops, because this procedure has the best statistical properties and provides more-accurate results.
Resende et al. [19] and Okeke et al. [89] stated that one of the main advantages of using multivariate models is higher selection accuracy. Okeke et al. [89] also reported that multi-environment models were useful for understanding G×E interactions. Higher Ac prog were observed when obtained via MTME for all traits, with SY standing out with 42.29 and 21.37% increases in the frequentist and Bayesian models, respectively. Greater accuracy and efficiency of multiple-trait models were also reported by Viana et al. [15] in selection among and within half-sib families.
However, it must be stressed that the BSTME model obtained superior accuracy in comparison with the FSTME model for DM and SY, despite the similar broad-sense heritability values. This is explained by the use of the estimator of selection accuracy. In this regard, Resende et al. [19] described that when Bayesian accuracy is higher than frequentist accuracy, the distribution of the parameters attributed to Bayesian approach were probably more adequate than those associated with the traditional model. The opposite can be considered true for the SW variable, for which the frequentist models obtained better results due to a better adjustment of the normal distribution of the parameters attributed to the data. These conclusions are also valid for the MTME models, which exhibited different obtained accuracies; however, for the SY variable, the FMTME model showed the best fit according to the mean accuracy of the progeny ( Table 3).
As can be seen in Tables 7 and 8, desirable gains are obtained in selecting the best progeny for the DM, SW, and SY variables for all models based on the AGI. However, the high positive correlation between the DM and SY traits can favor the selection of high-yielding and latecycle progeny. In this case, selection indices can help breeders select progeny that exhibit gains for both traits simultaneously [90]. Although similar gains were found for both approaches employed, it should be stressed that there was a slight increase in gain (%) when the 30 best progeny were selected for the DM and SY traits using Bayesian approach. This is a desirable factor that should be taken into account by breeders.
According to Silva et al. [1], individual plants or progeny in the F 2:3 or F 2:4 generations can be selected aiming at the adoption of the recurrent selection method described by Hallauer et al. [91], by exploiting genetic variability among and within progeny. Silva et al. [1] stated that the F 2:4 generation is suitable for selection, since 87.5% (1.75) of the total additive genetic variance (2s 2 a ) that will be available in F 1 is already available in F 2:4 . Thus, progeny selection in F 2:4 through a more precise method is relevant.
Early in the generation of the base populations of the soybean breeding programs, many populations are commonly obtained at the expense of the number of progeny to be evaluated; i.e., the evaluation of future lines, be them for high or low heritability, is based on samples with a finite (small) number of progeny. Thus, Bayes' theorem is recommended for those situations, as it gives precise solutions to the problem of finite-size samples, because for each data setlarge or small-there is an exact posterior distribution to draw inferences.
The MTME models provided better results than the single-trait models using frequentist and Bayesian approach. Therefore, the former procedures can be efficiently applied in the genetic selection of segregating soybean progeny. However, it is necessary to use an adequate statistical tool that provides algorithms and routines to efficiently perform the analyses. Though not necessarily easy, the use of Bayesian inference in quantitative genetics in the breeding of crop species [69,72] is a tendency in breeding programs [5].

Conclusion remarks
For our data set, the average BMTME processing time using an Intel(R) i7-5500U (2.4 GHz) processor with 8 GB of RAM was 1 h 40 min and 35 s, corresponding to approximately 0.006 s for each MCMC iteration. Silva et al. [72] considered this performance plausible, but pointed out that improvements can be obtained using the conditional decompositions proposed by Hallander et al. [92]. For the same purpose, in addition to improving the prior information, Montesinos-López et al. [13] proposed a Bayesian model for analyzing multiple traits and multiple environments for the whole-genome prediction model. The authors also developed an Rsoftware package that offers specialized and optimized routines to efficiently perform the analyses under the proposed model. By contrast, the FMTME model took approximately 14 s to converge. Despite the considerable difference in processing time of the analysis and output size of the results (around 1.03 GB) due to the high number of interactions adopted, the Bayesian model showed to be efficient for the proposed objective. Furthermore, it provided additional results to those obtained by the frequentist approach, with noteworthy credibility intervals.
The Bayesian models have desirable potentials when using informative prior distributions, providing parameters with lower standard deviations and/or possible genetic gains. However, the quality of the informative prior may have questionable origins and may not generate considerable advantages. Silva et al. [72] showed it can be advantageous to implement a Bayesian framework for mixed-model analysis in the breeding of crop species using informative priors. However, for potential future studies in plant breeding, the implementation of informative prior fitted to MTME models can be the next step to be assessed.