Identification of Key Genes Related to Dormancy Control in Prunus Species by Meta-Analysis of RNAseq Data

Bud dormancy is a genotype-dependent mechanism observed in Prunus species in which bud growth is inhibited, and the accumulation of a specific amount of chilling (endodormancy) and heat (ecodormancy) is necessary to resume growth and reach flowering. We analyzed publicly available transcriptome data from fifteen cultivars of four Prunus species (almond, apricot, peach, and sweet cherry) sampled at endo- and ecodormancy points to identify conserved genes and pathways associated with dormancy control in the genus. A total of 13,018 genes were differentially expressed during dormancy transitions, of which 139 and 223 were of interest because their expression profiles correlated with endo- and ecodormancy, respectively, in at least one cultivar of each species. The endodormancy-related genes comprised transcripts mainly overexpressed during chilling accumulation and were associated with abiotic stresses, cell wall modifications, and hormone regulation. The ecodormancy-related genes, upregulated after chilling fulfillment, were primarily involved in the genetic control of carbohydrate regulation, hormone biosynthesis, and pollen development. Additionally, the integrated co-expression network of differentially expressed genes in the four species showed clusters of co-expressed genes correlated to dormancy stages and genes of breeding interest overlapping with quantitative trait loci for bloom time and chilling and heat requirements.


Introduction
Bud dormancy is a critical stage in the phenology and reproductive cycle of Prunus species, as the survival and production of these crops depend on synchronization between growth and dormant periods with seasonal temperatures [1]. Cultivars must be adapted to temperatures in a particular growing area to avoid inadequate flowering and fruit set and therefore ensure fruit production. Thus, in Prunus crops, flowering occurs after a dormancy period during the autumn and winter, in which meristem growth is inactive to avoid unfavorable environmental conditions [2,3]. This dormant period, common in all Prunus species, is divided into two main stages [4]: endodormancy, in which bud growth is inhibited and specific amounts of chilling temperatures are required to resume growth, and ecodormancy, which refers to a period of heat accumulation needed to bud burst after chill satisfaction [5]. The transition between these stages is irreversible and is determined by the fulfillment of the chilling requirement, which is a genotype-dependent trait. This requirement is calculated by modeling temperatures in a particular environment using different approaches to estimating the chilling accumulation needed to break dormancy. The simplest model, called the 'chill hour' model, counts the number of hours below 7.2 • C. In contrast, a more advanced 'dynamic model' reports chill accumulation as chill portions that account for temperature fluctuations and is a more accurate way of estimating chill accumulation, especially in warm regions [6]. Within Prunus species, chilling requirement and bloom times (Table 1). For all of them, sample points at the same dormancy stage (Endo1, Endo2, and Eco1) were selected, and a total of 136 transcriptomes corresponding to these three stages were used in analyses (Table S1). Pairwise comparisons of these transcriptomic data between endodormancy initiation (Endo1), endodormancy progression (Endo2), and ecodormancy (Eco1) were performed to identify genes associated with these dormancy stages. Consequently, 13,018 of 47,089 annotated genes in the peach genome v2.0.a1 were differentially expressed (FC ≥ 2; p-value ≥ 0.05) in at least one of the analyzed cultivars (Table S2). Of these, 3603, 4595, 6246, and 6315 genes were differentially expressed in almond, apricot, peach, and sweet cherry cultivars, respectively. The projection of these DEGs into a two-dimensional space (principal component analysis, PCA) revealed three sample clusters, each corresponding to a different species, except for almond and apricot samples that were grouped together. However, the samples were not grouped according to a dormancy stage ( Figure S1).

