RNA-seq analysis reveals different transcriptomic responses to GA3 in early- and mid-season varieties before ripening initiation in sweet cherry (Prunus avium L.) fruits

Gibberellin (GA) negatively affects color evolution and other ripening-related processes in non-climacteric fruits. The bioactive GA, gibberellic acid (GA 3 ), is commonly applied at the light green-to-straw yellow transition to increase rmness and delay ripening in sweet cherry (Prunus avium L.), though causing different effects depending on the variety. Recently, we reported that GA 3 delayed the IAD parameter (a ripening index) in a mid-season variety, whereas GA 3 did not delay IAD but reduced it at ripeness in an early-season variety. To further explore this contrasting behavior between varieties, we analyzed the transcriptomic responses to GA 3 applied on two sweet cherry varieties with different maturity time phenotypes. At harvest, GA 3 produced fruits with less color in both varieties. Similar to our previous report, GA 3 delayed fruit color initiation and IAD only in the mid-season variety, and reduced IAD at harvest only in the early-season variety. RNA-seq analysis of control- and GA 3 -treated fruits of both varieties revealed that ripening related categories were overrepresented in the early-season variety, including 'photosynthesis' and 'auxin response.' GA 3 also changed the expression of carotenoid and abscisic acid (ABA) biosynthetic genes in this variety. In contrast, overrepresented categories in the mid-season variety were related mainly to metabolic processes. In this variety, some PP2Cs putative genes were positively regulated by GA 3 , which are negative regulators of ABA responses. This was accompanied with downregulation of MYB44-like genes, putative repressors of PP2Cs expression. These results show that GA 3 differentially modulates the transcriptome at the onset of ripening in a variety-dependent manner and suggest that GA 3 impairs ripening through the modication of ripening associated gene expression only in the early-season variety; whereas in the mid-season variety, possibly a control of the ripening timing occurs through the PP2C gene expression regulation. This work contributes to the understanding of the role of GA in non-climacteric fruit ripening.


Introduction
Fruit ripening is a complex process that involves changes in the cell wall composition, accumulation of sugars and pigments in the peel (exocarp), cell enlargement, and the decrease in organic acid content. In non-climacteric fruits, which do not depend on an ethylene rise for ripening initiation, many hormones control this process, with abscisic acid (ABA) being the most important 1,2 . Exogenous application of ABA triggers the ripening process in non-climacteric species such as grapevine 3 (Vitis vinifera), strawberry 2 (Fragaria spp.) and sweet cherry 4 (Prunus avium). Other hormones are involved in this process, but their participation is poorly characterized. ABA's most characteristic feature is its sharp increase a few days before the anthocyanin accumulation 3,5 . However, ABA is not the only hormone that is present during this period. In sweet cherry, auxin, gibberellin (GAs), and cytokinins are present 5 . These hormones may modulate the ripening process together with ABA. For example, underripe grapevine berries have higher seed content and higher auxin levels in veraison bunches 6 , and NAA (α-naphthalene acetic acid) application delays ripening in this non-climacteric species 7 . In sweet cherry, indole-3-acetic acid (IAA) downregulates the ABA signaling pathway genes 8 .
Regarding GA, mainly GA 4 is detected after the grapevine berry set 9 . In strawberry, GA 4 decreases from the white to the red stage of the fruits 10 . In sweet cherry, GA 4 signi cantly and inversely correlates with ripening parameters 5 . Taking these results together, GA possibly exerts a negative role in the ripening process. At the molecular level, at veraison, some GAs oxidase genes change their expression compared to the previous stage in grapevine fruits 11 . GA 3 treatment at veraison delays ripening and affects the transcript levels of putative PP2Cs, which are negative regulators of ABA signaling 12 . In sweet cherry, GA 3 application delays ABA accumulation and reduces anthocyanin levels 13 . Additionally, GA 3 delays fruit size increase, ripening, and color development when applied at the degreening stage 14,15,16 . Choi et al. 14 reported that the mid-season varieties had a delayed fruit size increase and polygalacturonase and Cx-cellulase activities after the GA 3 treatment, whereas GA 3 did not affect fruit growth in the early-season varieties. Kuhn et al. 17 reported that the Index of Absorbance Difference (IAD), a ripening index, was delayed after the GA 3 treatment only in a mid-season variety, whereas the treatment reduced the IAD only at harvest in an early-season variety; interestingly, both varieties presented modi ed expression of ABA pathway related genes.
Despite the available information on the participation of GA in the fruit ripening process, the effect of this hormone has been barely explored at the whole gene expression level during the sweet cherry fruit ripening.
There are a few transcriptomic analyses in sweet cherry. These works explore the transcriptomic features underlying fruit development 18 , fruit coloring 19 , and light-dependent anthocyanin accumulation in sweet cherry fruits 20 . Guo and co-workers 20 found that some ABA and GA pathway genes co-expressed with light-regulated genes and suggested that both hormones may play crucial roles in light-dependent anthocyanin biosynthesis. Therefore, further exploration of the interplay of GA with ABA and other hormones during the ripening process of non-climacteric fruits is required. In grapevine, some transcriptomic studies of the effect of GA 3 treatment during the initial fruit development have been performed 21 . In loquat (Eriobotrya japonica), GA 3 applied at bloom changed the expression of ABA and auxin biosynthesisrelated genes 22 . Finally, RNA-seq performed in chinese sour cherry (Cerasus pseudocerasus) identi ed transcriptomic changes in response to GA 3 -induced parthenocarpy 23 . These studies focus on initial fruit development in which the GA participation is quite clear; however the GA's role in the ripening process at the transcriptomic level has not been addressed to date.
Here we analyze the effect of GA 3 at the physiological and global gene expression level in early-and mid-season sweet cherry varieties. This work aims to contribute to the understanding of the participation of GA in the fruit ripening process in non-climacteric fruits.

