VcRR2 regulates chilling-mediated flowering through expression of hormone genes in a transgenic blueberry mutant

The molecular mechanism underlying dormancy release and the induction of flowering remains poorly understood in woody plants. Mu-legacy is a valuable blueberry mutant, in which a transgene insertion caused increased expression of a RESPONSE REGULATOR 2-like gene (VcRR2). Mu-legacy plants, compared with nontransgenic ‘Legacy’ plants, show dwarfing, promotion of flower bud formation, and can flower under nonchilling conditions. We conducted transcriptomic comparisons in leaves, chilled and nonchilled flowering buds, and late-pink buds, and analyzed a total of 41 metabolites of six groups of hormones in leaf tissues of both Mu-legacy and ‘Legacy’ plants. These analyses uncovered that increased VcRR2 expression promotes the expression of a homolog of Arabidopsis thaliana ENT-COPALYL DIPHOSPHATE SYNTHETASE 1 (VcGA1), which induces new homeostasis of hormones, including increased gibberellin 4 (GA4) levels in Mu-legacy leaves. Consequently, increased expression of VcRR2 and VcGA1, which function in cytokinin responses and gibberellin synthesis, respectively, initiated the reduction in plant height and the enhancement of flower bud formation of the Mu-legacy plants through interactions of multiple approaches. In nonchilled flower buds, 29 differentially expressed transcripts of 17 genes of five groups of hormones were identified in transcriptome comparisons between Mu-legacy and ‘Legacy’ plants, of which 22 were chilling responsive. Thus, these analyses suggest that increased expression of VcRR2 was collectively responsible for promoting flower bud formation in highbush blueberry under nonchilling conditions. We report here for the first time the importance of VcRR2 to induce a suite of downstream hormones that promote flowering in woody plants.