Endodormancy-Related Genes
To better understand whether specific genes could be related to a particular dormancy stage, we further investigated genes showing differential expression during endodormancy. The endodormancy-related genes comprised 6860 (53% of dormancy DEGs) transcripts that were differentially expressed between the beginning and progression of endodormancy (Endo1 vs. Endo2). Endodormancy DEGs in each cultivar ranged from 472 ('Cristobalina') to 2375 ('Royal Dawn'), with most of the cultivars having 650-1000 DEGs and displaying a higher proportion of upregulated genes, except the two apricots and the sweet cherry cultivar Cristobalina ( Figure 1A). Of all DEGs at endodormancy, only 55, 91, 107, and 381 genes were differentially expressed in all sweet cherry, peach, almond, and apricot cultivars, respectively ( Figure 1B). Even though no common DEGs were observed in all evaluated cultivars, 139 genes were of particular interest (Table S3), as they were differentially expressed in at least one cultivar of each species ( Figure 1C). GO enrichment analysis of these genes revealed terms involved in biological processes related to nitrogen compound metabolic process, external encapsulating structure organization, cell wall metabolic process, and molecular functions associated with transcription regulator activity, DNA-binding transcription factor activity, and xyloglucan:xyloglucosyl transferase activity ( Figure 1D). Furthermore, hierarchical clustering of these genes revealed two main groups based on their expression profiles (clusters 1 and 2; Figures 1E and S2; Table S3). Cluster 1 corresponds to genes (N = 107) showing maximum expression levels at the beginning of endodormancy (Endo1) and downregulation through the progression of endodormancy (Endo2) and ecodormancy (Eco1) in almond, peach, and sweet cherry. In the two apricot cultivars, a distinctive expression profile in these genes, with maximum expression during Endo2 and Eco1, was observed ( Figures 1E and S2). The other cluster (cluster 2; N = 32) grouped genes only differentially expressed during Endo1 and Endo2 in sweet cherry and peach cultivars and downregulated in the four species at ecodormancy (Figures 1E and S2).
Genes responsible for hormone signaling were also observed to have key functions during endodormancy. Several genes highly transcribed at the beginning of endodormancy in the four Prunus species were annotated with functions related to hormone pathways such as allene oxide synthase 1 (Prupe.1G386300), probable 9-cis-epoxycarotenoid dioxygenase (Prupe.4G082000), abscisic acid 8 -hydrolase 1 (Prupe.5G013100), isoflavone reductase (Prupe.5G022800), and NDR1/HIN1-like protein (Prupe.7G097200) (Table S3). A dormancyassociated MADS-box gene Prupe.1G531700 is also a member of this set of genes with a maximum expression profile at the beginning of endodormancy (Endo1) and downregulation during Endo2 and Eco1 in the four species.

Ecodormancy-Related Genes
The ecodormancy-related genes comprised DEGs between the transition from endodormancy (Endo2) to ecodormancy (Eco1). The 7506 DEGs (58% of dormancy DEGs) displayed a similar expression pattern, with sweet cherry 'Kordia' (31) and apricot 'Palsteyn' (2087) having the most contrasting number of DEGs in this transition ( Figure 2A). However, as foreseen in endodormancy-related genes, none of the 7506 genes related to ecodormancy were differentially expressed in the fifteen cultivars at the same time, with only 16, 52, 75, and 316 genes commonly expressed in sweet cherry, peach, almond, and apricot cultivars, respectively ( Figure 2B). Of these, 223 DEGs had key functions in the transition from endoto ecodormancy in Prunus species, as they were differentially expressed in at least one cultivar of each species ( Figure 2C). GO analysis of these genes revealed highly enriched categories associated with biological processes, including nitrogen compound metabolic, macromolecule metabolic, gametophyte and pollen development, and carbohydrate transport ( Figure 2D). In addition, the hierarchical clustering of these genes highlighted two main groups based on their expression profile ( Figures 2E and S2; Table S4). Cluster 1 corresponds to genes (N = 145) downregulated during endodormancy (Endo1 and Endo2) and showing maximum expression levels after chilling fulfillment (Eco1) in all cultivars except sweet cherry cultivars Burlat, Kordia, and Royal Dawn. Cluster 2 corresponds to genes (N = 78) downregulated during ecodormancy (Eco1) that were displaying high upregulation at the end of endodormancy (Endo2) (Figures 2E and S2; Table S4).