Results
'Celeste' and 'Bing' are described as early-and mid-season varieties 24 . Both varieties were characterized for fruit growth and phenology in two seasons ( Fig. 1A and Fig. S1A and Tables S1 to S4). In the 2017-2018 season, the slow growth period was more prolonged in the mid-season variety, in which growth resumption occurred from 39 DAFB (Fig. 1A). On the other hand, pink color initiation started at 40 DAFB in the early-season variety (Fig. 1B), whereas in the mid-season variety, it started at 50 DAFB ( Fig. 1A and Table S2). GA 4 and GA 1 were detected in the fruits of both varieties from 34 to 44 DAFB (Fig. S1B).
Several GAs are present at the time of fruit color initiation in sweet cherry and other non-climacteric fruits 5,10 ; therefore, we aimed to investigate the role of GA in ripening triggering by altering GA homeostasis. For this, GA 3 treatment was performed when both varieties transitioned from light green to straw yellow fruit color (35 DAFB the 2017-2018 season and 34 DAFB the 2018-2019 season). At the time of GA 3 application, both varieties were phenologically similar ( Fig. 1; Tables S1 to S4). GA 3 affected ripening parameters at fully ripeness in both varieties, in the 2017-2018 seasons (Table 1), including rmness in both varieties, early-season variety fruit weight, and mid-season variety SSC. Similar results were observed in the 2018-2019 season (Table S5). To further characterize the effect of GA 3 , IAD was measured from the date of the GA 3 treatment in both varieties ( Fig. 3).
IAD is a non-destructive ripening index utilized in sweet cherry, correlated with the anthocyanin content 25 . It reveals the presence of phenolic compounds that act as a screen for chlorophyll absorbance by measuring absorbance differences between 560 nm, 640 nm, and the 750 nm reference 26 . After the GA 3 treatment, the early-season variety had a lower IAD value only at harvest in both seasons ( Fig. 3A and Fig. S3A), whereas the mid-season variety had a delayed increase in the IAD values ( Fig. 3B and Fig. S3B). This delay observed in the mid-season variety was accompanied by color differences between control-and GA 3 -treated fruits as soon as 15 days after the treatment (50 DAFB), whereas control-and GA 3 -treated fruits of the early-season variety had no differences in color at the same date ( Fig. 3C).
As both varieties were differentially affected by GA 3, we explored the transcriptomic effect produced by GA 3 treatment through RNA-seq analysis. For this, we analyzed fruit samples collected minutes before the GA 3 treatment in the 2017-2018 season (T0, 35 DAFB), and control-and GA 3 -treated fruit samples collected four days after the GA 3 treatment (CT4 and GT4, respectively, 39 DAFB). In total, 18 samples were sequenced, nine per variety (three points: T0, CT4 and GT4, and three replicates) that yielded a total of 859 million ltered reads (84 Gb of data), with an average length of 98 bp, and a similar quantity of mapped reads to the P. avium transcriptome reference that was ca. 68% of the total ltered read count (Table S6).
Clustering of the 100 genes having the highest expression levels revealed that samples behaved as expected, i.e., replicates grouped, except for one replicate of T0 of early-season variety (Sample_ID: CT0_3) that was excluded from subsequent analysis (Fig. S4).
DEGs analysis was performed on the genes that changed between T0 and CT4 (CT4-T0) and between T0 and GT4 (GT4-T0) with a cutoff value of two-fold change and FDR < 0.05. Up and downregulated genes were identi ed in both varieties (Fig. 4). In the early-season variety, 381 and 808 genes changed only in GT4-T0 and CT4-T0, respectively ( Fig. 4A). In mid-season variety, 249 and 495 genes changed only in GT4-T0 and CT4-T0, respectively (Fig. 4B). There were 482 and 729 genes in the early-season variety, and 563 and 314 genes in the mid-season variety, that changed independently of the GA 3 treatment (they are in the intercepts of the diagrams). In contrast, the 381 GT4-T0 genes plus the 808 CT4-T0 genes in the early-season variety, and the 249 GT4-T0 genes plus the 495 CT4-T0 genes in the mid-season variety, were GA 3 -modulated. Therefore, they were used in the subsequent GO enrichment analysis focused on the Biological Process categories.
In the early-season variety, 'response to stimulus', 'response to hormone' and 'response to auxin' were some of the overrepresented GO terms (FDR < 0.05) in the GT4-T0 comparison. On the contrary, 'protein phosphorylation' and 'photosynthesis' were some of the overrepresented GO terms (FDR < 0.05) in the CT4-T0 comparison ( Table 2). *GT4-T0 up and CT4-T0 down comparisons of mid-season variety had more overrepresented categories than GT4-T0 total and CT4-T0 total, respectively; therefore, they were included as these comparisons are more informative.
Regarding the mid-season variety, the GO term 'carbohydrate metabolic process' was overrepresented (FDR < 0.05) in the GT4-T0 comparison with upregulated genes, whereas 'oxidation-reduction process' and 'response to stress', among others, were overrepresented GO terms (FDR < 0.05) in the CT4-T0 comparison that includes downregulated genes ( Table 2). In general, GO terms were speci c to a particular comparison, except for the GO category 'oxidation-reduction process' in the early-season variety (Fig. 5A) and 'carbohydrate metabolic process' in the mid-season variety (Fig. 5B), with downregulated genes in CT4-T0 and upregulated genes in the GT4-T0. The GO term 'metabolic process' had the highest number of DEGs modulated by GA 3 in the mid-season variety (Fig. 5B). In the early-season variety, this GO term and 'cellular process' were the most abundant (Fig. 5A).
CT4-T0 up and down comparisons contain GA-modulated genes that up or downregulated from T0 to T4, but they no longer changed after the GA 3 treatment. In the early-season variety, the CT4-T0 down comparison contained DEGs in 'photosynthesis', 'protein phosphorylation', and 'transmembrane transport' GO categories, among others (Fig. 5A). In contrast, in the mid-season variety, it contained DEGs in the 'response to stress' and 'plant-type secondary cell wall biogenesis' GO categories, among others (Fig. 5B). On the other hand, the 'carotenoid metabolism' GO category contained GA-modulated genes that upregulated from T0 to T4 in the early-season variety.
GT4-T0 up and down comparisons contain GA-modulated genes that did not change from T0 to T4 but are up or After that, we investigated the variety-speci c genes that respond to GA 3 treatment (Fig. 6). In the GT4-T0 comparison of the early-season variety, 436 and 652 genes were up and downregulated, respectively. On the other hand, 477 and 145 genes were up and downregulated in the mid-season variety.
As shown in Table 3, the GO terms 'response to auxin', 'response to stress', and 'DNA packaging', among others, were overrepresented (FDR < 0.01), with genes that exclusively changed in the early-season variety. On the other hand, in the mid-season variety, 'regulation of transcription, DNA-templated' and 'regulation of nucleic acid-templated transcription' were some of the overrepresented GO terms (FDR < 0.01) in the upregulated DEGs that were exclusive of this variety In Fig. 7, the enriched GO terms that had more DEGs unique to early-season variety were 'response to stimulus', 'oxidation-reduction' process', 'cellular component organization or biogenesis', 'chromosome organization', 'cell cycle', 'response to hormone', 'response to auxin', among others. In contrast, the GO terms that had more DEGs unique to midseason variety were 'biological regulation', 'aromatic compound biosynthetic process', 'regulation of transcription, DNA-templated', and 'regulation of nucleic acid-templated transcription', among others (Fig. 7).
As 'response to hormone' was an overrepresented GO, a directed search of hormone-related genes was performed in the CT4-T0 and GT4-T0 comparisons of each variety (Table 4). Table 4 Differential gene expression as fold change of genes related to hormone pathways in pairwise comparisons CT4-T0 and GT4-T0 of each variety. In gray, comparisons having genes whose absolute fold change value was at least two. Several gibberellin 2-beta-dioxygenase related genes were downregulated in the CT4-T0 and GT4-T0 comparisons only in the early-season variety. Regarding the ABA biosynthetic pathway, two abscisic acid 8'-hydroxylase related genes were downregulated in the early-season variety. Concerning the ABA response, several abscisic stress ripeningrelated genes were upregulated in the CT4-T0 comparison only in the early-season variety, which was less marked in the GT4-T0 comparison. Regarding the auxin response, several genes that encode putative auxin AuxIAA response repressors were induced in the CT4-T0 comparison of the early-season variety, and they were less upregulated in the GT4-T0 comparison. Some genes coding for auxin-induced protein 15A and 15A-like were upregulated in CT4-T0 and were upregulated even more in GT4-T0, whereas other groups of auxin-induced protein 15A-like were downregulated in the mid-season variety. In contrast, several ethylene-responsive transcription factor genes were more upregulated in the mid-season variety than in the early-season variety.
Given that 'negative regulation of gene expression' and 'DNA packaging' and 'nucleosome assembly' GO terms were overrepresented, a directed search of genes related to these categories was performed (Table 5). Many genes coding for MYB44-like transcription factor were upregulated in CT4-T0, but they were less induced in the GT4-T0 comparison in the mid-season variety. In contrast, other transcription factor genes were downregulated in the early-season variety, including a gene coding for a putative tannin-related R2R3 MYB. Finally, both varieties had downregulated genes encoding putative histone methyl transferases and a DNA (cytosine-5)-methyltransferase 1-like in the CT4-T0, which were slightly less expressed in the GT4-T0 comparison.