Introduction
Deciduous woody plants have evolved to be winter dormancy-dependent for survival under freezing temperatures and require accumulation of a minimum level of chilling (chilling requirement) to allow resumption of normal growing and flowering in coming seasons [1][2][3][4] . Photoperiod and temperature are the major environmental factors in inducing and maintaining dormancy 1,3,5 .
Insufficient chill can reduce fruit production by reducing bud break and lessening flower quality 6 . With climate changes, declining chilling has been occurring for decades and is predicted to continue to do so; and consequently, insufficient winter chill has been recognized as a potential limiting factor on fruit production for both wild and cultivated species [7][8][9] .
Flowering is a prerequisite for fruiting and plays a significant role in the life cycle of angiosperms 4,10 . Vernalization/chilling is necessary for normal flowering of winter-annual and deciduous perennial plants, and several critical genes controlling the process have well been characterized in annual species 4 . For example, FLOW-ERING LOCUS C (FLC), a MADS-BOX gene, plays a central role in the vernalization pathway of winter ecotypes of annual A. thaliana and is a critical component of the vernalization regulatory loop of FRIGIDA (FRI)-FLC-FT (FLOWERING LOCUS T), which seems being conserved in the family Brassicaceae [10][11][12] . Similarly, in vernalization-requiring wheat and barley, a VRN1-VRN2-TaFT (VRN3) loop regulates vernalizationmediated flowering 13 . In parallel, a cluster of six MADSbox transcription factors identified through forward genetics, analogous to FLC in A. thaliana and VRN2 in cereals and termed as DORMANCY-ASSOCIATED MADS-BOX (DAM), are involved in dormancy regulation in peach (Prunus persica) 14,15 . These DAM genes show high similarities to A. thaliana AGAMOUS-LIKE 24 (AGL24) and SHORT VEGETATIVE PHASE (SVP) genes [16][17][18] . Functional analysis of the DAMs through reverse genetics has not been reported in peach. To date, in woody plants, there remains a lack of understanding of the physiological, molecular, and genetical basis of chilling-regulated flowering 7 .
Highbush blueberry (Vaccinium corymbosum L.) is a major cultivated Vaccinium crop in the family Ericaceae (Syn. Heath) containing ∼450 species 19,20 . As a perennial shrub, blueberry floral bud initiation often starts before endodormancy 21 . Enough chilling accumulation during endodormancy, ranging from 150 to 1500 chill units, is required for breaking dormancy, and insufficient chilling prevents normal bud break and leads to reduced blueberry production 20 . Compared with annual and other woody species, the knowledge of genetic and molecular control of chilling requirements for flowering is very limited in blueberry. Our recent comparative transcriptome analysis of nonchilled, chilled, and late-pink bud revealed that the orthologs of many well-known flowering pathway genes in annual species, such as FT, PROTEIN FD (FD), TERMINAL FLOWER 1 (TFL1), and LEAFY (LFY), MADS-box genes, and hormone genes, were involved in chilling-mediated blueberry bud break in blueberry 22 . However, like what was noted in other woody plants, functional FLC orthologs were not found in blueberries [22][23][24] .
Hormone genes also played roles in blueberry flowering 23,24 . Overexpression of a blueberry FT (VcFT: hereafter, Vc included in front of any gene refers to the ortholog of the gene in blueberry) induced differential expression of a group of hormone pathway genes in transgenic blueberry leaves 23 . Recently, we identified a unique transgenic blueberry event (line) in the background of the cultivar 'Legacy' and this line, named as Mu-legacy, was found to be a valuable genetic stock for studying chilling-induced flowering 25 . Unlike nontransgenic 'Legacy' and other transgenic events, Mulegacy and its T 1 transgenic plants had early flower bud formation and were able to flower under nonchilling conditions. Our initial characterization of this mutant, compared with the nontransgenic control and other transgenic events, revealed several hormone-related pathway genes DE and concluded that hormone changes were likely responsible for the phenotypic changes in Mulegacy 25 . This was consistent with our previous work that hormones, at least at transcript levels, might function as physiological signals and play a significant role in regulating blueberry flowering, especially chilling-mediated flowering 22,23,25,26 .
To reveal the potential roles of VcRR2 and hormone(s) in chilling-mediated blueberry flowering, further analysis of the Mu-legacy was conducted in this study through hormone measurements and RNA sequencing. We revealed profiles of both hormones and DE transcripts (DETs) in leaf and bud tissues of the Mu-legacy plants. We found that the decreased auxin and the increased gibberellin 4 (GA4) contents in Mu-legacy leaves likely contributed to plant dwarfing, and promoted flower bud formation and nonchilled flower bud break of the Mulegacy plants. The Mu-legacy plants provided evidence to show that an increased expression of VcRR2 played a critical role in regulating blueberry flower bud formation and bud break mainly through hormone pathway genes but not flowering pathway genes. The results provided further evidence that hormones may play critical roles in chilling-mediated flowering in blueberry and possible other woody plants.

Flowering behaviors of Mu-legacy plants
We exanimated the flowering behaviors of nonchilled 6year-old 'Legacy', Legacy-VcDDF1-OX, Mu-legacy, and Mu-legacy-T 1 plants. Under the greenhouse conditions (the lowest temperature >20°C, natural light condition), all the Mu-legacy and the Mu-legacy-T 1 plants, in contrast to the 'Legacy' plants, showed a new break of some flower buds in mid-October of 2017 (Fig. 1); leaf tissues at this developmental stage for both the 'Legacy' and the Mu-legacy plants were collected for transcriptome and hormone analyses. The results are consistent with our previous observation 25 . The continuous flowering of the Mu-legacy plants lasted for about 9 months (from October to June) 25 , suggesting that the transgenes and insertion position in the Mu-legacy are responsible for the reduced plant size, the promoted flower bud formation, and the reduced chilling requirement for flowering. Interestingly, the observed plant dwarfing is similar to a typical phenotype of GA deficiency in Arabidopsis 27 , while the promoted flower bud formation and reduced chilling requirement for flowering are similar to the phenotypic changes of early flowering and seed dormancy induced by GA-overproduction phenotypes in Arabidopsis 28 .

