Parameter estimation and selection efficiency under Bayesian and frequentist approaches in peach trials

Identification of stable and high-yielding genotypes is a real challenge in peach breeding, since genotype-by-environment interaction (GE) masks the performance of the materials. The aim of this work was to evaluate the effectiveness of parameter estimation and genotype selection solving the linear mixed models (LMM) under frequentist and Bayesian approaches. Fruit yield of 308 peach genotypes were assessed under different seasons and replication numbers arranged in a completely randomized design. Under the frequentist framework the restricted maximum likelihood method to estimate variance component and genotypic prediction was used. Different models considering environment, genotype and GE effects according to the likelihood ratio test and Akaike information criteria were compared. In the Bayesian approach, the mean and the variance components were assumed to be random variables having a priori non-informative distributions with known parameters. According the deviance information criteria the most suitable Bayesian model was selected. The full model was the most appropriate to calculate parameters and genotypic predictions, which were very similar in both approaches. Due to imbalance data, Cullis’s method was the most appropriate to estimate heritability. It was calculated at 0.80, and selecting above 5% of the genotypes, the realized gain of 14.80 kg tree1 was attained. Genotypic frequentist and Bayesian predictions showed a positive correlation (r = 0.9991; P = 0.0001). Since the Bayesian method incorporates the credible interval for genetic parameters, genotypic Bayesian prediction would be a more useful tool than the frequentist approach and allowed the selection of 17 high-yielding and stable genotypes.