Discussion
Physiological differences between early-and mid-season varieties Therefore, the varieties utilized in this work differed mainly in their fruit developmental length and, more speci cally, in the duration of the slow-growth period that occurs during the light green-to-straw yellow transition, which was longer in the mid-season variety (Fig. 1A), similar to a previous report 14 . These results imply that the ripening processes, including color initiation, begin earlier in the early-season varieties. In fact, in the early-season variety, pink color initiates at 40 DAFB (Fig. 1B), whereas in the mid-season variety, it occurred from 50 DAFB onwards, several days after growth resumption (Fig. 1A).
GA 1 and GA 4 have been reported to be present at the onset of ripening in sweet cherry fruits, and GA 4 negatively correlates with ripening parameters 5 . We found that the color initiation differences between both varieties were accompanied by differences in the GA 4 and the GA 1 content. For instance, the early-season variety, which colors rst, had lower GA 4 levels than the mid-season variety at 38 and 44 DAFB (Fig. S1B); and GA 1 levels were higher in the midseason variety at 38 DAFB (Fig. S1C); therefore, possibly low levels of GAs are needed for the ripening triggering. In strawberry, another non-climacteric fruit, color change coincides with the GA 4 decrease between the white and the red stage 10 . In the non-climacteric sweet pepper (Capsicum annuum L.), silencing of DNA methylase gene CaMET1-like1 caused repressed expression of GA biosynthesis genes and led to premature ripening 27 . Taking this evidence together, GA could be a negative regulator of ripening, explaining the differences in the timing of the color development initiation between contrasting sweet cherry varieties.