Contents of major hormones in Mu-legacy leaves
Compared with 'Legacy', Mu-legacy carries a transgenic, overexpressed blueberry DWARF AND DELAYED FLOWERING 1 (DDF1) (VcDDF1), as well as a mutated and overexpressed blueberry RESPONSE REGULATOR 2like gene (VcRR2), which accidently resulted from the position effect of the introduction of the transgene VcDDF1 25 . In A. thaliana, the DDF1 activates GA2OX7 by binding to the DRE-like motifs (GCCGAC and ATC-GAC) of the GA2OX7's promoter 29 . The GA2OX7 catalyzes the decay of active GAs [29][30][31][32][33] . On the other hand, ARABIDOPSIS RESPONSE REGULATOR 2 gene (ARR2) is responsible for cytokinin responses in the cytokinin signaling pathway 34,35 . Since both VcDDF1 and VcRR2 are related to hormone signaling, we expect that overexpression of these two genes, especially the VcDDF1, would change the hormone profiles of the Mu-legacy plants.
A total of 41 metabolites of six groups of hormones were analyzed in and compared between Mu-legacy and 'Legacy' plants ( Fig. 2, Table 1). Of the 14 gibberellins (GAs) that were measured, GA19 was detected in the leaf samples of both Mu-legacy and 'Legacy' plants, and the Mu-legacy leaves, compared with 'Legacy' leaves, showed much decreased GA19 (Fig. 2a). GA20 was detected high in the Mu-legacy sample [74.5 pmol/g dry weight (DW)] and low in two 'Legacy' samples (18.2 and 25.4 pmol/g DW, respectively) ( Table 1). More importantly, a low level of GA4, which is an active GA form 36 , was present in two Mu-legacy samples but absent in all 'Legacy' samples ( Table 1), indicating that the active GA4 was likely responsible for the promotion of flower bud formation as well as the nonchilled flower bud break (Fig. 1a). Seven ABA and ABA metabolites analyzed were all detected. The Mu-legacy leaves had a higher total ABA content (7.8 vs. 4.7 nmol/g DW) than the 'Legacy' leaves ( Fig. 2b). Of the six auxin and auxin metabolites, only IAA was detected, and the Mu-legacy leaves showed a significantly (p = 0.01) lower content than the 'Legacy' leaves (Fig. 2c). The reduced IAA level likely contributed to the reduced size of the Mu-legacy plants (Fig. 1a). Of the ten cytokinin and cytokinin metabolites measured, four were detected in 'Legacy' and Mu-legacy samples; measurable isopentenyladenine (iP) was only found in the 'Legacy' leaves ( Fig. 2d, Table 1). The overall levels of the detected cytokinin metabolites were higher in the 'Legacy' leaves ( Fig. 2d). Salicylic acid (SA), jasmonic acid (JA), methyl jasmonate, and jasmonoyl isoleucine were measured and all detected in both the 'Legacy' and the Mu-legacy leaves. For the contents of JA and JA metabolites, there was no statistical difference between the 'Legacy' and Mu-legacy samples (Fig. 2e). The SA content was slightly higher in the Mu-legacy leaves, although not significant (Fig. 2f). These results showed that overexpressed VcDDF1 and VcRR2 induced content changes of some hormones in the Mu-legacy leaves.

