Enhanced Bayesian model for multienvironmental selection of winter hybrids maize: assessing grain yield using ‘ProbBreed’

Abstract

Background

Crossover interactions stemming from phenotypic plasticity complicate selection decisions when evaluating hybrid maize with superior grain yield and consistent performance. Consequently, a two-year, region-wide investigation of 45 hybrids maize across Nepal was performed with the aim of disclosing both site and wide adapted hybrids. Utilizing an innovative “ProbBreed” package, based on Bayesian probability analysis of randomized complete block designs with three replicated trials at each station, this study substantively streamlines hybrids maize selection.

Results

This finding revealed substantial genetic, environmental, and interactive influences on grain yield (p < 0.05). Among the hybrids, DKC9149 (8.8 tons/ha) emerged as the elite with probability coefficient of (0.39), followed by NK6607(0.35 & 8.6 tons/ha). Joint probability analysis identified RMH1899 super (0.23 & 8.3 tons/ha), followed by RMH 666 (0.15 & 8.4 tons/ha) and Uttam 121 (0.11 & 8.6 tons/ha), all of which accounted for overall environmental conditions. Additionally, over the years, DKC 9149, NK 6607(0.18 & 8.6 tons/ha), GK 3254(0.18 & 8.5 tons/ha), Shann 111(0.12 & 8.4 tons/ha), Sweety 1(0.13 & 8.4 tons/ha), and ADV 756(0.10 & 8.2 tons/ha) consistently demonstrated superior performance and stability. Delving with site specific recommendations include Nepalgunj: RMH 9999(8.5 tons/ha), NK 6607(8.6 tons/ha); Parwanipur: DKC 9149, MM 2033(8.5 tons/ha); Rampur: ADV 756, DKC 9149, MM 2929(8.6 tons/ha); and Tarahara: GK 3254(8.5 tons/ha), NK 6607(8.6 tons/ha), Uttam 121.

Conclusion

Thus, Selected hybrids are predicted to outperform within the recommended domain. Over and above, integrating genomic information into Bayesian models expected to enhance prediction accuracy and expedite breeding progress.

Background

Selecting elite hybrids with superior grain yield potential and stability across a spectrum of environmental conditions necessitates the optimization of crop breeding strategies, which is a critical priority for crop breeding and improvement [1]. Although classical breeding and crop selection tactics have proven efficacious, often necessitate substantial exacting and exhibiting constraints in their ability to comprehensively elucidate the nuanced interplay between genetic (GEN) and environmental factors [2]. Neglecting genotype‒environment (GEI) interactions in multi environment trials (METs) increases the risk of erroneous cultivar recommendations for farmers, which exacerbates the yield gap between the research stations and farmers’ fields [3, 4]. While traditional statistical models have been employed to investigate the GEI and cultivar stability, each approach has inherent limitations. For example, Wricke’s ecovalence, a metric introduced by Wricke in 1962, reviewed by [5], quantifies the stability of a genotype by assessing its interaction with environmental conditions. However, its assumption of homogeneous variance across environments apt to inflated type I error rates, particularly when random GEI are included, potentially compromising the sensitivity of the analysis [6, 7]. Moreover, Finlay‒Wilkinson regression, a statistical method devised by Finlay and Wilkinson in 1963 [6], helps breeders evaluate the stability and adaptability of myriad genotype. By fitting a regression model, this approach characterizes the sensitivity of each genotype to environmental effects. However, its treatment of both environmental means and genotype effects as fixed can lead to substantial sampling variances in estimates, especially in incomplete designs, and unaccounted uncertainty compromise the accuracy of the results [6, 8]. In addition, as proposed by Eberhart and Russell (1966), the limited responsiveness model offers the potential to identify genotypes that exhibit superior performance in specific environments but may not possess widespread adaptability [9]. Further, the assumption of linearity in statistical models posits a linear relationship between genotype and environmental factors; this assumption may not accurately represent the more convoluted interactions observed in empirical data [10]. Furthermore, as a model designed by Shukla (1972) for stability analysis, the intricate mathematical underpinnings of this model are strongly likely to hinder its widespread adoption, demanding a greater degree of statistical proficiency from users than more straightforward models such as Eberhart and Russell [11]. Despite, although akin to Eberhart and Russell, this model may not fully account for nonlinear interactions or multifaceted environmental factors influencing genotype performance [12]. Similarly, the additive main effects multiplicative interaction model (AMMI) proposed by Gauch et al.l (1990) presupposes a normal distribution of residuals. Departures from this assumption likely lead to biased results [13]. What’s more, AMMI relies on the assumption of homoscedasticity, where the variance of residuals remains consistent across all genotypes and environments. Violations of this assumption may compromise the model’s ability to accurately characterize interaction effects [14]. Similarly, the genotype‒genotype‒environment (GGE) biplot method introduced by Yan et al. (2000) may exhibit limitations in intricate scenarios. While GGE biplots are valuable tools in many contexts, their effectiveness can be compromised when dealing with datasets characterized by high dimensionality or particularly complex interactions [15]. Thus, recent advancements in probability theory coupled with cutting-edge software packages, developed by Saulo team [16]—a novel Bayesian method that employs the posterior distribution to obtain (HMC)Hamiltonian Monte Carlo estimates of performance and stability probabilities—offer a more streamlined decision-making process for selecting suitable genotypes across diverse environments, nevertheless prior knowledge of coding is needed. Predictive breeding programs routinely assess experimental genotypes through METs [17, 18]. The phenotypic manifestation of quantitative traits in METs is shaped by GEIs, which complicates selection owing to crossover interactions. Thus, the ability of deep learning to extract complex patterns from large datasets highly likely complements the strengths of Bayesian methods in handling uncertainty (unknown, unknown) and incorporating prior knowledge(risk) [19]. Besides, by furnishing probability distributions for parameters, Bayesian statistics facilitate a more nuanced comprehension of uncertainty. This approach empowers researchers to articulate the likelihood of a parameter residing within specific intervals, thereby offering a richer and more informative representation of inferential outcomes [20]. Also, capability to effectively accommodate unobserved or missing data by directly integrating the modeling of these uncertainties and credible interval [21].Thus, By harnessing this model, this study aimed to perform a probabilistic assessment of diverse winter hybrid maize varieties for predicting and informing recommendations on the basis of prime performance and stability across both general and specific environments.