Differential effect of GA 3 on fruit quality parameters between earlyand mid-season varieties
To further investigate the GA's role in ripening triggering, we altered GA homeostasis by treating fruits with GA 3 , a GA commonly used in agronomic practices. The treatment increased rmness in both varieties. The rmness effect of GA 3 may be associated with changes in cell wall composition. Kondo and Danjo 13 found that during sweet cherry fruit ripening, cell wall-bound neutral sugars -galactose, arabinose, among others -were reduced, which is associated with fruiting softening, and that GA 3 treatment prevented these decrease, supporting the role of GA in fruit rmness maintenance.
It is commonly accepted that early-season varieties are unresponsive to GA 3 ; this is based on some reports showing no effect of GA 3 on some ripening related parameters. For instance, Choi et al. 14 found that GA 3 did not affect fruit size and rmness in early-ripening genotypes. In contrast, we found that GA 3 impaired color, reduced IAD and increased fruit weight at harvest in the two seasons evaluated (Fig. 3A and Fig. S3A; Table 1 and Table S5), similar to our recent report in other early-season variety, where IAD was lower in GA 3 -treated fruits compared with control fruits 17 . In the early-season variety Satohnishiki, GA 3 treatment decreased anthocyanin fruit content and modi ed the sugar accumulation pattern 13 .
On the other hand, it is usually a rmed that GA 3 affects ripening-related parameters in mid-or late-season varieties; this is because of a delay in softening, fruit size increase, and soluble solids accumulation 14,15,16 . However, since the fruits ripen normally and only need more time to acquire the desired features, it would be correct to refer to a delay.
The idea of a ripening impairment arises from studies that analyzed only the harvest point and found less color in the GA 3 -treated fruits 16,28 . Therefore, the monitoring of fruit-ripening parameters is crucial for a better understanding of the physiological effect of GA 3 on late-or mid-season varieties. For this, we used the IAD ripening index 26 , already utilized in the context of sweet cherry fruit ripening 25 . We observed a color difference as soon as 15 days after the GA 3 treatment (Fig. 3C), accompanied by an IAD delay in the GA 3 -treated fruits of the mid-season variety (Fig. 3B). This delay was also found in other mid-season variety 17 , which supports the idea that mid-season varieties respond differently to GA 3 treatment than early-season varieties. GA 3 -treated fruits of the mid-season variety had less color, anthocyanin content, and IAD value at harvest ( Fig. 2 and S2). However, the IAD value of GA 3 -treated fruits reached similar control values four or ve days later ( Fig. 3 and S3), suggesting that there is no ripening impairment in the midseason variety; rather, a normal ripening process occurs but is delayed. A previous report shown that late-and midseason varieties had delayed color development after the GA 3 treatment, and maturity was two to three days delayed in the GA 3 -treated fruits 16 .