Individual and Integrated Prunus Co-Expression Networks
Co-expression networks were independently constructed in almond, apricot, peach, and sweet cherry using DEGs in each species to identify genes displaying similar expression profiles during dormancy stages ( Figure 3A). From WGCNA, nine modules of co-expressed genes in almond and apricot and eleven in peach and sweet cherry, containing 19 (Pdu_ME9) to 1487 genes were found (Par_ME1) ( Figure 3B). Moreover, the correlation between these modules and dormancy phases revealed co-expressed genes highly associated with particular dormancy stages ( Figure 3C). Modules Pdu_ME1, Pdu_ME4, Par_ME6, Ppe_ME6, Ppe_ME11, and Pav_ME8 were positively correlated to Endo1, in which GO analysis revealed significant enrichment in the regulation of the metabolic and cellular process, lipid metabolism, and cell wall regulation of co-expressed genes (Table S5). In addition, five modules (Par_ME3, Ppe_ME3, Ppe_ME8, Pav_ME7, and Pav_ME10) showed a positive correlation with the progression of endodormancy (Endo2) in apricot, peach, and sweet cherry. However, none of the almond modules presented a statistically significant correlation to this stage ( Figure 3C). Enrichment analyses of Endo2-related modules revealed genes associated with responses to abiotic stimulus and secondary metabolic processes (Table S5). Finally, seven modules of co-expressed genes displayed positive (Pdu_ME2, Pdu_ME5, Par_ME7, Pav_ME1) and negative (Ppe_ME1, Par_ME8, Pav_ME8) correlations to ecodormancy (Eco1), with genes in positively correlated modules enriched in pathways related to carbohydrate metabolism, hormone responses, and responses to a stimulus (Table S5). The independent co-expression networks for each species provided important information about co-expressed genes that might be relevant to certain dormancy stages, with various endo-and ecodormancy correlated modules sharing similar GO terms between species. Thus, to explore common expression patterns in the Prunus genus, we integrated the apricot, almond, peach, and sweet cherry datasets to perform a coexpression network. The 13,018 genes differentially expressed in at least one of the analyzed cultivars during the dormancy stage were considered for network construction and clustered into 25 modules of genes with similar expression profiles ( Figure 4A) that contained from 122 (ME11) to 1344 (ME1) genes ( Figure 4B). Of these, only ME11, which includes 122 genes, showed a significant correlation with Endo1, Endo2, and Eco1 ( Figure  4C). A similar expression profile in cultivars of the four species (upregulation at Endo1 and progressively downregulation through Endo2 and Eco1; Figure 5) was observed for genes in this module. Enrichment analysis of ME11 genes did not reveal any significant GO term. An additional module, ME12, was also related to both endodormant stages (Endo1 and Endo2; Figure 4C). It contains 140 genes mainly upregulated in sweet cherry cultivars during Endo1 and further downregulated during Endo2 and Eco1, with GO terms related to response to oxidative stress, response to abiotic stimulus, and carbohydrate catabolic process (Table S6). The ME22 was likewise correlated with Endo2 The independent co-expression networks for each species provided important information about co-expressed genes that might be relevant to certain dormancy stages, with various endo-and ecodormancy correlated modules sharing similar GO terms between species. Thus, to explore common expression patterns in the Prunus genus, we integrated the apricot, almond, peach, and sweet cherry datasets to perform a co-expression network. The 13,018 genes differentially expressed in at least one of the analyzed cultivars during the dormancy stage were considered for network construction and clustered into 25 modules of genes with similar expression profiles ( Figure 4A) that contained from 122 (ME11) to 1344 (ME1) genes ( Figure 4B). Of these, only ME11, which includes 122 genes, showed a significant correlation with Endo1, Endo2, and Eco1 ( Figure 4C). A similar expression profile in cultivars of the four species (upregulation at Endo1 and progressively downregulation through Endo2 and Eco1; Figure 5) was observed for genes in this module. Enrichment analysis of ME11 genes did not reveal any significant GO term. An additional module, ME12, was also related to both endodormant stages (Endo1 and Endo2; Figure 4C). It contains 140 genes mainly upregulated in sweet cherry cultivars during Endo1 and further downregulated during Endo2 and Eco1, with GO terms related to response to oxidative stress, response to abiotic stimulus, and carbohydrate catabolic process (Table S6). The ME22 was likewise correlated with Endo2 and Eco1. It displayed a common expression pro-file in all cultivars, downregulated during the last stage of endodormancy and upregulated during ecodormancy ( Figure 5). GO analysis in this module exhibited terms mainly related to cell cycle control (Table S6). In addition to modules that showed correlations with various dormancy stages simultaneously, we observed modules significantly associated with a particular dormancy stage ( Figure 4C). This is the case for ME10, ME13, and ME21, which were positively correlated to Endo1. However, the evaluation of the expression profiles of these modules revealed genes only differentially expressed in peach, sweet cherry, and almond cultivars in each module ( Figure 5). Finally, three modules (ME2, ME3, and ME14) were positively correlated with ecodormancy, consisting of genes over-expressed in this stage compared with endodormancy. ME2 and ME14 appeared to be specific to cultivars of some species (peach and apricot for ME2; apricot for ME14), and only ME3 was observed to be correlated to ecodormancy in all four Prunus species ( Figure 5). The GO analyses of these modules revealed terms related to the nitrogen and macromolecule metabolic process (ME2), transmembrane transport, small molecule and carbohydrate metabolic process (M14), and hormone biosynthesis processes (ME3) ( Table S6). and Eco1. It displayed a common expression profile in all cultivars, downregulated during the last stage of endodormancy and upregulated during ecodormancy ( Figure 5). GO analysis in this module exhibited terms mainly related to cell cycle control (Table S6). In addition to modules that showed correlations with various dormancy stages simultaneously, we observed modules significantly associated with a particular dormancy stage ( Figure 4C). This is the case for ME10, ME13, and ME21, which were positively correlated to Endo1. However, the evaluation of the expression profiles of these modules revealed genes only differentially expressed in peach, sweet cherry, and almond cultivars in each module ( Figure 5). Finally, three modules (ME2, ME3, and ME14) were positively correlated with ecodormancy, consisting of genes over-expressed in this stage compared with endodormancy. ME2 and ME14 appeared to be specific to cultivars of some species (peach and apricot for ME2; apricot for ME14), and only ME3 was observed to be correlated to ecodormancy in all four Prunus species ( Figure 5). The GO analyses of these modules revealed terms related to the nitrogen and macromolecule metabolic process (ME2), transmembrane transport, small molecule and carbohydrate metabolic process (M14), and hormone biosynthesis processes (ME3) ( Table S6).