Differentially expressed genes in Mu-legacy leaves
Previous comparisons of transcriptomic profiles of leaves between 'Legacy' and Mu-legacy plants have revealed that the increased expression of VcRR2 likely co-functioned with the overexpressed VcDDF1 for the promoted flowering and dwarfing of the Mu-legacy and Mulegacy-T 1 plants 25 . In this study, we identified 2543 DE transcripts (DETs) when transcriptomic profiles were compared between 'Legacy' and the Mu-legacy leaves and these DETs were much more than the 260 DETs identified in the previous transcriptomic comparison study [25]. Such a large difference likely reflects some variation in the biological samples (e.g., stages and sampling) and sequencing coverage. Nevertheless, among the 260 DETs identified from the previous study, 56 were overlapped  Tables 1 and 2).
Only two DE genes (DEGs) of the flowering pathway, including orthologues of Arabidopsis FRUITFULL (FUL) and Arabidopsis RAV1, showed up-and down-regulation, respectively; both of which had repressive effects on flowering bud formation and flowering in the 2018 RNAseq data. Although this result is inconsistent with the previously identified three DE flowering pathway genes (i.e., VcTCP8, VcARP6, and VcNFYC1) 25 , it still supports the previous conclusion that no orthologs of the major Arabidopsis vernalization pathway genes, e.g., VcFT and VcSOC1, showed differential expression. RAV1 encodes an AP2/B3 domain and is a negative regulator of flowering by repressing both FT and GA 37 . If the function of the RAV1 in Arabidopsis is conserved in blueberries, the upregulated VcRAV1 would repress flowering of the Mulegacy plants. FUL is redundant with APETALA1 (AP1) and CAULIFLOWER (CAL) in promoting to the transition to reproductive meristems 38 . The decreased expression of the VcFUL would repress flower bud formation, which is inconsistent with the promoted flower bud formation of the Mu-legacy plants. Overall, the two DE flowering pathway genes were not likely responsible for the promoted flowering in the Mu-legacy plants.
We identified seven DEGs related to the biosynthetic pathways of GA (2 DEGs), SA (3), JA (1), and cytokinin (1) but none for the ABA and auxin pathways ( Table 2, Figs. 2, 3). VcCPS and VcGA20OX1 were the two GA related genes. We found that VcCPS expression increased while VcGA20OX1 expression decreased. Ent-COPALYL DIPHOSPHATE SYNTHASE (CPS) catalyzes the first step of GA biosynthesis 39 . The increased VcCPS could lead to an increase in bioactive GAs, this matched the result of the presence of the GA4 in the Mu-legacy leaves but not in the 'Legacy' leaves 36 . GA20OX1 catalyzes a series of intermediate oxidation reactions during the biosynthesis of GA 36 , the slightly reduced expression of VcGA20OX1 should have reduced GA production. Whether the increased VcCPS and decreased VcGA20OX1 expression were part of the cause for the reduced GA19 and increased GA4 content is unknown. CYTOKININOXIDASE6 (CKX6) catalyzes the degradation of cytokinins in Arabidopsis 40 , and we found that VcCKX6 expressed reduced in this study. The reduced VcCKX6 expression explained well the decreased content of cytokinins observed 41 . Similarly, we observed an increased expression of cytochrome P450 94C1 (CYP94C1), which is involved in the oxidation of jasmonoyl-L-isoleucine (JA-Ile) in Arabidopsis 42 . However, the increased VcCYP94C1 did not lead to an increase of JA-Ile accumulation in this study. On the other hand, we found three SA-related DEGs in this study: SALICY-LATE/BENZOATE CARBOXYL METHYLTRANSFERASE (BSMT1), SGT1, and METHYLESTERASE 17 (MES17). Decreased expressions of these genes were found to lead to accumulation of SA in Arabidopsis 43 , accordingly there was no surprise that the decreased VcBSMT1, VcSGT1, and VcMES17 contributed to the increased SA in the Mulegacy leaves. We observed the same association of increased expression of these genes (i.e., VcBSMT1, VcSGT1, and VcMES17) with increased accumulation of SA in the Mu-legacy leaves.
Because VcRR2 could enhance cytokinin responses 35,44 , we searched for the DETs of the orthologues of the other ARR genes. Two DETs of VcRR9 showed downregulation, which could potentially reduce cytokinin responses (Table 2, Fig. 3). Since cytokinins are often positively correlated to the IAA level 45 , the reduced cytokinin responses supported the reduced IAA content (Fig. 2c).  A typical Legacy-VcDDF1-OX transgenic event, as expected, had the VcDDF1 gene overexpressed in its leaves but did not show altered flowering behavior The VcDDF1 gene was also overexpressed in the leaves of Mulegacy Mu-legacy, which, however, was promoted to flower early. Interestingly, VcRR2 and VcCPS were also overexpressed in Mu-legacy, but not in Legacy-VcDDF1-OX, suggesting that VcRR2 and VcCPS were likely the causes for flowering promotion in Mu-legacy and had no direct connections with the overexpression of the VcDDF1 gene in the background (Table 2). This conclusion was also confirmed in the Mu-legacy-T 1 plants derived from self-pollinated Mu-legacy. The Mu-legacy-T 1 plants showed a similar flowering behavior as the Mulegacy. Very fortunately, in the DEGs between the 'Legacy' and Mu-legacy-T 1 leaves, the VcDDF1 did not showed differential expression but the VcRR2 and the VcCPS were both upregulated (Table 2, Fig. 3). This suggests that the increased VcRR2 and VcCPS promoted early flowering of the Mu-legacy-T 1 plants under nonchilling conditions, most likely, through the gibberellin pathway genes.    additional DETs, which are the orthologs of 17 A. thaliana genes. These genes belong to the biosynthesis pathway genes of five groups of hormones (Table 2, Fig. 3) and are likely involved in the dormancy release of the NB of the Mu-legacy plants.
Eight DETs representing three genes (ABA2-4) in the pathway of ABA synthesis and the A. THALIANA BETA-GLUCOSIDASE 1 gene, which is a positive regulator of ABA-activated signaling pathway, were upregulated. These upregulated DETs could increase ABA content and enhance ABA signaling, thus potentially resulting in delaying flowering, which is contrary to the early flowering of the Mu-legacy (Table 3, Supplementary Table 2). Auxin interacts with ABA in controlling A. thaliana seed dormancy/germination 46 . Eight DETs of four auxin biosynthetic genes exhibited a potential for an increase in auxin biosynthesis (Table 3). Brassinosteroids (BRs) play a critical opposite role to ABA during seed germination and deficiency of BRs is often associated with dwarf phenotypes 47 . Five DETs of three genes in the pathway of BRs biosynthesis, including four upregulated and one downregulated DETs (Table 3), most likely inactivate BRs 48,49 . The reduced active BRs, working with others, can consequently lead to the dwarfing and delayed flowering in A. thaliana 50 . Cytokinins are multifunctional in plants. Four DETs representing three DE cytokinin genes similar to the CYP735A1, UGT85A1, and APT5 genes in Arabidopsis were identified (Table 3), the overall impact of the DETs can lead to cytokinin-deficient plants due to their potential to reduce active cytokinins through enhanced O-glucosylation [51][52][53] . JA is involved in chilling-dependent dormancy release through its interaction with ABA 54 . Four DETs of four genes in the JA biosynthetic pathway were identified (Table 3), and they all showed an upregulation that could either increase or reduce the production of JA 55,56 . It is interesting that unlike in the leaf tissues we did not detect any DEGs of the GA pathway in the transcriptome comparison between the nonchilled buds of the Mu-legacy and the 'Legacy', suggesting that the promoted flowering of the nonchilled Mu-legacy buds had little to do with the expression the GA pathway genes in this case. Apparently, no convincing evidence shows that any of the individual hormone(s) or the DEG of these hormones is responsible for the promoted flowering of the no-chilled Mu-legacy buds.
Since hormone genes are involved in blueberry flower bud dormancy release 22 , we further looked into the expression of the 29 DETs of 17 hormone genes identified in the NB of Mu-legacy plants in CB and late-pink buds (LPB) of both Mu-legacy and 'Legacy' plants (Table 3). In the comparison of CB and NB of 'Legacy', 23 transcripts showed differential expressions, of which 22 showed a consistent up-or down-regulation with those in the 29 DETs resulting from the comparison between NB of Mu-legacy and 'Legacy' (Table 3), supporting that the 29 DETs are responsible for the promoted flowering of the "Mulegacy" plants under nonchilling conditions. In the comparison of CB and LPB of 'Legacy', 18 of the 29 transcripts exhibited differential expression during flower bud break (Table 3). On the other hand, according to the low percentage of bud break and the reduced number of flowers per bud 25 , the 29 DETs did not initiate the normal flowering of nonchilled Mu-legacy plants, most likely because more chilling-driven DETs of other hormone genes or the flowering pathway genes are needed to induce normal flowering (Supplementary Table 2). In addition, two comparisons for Mu-legacy plants, including chilled buds versus nonchilled buds and LPB of the Mu-legacy plants, respectively, revealed 19 and 24 DETs of the 29 transcripts of 17 hormone genes were inducers of the Mulegacy plant flowering under chilling conditions (Table 3). This provides further evidence that the 29 DETs are responsible for flowering of the "Mu-legacy" plants under nonchilling conditions. The 17 DETs identified in the comparison of chilled 'Legacy' buds and nonchilled "Mulegacy" buds indicated that insufficient changes of these transcripts could also be responsible for the unusual flowering behaviors of the nonchilled Mu-legacy plants.  Table 3).
Only six DETs of four flowering pathway genes, including VcFD, VcTFL1, VcARP6, and VcDOF5.3, were detected and were all repressed in the NB of the Mulegacy plants (compared with the nonchilled 'Legacy' buds) 25 . We analyzed expression of the six transcripts in the fully vernalized flower buds and the LPB of the Mulegacy plants. None of the six transcripts showed 'Legacy' and Mu-legacy plants, respectively, were also included. #N/A: No differential expression differential expressions in the comparison of the CB with those of nonchilled CB. In contrast, five of the six transcripts showed differential expressions in the comparison of the CB with the LPB, of which expression of the VcARP6 was repressed and expression of the other three genes were upregulated (Supplementary Table 4). This indicates that the four genes were involved in flowering of the CB of the Mu-legacy plants. In the 'Legacy' plants, we further looked into expression of the six transcripts in both the CB and the LPB. Five transcripts of the four genes were all repressed in the CB (compared with the NB). The result was similar to the six DETs identified in the comparison between the NB of the Mu-legacy plants and the NB of the 'Legacy' plants (Supplementary Table  4). This suggests that the four genes were chilling responsive in the 'Legacy' flower buds and their expression changes in the NB of the Mu-legacy plants were most likely responsible, at least in part, for the flowering of the nonchilled Mu-legacy plants. All of the six transcripts were differentially expressed in the LPB (compared with the CB) of the 'Legacy' plants; the result is similar to that of the comparison of the CB with those of LPB of the Mulegacy plants (Supplementary Table 4).