GA ripening control differs between early-and mid-season varieties at the transcriptomic level
The differential physiological response to exogenous GA 3 between early-and mid-season varieties, and the advanced color initiation in the early-season variety, which is accompanied by less GAs levels, led us to hypothesize that the ripening process of these contrasting varieties is very different at the molecular level, and possibly it may be differentially modulated by GA. Therefore, we studied the effect of GA at the transcriptomic level by altering the GA homeostasis with exogenous application of GA 3 . This treatment was applied on the same date in both varieties in the 2017-2018 season, when they were at a similar physiological stage, as revealed by the IAD parameter ( Fig. 3A and B); allowing the transcriptomic data to be comparable between varieties. The date selected for the treatment was 35 DAFB in the 2017-2018 season, ve and 15 days before early-and mid-season variety color initiation, respectively, to characterize the molecular events modulated by GA at the onset of ripening and the early differences between both varieties.
To characterize the effect of GA 3 , we considered CT4-T0 and GT4-T0 comparisons, which contain GA-modulated genes since they exclude the genes that changed independently of the GA 3 treatment (the intercepts of the diagrams; Fig. 4). As mentioned, CT4-T0 up and down comparisons contain genes that changed from T0 to T4 but did not change in the GT4-T0 comparison, i.e., GA-modulated genes since they no longer changed after the GA 3 treatment.
Therefore, it can be deduced that GA 3 avoids their up or downregulation.
GT4-T0 up and down comparisons also contain GA-modulated, since their gene expression did not change from T0 to T4 in the control fruits, but the GA 3 treatment up or downregulated them.
The GO terms 'protein phosphorylation' and 'photosynthesis' were overrepresented in the CT4-T0 comparison of the early-season variety (Table 2), with genes downregulated in this comparison (Fig. 5A). This nding means that these genes downregulated from T0 to T4 in the control fruits, but they no longer changed after the GA 3 application.
Therefore, it is deduced that GA 3 impedes this downregulation, which was expected since the photosynthetic rate and capacity are usually positively controlled by GAs 29,30 .
On the other hand, we observed that 'auxin response' was an overrepresented GO category in the GT4-T0 comparison, which means that the transcript levels of the genes in this category only changed after the GA 3 -treatment. The genes in the overrepresented GO category 'auxin response' upregulated in the GT4-T0 comparison (Fig. 5A), and some of them were exclusive of the early-season variety (Fig. 7). Several putative genes of the Aux/IAA gene family coding for auxin-responsive proteins IAA29, IAA13, IAA6, among others, were upregulated from T0 to T4. These are mainly repressors of auxin responses in fruits 31 , therefore in the early-season variety, the auxin responses should turn off as ripening is initiating, since auxin usually antagonizes the ABA and ethylene promoting effect 6,7 . In this regard, auxin response should be activated by GA, since several genes encoding putative Aux/IAA repressors were downregulated in the GT4-T0 comparison. Reduced expression of these repressors should activate the auxin response; therefore, GA may be a positive regulator of the auxin pathway.
Interestingly, several genes encoding auxin-induced protein 15A and 15A-like were more upregulated in GT4-T0 comparison of the early-season variety than in the CT4-T0 comparison, whereas other different genes encoding a protein 15A-like were downregulated in the mid-season variety. These results demonstrate that there is an antagonistic IAA response between both varieties. It also suggests a gene specialization, where only early-season auxin-induced protein 15A and 15A-like genes might be GA-responsive.
Carotenoids are precursors of ABA, a key regulatory hormone that triggers ripening in non-climacteric fruits 2,32,33 .
Several genes encoding putative ASR (abscisic acid stress ripening proteins) were upregulated only in the earlyseason variety. The ASR genes encode transcription factors induced by ABA and abiotic stress in several plant tissues, including fruits 34 . Additionally, two genes encoding an ABA degrading enzyme, abscisic acid 8'-hydroxylase 1 (CYP707A1), were more downregulated in the CT4-T0 than in the GT4-T0 comparison; hence, it is deduced that GA is a positive regulator of these genes. The silencing of PavCYP707A2 in sweet cherry led to a delay in the transcript accumulation of several genes of the ABA and anthocyanin synthetic pathways 35 , supporting its role as a negative regulator of ripening. GA 3 enhances its expression according to our results (Table 5), in agree to our recent report showing that GA 3 increases the transcript abundance of PavCYP707A2 ve days after the treatment, but only in the early-season variety 17 .
In the mid-season variety, the GO term 'carbohydrate metabolic process' was overrepresented in both the GT4-T0 and the CT4-T0 comparisons (Table 2). Interestingly, this category's genes were downregulated in the CT4-T0 comparison, whereas they were upregulated in the GT4-T0 comparison (Fig. 5B), implying that GA strongly activates these genes. Several genes encoding putative xyloglucan endotransglucosylase/hydrolases were strongly upregulated by GA 3 in the mid-season variety (data not shown). These enzymes are cell wall-modifying proteins and loosen cell walls hence conferring cell wall extensibility. Interestingly, transient size increase was observed only in the mid-season variety four days after the GA 3 treatment (Fig. S5).
Another overrepresented GO term in the mid-season variety was 'response to stress' ( Table 2). Genes in this category were downregulated in CT4-T0 but not in GT4-T0. Therefore, it is deduced that GA prevents their downregulation. A gene coding for a putative ethylene biosynthetic gene, 1-aminocyclopropane-1-carboxylate oxidase homolog 1-like, was also downregulated in CT4-T0 but not in GT4-T0. Possibly GA has a positive effect on ethylene synthesis, and this hormone, in turn, maintains activated some stress responses. Ethylene seems to play a promoting role in nonclimacteric fruit ripening 36 ; therefore, it is intriguing that GA, a hormone that antagonizes ripening, increases ethylene levels. This supports the idea that GA, in late-or mid-season varieties, does not negatively affect ripening. Additionally, there is a more increased ethylene response in the mid-season variety; hence, ethylene pathway might be signi cant in this variety.
The GO terms 'regulation of transcription, DNA-templated' and 'regulation of nucleic acid-templated transcription' were overrepresented in the mid-season variety, and they included genes that exclusively changed in this variety (Table 3; Fig. 7). Several genes coding for MYB44-like transcription factors were upregulated only in the mid-season variety (Table 5). Additionally, they were GA-downregulated as they were less upregulated in the GT4-T0 comparison. Also, a gene coding for an MYB44 transcription factor was present in both varieties, but it was more strongly expressed in the mid-season variety. MYB44 proteins suppress the expression of PP2C genes encoding type 2C protein phosphatases 37 . PP2Cs are key proteins of the ABA regulatory network and are negative regulators of subfamily 2 of SNF1-related kinases (SnRK2s), therefore negative regulators of ABA responses 38 . PavSnRK2s and PavPP2Cs are expressed during sweet cherry fruit development 8 . Since MYB44 related genes were upregulated in the mid-season variety, they could positively regulate the ABA response.
Interestingly, in the GT4-T0 comparison, these genes were less upregulated than in the CT4-T0 comparison, suggesting that GA downregulates the expression of MYB44 genes. This could be a mechanism for GA to antagonize the ABA pathway. MYB44 transduce environmental signals and mediate stress responses 39 . It also mediates the interaction of hormones; for instance, it regulates the antagonistic interaction between salicylic acid and jasmonic acid pathways 39 , besides the modulation of the ABA signaling. Therefore, it seems to be a crucial regulatory element for converging pathways.