Candidate Genes in Major Prunus Bloom Time and Dormancy-Related QTLs
Ten regions of the Prunus genome with significant bloom time, chilling, and heat requirement QTLs in almond, apricot, peach, and sweet cherry were chosen to search for candidate genes ( Figure 6; Table 2). These regions were distributed along the Prunus genome, except chromosome (chr.) 3. Of the 134 DEGs, 49 showed annotations and functional predictions associated with dormancy control (Figure 6; Tables 2 and S7). Three of these genome regions, on chrs. 1 (43-43.5 Mbp), 4 (1.5-7 Mbp), and 7 (11)(12)(13)(14)(15)(16)(17)(18)(19), appear to be involved in dormancy control within the Prunus genus, as evidenced by QTL overlap in various species, and percentages of phenotypic variance explained. Within the interval on chr. 1 , in which dormancy-related QTLs were reported for almond, apricot, peach, and sweet cherry studies, three DEGs-Prupe.1G531100, Prupe.1G531600, and Prupe.1G531700, annotated as MADS-box-were found in this region (Table 2). These genes show maximum expression at the beginning of endodormancy and downregulation during ecodormancy in all cultivars ( Figure 6).
Finally, in the region on chr. 7 (11)(12)(13)(14)(15)(16)(17)(18)(19), in which highly significant flowering QTLs were detected for all species but sweet cherry, three of the ten DEGs were mapped with annotations related to dormancy regulation (Prupe.7G142500; Prupe.7G161100; Prupe.7G168200) ( Table 2). The Prupe.7G161100, annotated as cold shock protein CS66, showed a similar expression pattern in all cultivars, overexpressed at the final stage of endodormancy (Endo2) and downregulated during ecodormancy ( Figure 6). For the other QTL regions associated with particular species, detailed information on DEGs and candidate genes related to dormancy control is provided in Tables S7 and 2, respectively. Table 2. Genes differentially expressed in Prunus cultivars with annotations related to dormancy control within main bloom time, chilling and, heat requirement QTLs previously found in almond, apricot, peach, and sweet cherry. Gene ID, gene position, and QTL region based on peach genome v2.a.01.

Region
Gene