DE hormone genes are involved in dormancy release in fully chilled flower buds of Mu-legacy plants
In the 'Legacy' plants, 703 and 804 DETs of the genes of four hormones (i.e., ABA, ethylene, auxin, and GA) and ARR genes, were identified in CB (compared with NB) and LPB (compared with CB), respectively 22 . In the Mulegacy plants, 565 and 276 DETs of nine groups of hormone genes were identified in CB (compared with NB) and LPB (compared with CB), respectively (Supplementary Table 2).

Discussion
Cytokinins and gibberellins promote flower bud formation in the Mu-legacy plants Leaves are the primary sources of florigen signals in inducing flowering 57 . For blueberries, overexpression of VcFT (VcFT-OX) is able to stimulate continuous flower bud formation of transgenic "Aurora" through both the DE flowering pathway genes and the DE hormone genes (compared with nontransgenic leaves) 23,24,58 , which was associated with changes in the contents of different hormones (unpublished data). Apparently, expression of hormone genes and hormones abundances are responsive to the VcFT-OX, although it is impossible to conclude whether or not these hormones are part of the florigenic signals because of the well-documented leaves-to-buds transport of FT proteins in many herbaceous plants [59][60][61][62][63][64][65][66][67] . In another case, SOC1 is a major integrator downstream of FT 10 . Overexpression of the K-domain of the VcSOC1 (VcSOC1K-OX) in "Aurora" leaves enhanced flower bud formation through a group of DE hormone genes and DE flowering pathway genes (compared with nontransgenic leaves) where the VcFT did not show differential expression 26 . Thus, the promoted flower bud formation is not likely, at least at transcript levels, entirely through the enhanced VcFT production and transport. In fact, phenotypic changes in both the VcFT-OX and the VcSOC1K-OX plants suggest that hormones could have played essential roles in flower bud formation and dormancy release 23,26 .
The Mu-legacy plants are phenotypically similar to both VcFT-OX and VcSOC1K-OX plants and showed promoted flower bud formation and flowering under nonchilling conditions 25 . However, the major difference is that only two DE flowering pathway genes (compared with nontransgenic 'Legacy' leaves), in contrast to the many more identified in the VcFT-OX and the VcSOC1K-OX plants 24,26 , were detected. More interestingly, the two DE flowering pathway genes, including the upregulated VcRAV1 and the downregulated VcFUL (Table 2, Fig. 3), have repressive effects on flowering bud formation and flowering. This provides evidence that the signals promoting flower bud formation in the Mu-legacy leaves are independent of expression of VcFT and flowering pathway genes. Therefore, we believe that the seven DEGs related to biosynthetic pathways of GA, SA, JA and cytokinin, especially the contents of cytokinins and GAs, were most likely responsible for the promoted flower bud formation in the Mu-legacy plants ( Table 2, Figs. 2, 3).

