Abstract
The objective of this study was to evaluate the adaptability and phenotypic stability of watermelon hybrids across diverse environments in the semiarid region of Rio Grande do Norte, Brazil, using Toler’s simple linear and piecewise linear models under both frequentist and Bayesian frameworks. The Bayesian analysis employed noninformative uniform priors and yielded classifications consistent with the frequentist results, reinforcing the robustness of the classification results. Genotypes in group A (WM-5, WM-6, WM-9, and WM-10) exhibited low sensitivity to adverse environmental conditions and high responsiveness in favorable environments, along with high phenotypic stability. The results demonstrate that Bayesian approaches are a reliable complement to conventional analyses, thereby supporting breeding strategies and cultivar recommendations for semiarid regions while enhancing the robustness of frequentist inference through explicit quantification of uncertainty.
Keywords:
Bayesian inference; Citrullus lanatus; plant breeding
INTRODUCTION
Watermelon [Citrullus lanatus (Thunb.) Matsum. & Nakai] is one of the most widely cultivated and consumed fruits worldwide, with significant economic and social importance. In Brazil, watermelon is cultivated in nearly all regions, with the Northeast region standing out as the primary production area, contributing substantially to both domestic consumption and exports. In 2024, national watermelon production exceeded 1.7 million Mg, positioning Brazil among the leading global producers (IBGE 2024). Export-oriented watermelon production in the Mossoró-Assu agricultural pole is conducted primarily by companies employing advanced technology and single-cross hybrid seeds.
To meet farm requirements, seed companies continually develop new watermelon cultivars and release new hybrids annually. However, evaluating new hybrids for high yield and fruit quality is essential before recommending cultivars for export (Aragão et al. 2015). Watermelon hybrids in Rio Grande do Norte are grown under diverse environmental conditions, resulting in significant genotype × environment (G × E) interaction. This interaction influences phenotypic expression, such that cultivars with superior performance in one environment may not exhibit equivalent performance in another. This complicates cultivar recommendation by breeders (Lima et al. 2025).
Several approaches can be used to reduce the effects of G × E interaction, including the selection of cultivars with superior adaptability and phenotypic stability (Nunes et al. 2006, Nunes et al. 2011). In plant breeding, adaptability and stability are fundamental concepts. Adaptability refers to the capacity of a genotype to respond positively to improvements in environmental conditions, whereas stability refers to the consistency of its performance across varying environments (Cruz et al. 2012, Oliveira et al. 2019). Genotypes with high stability are less influenced by environmental variability and are particularly suitable for heterogeneous or unpredictable conditions. Stability may be static, defined as consistent performance across environments, or dynamic, defined as performance that changes proportionally with environmental quality (Lin et al. 1986, Becker and Leon 1988). Evaluation of both traits together provides a more comprehensive assessment of genotype performance in multi-environment trials. Thus, statistical models capable of accurately characterizing genotypic responses to environmental variation are required for reliable cultivar recommendation (Rosa et al. 2023).
The model proposed by Toler and Burrows (1998) advances the analysis of G × E interaction through a piecewise bilinear regression approach that estimates distinct genotype responses in unfavorable and favorable environments. This model yields parameters that characterize mean performance, sensitivity under unfavorable conditions, and responsiveness in favorable environments, providing a robust method for evaluating G × E interaction relative to conventional linear regression models (Kvitschal et al. 2009, Rosa et al. 2023). However, application of this model to various crops, including watermelon, remains limited, indicating potential for further investigation of its benefits in genotype selection.
The Bayesian approach to the Toler model offers a promising alternative, allowing the incorporation of prior information, estimation of posterior parameter distributions, and more accurate quantification of uncertainty (Rosa et al. 2023, Lima et al. 2025). Coherent and informative priors can improve parameter estimation; in contrast, poorly specified priors may bias results, whereas vague or noninformative priors generally produce estimates similar to those obtained using the conventional frequentist approach. Additionally, the Bayesian framework provides more precise estimates of standard deviations and yields probabilistic measures of uncertainty, such as highest posterior density intervals (HPD intervals), also referred to as highest density intervals (HDI). Furthermore, the Bayesian approach is particularly advantageous in situations with few environments or limited sample sizes, which are common in plant breeding trials.
Thus, the objective of this study was to apply Bayesian estimation to assess adaptability and stability parameters of watermelon genotypes using the Toler model, aiming to identify genotypes that combine high mean performance with consistent responses in both unfavorable and favorable environments.
MATERIAL AND METHODS
Experimental details
Fifty watermelon hybrids were evaluated in four municipalities representative of the Mossoró-Assú agricultural pole (Mossoró-Assú Agropolo) in Rio Grande do Norte, Brazil: Mossoró, Baraúna, Assú, and Ipanguaçu. Evaluations were conducted in these four municipalities over three consecutive years (2017, 2018, and 2019), totaling 12 environments.
In all environments, the crop was irrigated by drip irrigation with fertigation, with row spacing of 3.0 m and dripper spacing of 0.4 m. Fertilizers were applied according to soil analysis recommendations for each site. Other crop management practices followed the state recommendations for watermelon cultivation (Silva et al. 2008). The experiments were conducted in a randomized complete block design with three replications. Each plot consisted of three rows, each 4.8 m in length, with one plant per planting hole. The evaluation area comprised the central row, excluding the border plants at each end.
The simple and complex components of the genotype × environment (G × E) interaction were estimated using the method proposed by Cruz and Castoldi (1991). Decomposition of the hybrid × environment interaction was performed using the GENES software (Cruz 1997).
Toler’s model
This method assesses genotype adaptability and stability across an environmental gradient (), extending the linear regression model of Eberhart and Russell (1966) to a piecewise linear model with two slopes ( and ) connected at a breakpoint. In this study, the version described by Nascimento et al. (2020) was used:
where is the yield of genotype in environment , is the intercept representing mean performance, and are the slopes before and after the breakpoint at , is an indicator variable (1 if and 0 if ), and is the residual error.
After model fitting, hybrids were classified according to their response patterns to environmental variation. When the hypothesis is rejected, two response patterns are defined: group A, if , in which the genotype exhibits a convex response (doubly desirable); or group E, if , in which the genotype exhibits a concave response (undesirable).
If the hypothesis is not rejected, the response is linear, and a simplified Eberhart and Russell (1966) model is fitted:
where is the yield of genotype in environment , is the mean response in the average environment (μ′ⱼ = 0), is the common slope, and is the residual error.
Within this linear framework, three classifications are possible: group B, if , indicating responsiveness in favorable environments; group C, if , indicating stable performance across environments; and group D, if , indicating responsiveness in unfavorable environments. In addition to evaluating genotype adaptability through the slope parameters and , assessment of phenotypic stability is necessary, as it quantifies the consistency of performance across environments.
In this study, phenotypic stability was estimated using the variance of deviations (, calculated as , where is the mean square error of the fitted model, is the residual mean square from the analysis of variance (ANOVA), and is the number of replications. This variance quantifies the unexplained variation in genotypic responses after fitting the linear regression model. Occasionally, this estimate may be negative and therefore uninterpretable; in such cases, was assigned, indicating maximum phenotypic stability, following Nascimento et al. (2020) and de Lima et al. (2025).
Bayesian inference
Bayesian inference, based on Bayes' theorem, was used to estimate parameters by integrating prior information with observed data, resulting in robust posterior distributions. Initially, a frequentist analysis was conducted using piecewise regression with a breaking point at zero, fitted individually for each genotype and compared with the simple linear model using t-tests and F-tests to characterize the expected behavior of the parameters. These data were used as prior information and incorporated into the Bayesian inference. The likelihood function was specified from the model residuals, assuming a normal distribution with the mean adjusted to the model predictions and variance estimated from the residual standard deviation. Based on these results, uniform priors were specified as for and , and for ; these priors were applied to all genotypes, as derived from the frequentist analysis, and bounds were also defined for the environmental quality index ( to ensure consistent posterior estimates.
A total of 3,000 posterior samples per chain were retained for inference (draws = 3,000) across four independent chains (chains = 4). Thus, with 3,000 post-burn-in samples from each chain, the total posterior sample size was 12,000. The number of tuning iterations was set to 1,000 (tune = 1,000), such that the first 1,000 iterations of each chain were discarded as burn-in, a standard practice in Markov chain Monte Carlo (MCMC) sampling.
Initial parameter estimates, used as prior information, were obtained using the sm.OLS function from the statsmodels package in Python (Seabold and Perktold 2010). Bayesian regression was implemented using the PyMC package (Abril-Pla et al. 2023). In contrast to the ordinary least squares (OLS) method, which yields point estimates, the Bayesian approach explicitly models parameter uncertainty and provides full posterior distributions. MCMC sampling was used to fit the model to the data and generate posterior distributions that quantify parameter uncertainty.
Because Toler’s model can yield extreme values for and in unfavorable environments, wide uniform priors (−1000, 1000) were initially specified. Sensitivity analyses indicated that narrower uniform priors produced essentially identical posterior estimates while substantially reducing computational cost. Therefore, the final prior bounds were defined as the narrowest interval that ensured MCMC convergence without affecting posterior inference. Convergence was considered achieved when the maximum absolute change in posterior parameter estimates across iterations fell below 10−4.
The original decision rules proposed by Toler were retained in the Bayesian framework; however, genotype classification relied on the 95% highest posterior density intervals (HPD intervals) of the posterior distributions rather than on point estimates. Specifically, genotypes were assigned to group A or E when the HPD intervals for both β₁ and β₂ excluded the value 1 and satisfied the conditions (convex response) or (concave response), respectively. When the HPD intervals for and overlapped and provided no evidence against the hypothesis , the simplified linear model was adopted , and groups B, C, and D were assigned according to whether the HPD interval for the common slope lay entirely above 1, included 1, or lay entirely below 1. This approach provides a probabilistic and uncertainty-aware classification of genotypes into Toler’s groups.
RESULTS AND DISCUSSION
Analysis of variance showed a significant genotype × environment (G × E) interaction effect (p < 0.01), whereas the main effects of genotype and environment were not significant (Table 1). Similarly, Silva et al. (2008) reported a significant G × E interaction when evaluating seven cultivars across six environments in Rio Grande do Norte. Dia et al. (2016a) also reported a significant G × E interaction when evaluating 40 watermelon genotypes across 20 locations in the United States. In contrast, Correa et al. (2020) observed no significant interaction when evaluating 20 cultivars in two locations in Texas and attributed this result to similar environmental conditions during plant growth and comparable field management practices across environments.
Analysis of variance for fruit yield of watermelon genotypes evaluated across 12 environments in Rio Grande do Norte, Brazil
The presence of a significant G × E interaction indicates inconsistent genotype performance across evaluation environments. The G × E interaction results from two components. The first component, termed the simple or scale component, results from differences in variance among genotypes across environments, whereas the second, termed the complex component, is associated with low genetic correlations among genotypes across environments (Cruz and Castoldi 1991, Olivoto et al. 2019). In the present study, the complex component predominated (64.56%), with estimated genetic correlations below 0.40. Similarly, Dia et al. (2016b) identified a predominance of crossover G × E interaction for fruit yield when defining mega-environments for evaluating watermelon genotypes in the United States. The predominance of crossover G × E interaction complicates genotype recommendation because it leads to rank changes across environments (Cavalcante et al. 2025).
In the frequentist analysis using Toler's model (Table 2), most genotypes were classified in group C, characterized by a linear response to the environmental quality index without deviation from the expected response (Bradshaw 1965, Toler and Burrows 1998). These genotypes respond predictably to environmental variation but do not exhibit pronounced adaptive plasticity. In contrast, genotypes WM-5, WM-6, WM-9, and WM-10 were classified in group A, exhibiting a convex response pattern considered doubly desirable, with low sensitivity in unfavorable environments and high responsiveness in favorable environments. These genotypes represent promising material for breeding programs because they combine phenotypic stability under limiting conditions with high yield potential in favorable environments. Conversely, genotypes WM-2, WM-11, WM-13, and WM-15 were assigned to group E, exhibiting a concave response pattern considered undesirable due to high sensitivity in unfavorable environments and limited responsiveness in favorable environments (Toler and Burrows 1998).
Figure 1 illustrates the response functions (concave, convex, or linear) of genotypes as determined by the model under the Bayesian framework. Considering the environmental quality index (), a positive indicates an environment with above-average yield, whereas a negative indicates a below-average-yield environment.
Regression plot based on the Toler model applied to watermelon data. The x-axis represents the environmental quality index (μⱼ), and the y-axis represents the predicted genotype mean yield. The colored lines represent the fitted responses for each genotype obtained using Bayesian inference. The letters preceding genotype identifiers indicate the classification of each genotype based on the estimated model parameters. The black line represents the overall mean yield across environments and genotypes. The gray dashed line indicates the threshold separating unfavorable unfavorable (μⱼ < 0) and favorable (μⱼ > 0) environments.
The Bayesian analysis confirmed the classifications from the frequentist analysis, reinforcing the robustness of the results (Table 3). In addition to replicating classifications into groups A, C, and E, the Bayesian analysis provided 95% highest posterior density (HPD) intervals for the parameters of each genotype, allowing quantification of the uncertainty associated with the estimates. In all cases, the posterior distributions were consistent with the frequentist interpretations while providing greater depth in parameter estimation, as shown by the distinct separation between and in genotypes classified in group A and the corresponding differences observed in group E. The coefficient of determination () ranged from 0.10 to 0.86 and was higher for genotypes in the extreme groups (A and E), indicating a better fit of Toler's model to the responses of these genotypes. Similar patterns have been reported in recent studies using comparable approaches (Oliveira et al. 2020, Partelli et al. 2022, Lima et al. 2025, Silva et al. 2025). Nascimento et al. (2020), however, reported that informative priors can narrow credible intervals compared with non-informative priors.
The preference for the Bayesian framework with this dataset is supported not only by its agreement with the frequentist Toler classifications but also by its capacity to express results as probabilistic statements about genotype performance. In contrast to the frequentist approach, which provides single-point estimates for α i , β 1i , and β 2i , the Bayesian approach yields full posterior distributions, facilitating transparent and coherent quantification of uncertainty. This capability is particularly valuable in multi-environment trials involving few environments, in which uncertainty can substantially affect genotype rankings and recommendations. Furthermore, Bayesian evidence [ and assessment of posterior chain stability provide complementary diagnostics unavailable in the frequentist implementation of Toler’s model. By incorporating these features, the Bayesian analysis not only replicates the conventional classification but also enhances it through explicit quantification of uncertainty, thereby increasing confidence in the identified adaptability patterns and supporting recommendations for contrasting environments.
Figure 2 shows trace plots from Markov chain Monte Carlo (MCMC) sampling for the parameters of genotype WM-1, together with the corresponding marginal posterior distributions, illustrating chain convergence. All other watermelon genotypes exhibited similar chain patterns.
Trace plots of the parameters α (intercept), β 1 (regression coefficient for unfavorable environments), and β 2 (regression coefficient for favorable environments), obtained by Bayesian sampling for watermelon genotype WM-1. Each line represents iterations of a Markov chain for a given parameter, and the histograms show the marginal posterior distributions. The stable traces and good mixing among chains indicate satisfactory convergence. All genotypes exhibited similar convergence patterns.
The distribution of genotypes across yield groups and adaptability patterns reinforces the importance of joint analysis of yield, adaptability, and phenotypic stability, consistent with recent studies on watermelon and other crops (Dia et al. 2016a, Dia et al. 2016b, Lima et al. 2025). Despite the high overall mean yield (36.81 Mg ha−1), only a few genotypes exhibited the doubly desirable response (Group A). This pattern is consistent with previous reports indicating that the combination of high yield, stability, and wide adaptability is uncommon and remains one of the main challenges in plant breeding (Dutra Filho et al. 2021).
Overall, the results obtained for watermelon align with those previously reported for melon, in which a predominance of genotypes classified in group C and few genotypes in the extreme groups (A and E) were observed (Rosa et al. 2023, Lima et al. 2025). The agreement between the frequentist and Bayesian analyses supports the applicability of Toler's approach. The primary advantage of the Bayesian method lies in its capacity to quantify uncertainty in parameter estimates and to provide additional diagnostics, such as credible intervals and Bayesian evidence (Van de Schoot et al. 2021). These features improve interpretation of the results and increase the reliability of recommendations for contrasting environments.
Moreover, the Bayesian framework allowed a more informative assessment of genotypic performance by providing direct access to posterior distributions for genotypic parameters and associated uncertainty. These results, detailed in Table 3, include posterior means and 95% HPD intervals for the primary parameters (α i , β 1i , and β 2i ), contrasts between β 1i and β 2i , the common slope (β common ), Bayesian evidence [, and genotype classifications. This information supports clearer identification of high-yielding and widely adaptable genotypes, as the HPD intervals indicate the probability that each parameter falls within a given range. For example, the non-overlapping HPD intervals of β 1i and β 2i in group A genotypes provide strong support for the convex response pattern (doubly desirable), thereby improving the reliability of genotype recommendations for contrasting environments.
Thus, the significant G × E interaction, combined with the identification of distinct genotype groups (A, C, and E), supports the application of Toler's model for evaluating watermelon genotypes. This approach distinguishes genotypes with wide adaptability and phenotypic stability (group A) from those with average and predictable responses (group C) and from those with undesirable responses (group E), thereby supplying essential information for breeding programs and cultivar recommendations across diverse growing environments.
CONCLUSION
The robustness of genotype classifications using Toler's model was confirmed by the agreement between frequentist and Bayesian analyses. The Bayesian approach not only validated the identification of group A as the most promising genotypes but also quantified the associated uncertainty through posterior distributions, thus increasing confidence in recommendations for contrasting environments.
Genotypes in group A, particularly WM-5 and WM-6, emerge as the most promising for recommendation across contrasting environments, exhibiting high responsiveness in favorable environments and low sensitivity in unfavorable environments.
ACKNOWLEDGEMENTS
The authors thank the Postgraduate Program in Plant Science (PPGFito) and the Federal Rural University of the Semi-Arid Region (UFERSA) for providing the facilities and support necessary to conduct this experiment. This study was funded by the Coordination for the Improvement of Higher Education Personnel (CAPES) and the National Council for Scientific and Technological Development (CNPq), Brazil. M.M.F. de Lima acknowledges a fellowship from CNPq (process no. 152681/2024-8).
Data Availability Statement
The datasets generated and/or analyzed during this research are freely accessible in a repository on GitHub, available at https://github.com/MaxsuelMFLima/Adaptability-Stability-in-Watermelon.
REFERENCES
- Abril-Pla O, Andreani V, Carroll C, Dong L, Fonnesbeck CJ, Kochurov M, Zinkov R2023 PyMC: a modern, and comprehensive probabilistic programming framework in PythonPeerJ Computer Science 9:e1516
- Aragão FASD, Nunes GHDSS, Queiróz MAD2015 Genotype x environment interaction of melon families based on fruit quality traitsCrop Breeding and Applied Biotechnology 15:79-86
- Becker HC, Leon JI1988 Stability analysis in plant breedingPlant Breeding 101:1-23
- Bradshaw AD1965 Evolutionary significance of phenotypic plasticity in plantsAdvances in Genetics 13:115-155
- Cavalcante JG, Martins AF, Silva EMD, Moreira SDA, Nunes EWLP, Negreiros AMP, Melo SB, Nunes GHDSS2025 Impact of environmental covariates on genotype-environment interactions in a semi-arid region of Rio Grande do Norte State, BrazilActa Scientiarum Agronomy 47:e70904
- Correa E, Malla S, Crosby KM, Avila CA2020 Evaluation of genotypes and association of traits in watermelon across two southern Texas locationsHorticulturae 6:67
- Cruz CD1997 Programa GENES - Aplicativo computacional em genética e estatística. Editora UFV, Viçosa MG, pp. 442
- Cruz CD, Castoldi FL1991 Decomposição da interação genótipos x ambientes em partes simples e complexaRevista Ceres 38:219
- Cruz CD, Regazzi AJ, Carneiro PCS2012 Modelos biométricos aplicados ao melhoramento genético. Editora UFV, Viçosa, 514p.
- Dia M, Wehner TC, Hassell R, Price DS, Boyhan GE, Olson S, Juarez B2016a Value of locations for representing mega‐environments and for discriminating yield of watermelon in the USCrop Science 56:1726-1735
- Dia M, Wehner TC, Perkins-Veazie P, Hassell R, Price DS, Boyhan GE, Juarez B2016b Stability of fruit quality traits in diverse watermelon cultivars tested in multiple environmentsHorticulture Research 3:16066
- Dutra Filho JA, Calsa Júnior T, Simões Neto DE, Souto LS, Souza ADS, Luna RG, Costa MLL2021 Genetic divergence for adaptability and stability in sugarcane: Proposal for a more accurate evaluationPLOS One 16:e0254413
- Eberhart ST, Russell W1966 Stability parameters for comparing varieties 1Crop Science 6:36-40
-
IBGE - Instituto Brasileiro de Geografia e Estatística2024 Produção agrícola municipal. IBGE, Rio de Janeiro. Available at <Available at https://sidra.ibge.gov.br/pesquisa/pam/tabelas >. Accessed on September 14, 2025.
» https://sidra.ibge.gov.br/pesquisa/pam/tabelas - Kvitschal MV, Vidigal Filho PS, Scapim CA, Gonçalves-Vidigal MC, Sagrilo E, Pequeno MG, Rimoldi F2009 Comparison of methods for phenotypic stability analysis of cassava (Manihot esculenta Crantz) genotypes for yield and storage root dry matter contentBrazilian Archives of Biology and Technology 52:163-175
- Lima MMF, Martins AF, Silva AS, Alves CS, Nunes EWLP, Silva EM, Nunes GHS2025 Adaptability and stability of melon genotypes through linear regression and segmented linear regression models: a bayesian approachEuphytica 221:80
- Lin CS, Binns MR, Lefkovitch LP1986 Stability analysis: where do we stand? Crop Science 26:894-900
- Nascimento M, Nascimento ACC, Silva FF, Teodoro PE, Azevedo CF, Oliveira TRA, Carvalho LP2020 Bayesian segmented regression model for adaptability and stability evaluation of cotton genotypesEuphytica 216:30
- Nunes GH, Santos Júnior H, Grangeiro LC, Bezerra Neto F, Dias CT, Dantas MS2011 Phenotypic stability of hybrids of Galia melon in Rio Grande do Norte state, BrazilAnais da Academia Brasileira de Ciências 83:1421-1434
- Nunes GHDSS, Madeiros AEDSS, Grangeiro LC, Santos GM, Sales Júnior R2006 Estabilidade fenotípica de híbridos de melão amarelo avaliados no Pólo Agrícola Mossoró-AssuPesquisa Agropecuária Brasileira 41:1369-1376
- Oliveira LADA, Cardoso EDA, Ricarte ADO, Martins AF, Costa JMD, Nunes GHDSS2019 Stability, adaptability and shelf life of Cantaloupe melon hybrids. Revista Brasileira de Fruticultura 41: e-418.
- Oliveira TRA, Carvalho HWL, Nascimento M, Costa EFN, Oliveira GHF, Gravina GA, Amaral Júnior AT, Carvalho Filho JLS2020 Adaptability and stability evaluation of maize hybrids using Bayesian segmented regression modelsPLOS One 15:e0236571
- Olivoto T, Lúcio AD, Silva JA, Marchioro VS, Souza VQ, Jost E2019 Mean performance and stability in multi‐environment trials I: combining features of AMMI and BLUP techniquesAgronomy Journal 111:2949-2960
- Partelli FL, Silva FA, Covre AM, Oliosi G, Correa CCG, Viana AP2022 Adaptability and stability of Coffea canephora to dynamic environments using the Bayesian approachScientific Reports 12:11608
- Rosa JC, Uhdre RS, Faria MV, Pinto RJB, Rossi RM, Scapim CA2023 A new Bayesian approach to the Toler model for evaluating the adaptability and stability of genotypesCrop Breeding and Applied Biotechnology 23:e440523111
- Seabold S, Perktold J2010 Statsmodels: econometric and statistical modeling with pythonSciPy 7:92-96
- Silva EVP, Davide LMC, Gianlup C, Oliveira WJS, Oliveira LA, Silva AQ, Khan S2025 Assessing the stability and adaptability of maize hybrid yield with the Bayesian AMMI modelEuphytica 221:43
- Silva JR, Sousa Nunes GHD, Negreiros MZ, Torres JF, Dantas MSSM2008 Interação genótipo x ambiente em melancia no estado do Rio Grande do NorteRevista Caatinga 21:95-100
- Toler JE, Burrows PM1998 Genotypic performance over environmental arrays: a non-linear grouping protocolJournal of Applied Statistics 25:131-143
- Van de Schoot R, Depaoli S, King R, Kramer B, Märtens K, Tadesse MG, Yau C2021 Bayesian statistics and modellingNature Reviews Methods Primers 1:1
Edited by
-
SCIENTIFIC EDITOR:
Luiz Antônio dos Santos Dias




