Dormancy in fruit species such as cherry (Prunus avium L.) is an essential process that synchronizes flower bud growth and reproduction with seasonal cycles (Considine and Considine, 2016). The bud dormancy has been classified into three widely accepted phases: paradormancy, endodormancy, and ecodormancy, which respond to physiological and environmental signals, allowing buds to survive adverse winter conditions and resume development in spring (Cline and Deppong, 1999; Considine and Considine, 2016). The dormancy has generally been assessed using destructive phenological assays (e.g., forcing tests; Baumgarten et al., 2021) and agroclimatic models (chill hours, chill units, or chill portions; Luedeling et al., 2011; Baumgarten et al., 2021). However, these approaches have significant limitations in terms of accuracy and interannual application due to dynamically changing environmental conditions (Alonso-González et al., 2020). Therefore, several studies have sought molecular biomarkers that allow more accurate, rapid, and noninvasive monitoring of bud physiological status.

During the past decade, progress has been made in the study of the mechanisms associated with dormancy stages in various species, improving the general understanding of dormancy transitions in different species (Rothkegel et al., 2017; Yang et al., 2021; Rothkegel et al., 2020; Ding et al., 2024). These studies have revealed that dormancy phases are regulated by DORMANCY-ASSOCIATED MADS-box (DAM) genes, phytohormones, carbohydrates, temperature, photoperiod, reactive oxygen species, water deprivation, cold acclimation, and epigenetic regulation. Particularly, DNA methylation has emerged as one of the key epigenetic mechanisms involved in plant development and environmental responses. In this sense, dynamic changes in genome methylation patterns throughout the dormancy cycle, suggesting a functional role for this epigenetic mark in determining bud fate (Rothkegel et al., 2020; Yang et al., 2021; Kumar et al., 2016). In this sense, Rothkegel et al. (2017) indicated that in Prunus avium L, genome methylation patterns are one of the mechanisms that regulate the MADS-box genes controlling bud dormancy. Similarly, Kumar et al. (2016) showed that DNA methylation has been linked to chilling acquisition during dormancy in Malus domestica. Additionally, high levels of DNA methylation have been observed during the induction of endodormant floral buds in blueberry compared to those in the ecodormant stage (Li et al., 2015), which indicates a dynamic regulation throughout the dormancy stages.

The development of high-throughput sequencing technologies, such as whole-genome bisulfite sequencing (WGBS), enables the precise examination of methylation status at individual CpG sites with high resolution (Zhou et al., 2019). However, this approach generates highly complex and high-dimensional datasets, where the number of features (methylated sites) far exceeds the number of available samples. This problem, coupled with biological variability between cultivars and climatic heterogeneity between seasons, poses significant challenges for conventional statistical analyses and hampers biomarker discovery and straightforward biological interpretation (Zhang et al., 2018).

Artificial intelligence (AI)-based approaches have shown great promise for the analysis of complex and high-dimensional molecular data (Bzdok et al., 2018). However, one of the main limitations of these algorithms in biology is their interpretability and identification of biologically relevant variables. In this sense, feature selection is essential to eliminate noisy features, improve model performance, and optimize interpretability (Lundberg and Lee, 2017). Ensemble learning algorithms, such as Random Forest (RF) and eXtreme Gradient Boosting (XGBoost), have demonstrated exceptional performance in genomic studies due to their ability to handle noisy, imbalanced, and high-dimensional data while providing internal feature importance metrics (Raihan et al., 2023; Peng and Yu, 2024; Chu et al., 2024; Li et al., 2024; Zhang et al., 2024). Furthermore, the SHapley Additive exPlanations (SHAP) method offers a game-theoretic interpretation of each variable’s individual contribution to model predictions (Lundberg and Lee, 2017). These models (RF, XGB, and SHAP) have been successfully used to identify important variables in predictive models for human diseases (Raihan et al., 2023), agronomic traits (Zhang et al., 2024), precision agriculture (Li et al., 2024), and seed viability (Chu et al., 2024).