Discussion
Recent literature on RNAseq-related studies reporting response to dormancy control in Prunus species has revealed a significant and unmanageable number of genes associated with genetic regulation of this period [28][29][30][31][32][33][34][35][36][37][38][39][40][41]. Most studies focused on understanding dormancy at the species level, while only two attempted to combine transcriptome data of various Prunus species to identify genes with similar expression patterns during dormancy transitions [30,53]. In this study, an extensive dataset comprising almond, apricot, peach, and sweet cherry cultivars was analyzed through a meta-analysis of RNAseq data to search for genes related to length and transition between endo-and ecodormancy phases. During this transition, these four Prunus species have conserved similar pathways, and their activation is triggered at comparable dormancy points. Our results suggest that although several pathways were activated during this period, dormancy is mainly regulated by a complex interaction of hormones, carbohydrate balance, and variations in cell wall structure, all in response to internal and external stimuli such as abiotic stresses.
In temperate fruit species, several hormones such as abscisic acid, auxins, gibberellins, ethylene, cytokinins, and jasmonates have been widely reported to control development and growth stages triggered by environmental factors [48]. Of these, abscisic acid is considered the primary hormone controlling dormancy activation and growth resumption, as its concentration was the highest at endodormancy induction and progressively decreased towards chilling accumulation in peach [58] and sweet cherry [59]. In this meta-analysis study, we found abscisic acid-related genes such as abscisic acid 8 -hydroxylase (Prupe.5G013100), phosphatase 2C protein (Prupe.5G118100 and Prupe.8G139700), and 9-cis epoxycarotenoid dioxygenases (NCED; Prupe.4G082000) playing a central role on dormancy regulation in all Prunus species. Among these, Prupe.4G082000 had high expression at the beginning of endodormancy and progressively decreased during chilling accumulation, as previously observed in peach, sweet cherry, pear, and grapevine [33,49,58,60,61]. This gene overlapped with the main bloom time and chilling requirement QTLs on chromosome four in all Prunus species [8,10,11,[18][19][20]22,56,62] and appears to be a key gene related to dormancy control in Prunus species. The Prupe.4G82000 is associated with abscisic acid accumulation [63][64][65], and together with CYP707A genes that are associated with abscisic acid catabolism, such as Prupe.5G013100 [61,66], regulate the balance of this hormone during dormancy. These two genes, NCED and CYP707A, had similar expression profiles in almond, peach, and sweet cherry cultivars of this study, which agrees with the hypothesis that abscisic acid catabolism is activated at the same time as metabolism to support homeostasis [67]. However, high expression of Prupe.5G013100 was observed in the two apricot cultivars included in this study, as reported for other CYP707A genes in some peach and Japanese pear cultivars during dormancy release [58,61,68], revealing slight differences within genotypes. In addition to abscisic acid, maintaining gibberellic acid equilibrium is another crucial process regulating dormancy [40,48]. This hormone is inversely correlated with abscisic acid, exhibiting low concentration during the early stages of endodormancy and progressively increasing until reaching the concentrations during dormancy release, as observed in sweet cherry and Japanese apricot [69][70][71]. The pattern of gibberellic acid concentration agrees with the expression profile of Prupe.4G080700.1, a GA20x gene associated with gibberellic acid catabolism. GA20x appeared essential during the transition from endodormancy to ecodormancy, as its expression was inhibited during chilling accumulation in the four Prunus species, as well as in rose [72] and Japanese apricot [73]. Thus, gene expression and hormone content suggest that high concentrations of abscisic acid inhibit dormancy release, while gibberellic acid promotes the transition from endodormancy to ecodormancy. Interaction of these three genes (NCED, CYP707A, and GA20x) is crucial to maintain an adequate hormone balance and, therefore, dormancy control since alteration in some of these genes, such as the overexpression of GA20x, resulted in delayed bud burst [74] and different hormone balance [75].
Flower bud protection from freezing temperatures during the winter is another critical phase during dormancy. In this study, several genes associated with cold response have been observed, mainly upregulated during the early stages of chill accumulation, ensuring flowering took place under optimal temperatures. In particular, we observed C-repeat/dehydration-responsive element-binding factors (CBF/DREB) playing a pivotal role in cold regulation, especially during the induction of endodormancy in peach, almond, and sweet cherry. In apricot cultivars, maximum expressions were reported at the end of endodormancy [30]. These genes control large parts of the cold-induced changes and have been implicated in endodormancy control in various perennial species [76]. For example, overexpression of PpeCBF1 in apple caused a critical alteration in cold acclimation and dormancy [77]. Moreover, CBF/DREBs have been reported to control dormancy by binding to DAM genes and regulating their expression in apple, Japanese apricot, and pear [43,[77][78][79]. Similarly, other transcription factors also appear to be involved in cold responses, such as MYB14 (Prupe.3G227800), which showed maximum expression at the induction of endodormancy and was progressively suppressed during chilling accumulation, as previously reported in sweet cherry [33], Arabidopsis [80] and grapevine [81]. This suggests these transcription factors suppress genes subsequently activated during ecodormancy. Another cold-related gene, Prupe.2G294400, was of particular interest because of its differential expression during the transition from endodormancy to ecodormancy and overlap within a major QTL at the bottom region of chromosome 2 in sweet cherry (bloom time; [22]) and almond (heat requirements; [10]). This gene is predicted to encode a low-temperature-induced 65kDa protein, homolog to CAP160 protein in spinach that increases freezing tolerance when overexpressed [82][83][84], making it a suitable candidate associated with the phenotype variation of this QTL. In addition, two genes, Prupe.2G139000 and Prupe.2G139500, also associated with cold hardiness in peach [85], were reported here in the four species, showing maximum expression at the end of endodormancy and downregulation during deacclimation. This expression pattern is similar to that previously observed in peach and supports the hypothesis that these genes increase cold hardiness in tree shoots [85]. Another DEG, Prupe.7G161100, appeared essential in the transition from endodormancy to ecodormancy in Prunus. Prupe.7G161100 is a cold shock protein (CSP) co-located within the main bloom time QTL of Prunus species on chromosome seven [8,10,11,18,19,22] and highly expressed before dormancy release and downregulated during ecodormancy in the four species. This expression pattern correlates with that observed in Arabidopsis, in which distinctive expression patterns of CSP genes were associated with increasing cold tolerance [86] and negative regulation of bloom time [87]. Additionally, AtAGL15, a key gene controlling flowering in Arabidopsis, has been observed to interact with the promoter regions of two CSPs, affecting their expression [88], which supports the role of Prupe.7G161100 in the dormancy control of Prunus species.
In addition to hormone and cold-related genes, GO analysis revealed a large number of DEGs associated with cell wall modifications and plasmodesmata in all Prunus cultivars. These genes appear to be involved in critical structural mechanisms that regulate long-distance transport during dormancy onset and release through alterations in cell membranes [89,90]. For example, in this study, we found that genes involved in cell wall organization, such as those encoding for non-specific lipid transfer proteins, were upregulated during ecodormancy and facilitated lipid transport between membranes [91]. This agrees with the observation of lipid content modification throughout dormancy and its potential effect associated with cytoskeleton modification [89]. Similarly, the accumulation of callose during endodormancy establishment has been observed to be essential to block plasmodesmata and regulate dormancy [92] through the expression of CALLOSE SYN-THASE 1. This gene is activated by abscisic acid and represses glycosyl-hydrolase genes that, after chilling fulfillment, hydrolyze callose to resume growth [93,94]. This model was confirmed for all Prunus species in this study, with the upregulation of abscisic acid genes during endodormancy and the later activation of glucan endo-1,3-beta-glucosidases genes in ecodormancy that progressively degrade callose, allowing metabolite transport such as carbohydrates through flower buds [40]. These carbohydrates were essential for transition from endodormancy to ecodormancy, as an energy source, and as sugar-related signals [95,96]. In Prunus crops, the accumulation of soluble sugars increases during the induction of cold conditions and then declines after the chilling requirement is satisfied [59,71,97], which was also suggested as a mechanism to increase protection against low temperatures [40,98]. After chilling fulfillment, the activation of genes related to sugar transport in ecodormancy for cultivars of the four species (Prupe.1G220700, Prupe.3G283400, Prupe.5G083900, Prupe.5G090900, Prupe.5G146500 and Prupe.8G017500) was observed in this study, suggesting reserve mobilization as an essential step in resuming bud growth.
Even though both up-and downregulated genes during endo-and ecodormancy were reported in the above discussed pathways, in this meta-analysis, we observed a large number of genes associated with flower development, whose activation was initiated during ecodormancy in the four species, and also reported in peach and apricot [30,99]. For example, it includes genes such as Prupe.1G192200 that encodes for a strictosidine synthase, homolog to OsSTRL2 in rice, and is observed to play a key role in pollen formation and anther development [100]; the MYB80 transcription factor Prupe.6G032800, which is essential for pollen development [101]; or the ABORTED MICROSPORES gene (Prupe.8G166900), whose function was associated with pollen wall formation and male sterility [102,103]. This same expression pattern, with pollen development-related genes activated after chilling satisfaction, was observed in all studied species. It is supported by microscopy observations in sweet cherry [104] and apricot [105], where stamen and pollen development were reported after endodormancy release (ecodormancy).

Selection of RNAseq Studies and Sample Points Normalization
Transcriptomic studies reporting response to dormancy release in Prunus species were identified in the literature by searching the keywords 'Prunus', 'dormancy', 'flower buds', 'RNAseq', and 'transcriptomic' in the NCBI Sequence Read Archive (SRA) database. RNAseq datasets of these studies were selected and retrieved from NCBI SRA if they fit the following criteria: (i) RNAseq dataset was publicly available; (ii) a Prunus species was evaluated; (iii) transcriptome data from flower buds were sampled during dormancy stage; (iv) dormancy status of each sample point was known or could be inferred; and (v) Illumina technology was used for sequencing. Following these requirements, six studies from four Prunus species (almond, apricot, peach, and sweet cherry) were selected (Table 1), and transcriptomic datasets from fifteen unique cultivars were downloaded (Table S1). For the comparative study, sample points at three common dormancy stages were identified in all cultivars according to their phenological status (Endo1, Endo2, and Eco1). Endo1 corresponds to an initial endodormancy point in which 0-20% of the chilling requirement of each cultivar was fulfilled. Similarly, Endo2 represents an endodormancy point at which more than 80% of the chilling requirement was fulfilled but heat accumulation had not been started. Finally, Eco1 refers to an ecodormancy point at which every cultivar had fulfilled its chilling requirement, but flowering had not been reached. The chilling requirements of cultivars considered in this work were calculated in each independent study by using chill hour ('Cristobalina', 'Garnet', 'Regina', 'Burlat', 'Royal Dawn', 'Kordia', 'A209 , 'A340 , 'A318 , and 'A323 ) and chill portion ('Desmayo largueta', 'Penta', and 'Tardona') models, except for the two apricot cultivars (Palsteyn and Bergeron) whose chilling requirements were predicted based on correlations between bloom dates and chill accumulation.

Read Alignment and Differentially Expressed Gene Analysis
Raw transcriptomic sequences of selected studies were retrieved from NCBI SRA (https:www.ncbi.nlm.nih.gov/sra, [106]) according to their accession numbers, and sequence quality was assessed through FastQC to spot potential problems in raw sequencing datasets. Low-quality reads and adaptors from Illumina sequencing were removed using Trimmomatic software v.0.38 [107], and clean reads were mapped to the peach genome v2.0.a1 [57] using the Bowtie2 short-read aligner [108]. Peach genome v2.0.a1 was used as a reference for all species because of its assembly and annotation quality to allow a more straightforward comparison between studies. Prunus genus shares significant collinearity and synteny; thus, mapping sequences from other species to the peach genome is possible [109]. For each gene in a sample, transcript abundance was quantified using RSEM [110], which counts the total number of raw reads mapped per gene. Additionally, we carried out a differential expression analysis on these cultivars using edgeR [111]. Significant differentially expressed genes (DEGs) were identified with a fold change (FC) ≥ 2 and a probability of differential expression higher than 0.95. Graphical visualization of gene expression was conducted by heatmap.2 function in the gplots package implemented in R [112].
To allow expression pattern representation, we performed data normalization using a z-score for each gene using the following formula:

GO Enrichment and Expression Profile Analysis
Gene ontology enrichment (GO) terms were retrieved from the http://www.geneontology. org database [113,114] that provides the functions of the gene products. Enrichment analyses were conducted on GO terms for biological processes, cellular components, and molecular functions using Fisher's exact test. Only GO terms with a false discovery rate (FDR) lower than 0.05 were considered.

Gene Co-Expression Network Analysis
Weighted gene co-expression network analysis was performed using WGCNA v.1.70-3 package implemented in R [115]. Co-expression networks were independently constructed for cultivars of each species (almond, apricot, peach, and sweet cherry) and combined for all Prunus cultivars. Only DEGs during dormancy were included for network construction, considering a minimum module size of 10 and a merge cut height value of 0.25 as parameters. The first principal component of the expression matrix of the corresponding modules was used to define module eigengene (ME), and Pearson's correlation coefficient at the significant level of p < 0.05 was used to determine correlations between modules and dormancy stages (Endo1, Endo2, and Eco1), where modules were denoted by their ME and dormancy stage by a numeric vector ('1 for the considered stage and '0 for the others).

Candidate Gene Analysis
To identify candidate genes related to bloom time, chilling, and heat requirements in Prunus species, we searched for DEGs that co-localized within QTLs previously reported for these traits (Genome Database for Rosaceae; [116]). Thus, genome regions with overlapping major QTLs for almond, apricot, peach, and sweet cherry were investigated, and their predicted physical position in the peach genome v2.0.a1 were obtained. Genes commonly differentially expressed during dormancy in at least one cultivar of each species and mapped within these QTL regions were considered for further analysis. For these genes, literature was searched for evidence (annotation and predicted function) of involvement in dormancy control (cold response, cell wall modification, oxidative process, carbohydrate regulation, and phytohormone biosynthesis).

Conclusions
With this work, we aimed to provide a general overview of the main transcriptome changes in flower buds during dormancy in four main Prunus crops. We reported common pathways and genes activated at similar dormancy points in the four species through the meta-analysis of publicly available RNAseq data from fifteen almond, apricot, peach, and sweet cherry cultivars. The results revealed the presence of a similar genetic control of this dormant period in the genus, with genes associated with cold response, hormone biosynthesis, cell wall organization, and sugar content having pivotal roles in the control of dormancy. Additionally, the overlap of DEGs within previously reported main QTLs for chilling requirement and bloom time provides further support for candidate genes associated with the phenotypic variation of these QTLs. Prupe.4G082000, an NCED gene related to abscisic acid biosynthesis, and the dormancy-associated MADS-box genes located in the interval of the main bloom time QTLs on chr. 4 and 1, respectively, appear to be targeted genes for increased resiliency to climate change in Prunus breeding because of their correlations to endo-and ecodormancy transitions in the genus.
Supplementary Materials: The following supporting information can be downloaded at: https: //www.mdpi.com/article/10.3390/plants11192469/s1, Figure S1: Principal component analysis (PCA) of transcriptomic data using the 13,018 differentially expressed genes in the four Prunus species. Figure S2: Expression patterns in the four Prunus species of gene cluster related to endodormancy (left) and ecodormancy (right); Table S1: List of samples available from each study and metadata information. Table S2: List of differentially expressed genes. The trimmed mean of M-values (TMM) for each sample and gene is shown. Table S3: Differentially expressed genes in at least one cultivar of each species (almond, apricot, peach, and sweet cherry) during endodormancy (Endo1 vs. Endo2). Gene cluster based on expression profile and annotation is also provided. Table S4: Differentially expressed genes in at least one cultivar of each species (almond, apricot, peach, and sweet cherry) during ecodormancy (Endo2 vs. Eco1). Gene cluster based on expression profile and annotation is also provided. Table S5: Significant gene ontology terms (biological process, molecular function, and cellular component) of gene modules correlated with endodormancy (Endo1 and Endo2) and ecodormancy (Eco1) in the individual co-expression network for each species. Table S6: Significant gene ontology terms (biological process, molecular function, and cellular component) of gene modules correlated with endodormancy (Endo1 and Endo2) and ecodormancy (Eco1) in the integrated coexpression network for individuals of all the species. Table S7: Genes differentially expressed in Prunus cultivars within main bloom time, chilling, and heat requirements QTLs previously found in almond, apricot, peach, and sweet cherry. Genes highlighted in red showed annotation related to dormancy control.