On the basis of the scholar search engines—google scholar, WorldCat, ScienceDirect, DOAJ, AGORA, Bielefeld Academic Search Engine, PubMed, JSTOR, AGRIS— using keywords of “Bayesian methods”, “ProbBreed model”, “Statistical modeling”, “Bayesian probability”, only a few papers have been documented; thus, this study delves into 3 major aspects: (1) probabilistic assessment, (2) integration of GEIs, and (3) risk-based recommendations (known, unknown), which are widely applicable for all types of crops worldwide and appear underexploited as well as unparalleled. Thus, compared with traditional breeding methods, Bayesian probability methods do not significantly improve the accuracy of performance predictions, as indicated by the null hypotheses.

Method

Plant materials, experimental site & design

In a two-year experiment conducted during the year 2021–2023 across Nepal, 45 maize hybrids were evaluated at four research stations representing diverse climates with the primary goal of optimizing grain yield. Considering the winter maize cultivation area of Nepal where cold waves is more prevalent this study performed the four research stations leveraging randomized block design to evaluate the hybrids. A planting density of 66,666 plants per hectare was established by adopting a spacing of 60 cm between individual plants and 25 cm between rows with fivemeter row length/plot. Sowing was performed between November 13th and 23rd, progressing from east to west. Anthesis occurs 110–115 days after sowing (DAS) and harvesting is performed 1140 DAS. To capture the variability across climatic zones, four research stations were strategically selected: the eastern humid (Tarahara), the central intermediate zone (Bara, Pawanipur), the maize research hub (Rampur-Chitwan), and the far-western drylands (Nepalgunj). The experimental sites are shown in (Fig. 1). The code associated with this study map is deposited in the gitHub library of Bikas Basnet [22] (https://github.com/Bikasbasnettest/Loaction-map-in-R).

Fig. 1

Location map depicting the study area of 4 research station in winter season in Nepal

Statistical analysis

This study analyzed the data via a probabilistic approach, employing the “ProbBreed” package for Bayesian probability modeling [23]. This analysis was conducted in the R programming language, specifically version 4.3.1 (2023–06–16 ucrt) [24]. The code used in this study was archived on Bikas Basnet’s GitHub account [25](https://github.com/Bikasbasnettest/Bayesian-Model-in-R-for-study-of-MET).

Probability of superior performance

$$Pr({g_j} \in \Omega |y) = 1/S\sum\limits_{s = 1}^s {I(g_j^{(s)} \in \Omega k|y)}$$

where, J genotypes (j = 1,2, …………. j) were assessed in K environments (k = 1,2, …….k), with the observed phenotypes y(yield/ha). A subset of the high-performance selected genotypes, Ω, was determined according to the intensity of selection.

Probability of superior stability

$$\:\text{P}\text{r}\left[\text{v}\text{a}\text{r}\right({ge}_{jk}\:)\in\:{\Omega}|\text{y}]=\frac{1}{S}\sum\limits_{s=1}^{s}I\left[\text{v}\text{a}\text{r}\right({ge}_{jk}^{\left(s\right)})\in\:{\Omega\:}|\text{y}]$$

where $\:I\left[\text{v}\text{a}\text{r}\right({ge}_{jk}^{\left(s\right)})\in\:{\Omega\:}\:|\text{y}]$ indicates whether $\:\text{v}\text{a}\text{r}({ge}_{jk}^{\left(s\right)}$ exists in Ω [1] or not (0).

Joint probability of superior performance and stability

$$\begin{gathered} Pr[({g_j}) \in \Omega \cap \operatorname{var} (g{e_{jk}}) \in \Omega ] = \Pr ({g_j} \in \Omega )\, \times \,\Pr [\operatorname{var} (g{e_{jk}}) \in \Omega ] \hfill \\\hfill \\ \end{gathered}$$

Probability of superior performance within an environment

$$Pr[var({g_{jk}}) \in {\Omega _k}|y] = 1/S\sum\limits_{s = 1}^s {I[var(g_{jk}^{(s)}) \in {\Omega _k}|y]}$$

Results

Estimation of variance components and goodness of fit of the Bayesian P value

Bayesian analysis revealed that not only genotypic (GENs) factors but also environmental factors (ENVs) notably influenced maize grain yield/ha. However, while the environmental component (ENV) has the highest variance (3.861) and standard deviation (14.008), which plays a dominant role in yield variability, GENs, with a variance of 0.109, also contributed significantly to the substantial fluctuations in yield over time, as evidenced by the high variance associated with the year factor (1854.022), underscoring the critical influence of environmental conditions and management practices (Table 1). Moreover, GEN: ENV (40.070) and GEN: YEAR (0.335) elucidate the intricate relationships between genetic traits and environmental factors, stressing the need for nuanced predictive breeding strategies.

Table 1 Appraised components of the declared effect with their respective standard deviation (SD), naive standard error (naive SE), and inferior and superior high posterior density intervals [HPD (0.05) and HPD (0.95)] for grain yield/ha

Following the goodness of fit via Bayesian analysis, the p values for the median (0.344), mean (0.505), and standard deviation (0.545) were all greater than the conventional threshold of 0.05, indicating a lack of significant deviations from the expected distributions. Furthermore, the lower WAIC2 value (2903.349) emphasizes a strong fit to the data, as evidenced by the similarity between the sample distribution and empirical distributions. After the best model was fit, analysis revealed that the homogeneous residual variance was smaller than the heterogeneous residual variance (Table 2). Therefore, a smaller value is selected for further analysis, resulting in a more concise and precise model. Additionally, the effective parameters (94.285) imply a moderately complex model, whereas the mean R-hat value of 1.007 indicates vague convergence of the Markov chain Monte Carlo (MCMC) algorithm[14]. Finally, the effective sample size (1.158) confirmed that the provided data were sufficient to estimate the model parameters accurately. (Fig. 2-3) illustrates the posterior effects and goodness-of-fit diagnostics, and most of the overlapped portions indicate that the model is best fitted.

Table 2 Goodness-of-fit parameters: bayesian “P values” of test statistics [maximum, minimum, median, mean, and standard deviation], the effective number of parameters, WAIC2(widely Applicable Information Criterion), potential scale reduction factor (𝑅), and effective sample size for grain yield/ha

Fig. 2

Graphics depicting the posterior effect distribution leveraging density and histogram

Fig. 3

Graphics illustrating distribution of empirical vs. sample density plot yield/ha among the hybrid’s maize

HPD of the posterior genotypic main effects and probability of superior performance across the environment

On the basis of a comparative study of the genotypic main effect and probability of superior performance (Fig. 4), at 20% selection intensity to boost average grain yield. DKC9149 has a 0.53 probability of superior performance, revealing 47% poor performance at the given selection intensity follows the same interpretation for other hybrids, and NK6607 are the vertex hybrids with the superior coefficient value of the genotypic main effect. DKC 9149: Positive mean (0.1125), with HPD intervals suggesting strong potential. MM 2929 also shows promising performance, with a mean of 0.0758. The top-performing hybrid, DKC 9149, has a probability of 0.397, indicating a strong likelihood of achieving superior yields. NK 6607 and Uttam 121, with probabilities of 0.374 and 0.359, respectively, closely follow, highlighting a competitive edge in performance. As the probabilities decrease, genotypes such as RMH 666 and DKC 9144 also show promise but with diminishing advantages, followed by other genotypes. The lower probabilities for hybrids such as Delta 9999 and RH 10 suggest limited performance potential. Moreover, high variability is observed in the HPD intervals, particularly in hybrids such as MM 9440, which has a very low mean (-0.0959) and wide uncertainty.

Fig. 4

Graphics depicting posterior effects (95% credible intervals) of hybrids, along with probability of hybrid being superior

Pairwise probability of superior performance

To directly compare a selected hybrid with other promising hybrids, this study computed the pairwise probability of superior performance. The green, yellow, and red color legends represent low, moderate, and high probabilities, respectively, of one genotype outstripping another. Over and above, the diagonal line from the top left corner to the bottom right corner, representing comparisons of a genotype with itself, is always green (0 probability)(Fig. 5), whereas off-diagonal cells reveal the probability of the genotype on the x-axis (Compared) outperforming the genotype on the y-axis at interaction cell (references weak performer). In this empirical study, NK-6607 consistently outperformed the other hybrids in the red band (RML95/RML96, Rajkumar, RH-10, MM9440, MM9488, Delta3333, MM2424, MM2122, etc.). Similarly, DKC-9149 was superior to RML95/RML96, Rajkumar, and Delta3333 with the chance of close to 100%. Likewise, RML95/RML96< Rajkumar(65%)< Delta 333(63%)< Delta-9999(53%) superior performance compared with the y- axis cell interacted hybrids, as indicated by the predominance of red bands in their rows.

Fig. 5

Heatmap depicts the probability of GENs in the x-axis being superior to those in the y-axis

Probability of superior stability

The hybrid with a low G*E interaction variance was considered stable. When the genotype was assessed by location factors, this lollipop plot revealed Rajkumar and Delta 90V90 as the most stable hybrids. Moreover, three hybrids, RMH 1899 Super, Delta 3333, and MM 2424, demonstrated comparable stabilities, with a probability stability coefficient of 0.61. However, a significant number of hybrids, including DKC 9144, GK 3226, RMH 567, DKC 9149, GK 3255, Rampur Hybrid-6, RH-10, GK 3156, and MM 2050, exhibited poor stability because their probability coefficient values fell below 0.10 (Fig. 6).

After examining the genotype-by-year factors, TMMH 812, ADV756, Delta 9999, and RH 10 emerged as the top stable hybrids, with probability coefficients that surpassed 0.50. Nevertheless, more than 20 hybrids displayed low stability, with coefficient values less than 0.10, and decreased to zero in cases such as GK 3156/57 and MM 2424.

Fig. 6

Probabilities of Superior Stability for Selected Hybrids by Year and Environment: A Lollipop Plot Analysis

Pairwise probability of superior stability

The pairwise probability of superior performance, as illustrated in (Fig. 7), assumes critical importance during the concluding stages of a breeding program. Because it precedes recommendation and allows the identification of one or more genotypes with performance superior to the variety released commercially. Considering GEI, the x-axis hybrids presented a lower GEI variance than the y-axis hybrids did. Darker blue hues indicate a greater probability of the x-axis genotypes having lower GEI variance, whereas redder hues suggest the opposite. While DKC 9149 exhibits negligible GEI variance compared with RML 95/96, Rajkumar, Delta 333, Delta 5555, NK 7660, Delta 90V90, MM2424, Maan 131, GK 3155, suggesting greater stability across the environment despite the grain yield. A comparison between DKC 9149 and MM 2033 (80% higher variance of GEI) reveals a higher probability of lower GEI variance for DKC 9149, indicating its potential for superior stability across the environment. Conversely, RH 10, Hero no 1, NK7884, and GK 3255 exhibit greater fluctuations in performance, rendering them less suitable for selection because of their pronounced GEI interactions across years and environments. Thus, selected hybrids are more likely to perform well under various environmental conditions. Assuming Rampur Hybrid-6 to be a standard check hybrid, it exhibits a higher probability of variance in GEI when compared to hybrids such as DKC9149, NK6607, and MM2033. However, it demonstrates a lower variance of GEI when contrasted with hybrids RMH9999, Uttam 121, and DKC144 across the different years of evaluation.

Fig. 7

Graphics depicting the pairwise probability of the superior stability among the hybrid’s maize

Joint probability of superior performance and superior stability

Understanding the joint probability of superior performance and stability is crucial for selecting hybrid maize varieties that are both high yielding and adaptable to varying environmental conditions. Breeders aim to develop stable, high-yielding genotypes to ensure the success of recommended varieties. Farmers prioritize minimizing the risk of crop failure and rely on breeders to achieve this goal. Given the joint ability of superior performance and stability analysis, the highest value of coefficient RMH1899 (0.23) was followed by RMH 666 (0.15) and Utam 121 (0.11) in terms of the environement, whereas DKC 9149 (0.19) was followed by NK 6607 (0.19), GK 3254 (0.16), Shann 111 (0.14), Sweety 1 (0.11), ADV 756 (0.11) (Fig. 8)etc.

Fig. 8

Circular Graphics with lollipop plot illustrating the joint probability of superior performance and stability

Probability of superior performance within the environment

A recent analysis of hybrid maize performance in Nepal revealed that several varieties consistently outperformed others in multiple locations. On the basis of the coefficient values, hybrids exhibiting superior performance across different environments in Nepalgunj, Rampur, Parwanipur, and Tarahara were identified. In Nepalgunj, RMH 9999 (0.99), NK 6607 (0.886), and DKC 9144 (0.642) outperformed the others. Similarly, in Parwanipur, DKC 9149 [1], MM 2033 (0.989), and MM 2050 (0.903) demonstrated exceptional yields. Rampur favored ADV 756 (0.985), DKC 9149 (0.806), and MM 2929 (0.97), whereas Tarahara highlighted GK 3254 (0.948), NK 6607 (0.985), Uttam 121 (0.99), and MM 2929 (0.758) as top performers (Fig. 9). In essence, this study provides a foundation for more efficient, resilient, and sustainable maize production.

Fig. 9

Heatmap illustrating the probability of superior performance for different environments in Nepal

Discussion

Circular violin plots illustrate the mean yields of the winter hybrid maize. DKC-9149, with a remarkable 8.8 tons/ha, outperformed all the others, followed closely by Uttam 121 (8.6 tons/ha) and NK6607 (8.6 tons/ha). Conversely, RML95/RML96, Rajkumar, RH 10, etc., struggled, yielding only 5.8 tons per hectare (Fig. 10). On the other hand, Nepal’s average maize yield of 3.06 tons/ha, however, falls significantly short of the global average of 5.06 tons/ha [4]. This discrepancy is strongly attributed to the use of suboptimal data analysis methods in multiple environmental studies, leading to flawed recommendations and ultimately hindering crop performance. Nepal has adopted various traditional variety recommendation models, including AMMI (additive main effect multiplicative interaction) ANOVA (analysis of variance) for spring maize [26], weightage index method [27], Probit regression model [28]. Recently, two novel models, the best linear unbiased prediction (BLUP) and the multi trait stability index (MTSI) (linear mixed model), have been employed for the selection of cold wave-resilient as well as stable high performing hybrid maize varieties [3, 4], climate resilient mungbean genotypes selection [29, 30] but the application of Bayesian probability models has been groundbreaking in Nepal and many developing and developed countries, particularly owing to their recent availability. Bayesian models have been applied to assess probability-based uncertainty, risk estimation, and GEI in bean [31] and Tahiti acid lime cultivars, selection of inbreed maize for Nitrogen used efficiency heritable probability coefficient [32]. In this study, DKC 9141 presented the highest coefficient of probability, suggesting a greater likelihood of performance. This was followed by NK6607, MM2033, RMH9999, GK3254, MM2929, GK3157, Utaam 121, DKC9144, NK6702, Sweety 1, RMH666, and others. In contrast, hybrids with lower probability coefficients and HPDs (genotypic main effects), such as MM2122, Delta 3333, MM9440, RH10, Rajkumar, and RML95/RML96, demonstrated poorer performance (Fig. 4). Despite being considered an elite hybrid in the Godavari Valley region [33], Rajkumar was identified as a poor performer in this study and in previous linear mixed model analyses suggesting complex GEI interaction. Unveiling this disparity in performance among hybrids can be attributed to genetic factors, environmental interactions, and the efficacy of Bayesian modeling techniques in capturing these complexities. The insights gained from such analyses not only guide immediate breeding decisions but also contribute to long-term strategies aimed at enhancing crop resilience and productivity in changing climates. Similarly, joint probability analysis of superior stability and performance revealed DKC 9149, NK6607, GK 3254, Shann111, MM2033, Sweety 1, ADV 756, and GK3155 as the vertex hybrids, which exhibited exceptional probability coefficients and superior stability and performance (Fig. 8). Joint probability analysis highlights the significance of selecting hybrids that not only yield well but also exhibit stable performance across diverse environments, thereby maximizing agricultural productivity and sustainability.

To comprehensively evaluate yield and stability, a slew of researchers have employed MTSI analysis for interaction traits in various crops, including rapeseed [34], drought tolerance screening of wheat [35], elephant grass breeding for bioenergy production [36]and the Guinea yam [37], which was developed by [38] in the “metan” package of R [39]. Furthermore, screening and selecting Tahiti acid lime genotypes through the Markov chain Monte Carlo (McMC) algorithm also uses Bayesian interference to circumvent the flaw of cultivar recommendation [40]. In the context of location-specific genotype recommendations, a MET study, coupled with cutting edge Bayesian probability interference, revealed that in Nepalgunj, the RMH 9999, NK 6607, and DKC 9144 genotypes presented exceptional yields of 8.5, 8.6, and 8.18 tons/ha, respectively, outperforming other varieties. Similarly, in Parwanipur, the DKC 9149, MM 2033, and MM 2050 genotypes demonstrated superior performance, with yields of 8.81, 8.5, and 8.2 tons/ha, respectively. In Rampur, the ADV 756, DKC 9149, and MM 2929 genotypes yielded 8.2, 8.81, and 8.6 tons/ha, respectively. Tarahara favored the GK 3254, NK 6607, Uttam 121, and MM 2929 hybrids, each producing 8.5, 8.6, 8.6, and 8.6 tons/ha, respectively. Findings highlight the critical role of tailored hybrids selection in maximizing maize grain yields/ha across different regions in Nepal, supported by advanced statistical methods that enhance predictive accuracy and reliability in agricultural decision-making. Consequently, the chosen hybrids demonstrate promising potential for success.

Fig. 10

Comparative Mean Yield of Hybrid Maize Cultivars in Nepal (2021–2023)

Conclusion

In Bayesian analysis, evaluating convergence on the basis of samples from the posterior distribution is crucial for GEI analysis. This process aids in understanding crop performance, adaptability, and consistency across diverse environments, thereby revealing their potential for success or failure under unpredictable conditions. By employing probabilistic Bayesian models for winter hybrid maize, this study offers a nuanced assessment of hybrid performance, highlighting the likelihood of superior outcomes, pairwise performance as well as stability probabilities among hybrids both across and within research stations. On the basis of the analysis, DKC9149 and NK6607 emerged as top-performing hybrids, demonstrating exceptional performance and stability across various environments. In Nepalgunj, RMH 9999 and NK 6607 excelled, whereas DKC 9149 and MM 2033 stood out in Parwanipur. Rampur favored ADV 756, DKC 9149, and MM 2929, whereas Tarahara highlighted GK 3254, NK 6607, Uttam 121, and MM 2929. Therefore, we contend that the use of a Bayesian probabilistic model enhances cereal crop recommendations by providing clearer and more precise interpretations of elite genotype performance and persistence.