Coffee is one of the most globally beverages presenting importance in terms of its potential health, socioeconomic, and economic effects (Porto et al., 2019). These effects drive breeding programs worldwide to develop high-yielding, adaptable cultivars delivering superior bean quality (Barbosa et al., 2019). However, traditional selection methods based on phenotypic observations of the plants or their family history (pedigree) are expensive and time consuming, especially for perennial crops as coffee.

An alternative approach denoted genomic selection (GS) has been used as a successful tool in genetic improvement (Meuwissen et al., 2001). GS helps increase genetic gain per generation by allowing for earlier selection through improved prediction of the potential of individual plants (Daetwyler et al., 2013; Nascimento et al., 2019; Voss-Fels et al., 2019). In the field of coffee breeding, GS has been utilized with the dual aim of accelerating genetic gain through early selection and improving prediction accuracy (Sousa et al., 2019; Alkimim et al., 2020; Sousa et al., 2021; Coelho de Sousa et al., 2022).

Among several prediction models, Genomic Best Linear Unbiased Prediction (GBLUP) is the most widely used approach for genomic prediction due to its advantages (Zhang et al., 2021). Compared to other parametric methods, GBLUP allows to accurately estimate narrow-sense heritability (Li et al., 2019) and presents higher computational efficiency (Hernandez et al., 2020). GBLUP modeling is also flexible. It can be modified to incorporate additional genetic information beyond the typical single-nucleotide polymorphism (SNP) markers. Specifically, this modeling allows to account for non-additive genetic effects, environmental factors, and even genotype-by-environment interactions, enriching the analysis and potentially improving prediction accuracy (Jarquín et al., 2014).

In the Artificial Intelligence Era, the interest in semi- and non-parametric methods for GS is increasing (Larkin et al., 2019; Coelho de Sousa et al., 2022; Seyum et al., 2022). These approaches, such as Artificial Neural Networks and Decision Trees, do not require prior assumptions about the relationships between inputs (SNP markers) and the output (phenotypic observations), allowing great flexibility to handle complex non-additive effects, such as dominance and epistasis (McKinney et al., 2006; Abdollahi-Arpanahi et al., 2020; Coelho de Sousa et al., 2022). In general, despite their potential, these approaches do not outperform the traditional parametric methods (e.g., GBLUP, Bayesian Alphabet - Gianola et al., 2009) used to predict the genetic merit of individuals (Liang et al., 2021).

Aiming to enhance predictive ability (PA), Ensemble Learning (EL) combines predictions from multiple models (base learners) into a single prediction (meta-learner) (Mendes-Moreira et al., 2012; Ganaie et al., 2022; Mienye and Sun, 2022). This approach leverages the strengths of diverse models to potentially generate more robust results compared to relying on a single learner (Liang et al., 2021; Kalule et al., 2023). In the context of GS, EL has found application through methods such as Random Forest (RF) and Bagging (Bag) (Xu et al., 2019; Abdollahi-Arpanahi et al., 2020; Sousa et al., 2021; Costa et al., 2022). These methods, categorized as Homogenous Learning (HL), utilize a single framework to produce a single prediction value. Conversely, the Stacking Ensembles Learning (SEL) approach combines predictions from diverse methods, potentially outperforming HL (Mendes-Moreira et al., 2012). SEL has seen success in GS, improving PA in Chinese Simmental cattle, Dutch cattle, and pine (Liang et al., 2021), achieving higher accuracy than GBLUP for most evaluated traits.

Despite being interesting, EL approach arises some issues that needs to be considered. First, since the same individuals are used to fitting the model(s) in the EL approach, it is expected the existence of correlation between the predictions derived from the different methods. This well-known statistical problem is referred as multicollinearity (Montgomery et al., 2021) and causes high variability of the estimated effects. The second issue is related to which dataset should be used to fit the meta-learner. One option is to use directly the predicted values derived from the base learners. In this case, the simple mean or some regression model that accounts for multicollinearity problem (e.g., Ridge Regression) can be implemented to makeup a single prediction. An alternative option also could consider combining the predicted values with the genomic covariates (i.e., SNP markers|predicted values from the base learners) with the previous training data as new inputs. In this regard, in addition to the multicollinearity, the course of the dimensionality is another issue to consider mainly because these are covariables of different type.

Liang et al. (2021) used the predicted values derived from a multiple regression model as meta-learner. These authors obtained good results to improve PA compared to the conventional genomic prediction models on three different datasets. However, the use of an expanded training data augmented by SNP markers could be beneficial to further enhance the PA, and it emerges as an interesting approach. In this case, a model that addresses both multicollinearity and dimensionality problems should be used. One of the possible solutions can be considered to use a two-kernel GBLUP model as the meta-learner model. Another approach to evaluate is to consider only the predicted values provided by the best base learner models.

To date, no research has applied SEL to improve the prediction accuracy of important traits in coffee cultivars. This approach presents potential for coffee breeding, as it has been shown to outperform standard methods in other applications (Liang et al., 2021; Mohammed and Kora, 2023). By combining the strengths of multiple prediction models, SEL could lead to more reliable and accurate identification of valuable genetic traits in coffee plants, accelerating the development of superior coffee varieties.