In this study, DNA methylation profiles obtained via WGBS with supervised Machine Learning (ML) algorithms (RF and XGB) and SHAP interpretability analyses were integrated to identify informative cytosines and methylated regions capable of discriminating between different dormancy stages in sweet cherry floral buds. Classification models were constructed for two scenarios: (i) three classes (paradormancy, endodormancy, and ecodormancy stages) and (ii) two classes (endodormancy and ecodormancy stages). This integrative strategy not only improves classification accuracy of dormancy stages but also provides deeper insights into the epigenetic mechanisms governing this process, offering potential biomarkers for breeding programs and more precise tools for phenological management in perennial fruit species.

The accurate and timely assessment of dormancy stages is a critical challenge in perennial horticulture, directly impacting orchard management and yield optimization. In this study, high-resolution DNA methylation profiling and machine learning techniques were integrated to identify robust epigenetic markers for predicting dormancy in sweet cherry. In plants, methylation occurs in three different contexts that serve a variety of roles in mediating environmental signals with molecular regulation. mCG methylation is found throughout the genome, including gene bodies, promoters, and repetitive regions, while mCHG and mCHH methylation are particularly enriched in repetitive sequences, including transposable elements, and play a key role in silencing these elements to maintain genome stability.

The results showed that in the three-class scenario, when methylated cytosines data were used, the paradormancy stage was misclassified 100% of the time by the RF model and 78% of the time by XGBoost. In contrast, endodormancy and ecodormancy had misclassification rates of 0% and 31% with RF, and 25% and 50% with XGBoost, respectively. Similar results were observed in methylated regions, where paradormancy exhibited lower classification efficiency compared to endodormancy and ecodormancy. This pattern may be explained by the imbalance in the number of samples per class, which appears to correlate with the error rate, with paradormancy having the fewest samples (9) and endodormancy having the most (36). Due to this class imbalance, both models tended to misclassify paradormancy samples, often assigning them to the majority class, endodormancy. The poor predictive performance in classifying unbalanced datasets can be since machine learning algorithms are designed to maximize overall accuracy and may struggle to accurately classify the minority class (Imani et al., 2025). The results of this study showed that RF outperformed XGBoost in accuracy and efficiency across all classes, regardless of class imbalance. Similar findings were reported by Dube and Verster (2023), who observed that RF demonstrated robustness to class imbalance and outperformed several machine learning models, including naïve Bayes, XGBoost, and k-nearest neighbors, which struggled with imbalanced data.

Feature selection applied to methylated cytosines and regions improved classification accuracy by approximately 40% for RF and 34% for XGBoost, highlighting the effectiveness of this approach. Similar findings were reported by Hassan et al. (2023), where feature selection improved the accuracy of all classification models. This improvement is due to the effective selection of features, which ensures that the most relevant information is extracted from the dataset, enabling more accurate classification (Javidan et al., 2024). Bolón-Canedo and Alonso-Betanzos (2019) pointed out that different feature selection methods can yield varying subsets from the same dataset; thus, the results may lack stability. Therefore, combining multiple feature selection techniques, rather than relying on a single method, helps control variance, mitigates the limitations of individual approaches, and offers a more comprehensive view of feature relevance, which makes this approach successful (Bolón-Canedo and Alonso-Betanzos, 2019). Therefore, the approach presented in this study combines three algorithms for identifying relevant variables, since the combination provides a better approximation based on the idea that “two heads are better than one” (Bolón-Canedo and Alonso-Betanzos, 2019). Consistent with the classification results, RF outperformed XGBoost in accuracy following feature selection. Additionally, using methylated cytosine data proved more effective than using methylated regions, improving classification accuracy by 11% in RF and 19% in XGBoost. Similarly, Hassan et al. (2023) also showed that RF consistently outperformed XGBoost, both on the full dataset and after feature selection. Notably, the improvement in the classification after feature selection was also reflected in the t-SNE plots, since the clustering with the original data set, samples tended to cluster by cultivar rather than dormancy stage, whereas after feature selection, the clustering aligned more clearly with dormancy stages.