Responses of blueberry flowering pathway genes during dormancy release
Sufficient chilling hours and warm temperatures are two prerequisites to break blueberry dormancy for healthy flowering and growth 4 . Accordingly, blueberry flower buds often experience two major changes of their physiological statuses, including transition from nonchilled to chilled buds following full chilling and then from chilled buds to flowers upon exposure to warmer temperatures. In the previous transcriptome comparisons for 'Legacy' plants, the comparison of chilled and NB resulted in differential expression of 89% of the blueberry flowering pathway genes; and, the comparison of CB and late-pink flower revealed differential expression of 96% of flowering pathway genes 22 . Apparently, flowering pathway genes are highly involved in dormancy release of blueberry buds at transcript levels. Interestingly, only four (3.8%) flower pathway genes, including VcFD, VcTFL1, VcARP6, and VcDOF5.3, showed differential expression in the nonchilled buds of Mu-legacy, of which~50% were able to flower under nonchilling conditions 25 . Although this does not look reasonable for the Mu-legacy plants at transcript levels, it is possible at post-transcriptional or translational levels because little is known about what physiological signals are actually induced by these DE flowering pathway genes identified in dormancy release. One fact is that many DE hormone genes were also induced in 'Legacy' flower buds during dormancy release 22 . It remains unclear whether the DE flowering pathway genes or DE hormone genes determine the dormancy release. Thus, the value of the Mu-legacy flowering under nonchilling conditions associated with only 4 DE flowering pathway genes is that it demonstrates that other pathway genes, e.g., the 29 DE transcripts of 17 hormone genes (Table 3) 25 . In addition, for the plants of the same age, the Mu-legacy plants had more flower buds than the nontransgenic 'Legacy' plants 25 . Only two DEGs of the flowering pathway genes, including zero in the previous comparison 25 and two in this study, were detected in the two transcriptome comparisons between the leaf tissue of the Mu-legacy plants and that of the nontransgenic 'Legacy' plants. It was likely that the DE flowering pathway genes played a minor role in the promoted Mu-legacy flowering. Accordingly, if the florigenic signals for the promoted flower bud formation in the Mulegacy plants were from leaves, the DETs of hormone genes were more likely the major cause than the DETs of the flowering pathway genes (Fig. 4a). Based on the transcriptome data and the hormone profile in the leaves of the Mu-legacy plants, we believe that it was the expression of VcRR2 that enhanced plant response to cytokinins, which consequently promoted flower bud formation through the changes of either the content of hormones (e.g., GA4) or the expression of the hormone genes (e.g., VcCPS) ( Table 2, Fig. 4a).
In the NB of the Mu-legacy plants, only four DEGs of the flowering pathway genes, including VcFD, VcTFL1, VcARP6, and VcDOF5.3, were detected and were all repressed (compared with NB of the nontransgenic 'Legacy' plants) 25 . In addition, 29 DETs of 17 hormone genes of ABA, auxin, cytokinin, JA, and BR were found (Table 3). These 4 DEGs flowering pathway genes 17 DEGs of hormone genes were associated with mainly by the expression of VcRR2 and resulted in unusual flowering of the NB (Fig. 4b).
From nonchilled to CB and to LPB of both the 'Legacy' and the Mu-legacy plants, a great number of flowering pathway genes and hormone genes were involved (Supplementary Tables 2 and 3). For fully chilled plants, the interactions of these genes under warm conditions drove normal bud break and flowering (Fig. 4c).