In light of the mentioned points, the objective of this study was to use and evaluate the SEL to improve PA of important traits in Coffea Arabica. For that, the GBLUP, multivariate adaptive regression splines (MARS), Quantile Random Forest (QRF), and RF models were used as the base learner. Several approaches were considered as the meta-learner to construct the SEL framework. Specifically, the expanded- and non-expanded datasets were used for training. In addition, models that account for multicollinearity (Ridge Regression) and multicollinearity and dimensionality jointly (GBLUP) were also implemented.

In this study, we used the SEL approach to improve PA of four important traits in Coffea Arabica. Two of these traits are associated with the productivity (YL and NF) and the remaining two with disease resistance (LM and Cer). The population under study is comprise of 195 genotypes of Coffea Arabica genotyped for 5,970 SNP markers. We compared the PA of different approaches used in the Level 1 of the SEL to the results obtained with the base learners [GBLUP, MARS considering degrees equal to 1, 2, and 3, QRF considering nine quantiles (from 0.1 to 0.9, every 0.1) and RF]. Since the GBLUP is the most implemented prediction model (Zhang et al., 2021), their results were used as benchmark. The PA of the different approaches was assessed using a CV scheme repeated 10 times. The Spearman’s correlation and the agreement (based on the top 10%) coefficients between the GEBV values of the different models were also estimated. The genetic parameters were also estimated for the evaluated traits (YL, NF, LM, and Cer).

The heritability estimates for YL (0.55), NF (0.49), LM (0.30), and Cer (0.38) were consistent with those reported in the literature for this specie and same traits. Specifically, the heritability estimates varied between (0.1–0.74) [Alkimim et al. (2020) and Alemayehu, 2019], (0.30–0.55) [Gokavi et al., 2023 and Weldemichael et al., 2017], (0.30–0.51) [Chrigui et al., 2020 and Ferrão et al., 2023], and (0.09–0.61) Alkimim et al., 2021 and Ferrão et al., 2023] for the YL, NF, LM, and Cer, respectively. Although the Pearson correlation between YL and the disease resistance traits were not statistically significant, a significant and positive genetic correlation (0.52) was obtained between YL and NF.

The machine learning methods as base learners have been already used in genomic prediction (Long et al., 2011; Lenz et al., 2019; Montesinos-López et al., 2019; Coelho de Sousa et al., 2022; Costa et al., 2022). However, generally these methods do not outperform significantly the traditional genomic prediction approach based on parametric models such as GBLUP and Bayesian Alphabet (Liang et al., 2021).

Liang et al. (2021) used the SEL for improving PA in three real datasets on average by 7.70%, compared to GBLUP. The SEL uses predicted values from different machine learning implementations to obtain a single prediction value. These authors integrated/combined the results of three machine learning implementations (Support Vector Machine, Kernel Ridge Regression and Elastic Net) to compute the GEBVs.

In contrast to Liang et al. (2021), in our study, the GBLUP approach was used as one of the base learner methods for the SEL too. The GBLUP was considered since it is widely used for genome prediction (Zhang et al., 2021) due to its reduced computational demand and simplicity (Hernandez et al., 2020) compared to the other parametric methods (e.g., Bayesian Alphabet, Gianola et al., 2009). The MARS that allows automatically selecting and modeling nonlinear relationships and interaction effects of the predictor variables was also considered as base-learner method (Costa et al., 2022). In addition, the RF (James et al., 2023) and the QRF (Meinshausen, 2006) were also set as base learner methods. Specifically, the RF is a machine learning approach used to increase the predictive power and reduce the variance of the predicted values by averaging uncorrelated quantities (James et al., 2023). The QRF combines interesting characteristics from RF and Quantile Regression (QR) approaches. QR (Koenker and Bassett, 1978) allows fitting all the portions of the probability distribution of the trait, enabling a more complete picture of the conditional distribution than a single estimate of the center (Briollais and Durrieu, 2014; Nascimento et al., 2019).

Overall, for each evaluated trait (YL, NF, LM, and Cer), a different model presented the highest PA. These results show that there is not a single approach that outperforms the others in the evaluated data sets. Also, it could be case of a model performing better than the others in a given dataset but another model could perform better in a similar dataset (James et al., 2023). For example, the MARS models with 2° and 3° (model with interactions) presented higher PA for YL. These results are in line with those obtained by Coelho de Sousa et al. (2022). Using artificial neural networks to predict the genetic merit of genotypes of Coffea canephora these authors showed a larger dominance markers effect for YL when compared to the GBLUP additive dominant model based on additive marker effects.