Despite the overall increase in accuracy after feature selection, the classification error for paradormancy samples remained high, ranging from 33% (with methylated cytosines using RF) to 78% (with methylated regions using XGBoost). This highlights a significant challenge in accurately classifying the paradormancy stage compared to endodormancy and ecodormancy. These results suggest that neither RF nor XGBoost was able to overcome the issue of class imbalance, even after applying feature selection. It is worth noting that the SMOTE algorithm, which is designed to generate synthetic examples of the minority class and alleviate class imbalance (Fernández et al., 2018), was also tested (data not shown). However, its application did not lead to improved classification performance. As a result, in the second scenario, paradormancy was excluded from the dataset, and only endodormancy and ecodormancy were used to train and build the classification models.

In the second scenario, where only endodormancy and ecodormancy were considered, the RF model showed an increase in accuracy of 15% with methylated cytosines and 9% with methylated regions, compared to the three-class scenario. In contrast, XGBoost experienced a decline in accuracy of 24% and 18% for cytosines and regions, respectively. Consistent with the second scenario observations, both models showed improved accuracy when feature selection was applied, reinforcing the importance of identifying the most informative variables. However, it is important to note that the observed increase in accuracy does not necessarily reflect an improvement in the overall classification performance, since in the three-class scenario, much of the accuracy lost was tied to the misclassification of paradormancy samples. In both scenarios, RF consistently achieved the highest accuracy across all configurations, confirming its superior performance and greater stability under both balanced and imbalanced conditions. These findings further emphasize the limitations of XGBoost when dealing with skewed datasets, as well as the benefit of applying feature selection to enhance the model performance of dormancy stages classification based on epigenetic information.

Given that RF consistently outperformed XGBoost in accuracy and robustness across all scenarios, interpretability efforts based on SHAP focused on RF models trained with the informative feature sets of both cytosine-level and region-level methylation data. In the cytosine dataset, features such as chr_4_31092165 and chr_1_38152240 emerged as the most impactful, particularly for distinguishing the ecodormancy stage. This suggests that methylation at these loci could be tightly associated with epigenetic regulation mechanisms unique to ecodormancy. Other features, like chr_3_23416192 and chr_4_31050901, displayed stronger relevance for paradormancy, indicating that although paradormancy was difficult to classify due to data imbalance, certain epigenetic markers do exist that could potentially improve its prediction if more balanced data were available. Similarly, in the regional methylation dataset, regions like chr_3_9331371_9331488 and chr_1_27137350_27137417 had high SHAP values. Interestingly, some regions, such as chr_6_9727681_9727971, demonstrated relatively balanced contributions across all three dormancy stages. This may indicate shared or transitional roles in the broader regulation of dormancy, potentially pointing to core epigenetic signatures common to different dormancy stages. The individual SHAP plots for each dormancy stage ( Supplementary Figure S8 ) showed that features such as chr_1_38152240 contributed strongly and positively to the correct classification of ecodormancy samples, while exhibiting negative contributions in misclassified samples labeled as endodormancy. This highlights the high discriminative power of this feature in separating these two stages. In contrast, features like chr_3_8464601, although showing high SHAP values, presented similar contributions in both ecodormancy and endodormancy samples. This suggests that it may not be informative on its own but could contribute to the model’s predictions through interactions with other features, an effect often captured by ensemble models like SHAP (Štrumbelj and Kononenko, 2014).

There are several examples in fruit trees where methylation of individual cytosines plays key roles during dormancy. Kumar et al. (2016) examined four developmental stages in apple, including dormant buds, and showed, using bisulfite–sequenced MSAP fragments, that the number of methylated cytosines at single loci correlated with transcript levels of those genes (e.g., an acid phosphatase 1-like gene and a galactose oxidase-like gene). Prudencio et al. (2018) in almond found that over 90 % of cytosines are unmethylated, while approximately 1-1.3 % are fully methylated in CpG contexts; these polymorphic 5-mC sites were conserved across two years and between dormant and non–dormant buds. These studies demonstrate that methylation at even a handful of cytosines in specific genomic contexts can meaningfully modulate gene expression. Given their stability and stage specificity, single–cytosine marks hold considerable promise as robust biomarkers for dormancy transitions. Although the precise mechanism by which a single methylated cytosine influences gene expression remains unclear, we hypothesize that such effects occur when the cytosine resides within functionally relevant regions, particularly transcription factor binding sites or cis–regulatory modules.