Plant materials
Six 5-year-old plants for each of nontransgenic highbush blueberry 'Legacy' and transgenic Mu-legacy and 12 plants of three Legacy-VcDDF1-OX transgenic lines were grown in the courtyard between two greenhouses for phenotyping in 2017. The plants were moved into the greenhouse (heated for winter) in the October of 2017. In January of 2018, young leaf tissues, 5-10 g for each of the 'Legacy' and Mu-legacy plants, were harvested from multiple new shoots; half of these leaves were subjected to freeze drying immediately and another half were ground in liquid nitrogen and stored in a −80°C freezer. Six Mulegacy-T 1 plants, which are from one self-pollinated T 1 seed of Mu-legacy 25 , were grown for phenotyping in the greenhouse since 2016 and never exposed to chilling temperatures. All plants were developed from in vitro cultured shoots and grown under natural light conditions and a regular schedule of irrigation and fertilization using 0.2 g/L fertilizer (Nitrogen:Phosphorus:Potassium = 21:7:7) 70 .

RNA-seq and differential expression analysis
Total RNA of the same set of samples used in hormone analysis was each isolated from~0.5 g of ground tissues using a CTAB method 72 , followed by using the RNeasy Mini Kit for on-column DNase digestion and RNA purification (Qiagen, Valencia, CA, USA). The integrity of the RNA samples was assessed using the Agilent RNA 6000 Pico Kit (Agilent Technologies, Inc., Germany). All RNA samples submitted for RNA sequencing had an RNA quality score >8.0. Sequencing (150-bp pair-end reads) was conducted using the Illumina HiSeq4000 platform at the Research Technology Support Facility at Michigan State University (East Lansing, MI, USA). In total, 30-60 million reads were generated for each biological replicate. The FastQC program (www.bioinformatics.babraham.ac.uk/projects/fastqc/) was used to assess the quality of sequencing reads, with the per base quality scores ranging from 30 to 40.
The RNA-seq reads were analyzed using Trinity 73 and aligned to the transcriptome reference Reftrinity  24 , and the abundance for each of a single read was estimated using the Trinity command "align_and_estimate_abundance.pl". The Trinity command "run_DE_analysis.pl-method edgeR" was used to conduct a differential expression analysis 73 . The differentially expressed (DE) genes or isoforms with ta false discovery rate (FDR) value below 0.05 (p-value < 0.001) were used for further pathway analysis. Fragments per kilobase of transcript per million mapped reads (FPKM) were used to evaluate expression abundance. Most of the analyses were performed using the resources at the High Performance Computing Center of Michigan State University.
We focused our analyses on the genes of three pathways: hormone, sugar, and flowering. Synthesis pathway genes of nine phytohormones in Arabidopsis, including auxin, cytokinin, ABA, ethylene, gibberellin, brassinosteroid, jasmonic acid, salicylic acid, and strigolactones, were retrieved from the RIKEN Plant Hormone Research Network (http://hormones.psc.riken.jp/). Similarly, sugar synthesis pathway genes in Arabidopsis were used as references in analyzing sugar-related genes in this study (Supplementary Table 5). All of the selected Arabidopsis genes were used as queries to blast against the transcriptome reference Reftrinity 24 and the blueberry isoforms showing e-values less than e −20 were used for further various transcriptome comparisons. On the other hand, blueberry flowering pathway genes identified in our previous study 24 were referenced for analyzing the DETs related to flowering pathways in this study.
Some additional transcriptome comparisons were carried out by using data from our previous reports 22,25,74 , as listed in Supplementary Table 6. The previous transcriptome data involved in the comparisons had three biological replicates for each sample and an RNA quality score above 8.0. The data were collected by sequencing 100-bp pair-end reads using the Illumina HiSeq2500 platform at the Research Technology Support Facility at Michigan State University (East Lansing, MI, USA). All the raw sequencing data had been deposited in GenBank.