Another interesting result was obtained modeling the QRF where the distribution of the adjusted phenotypic values for YL and NF. These presented a higher concentration in the first quantiles (0.1–0.3) ( Supplementary Figure S2 ) and positive skewness ( Supplementary Figure S3 ). For YL, the best models were the QRF0.1 and QRF0.2, and QRF0.3 for NF ( Supplementary Figure S2 ). The distribution of the adjusted LM values presented tree major portions ( Supplementary Figure S4 ). The QRF modeling was able to distinguish these three different groups ( Figure 2 ). Finally, since the distribution of adjusted Cer phenotypic values did no present a specific pattern to highlight ( Supplementary Figure S5 ) all of the QRF models present similar PA ( Figure 2 ). A similar trend was shown by Nascimento et al. (2017). These authors found that the Quantile Regression approach outperform the traditional genomic prediction methods of not normal distributed traits.

An interesting approach to address the non-normality assumption is using multiple models to conduct the predictions, and then combine the predicted values to makeup a single prediction through the SEL approach. In general, the SEL outperforms the methods based on base learners only (Liang et al., 2021; Kandel et al, 2021; Kalule et al., 2023). In our study, the SEL approach outperformed all base learner methodologies ( Figures 2 , 3 ). However, it is important to emphasize that these results were observed by those SEL models that used only the predicted values provided from the base learners with PA higher than mean of all the base learners. Specifically, the Stacking Mean Best (SMB) presented the highest PA for all of the evaluated traits. The average of the predictions from several fitted models has been successfully implemented with Bagging and RF approaches (Breiman, 1996 and Breiman, 2001). The SEL approach allow to use several models to combine the predicted values, for example, XGBoost (Ghasemieh et al., 2023), Penalized methods (Kalule et al., 2023), Linear Regression (Liang et al., 2021). Similar to the single model approach, the performance of the different SEL implementation can vary from one data set to another with no one of these outperforming the others in all data sets. Thus, as it was recommended by these authors it is important to evaluate several models as meta-learners as well.

Regarding the MSE, as expected, the penalized models, showed larger values compared to the other evaluated methods. By design these methods induce bias aiming to reduce the variance of the estimations (Montgomery et al., 2021; Chan et al., 2022). However, these cannot guarantee the increasing of the PA compared to other methods. The SMB, which resulted to return the best results in terms of PA, also presented large values for the MSE. This can be a consequence that SMB-SEL uses predicted values derived from base learners that return large MSE values ( Supplementary Table S1 ).

Overall, the SEL models presented moderate to high Spearman’s correlation between them ( Figures 4 – 6 ). On the other hand, these were low to moderate between SEL approaches and the GBLUP model. Additionally, among the 10% of genotypes with the highest GEBVs for YL, NF, LM, and Cer, the agreement coefficient between the SEL and GBLUP models showed values varying from moderate to high, suggesting differences in the obtained classifications with these. In general, the cluster analysis of these results showed that the methods can be grouped into three distinct groups ( Supplementary Figures S6 - S13 ) with the GBLUP forming a group by itself.

Altogether, these results show that the use of SEL to predict the individual genetic merit of four important traits in Arabica Coffee is worth to investigate. The SEL approach showed higher estimates of PA compared with all evaluated base learning methods, in special to the traditional GBLUP method. In practice, SEL’s ability to combine methods with diverse characteristics facilitates a more comprehensive exploration of the relationships between variables leading to more accurate selection of breeding. This approach considers a wider range of factors and reduces the reliance on any single model’s limitations. However, evaluating phenotypes across multiple environments can pose challenges for SEL. Unlike GBLUP, which presents higher computational efficiency (Hernandez et al., 2020), many base learners in SEL methods are based on machine learning requiring significant computation time in certain scenarios. Studies have explored the use of single machine learning methods for multi-environment trials (METs). For example, Barreto et al. (2024) applied machine learning to predict hybrid performance in METs and achieved similar PA compared to GBLUP with non-additive effects. As highlighted by Montesinos López et al. (2022) in their study using RF for METs, training any machine learning model can be computationally demanding, especially for datasets where the training data sets are very large. The hyperparameter tuning for individual base learners within a SEL framework is a well-established approach to enhance model performance. However, it is important to acknowledge that SEL ensembles can achieve strong results even with default base learner parameters (Friedel et al., 2023). This aligns perfectly with the core principle of ensemble learning, that is, leveraging predictions from multiple models can outperform any single model. In our study, the high dimensionality of the data presented significant computational challenges for hyperparameter tuning. Additionally, the observed superiority of the SEL approach compared to traditional methods suggested that tuning might not be as critical for achieving good results.

The SEL method was able to predict the PA of important traits (YL, NF, leaf miner infesting, cercosporiosis resistance) in Coffea Arabica. In addition, SEL presented higher PA compared with those obtained for all base learner methods (GBLUP, MARS considering degrees equal to 1, 2, and 3, QRF considering nine quantiles, from 0.1 to 0.9, every 0.1 and RF).

The original contributions presented in the study are included in the article/ Supplementary Material . Further inquiries can be directed to the corresponding authors.

MN: Writing – review & editing, Writing – original draft, Validation, Software, Methodology, Investigation, Formal analysis, Conceptualization. AN: Writing – review & editing, Software, Methodology, Investigation. CA: Writing – review & editing, Methodology. AO: Writing – review & editing, Data curation. EC: Writing – review & editing, Investigation, Data curation. DJ: Writing – review & editing, Methodology, Investigation.