Differences in the number of relevant features were evident between the dormancy models. The 2-stages model generated the largest and most unique feature sets, resulting in 500 regions and 283 cytosine unique features. This genome-wide representation is consistent with the higher predictive power of 2-stage cytosine methylation data in Random Forest and XGBoost, suggesting that the additional loci due to the higher input samples, supplies informative data that improves the classification. By contrast, the 3-stage model contributed a leaner, yet distinct, repertoire of 182 unique regions and 58 unique cytosines, which may represent the early dormancy stage paradormancy signatures that the 2-stage model misses.

The reduced number of individual cytosine overlapping regions in the 2-stages model (4) may imply that single cytosines show more sensitivity to coverage noise or local sequence context, whereas regions provide a more robust, though less precise, relevant dormancy-associated loci. In order to address those questions, we searched for the genomic context of those relevant features from the 2-stage model to gather insights about the functional roles of those cytosines ( Figure 5B ). The genomic context frequencies rank in the same order across both feature types. From highest to lowest: TEs, promoter, gene bodies, downstream, and intergenic, indicating that cytosine and region features share similar genomic preferences. Since TEs were the highest genomic context with features, and they can be located in any of the above genomic contexts, we further study the genomic context of TEs with features and their corresponding class ( Figure 5 ). The number of TEs that remained unclassified is still high for P. avium Tieton annotation, highlighting the need for cultivar-specific transposable elements annotation.

Certain TE families exhibit preference for specific genomic features. A well known example is ONSEN, an LTR/Copia element in Arabidopsis thaliana that integrates preferentially within genes (Merkulov et al., 2025). Members of the LTR/Copia class often insert near gene-rich regions, unlike other TEs that target intergenic areas. Several Copia families including Copia87/ONSEN (Ito et al., 2011), COPIA37, TERESTRA, and ROMANIAT5 (Pietzenuk et al., 2016) have been well described as stress-responsive TEs, becoming active when temperature rises due to heat-responsive elements in their promoter regions. The LTR/Copia and Ty3, flanked by a promoter region, can undergo a replication cycle (Oberlin et al., 2017), with the later ability to confer neighboring genes the responsiveness to abiotic stressors (Xu et al., 2024). LTR/ty3-retrotransposons elements in P. avium were found predominantly in heterochromatic regions and LTR/Copia elements located in euchromatic domains, consistent with observations in almond (Alioto et al., 2020). TE activity is often suppressed through epigenetic mechanisms such as the RNA-directed DNA methylation (RdDM) pathway, which is responsive to environmental conditions. Nozawa et al. (2022) reported that ONSEN expression varies across Arabidopsis ecotypes depending on DNA methylation levels, with the Kyoto ecotype showing increased expression due to reduced methylation conferring release from transcriptional silencing. Although the functional dynamics of TE regulation during dormancy remain unclear in sweet cherry, we hypothesize that small RNAs involved in RdDM, previously characterized in P. avium (Rothkegel et al., 2017; 2020; Kuhn et al., 2025), may regulate TE transcription and potentially modulate the expression of adjacent genes involved in dormancy transitions.