Ripening differences between maturity time contrasting varieties: what is in common and what is different?
We observed that some ripening-related molecular features were more activated in the early-season variety, including genes related to photosynthesis, carotenoid, and ABA pathway. Photosynthesis is a negative biomarker of ripening 40 , and in sweet cherry, chlorophyll decrease is one of the rst events at the ripening initiation, even before sugar and anthocyanin content increase 33 . 'Photosynthesis' was an overrepresented category in the early-season variety (Table 2), and several photosynthesis-related genes were downregulated in the CT4-T0 comparison (Fig. 5A). This might be exclusive of the the early-season phenotype, as this GO term was not overrepresented in the mid-season variety.
Auxin counteracts the ripening processes 2,6 . Here we found that a probable indole-3-pyruvate monooxygenase YUCCA10 was upregulated in both varieties from T0 to CT4. This agrees with the increase in the IAA content detected in both varieties from 34 to 44 DAFB (Fig. S6A). Therefore, possibly this hormone prevents that the ABA ripening triggering occur too early. On the other hand, GA seems to control the IAA response, as discussed above, rather than the IAA biosynthesis, with differences between varieties in the expression of auxin induced genes.
Carotenoid and ABA biosynthesis is a positive biomarker of ripening 40 . In the early-season variety, the 'carotenoid metabolism' category contained genes upregulated in the CT4-T0 comparison (Fig. 5A). Phytoene synthase and NCED1 like genes were upregulated in this comparison in both varieties (Table 4), possibly involved in carotenoid production and ABA biosynthesis. Additionally, in the CT4-T0 comparison, a probable carotenoid cleavage dioxygenase 4 chloroplastic gene was upregulated in the early-season variety, which probably degrades carotenoids, since in chrysanthemum (Chrysanthemum morifolium) a possible orthologue of this gene contributes to white petal color 41 . Regarding the ABA biosynthesis, in sweet cherry, NCED1 transcript variations correlate with ABA increase during ripening 4 . The NCED1 upregulation in CT4-T0 implies that ABA increases in the early-season variety at the onset of ripening. We found that ABA content increases fromt 34 to 44 DAFB in the early-season variety, which occurred later in the mid-season variety (Fig. S6B).
Finally, phenylpropanoid early branches are negative biomarkers of non-climacteric ripening 40 . In this regard, we found that a tannin-related R2R3 MYB transcription factor was strongly repressed in the early-season variety, compared with the mid-season variety.
Chromosome reorganization, changes in histone methylation pro le and DNA hypomethylation are characteristic features of ripening initiation 42 . For instance, anthocyanin-de cient apple mutant had higher methylation levels in the promoter of a MdMYB gene 43 . Silencing of a methylase gene in the non-climacteric sweet pepper led to premature ripening 27 . We observed several downregulated genes in the CT4-T0 comparison encoding putative histone modi cation enzymes and DNA (cytosine-5)-methyltransferase 1-like (Table 5). These genes were downregulated in both varieties, suggesting that chromosome remodeling and DNA hypomethylation could be relevant at this stage of development. However, some were slightly more downregulated in the early-season varieties, suggesting that these processes could be advanced in this variety, which could lead to earlier ripening initiation in early-season varieties.
Additionally, there are some early-season variety-exclusive genes, including histone-lysine N-methyltransferase ATXR6, which is involved in the transcriptional repression in Arabidopsis 44 . On the other hand, some genes in the 'DNA packaging' GO categories were exclusive of the early-season variety ( Table 3), suggesting that though DNA related changes may occur in both varieties, they have variety-dependent particularities.
The GA pathway seemed to be downregulated earlier in the early-season variety. Supporting this idea, several gibberellin 2-beta-dioxygenase related genes were downregulated from T0 to T4 in the early-season variety (Table 4).
These genes encode putative enzymes involved in GA degradation, previously characterized in grapevine fruits 9 . With low GA levels, these genes downregulate as they act in a negative feedback for controlling GA levels. Therefore, possibly a low GA content may be required for the onset of ripening to occur in the early-season variety, as it was observed ( Fig. S1B and Fig. S1C).
GA 3 treatment upregulates the DNA (cytosine-5)-methyltransferase in the early-season variety, so possibly higher methylation levels occur when GA content is elevated, thus impairing the ripening processes. The effect of GA 3 on IAA and ABA pathways in the early-season variety, as discussed above, also supports that GA exerts a negative effect on the ripening processes.
Higher GA levels could be retarding the onset of ripening in the mid-season variety, possibly to reach full embryo development. The regulatory module controlling the timing of ripening initiation could involve GA repression of MYB44-like gene expression, a putative negative regulator of ABA signaling genes, PP2Cs. Less expression of MYB44like genes could lead to accumulation of negative PP2Cs transcripts. Therefore, more PP2Cs expression could imply that more ABA is needed to surpass this effect, and this could explain the ripening delay upon the GA 3 treatment. In the mid-maturing sweet cherry variety Lapins, GA 3 increases PavPP2C3 and PavPP2C4 trancripts and delays ripening 17 , which supports our ndings.
The summary of possible molecular interactions controlling ripening in early-and mid-season varieties is depicted in were selected from each tree, three control branches and three GA 3 -treated branches, following Usenik's methodology 16 . GA 3 (ProGibb ® 40% SG) was dissolved in water and applied to individual branches with a hand sprayer to run-off at a rate of 30 ppm when the fruits of each variety were in late Stage II of development, and fruit color transitioned from green to straw yellow 28 . Control branches were treated with water and protected from spraying with GA 3 , according to Usenik et al. 16 . The treatment was performed at the same hour of the day in both varieties and both seasons (12:00 GMT). Ten fruits per repeat (branch), i.e. 30 fruits per treatment, were sampled, with sample size according to Luo et al. 32 . Fruit samples were homogenous in color, size and form, and without any visible defect. The fruits were immediately frozen in liquid nitrogen after removal from the tree and stored at -80°C until used for RNA-seq analysis, anthocyanins estimations, and hormone quanti cation.