Introduction
The peach tree [Prunus persica (L.) Batsch] belongs to the Rosaceae family and is native to China (Hedrick 1917). It is one of the main fruit species grown in temperate and subtropical zones around the world, with relevant commercial, ornamental, and economic value. According to FAOSTAT (2020) the main peach and nectarine producing countries are China, Spain and Italy, while Argentina is the tenth largest producer in the world with 217,217 tons and the second in South America after Chile. The successful activity of many breeding programs worldwide has led to the release of more than 1000 new cultivars in the last century (Mas-Gómez et al. 2021). The main peach-producing areas in Argentina are Mendoza, Rio Negro (in the so-called "High Valley") and the northeastern of Buenos Aires province. San Pedro Agricultural Experimental Station of the National Institute of Agricultural Technology (INTA San Pedro) is located in this last region and has a wide peach and nectarine germplasm that is evaluated every season. Many of the peach varieties used in Argentina are initially selected in INTA San Pedro and the final evaluation is performed in each one of the producing regions. Despite the large number of cultivars obtained, exist a constant market and grower demand for new ones. In this way, breeders should consider different potential selection criteria such as performance, stability, adaptability and quality traits. The genetic basis of traits and the environment have a strong impact on phenotypic variation, from which useful information can be extracted. Such information includes determination of several parameters such as mean, variances, association between characters as well as heritability of each one of them. Genotypeby-environment interaction (GE) is a common phenomenon in multi-environment trial (MET). In peach, it can seriously affect stability production (Maulión et al. 2014(Maulión et al. , 2016a. The average fruit yield of genotypes evaluated in different environments could be a good parameter to select materials according to their performance and stability. However, phenotype is the result of the joint action between genes and environment. Depending on the nature of the GE, it could bring a serious problem for the breeder's work, because it can reduce the association between phenotypic and genotypic value. This is, GE may translate into inconsistent genotype behavior altering accession rank (crossover interactions-COI) or changes in the response without reordering (non-crossover interactions-NCOI) (Cruz and Regazzi 1997). Yield data from MET are frequently analyzed model in which the total variability is partitioned into the variance components due to genotypic, environmental and GE effects. Plant breeders have traditionally estimated variance and covariance components using a two-way analysis of variance (ANOVA) applying the moment-based estimators (Searle et al. 1992). Considering that the accuracy of genotypic value estimation is essential in a breeding program, it should be carried out by the most appropriate methods (Borges et al. 2010). The most popular approaches for estimating quantitative genetic parameters make use of the linear mixed models (LMM). They allow accounting for various confounding effects (Kruuk 2004;Wilson et al. 2010), achieving substantial progress in highlighting issues pertaining to fixed effects in quantitative genetic inferences (Wilson et al. 2010;Wolak et al. 2015). LMM are estimated by maximum likelihood or maximum restricted likelihood (REML), and also Markov-Chain Monte Carlo (MCMC). While REML is associated with the frequentist approach, the MCMC is linked with the Bayesian framework. Under frequentist LMM, direct genotypic prediction using the best linear unbiased prediction (BLUP) (Henderson 1975(Henderson , 1983White and Hodge 1988) is more reliable than selection based on phenotype mean only. In this case, genotype can be assumed as a fixed (best linear unbiased estimation-BLUE) or random effect (BLUP) (Searle et al. 1992). According to Piepho and Möhring (2006), the most reasonable is to consider that the genotypes under evaluation are a random sample of a hypothetical population in which selection was performed. The Bayesian method has increased its attractiveness to analyze data in several fields of knowledge through the MCMC procedure, since it is efficient even with an arbitrary number of random effects, and provides accurate estimates of the parameters of interest. Bayesian models have a number of advantages over frequentist ones (Bolker 2008;Clark 2005;Gelman et al. 2013). Bayesian parameter estimates are more accurate when the sample size is small; the interpretation of the results is easier and more direct since they indicate the probability that a parameter assumes a certain value. It is also possible to include measures of uncertainty, missing data and different levels of variability; and finally, it allows to Page 3 of 13 107 Vol.: (0123456789) specify the parameter distribution based on a priori information (prior). Several MCMC algorithms are used to estimate the posterior distributions of parameters such as Metropolis, Gibbs and Hamiltonian. The Hamiltonian Monte Carlo (HMC) is a generalization of the Metropolis algorithm that works faster, especially in high dimensions, eliminating the random walk behavior endemic to Random Walk Metropolis and the Gibbs samplers.
Peach yield is a character whose main characteristic is its COI (Maulión et al. 2014(Maulión et al. , 2016aAngelini et al. 2019), so that its statistical modeling is a real challenge. Measured and genetic parameters have been studied in reproductive characters (Souza et al. 1998a) and associated with peach fruits (Souza et al. 1998b) using LMM under the frequentist paradigm. However, this methodology was not used to model peach yield, or in the identification of high-yielding and stable genotypes. Moreover, to our knowledge, the Bayesian approach has not been used to estimate parameters in peach collection entries, so far. The objective of this work was to adjust frequentist and Bayesian LMM and compare their accuracy when genetic parameters are estimated, as well as to identify high-yielding and stable genotypes.

Dataset
Fruit yield data set was obtained from 18 season evaluations (SE) involving a total of 308 genotypes, the amount of them evaluated in each SE is shown (Table 1). They include peach and nectarines for fresh market consumption with different origin and pedigree: The replication number, the fruits yield per genotype by SE, as well as the name and origin of each accession have been reported (Supplementary File 1). Furthermore, the known pedigree of 117 genotypes used is indicated (Supplementary File 2).
Each SE corresponded to the evaluation of peach germplasm from INTA San Pedro located at 31° 41´12° S and 60° 47´32° W. Since all experiments were conducted in the same orchard but under different climatic conditions, each SE was considered different environment. Genotypes were arranged in a completely randomized design with a replication number between two and thirteen trees. Genotypes evaluated at least in three SE were included in the analysis. Each year the orchard received management recommended for the area; this included fungicide and insecticide sprays and pruning similar to the commercial orchards. To avoid overproduction and the consequent weakening of the plant, the amount of fruit was thinned 30 days after full blooming date, leaving one every 15-20 cm along the branch. Fruit yield was obtained for each replication by using a weighing machine and expressed in kg tree −1 .

Linear mixed model
The single-trait multi-environment statistical model is given by: where y ijk is the yield from the k-th replication of the i-th genotype in the j-th environment, is the overall mean, G i is the effect of the i-th genotype, E j is the effect of the j-th environment, GE ij is the interaction between the i-th genotype and the j-th environment, and ijk is the random error from the k-th replication of the i-th genotype in the j-th environment.
(1) Genotype (G), environment (E) and GE were assumed as random (Hinkelmann and Kempthorne 2005). Overall, apart from , the following random-effect assumptions were made:

Frequentist approach
The REML method was used to fit the mixed model in Eq. (1). To compare the fit of following models: M null (M N ): simply ; M W/G : only E; M W/GE : G + E and M full (M F ): E + G + GE, the likelihood ratio test (LRT) and Akaike information criteria (AIC) were used. The lowest value of AIC indicates the best model. LRT statistics approximately follows a chi-square distribution with degrees of freedom equal to the number of additional parameters in the more complex model. LRT and AIC have the following equations: where RL represents the restricted maximum likelihood; RL + and RL − are the maximum values of the restricted likelihood of the more complex and simpler models, respectively; d is the dimension of the model; n equals the number of observations valid for the REML estimate.

Bayesian approach
Mean and variance components were assumed to be random variables having a priori distributions, with known parameters. The most suitable prior information should be selected to fit the model. Several options for non-informative priors for scale parameters in hierarchical models were reviewed by Gelman (2006), and his recommendation includes the use of uniform, half-t, half-cauchy, and halfnormal families of distributions. The most suitable prior information should be selected to fit the model. Non-informative normal (0, 10 10 ) for , as well as eight standard deviation component distributions (P 1-8 ) were tested: P 1 : priors for the standard deviation components E , G , GE and follow inverse-gamma (10 −3 ,10 −3 ); P 2 : priors for the standard deviation components E , G , GE and follow gamma (0.5,0.5); P 3 : priors for the standard deviation components E , G , GE and follow log-normal (0,1); P 4 : priors for the standard deviation components E , G , GE and follow the positive part of the normal distribution denoted as half-normal (0,1) or equivalently as normal(0,1) + ; P 5 : priors for the standard deviation components E , G , GE and follow uniform (0,1000), P 6 : priors for the standard deviation components E , G , GE and follow the positive part of the normal distribution denoted as half-normal (0,10) or equivalently as normal(0,10) + ; P 7 : priors for the standard deviation components E , G , GE and follow the positive part of the cauchy distribution denoted as half-cauchy (0,5) or equivalently as cauchy (0,5) + ; P 8 : priors for the standard deviation components E , G , GE and follow the positive part of the cauchy distribution denoted as half-cauchy (0, 2.5*sd) or equivalently as cauchy (0, 2.5*sd) + , where sd is de standard deviation of the residuals of the frequentist model.
Model P 1 to P 8 were adjusted using four chains, each one comprising 4,000 iterations (the first 2,000 ones were discarded) applying the HMC algorithm. Gelman-Rubin (R) statistics, which compare the variation between and within chains (Gelman and Rubin 1992), were used, R̂ values lower than 1.1 were good evidence in favor of chain convergence. MCMC returns samples with some degree of autocorrelation, which varies with the model and algorithm selected. Sampling quality was evaluated using the effective sample size. It is a measure of the real sample size, discarding those values with correlation in successive simulations. Prediction accuracy of the Bayesian model was established through leave-one-out cross-validation, ê lpd loo . The difference of ê lpd loo between two models was expressed in a deviance information criteria (DIC) scale, multiplicating ê lpd loo by −2. A smaller value of DIC reflected a better suitability of the chosen priors and the most suitable model among M N , M W/G , M W/GE and M F . Genotypic prediction with Bayesian (G-BP) and frequentist (G-BLUP) models were obtained considering peach accessions unrelated, while the degree of similarity between both predictions was measured using Spearman correlation.
Stable and high-yielding genotypes were ranked and selected based on their G-BP and predicted performance ( ŷ i ) estimated as: ŷ i =̂ +Ĝ i +Ê +ĜE i (2), where ̂ is the overall mean estimate, Ĝ i is the genotypic prediction of i-th genotype; Ê is the mean E prediction and Ĝ E i is the mean GE prediction of the i-th genotype. Spearman's correlation to establish the association between ŷ i and G-BP of selected genotypes was computed.
Heritability and genetic gain due to selection According to Hanson (1963), heritability is the fraction of the selection differential expected to be gained when selection is practiced on a defined reference unit. Although some degree of disagreement in the estimated values is always expected, since the estimate of heritability depends on the variability of the population and the environment where the experiment was performed (Falconer and Mackay 1996), several heritability estimators can be evaluated in terms of response to selection. Two broad-sense heritability measures were estimated. First, the variance components estimated under frequentist and Bayesian methods can be calculated from: . However, Piepho and Möhring (2006) consider that this formula would not be suitable for computing response to selection, essentially, because the correlation between phenotype and response to selection differs among genotypes. That is, there is no linear relation between response to selection and selection differential, as in the trials with balanced data. Second, the estimator proposed by Cullis et al. (2006) is widely used to compute heritability when data consists of unbalanced and/or correlated genetic effects. In this case, the genetic term is fitted as a random effect by calculation of the best linear unbiased prediction (BLUP), and it uses the square of the standard error of the genetic estimates (across environments) to attain an approximation of the non-genetic variation. The equation is as follows: where v BLUP is the average variance of pair-wise differences between BLUPs.
The genetic gain (GG) based on mean over replications was calculated for selection intensity of p is C G h∕̂ , where C is a constant given by C = 1 p √ 2 e −z 2 p ∕2 and z p is the upper p quantile of standard normal distribution, h is the square root of the heritability, and ̂ is the mean response (Hinkelmann and Kempthorne 2005;Singh et al. 2012). C value is 2.06 for a 5% intensity of selection. The realized gain using h 2 B or h 2 BC is: 1-RG 1 = h 2 B .S and 2-RG 2 = h 2 BC .S, where S is the selection differential (Falconer and Mackay 1996).
Bayesian and frequentist analyzes were conducted with the R software (R Core Team, 2018). Bayesian inference was performed through the Stan language (Stan Development Team, 2017) using a script written in the lab available on request (angelini@cefobiconicet.gov.ar), while the frequentist procedure was implemented using lme4 package.

Selection models
Several Bayesian models with non-informative prior distribution according to DIC statistics were evaluated. According to Gelman et al. (2004), the prior set P 3 was the most appropriate since the lowest DIC values for each one was obtained (DIC = 32,808.70). In an ongoing breeding program over many years, information about distribution of estimated parameters is available. It can be included as informative prior to increase accuracy of inference estimated parameters . Since this is the first work of peach yield under Bayesian approaches, no prior information can be used. Therefore, genetic and nongenetic components were estimated using P 3 prior. Current trial estimates of variance components can be used as prior distributions for similar future trials. However, the posteriori distributions of the parameters naturally provide an update of the a priori distributions in the light of current datasets.
To evaluate the fit of the models M null (M N ): simply ; M W/G : only E; M W/GE : G + E and M full (M F ): E + G + GE, LTR and AIC measures were considered for the frequentist model, while the DIC was used for the Bayesian approach. Analyzing the measurements obtained with all the models (Table 2), it is evident that considering only the general average (M N ) was the least suitable, although it improved significantly after the environmental effect (E) was included (M W/G ). Since peach yield is sensitive to climatic variability because change in temperature and precipitation can alter blooming and fruit development period, a significant effect of E on the model is expected. Bayesian and frequentist environmental effects (E-BP and E-BLUP) were similar (Table 3). According to the magnitude of E-BP and E-BLUP the SE can be clustered (C) in C1, C2 and C3. The C1 grouped SE3, SE6, SE16 and SE17 whose high and positive predictions would favor performance. On the contrary, C2 clustered SE4, SE8, SE10, SE11 and SE15, which had a strong negative influence on yield, while C3 enclosed the remaining SE1, SE2, SE5, SE7, SE9, SE12, SE13, SE14 and SE18 with moderate positive or negative effects due to interval limit values. Maulión et al. (2014) have already indicated that the interaction between rain and heat accumulation during fruit development period was associated with performance instability. Later, a portion of yield variation was explained by heat accumulation during fruit development period, rainfall during floral bud endo-dormancy and rainfall from floral bud endo-to ecodormancy (Maulión et al. 2016a, b). The influence of climatic variables were reported previously by Erez and Couvillon (1987), and Lopez and DeJong (2007) among others, as essential for flowering and fruit development. Genotypes whose requirements do not match the climatic conditions of the region cannot be recommended for sustainable production. Angelini et al. (2019) showed that erratic behavior can be reduced by the release of cultivars adapted to a set of environments. Therefore, understanding of the nature of GE and, consequently, the matching of a genotype with appropriate locations to ensure a high and stable production is a major aim of any genetic program.
Inclusion of the G effect explains an important portion of data variation as shown by the lower AIC value and the significance of the LRT M W/GE model (Table 2). This genotypic variability is caused partly because accessions of peach germplasm at INTA San Pedro were originated in different regions around the world with different pedigrees and origin, making possible the selection and development of new superior cultivars or parents. Finally, the M F model had the lowest AIC, a significant LRT value and also the lowest DIC value (Table 2), representing the best model under frequentist and Bayesian frameworks to show the sensitivity to detect all variability sources. Since the presence of GE is a natural phenomenon in MET crops, including this effect in the final model to analyze agronomic data is a common fact, the full model in sorghum (Omer and Singh 2017), wheat (Mohammadi et al. 2015), cotton ), soybean (Volpato et al. 2019), among many others, was used. Although the GE characterization for main crops has been studied and reported extensively for many years, GE studies in peach started recently (Citadin et al. 2014;Maulión et al. 2014). Maulión et al. (2014) showed that peach yield is a character in which COI predominates over NCOI. Recently, Angelini et al. (2019) determined the most accurate tests about kind of crossover in a set of peach genotypes were present. Therefore, full model was the most proper to estimate the genetic parameters and predict genotypic values.

Estimates of variance components, heritability and genetic gain
The estimates of average, variance components, heritability and GG are shown (Table 4). Genetic and nongenetic variance components estimated under frequentist and Bayesian models were very similar. For all estimates, very narrow credible and confidence intervals for Bayesian and frequentist methodologies were obtained, respectively. Variance components are associated with random effects and their knowledge are extremely necessary in genetic breeding. Currently, that components are calculating using the REML method (Patterson and Thompson 1971), as well as the BLUP and genetic parameters. Variance component estimates via REML can be affected by missing data in any character measured on a single three (Cappa and Cantet 2006).
The authors proposed as an alternative the application of REML using the Bayesian approach through MCMC. Such procedures are used to determine the marginal posterior distributions of the parameters by algorithms converging to such marginal densities. The Bayesian procedure is an important tool in genetic evaluation since it considers the variability of all parameters in the model (Wright et al. 2000). As mentioned above, both approaches estimated parameter values and intervals of the same magnitude. Under these conditions there would be no difference between Bayesian and frequentist approaches. Bayesian is superior to frequentist inference when the posterior distribution of a variance component is bimodal (Mathew et al. 2012). Sorensen and Gianola (2002) reported that when the mixed-model parameters are assigned non-informative distributions, Bayesian and frequentist inferences should be equivalent. On the contrary, the Bayesian approach showed a large reduction of error in estimate measures when the prior distribution is informative .
Spearman coefficient correlation between Bayesian and frequentist genotypic predictions was 0.9997 (P = 0.0001). High correlation between the estimated predictions by REML/BLUP and independent chains (a variant of the methods of MCMC) were reported in maize (Mora and Arnhold 2006) cotton , and also between REML/BLUP and Gibbs algorithm in Eucalyptus (Soria et al. 1998;Mora and Perret 2007) and Pinus sylvestris L. (Waldmann and Ericsson 2006). Previously, Blasco (2001) demonstrated that BLUP can be considered a Bayesian estimator, when this was constructed using a uniform and normal a priori distribution for environmental and genetic effects, respectively. This study identified a high correlation between REML/BLUP and the HMC algorithm of genotypic predictions of estimates. Therefore, independently from different considerations on the REML/BLUP and Bayesian algorithms used, estimates of both methods appear to be consistent in plant breeding. According to G-BP and G-BLUP predictions the best 46 genotypes (about 15% of the selection index) were selected (Table 5). The high correlation (0.9997; P = 0.0001) indicates a strong linear association between the predictions, consequently, no discrepancy among the selected genotypes was identified. An advantage of Bayesian procedures is the possibility of construing credible intervals of genetic parameters, including the G-BP values, which are obtained directly from the a posteriori distribution (Wright et al. 2000).
Confidence intervals in the classical context can eventually be built by semi-parametric bootstrap techniques (Efron 1979). However, in the case of G-BLUP, Morris (2002) clarifies that bootstrap constantly underestimates data variation in finite samples. In this way, the G-BLUP will not be the best prediction. In the Bayesian perspective, the breeder may choose to select genotypes considering credibility limits obtained through a posteriori distribution for G-BP (Table 5). In this context, the genotypes that showed credibility intervals with negative values should not be selected because the range includes negative prediction values, indicating that the prediction values could be lower than the trial average. Consequently, the phenotypic efficiency would not be efficient and the genetic gain would be lower than the one obtained, 14.80 kg three −1 (Table 4), if only 17 from initial 46 genotypes, with positive values in their credible intervals, had not been selected (marked in gray color, Table 5).
Heritability is the key parameter to study genetic changes in a breeding population undergoing selection among alternative breeding procedures (Cockerham 1963;Hill 1971). In this work, using the Eq. (3), the heritability yield was estimated at 0.15 (Table 4), which is extremely low, even for a very complex trait such as peach yield. Based on the best 17 genotypes and a selection index around 5%, the genetic gain would be 0.13 kg tree −1 (Table 4). Minvielle (1990) mentioned that a high estimated value of heritability indicates that the correlation between phenotype and genotype of the two individuals is also high, and an efficient phenotypic selection can be made. Therefore, low heritability values estimated in this study suggest that simple selection methods cannot be implemented. According to Cullis et al. (2006), calculation of broad-sense heritability based on mentioned Eq. (3), is expected to perform not so well in the case of strong unbalanced and genetic correlation due to the pedigree, in particular when response to selection is computed from them. Here, response to selection was not based on the accession pedigree, but our experiment is unbalanced and correlation between genetic effects must be assumed as existent given the kinship of the modern peach germplasm (Aranzana et al. 2003). Therefore, the broad-sense heritability estimates according to Cullis et al. (2006), Eq. (4), was five times greater than calculated with Eq. (3), h 2 BC = 0.80, however, it represented an increase in GG C (0.13 to 0.22 kg tree −1 ). Both realized gains, RG 1 and RG 2 , were greater than GG and GGc, especially RG 2 , because a large differential selection (S) was combined with a higher heritability (Table 4). To achieve efficient phenotypic selection, heritability and genetic gain are key parameters in a breeding program, given that they lead to changes on the crop genetic base, through crossings, selection, and inclusion of promising genotypes along with the elimination of those with poor performance. In this context, the interest lies in the prediction of superior genotypes. Selection index close to 5% should lead to relatively rapid genetic progress for yield, 14.80 kg tree 1 . Selection based on performance should be considered in a broader context, e. g. flower density, flower per node, fruit set and size, among other traits, to achieve proper configuration traits to yield. At the same time, measurements and analysis of characteristics correlated (Nikolić et al. 2010;Cantín et al. 2010;Maulión et al. 2016a, b), simultaneous selection of characters  On the other hand, although accessions SP 10-32, M. Fortininer, SP 12-32, 55RA15 and Barceló were selected among the best, they did not show an outstanding GE-PB in any SE. Since none of the materials kept their performance prediction in all SE, the instability of the character due to its strong GE is confirmed once more. Because the genotype number evaluated in each environment was different, it is not possible to establish a rank based only on GE-BP to select the most stable and highly productive materials. Phenotypic predicted ( ŷ i ) value according to Eq. (2) for all 17 genotypes selected was calculated (Supplementary File 4). Coefficient correlation between the rank of G-BP and ŷ i showed a highly positive correlation: r = 0.85 (P = 0.0001). Small discrepancies between ŷ i and G-BP order indicates that simultaneous selection for high yield and stability could be made considering only the G-BP, without the calculation of ŷ i since it will simply delay the selection process. Previously, Maulión et al. (2014) showed that fruit yield was strongly correlated with superiority (Pi) (Lin and Binns 1988), reliability indexes (Ii) (Eskridge 1990), with non-parametric measurements for the genotype selection index (GSI) (Farshadfar 2008) and the rank sum (RS) (Farshadfar et al. 2011). However, Bayesian prediction can be used in unbalanced experiments or trials with missing data, but not the Pi, Ii, GSI and RS measurements. On the other hand, Bayesian approaches allow incorporation of more information than fruit yield as well as more accurate estimates than those obtained by other methods.

Conclusions
The outcome of this work can be summarized as follows: (1) Near all estimated parameters and genotypic predictions were very similar under frequentist and Bayesian methods, (2) Bayesian genotypic prediction is preferable when selecting high-yielding and stable genotypes, (3) A significant genetic gain was achieved when heritability was estimated according to the Cullis method and (4) Genotype-by-environment interaction as well as genotype rank variation across environments suggest that selection of superior genotypes should be performed in homogenous environments.
Acknowledgements The authors appreciate funding support from Consejo Nacional de Investigaciones Científicas y Técnicas (CONICET). PIP N° 0411.
Author contributions GDLC and GHV are corresponding author. GDLC and JA wrote the paper and conducted the statistical analysis. GHV coordinated trials, performed the field measurements, generated phenotypes and revised the paper. EBB, GSF and CFP wrote and edited different versions of the manuscript.