Mapping these loci along the eight chromosomes of P. avium ( Figure 6A ) revealed broad genome coverage, with particularly relevant regions on chromosome 1 (0–2 Mb), chromosome 4 (25–32 Mb), and chromosome 5 (10–15 Mb). Chromosome 4 (LG4, Figure 6B ) includes a well reported QTL hotspot where flowering date (FD), chilling requirement (CR), and heat requirement (HR) QTLs, identified in sweet cherry progenies with Regina as the progenitor (Regina × Lapins and Regina × Garnet), co-localize with maturity date (MD) QTLs (Castède et al., 2014; Calle and Wünsch, 2020; Calle et al, 2020). Our model identified a concentration of informative regions, with fewer cytosine features, suggesting a broader chromatin-level regulation rather than single-site CpG variation may underlie dormancy transitions. This is consistent with previous reports demonstrating that DNA methylation and chromatin remodeling regulate dormancy release in Prunus (Kuhn et al., 2025; Rothkegel et al., 2020; Zhu et al., 2020). Within chromosome 4, several biologically relevant candidate genes have been reported (Castède et al., 2014; Calle and Wünsch, 2020). These include EMF2, a MADS-box transcription factor involved in floral meristem identity (Yoshida et al., 2001); NUA (Nuclear Pore Anchor), which regulates nuclear-cytoplasmic transport and chromatin organization; NAC domain transcription factors associated with stress responses and developmental control; EXPA1 (Expansin A1), involved in cell wall loosening and bud outgrowth; bHLH (basic helix-loop-helix) transcription factors implicated in hormone and light signaling pathways; and WRKY transcription factors known for their role in stress-responsive gene regulation.

Additionally, a QTL on LG1 includes AGL24-like MADS-box genes, which are involved in floral transition and respond to environmental cues, like the DORMANCY-ASSOCIATED MADS-box (DAM) gene cluster (DAM1 to DAM6) (Castède et al., 2015), which also overlaps with regions selected by our model. These genes have been shown to play a central role in the repression of bud growth during dormancy and its release in response to chilling accumulation. Castède et al. (2014) identified this region as a key determinant of phenological variation in sweet cherry, particularly in high-chill cultivars such as Regina. Additionally, Calle and Wünsch (2020) demonstrated that this QTL on LG1 remains significant even in low-chill cultivars, supporting its broad importance across diverse genetic backgrounds. In our study, several methylation features selected by the model colocalized with this QTL region, suggesting that epigenetic regulation contributes to differential DAM genes expression among cultivars. Increasing evidence supports the role of DNA methylation in modulating DAM genes activity during dormancy transitions. Zhu et al. (2020) showed that changes in methylation levels at the promoters and gene bodies of DAM genes in peach were associated with their transcriptional downregulation during chilling accumulation. Specifically, DNA hypomethylation correlated with a reduction in DAM transcript levels, suggesting that chilling-induced epigenetic remodeling facilitates dormancy release. Similar patterns of dynamic methylation have been reported in other Prunus species, underscoring the conserved nature of this regulatory mechanism. Our results align with these findings and point to DNA methylation as a key regulatory layer influencing the dormancy-related function of DAM genes. These observations support a model in which genetic loci like DAM are subject to cultivar-specific epigenetic modulation, contributing to phenotypic diversity in dormancy behavior across sweet cherry germplasm. Altogether, these results show the relevance of the LG1 QTL and the DAM cluster as conserved regulators of dormancy, while also highlighting the potential role of epigenetic variability in shaping cultivar-specific responses to chilling.

This study presents a new integrating frame that combines the high-resolution DNA methylation profile with automatic learning approaches to classify dormancy stages in Prunus avium. The results showed that the methylation data at the cytosine level, when processed through the selection of features and interpreted through form values, provide highly informative epigenetic markers to distinguish endodormancy and ecodormancy stages with high precision. Among the proven models, RF constantly exceeded XGBoost in robustness, precision, and interpretability in balanced and unbalanced scenarios. It is important to note that some of the methylation features co-localize with related QTLs previously identified by other groups, such as chilling and heat requirement flowering date and maturity date. The results of this study contribute a valuable methodological and conceptual advance for epigenetic research in these traits. It provides a fundamental resource for future functional validation studies and prepares the scenario to develop predictive tools in horticultural management. When discovering key methylation markers and their genomic contexts, this study improves our understanding of dormancy regulation and underlines the importance of epigenetic mechanisms in the adaptation and resistance of perennial crops.