Fruit parameters
For fruit width, a caliper was used to measure the equatorial diameter at the fruits' widest point 14

Analysis of differentially expressed genes
Adaptors, low-quality bases, and short sequences trimming were performed using CLC Genomics Workbench version 11.0.1 following parameters: Q ≥ 20; no more than two ambiguities; nal read length ≥ 50 bp. Sequence mapping to the reference genome 49 was performed using CLC Genomics Workbench version 11.0.1 with the following parameters: similarity = 0.9; length fraction = 0.6; insertion/deletion cost = 3; mismatch cost = 3, and unspeci c match limit = 10.
Expression levels were normalized by calculating RPKM (Reads Per Kilobase Million) value. Transcript abundances were tted using a general linear model (GLM) and differential expression of treatments tested with the Wald test against control groups 50 . Differentially expressed genes (DEGs) were de ned as having an absolute fold-change value of at least two between T4 and T0, with an adjusted p-value using a false discovery rate (FDR) with at least a 95% signi cance 51 .
Functional annotation of P. avium reference transcripts 49 was performed by BLAST2GO software version 5.2.5 52 .
First, a BLASTx search was performed against the NCBI NR database 53 with an e-value cutoff of 1e-6 and HSP length cutoff of 33. Then, INTERPROSCAN analysis 54 was performed with BLAST2GO default parameters. Finally, the data from BLAST searches, INTERPROSCAN terms, enzyme classi cation codes (EC), and metabolic pathways (KEGG, Kyoto Encyclopedia of Genes and Genomes) were merged in gene ontology (GO) terms for a comprehensive functional range cover in the functional annotation. The BLAST2GO program defaults were used in all mapping and annotation steps, and the false discovery rate (FDR) cutoff was set to a 0.05% probability level. GO overrepresentation analysis was performed with the Fisher's Exact Test Enrichment Analysis using BLAST2GO tools and integrated FatiGO package 55 with default parameters.
Venn diagrams were constructed using the online tool VENNY 56  acetic acid solution containing internal standards (deuterium-labeled hormones; OlChemim Ltd., Olomouc, Czech Republic). The mix was shaken for one hour at 4ºC, and the extracted fraction was maintained at -20ºC overnight. The samples were centrifuged, and the supernatant was vacuum dried and then dissolved in 1% acetic acid. A reversephase column (OasisHLB) was used 57 , and the eluate was dried and dissolved in 5% acetonitrile − 1% acetic acid. An autosampler and reverse-phase UHPLC chromatography column, 2.6 µm Accucore RP-MS, 100 mm x 2.1 mm (ThermoFisher Scienti c, San Diego, CA, USA) were used. Then the hormones were separated using a gradient of acetonitrile (2%-55%) containing 0.05% acetic acid, at a rate of 400 µL/min over 22 min. ABA, IAA, GA 1

Statistical analysis
One-way ANOVA analysis with Tukey's post hoc test was used for assessing differences between control-and GA 3treated fruits in estimated anthocyanins and fruit parameters,where three replicates per date were used to establish the mean differences, using the InfoStat software 59 .