Transcriptomic changes reveal gene networks responding to the overexpression of a blueberry DWARF AND DELAYED FLOWERING 1 gene in transgenic blueberry plants

Background Constitutive expression of the CBF/DREB1 for increasing freezing tolerance in woody plants is often associated with other phenotypic changes including dwarf plant and delayed flowering. These phenotypic changes have been observed when Arabidopsis DWARF AND DELAYED FLOWERING 1 (DDF1) was overexpressed in A. thaliana plants. To date, the DDF1 orthologues have not been studied in woody plants. The aim of this study is to investigate transcriptomic responses to the overexpression of blueberry (Vaccinium corymbosum) DDF1 (herein, VcDDF1-OX). Results The VcDDF1-OX resulted in enhanced freezing tolerance in tetraploid blueberry plants and did not result in significant changes in plant size, chilling requirement, and flowering time. Comparative transcriptome analysis of transgenic ‘Legacy-VcDDF1-OX’ plants containing an overexpressed VcDDF1 with non-transgenic highbush blueberry ‘Legacy’ plants revealed the VcDDF1-OX derived differentially expressed (DE) genes and transcripts in the pathways of cold-response, plant flowering, DELLA proteins, and plant phytohormones. The increase in freezing tolerance was associated to the expression of cold-regulated genes (CORs) and the ethylene pathway genes. The unchanged plant size, dormancy and flowering were due to the minimal effect of the VcDDF1-OX on the expression of DELLA proteins, flowering pathway genes, and the other phytohormone genes related to plant growth and development. The DE genes in auxin and cytokinin pathways suggest that the VcDDF1-OX has also altered plant tolerance to drought and high salinity. Conclusion A DDF1 orthologue in blueberry functioned differently from the DDF1 reported in Arabidopsis. The overexpression of VcDDF1 or its orthologues is a new approach to increase freezing tolerance of deciduous woody plant species with no obvious effect on plant size and plant flowering time. Electronic supplementary material The online version of this article (doi:10.1186/s12870-017-1053-z) contains supplementary material, which is available to authorized users.

Overexpression of CBF2 in A. thaliana is mainly associated with delayed leaf senescence and extended plant longevity; additionally, overexpression of Muscadinia rotundifolia CBF2 gene in Muscadinia rotundifolia resulted in growth retardation, dwarfism, late flowering, and abiotic stress tolerance [40,49]. In this study, we showed a blueberry-derived CBF (BB-CBF), which was initially considered to be an orthologue of CBF2 that promoted freezing tolerance in A. thaliana [50], was more similar to A. thaliana DDF1. Overexpression of the BB-CBF (herein renamed as VcDDF1) enhanced cold tolerance in leaves and dormant buds but not in flower tissues of a southern highbush blueberry cultivar [51]. Regardless of whether BB-CBF is a DDF1 or CBF2 orthologue, further studies are needed to facilitate a better understanding of the CBF/DREB1-mediate gene networks in blueberry. Unlike A. thaliana, few studies have been conducted to investigate the overall impact of the overexpression of a CBF/DREB1 pathway gene on transcriptomic changes in woody plant species.
Comparative transcriptome analysis is a powerful tool used to identify differential gene expression caused by overexpression of a transgene [52]. For example, overexpression of blueberry FLOWERING LOCUS T (VcFT) in blueberry plants resulted in plant dwarfing and early flowering [53]. Transcriptome analysis of these transgenic plants revealed differentially expressed (DE) genes in flowering and phytohormone pathway genes that are involved in the phenotypic changes driven by VcFT-overexpression [54,55]. The aim of this study is to elucidate transcriptomic responses to the overexpression of VcDDF1 (herein, VcDDF1-OX) and predict overall performance of VcDDF1-OX transgenic blueberry plants. The analysis of DE genes focused on the pathways related to plant growth, flowering or freezing tolerance in blueberry such as plant flowering, CBF-mediated cold/ freezing tolerance, phytohormones and DELLA proteins [51,54,55].

VcDDF1 and VcDDF1-OX in blueberry
The VcDDF1 was initially designated as BB-CBF (Gen-Bank: FJ222601.1) due to its similarity to A. thaliana CBF2, and this reasoning is valid when DDF1 is not included in phylogenetic analysis [50]. However, in our recent transcriptome analysis of highbush blueberry using Trinity and Trinotate [56], the BB-CBF was annotated as DRE1E_ARATH (DDF1). Our designation of BB-CBF as VcDDF1 is the result of the phylogenetic analysis of A. thaliana CBF/DREB1 (i.e., CBF1, CBF2, CBF3, DDF1, and DDF2) and the blueberry-derived DRE1E_ARATH, DRE1A_ARATH, DRE1B_ARATH, and DRE1F_ORYSJ, which showed that BB-CBF is 52.5% similar to DDF1 compared 45.9% to CBF2 (Fig. 1a). The CBF2 orthologues in blueberry were then assigned to the other two gene contigs (c88132_g2 and c85919_g2 in Fig. 1a). It is interesting to note that VcDDF1 orthologues in many other woody plants are often annotated as DREB1 due to the conserved ERF/AP2 DNA-binding domains (Fig.  1b).
To investigate the effect of VcDDF1-OX at transcript levels in non-acclimated floral buds, comparative transcriptome analysis was conducted in non-transgenic 'Legacy' plants and plants of a representative transgenic event 'Legacy-VcDDF1-OX' with a single copy of transgenes (named as II7 in our previous report [51]). The 'Legacy-VcDDF1-OX' showed a 145-fold increase in the expression of the VcDDF1 in comparison to the nontransgenic 'Legacy' plants. The high VcDDF1 expression supported our previous observation that the 'Legacy-VcDDF1-OX' transgenic event showed high freezing tolerance in electrolyte leakage assays [51].
Effect of the VcDDF1-OX on plant freezing tolerance Constitutive expression of VcDDF1 resulted in increased freezing tolerance in detached tissues of A. thaliana and blueberry plants [50,51]. In this study, the VcDDF1-OX enhanced freezing tolerance in intact plants. The freezing tolerance in (45) four-year plants, one of non-transgenic 'Legacy' , two of 'Legacy-pCAMBIA' events, and 41 of 'Legacy-VcDDF1' transgenic events, was investigated. The 'Legacy-VcDDF1' transgenic plants showed a significantly higher survival rate (p = 0.000126) than those of the nontransgenic 'Legacy' and transgenic 'Legacy-pCAMBIA' controls ( Fig. 2a).
In the winter of 2015, we also investigated the freezing tolerance of 12 three-year plants for both the nontransgenic 'Legacy' and the 'Legacy-VcDDF1-OX' transgenic event. The 'Legacy-VcDDF1-OX' transgenic plants exhibited a higher plant survival rate (83.3%) than the non-transgenic plants (41.7%) (Fig. 2b). Applying a freezing shock of −12°C for 15 min resulted in visual differences between transgenic 'Legacy-VcDDF1-OX' and non-transgenic 'Legacy' plants during the plant recovery process. For non-transgenic plant, all leaves and over 90% of the buds showed dying symptoms and died in three weeks. In contrast, for transgenic plants, about 70% of the leaves had no survival leaf tissues, in three weeks; additionally, about 25% buds died. Overall, VcDDF1-OX enhanced freezing tolerance of the intact 'Legacy-VcDDF1-OX' plants (named as II7 in our previous report).
Effect of the VcDDF1-OX on plant growth and flowering The VcDDF1-OX did not alter the growth of transgenic blueberry plants. When four-year transgenic plants of 11 independent 'Legacy-VcDDF1' events, including the representative transgenic event 'Legacy-VcDDF1-OX' , were compared with those of non-transgenic 'Legacy' and transgenic control 'Legacy-pCAMBIA' plants, all plants looked similar in plant stature and appearance (Fig. 3a) and did not show any difference (P < 0.05) in plant height, the number of canes, or the number of flower buds (Fig. 3b, c). These results suggest that VcDDF1-OX has little phenotypic effect on blueberry plant growth and floral bud formation. Therefore, VcDDF1 do not share the designated role of DDF1 in inducing growth retardation and dwarfism.
In relation to both non-transgenic 'Legacy' plants and transgenic control 'Legacy-pCAMBIA' , delayed flowering was not found in any 'Legacy-VcDDF1' plants. For example, 'Legacy-VcDDF1-OX' and non-transgenic 'Legacy' plants did not show significant differences in the number of floral buds, the age of plant flowering, and the yearly flowering time. Moreover, VcDDF1-OX did not affect the chilling requirement of 'Legacy-VcDDF1-OX' plants ( Fig. 3d). Taken together, VcDDF1-OX is not associated with significant changes in plant growth and flowering of tetraploid blueberry plants unlike overexpression of DDF1 in A. thaliana [45,46].

Profile of differentially expressed (DE) genes induced by the VcDDF1-OX
To reveal the potential roles of the VcDDF1 at gene transcription levels, comparative transcriptome analysis was conducted between the 'Legacy-VcDDF1-OX' and nontransgenic 'Legacy' plants. The VcDDF1-OX in nonacclimated floral buds of the 'Legacy-VcDDF1-OX' plants resulted in 2463 DE genes and 3644 DE transcripts, of which 1668 DE genes were annotated. These DE genes were classified in 54 over-represented Gene Ontology (GO) terms (P < 0.05) in the analysis using the GOSlim_-Plant as the selected GO file and A. thaliana annotation as a reference (Fig. 4). Of the 27 over-represented GO terms in biological_process, two highly over-represented GO terms (i.e., GO:0006950-response to stress and GO:0009628-response to abiotic stimulus) revealed the potential function of the VcDDF1-OX in affecting plant freezing tolerance as well as other abiotic stresses. Additionally, two other highly over-represented GO terms (i.e., GO:0007275-multicellular organismal development and GO:0009791-post-embryonic development) suggested that VcDDF1-OX could affect plant growth and flower development (Fig. 4a). The VcDDF1-OX functions at cellular levels (GO:0005263: cell) and intracellular (GO:0005622: intracellular) levels through its catalytic activity (GO:0003824), DNA binding (GO:0003677), transcription factor activity (GO:0003700), and transcription regulator activity (GO:0030528) (Fig. 4b, c). Overall, the results suggest that VcDDF1 is a functional DREB1 transcription factor and the VcDDF1-OX has an impact on gene expression of multiple pathways in blueberry.
In the dormant buds of 'Legacy-VcDDF1-OX' plants, 11,162 DE transcripts of 725 VcCOR genes showed high similarities (e < −20) to 1085 COR genes of A. thaliana (Additional file 1: Table S1). Of these DE VcCOR genes, the up-regulated VcDDF1 and down-regulated VcCBF2  are the DE CBF/DREB1 genes of blueberry. These results suggest VcDDF1-OX regulated some VcCOR genes, which contributed to increase freezing tolerance in intact plants of the 'Legacy-VcDDF1-OX' (Fig. 2).

The responses of blueberry floral genes to the VcDDF1-OX
Whereas overexpression of DDF1 resulted in delaying A. thaliana plant flowering [45,46], VcDDF1-OX did not result in visible changes in flowering of tetraploid blueberry plants (Fig. 3c, d). To investigate the potential impact of VcDDF1-OX on blueberry flowering, we searched for DE floral genes in the dormant buds of 'Legacy-VcDDF1-OX' using the floral gene list of blueberry [54]. Twenty-one floral genes derived from 44 transcripts of 32 gene contigs showed differential expression, of which seven floral genes were up-regulated and 14 were down-regulated (Table 1). This suggests VcDDF1-OX affects flowering pathway genes. However, none of the 21 DE floral genes showed changes above four folds. The expression of blueberry SUPPRESSOR of OVEREXPRESSION OF CONSTANS 1 (VcSOC1) was reduced to 67.6% and the expression of blueberry CONSTANS-LIKE 5 (COL5)-like gene was repressed to 70.5% as much as non-transgenic 'Legacy' plants. This down-regulated expression of VcSOC1 and VcCOL5 is theoretically associated with delayed flowering [57]. However, the blueberry SHORT VEGETATIVE PHASE (SVP)(VcSVP) showed a decreased expression, which in contrast theoretically promotes plant flowering. In spite of these DE floral genes, the VcDDF1-OX was insufficient to promote significant changes in floral bud formation, chilling requirement and flowering time of the 'Legacy-VcDDF1-OX' plants ( Fig. 3c, d).

The responses of major phytohormone genes to the VcDDF1-OX
The VcDDF1-OX did not lead to dwarf 'Legacy-VcDDF1-OX' plants ( Fig. 3a, b), which is inconsistent with the designated function of DDF1 overexpression in causing dwarf A. thaliana plants [45,46]. To evaluate the potential effect of the VcDDF1-OX on plant growth, we identified the DE pathway genes of five major phytohormones [i.e., ABA, GA, auxin (IAA), cytokinin, and ethylene] in dormant buds of 'Legacy-VcDDF1-OX' plants. Except for the ABA pathway, DE transcript contigs were found for all other pathways, including 11 for IAA, 23 for GA, five for cytokinin, and 49 for ethylene (Additional file 2: Table S2). The GA pathway has six and eight transcript contigs shared with those in the IAA and ethylene pathways, respectively, indicating the interaction of these pathway genes (Fig. 5a). In the ethylene pathway, 85 out of 86 DE transcripts showed a less than 4-fold change; and only one DE transcript contig of an orthologue of ETHYLENE-INSENSI-TIVE5 (EIN5) was up-regulated to approximately ten fold. These DE phytohormone genes did not alter plant growth of the 'Legacy-VcDDF1-OX' plants ( Fig. 3a, b).

The responses of DELLA proteins to the VcDDF1-OX
We found 79 transcript contigs of 47 gene contigs in the blueberry transcriptome reference Reftrinity that show high similarities (e < −20) to the five DELLA protein genes of A. thaliana. Of the 79 transcript contigs, two DE transcript contigs of two genes are the RGL3 orthologues in the bud tissues of the 'Legacy-VcDDF1-OX' plants. One of them was repressed to 72.9% and another one was up-regulated to 143.8% (up-regulated by 43.8%). VcDDF1-OX poses little effect on the expression of DELLA protein genes in blueberry plants. This provides additional evidence to show the insignificant effect of VcDDF1-OX on blueberry plant growth and flowering (Fig. 3).

Confirmation of the expression of the selected DE transcripts
We designed six pairs of qRT-PCR primers, consisting of two pairs for GA and IAA pathways and one pair for    Table S3); only one DE transcript revealed by RNA-seq did not show significant difference (p < 0.05) in qRT-PCR analysis (transgenic 'Legacy-VcDDF1-OX versus non-transgenic 'Legacy' samples) but it showed an increase in a regular RT-PCR analysis (Additional file 4: Fig. S1). These results suggest that our RNA-seq data analysis appears reliable for identification of DE genes.

Discussion
Plant freezing tolerance depends on many factors, such as natural environment, plant species/genotypes, plant developmental stages, acclimation state, organs, and tissues [58]. For woody fruit crops, global warming poses concerns for its impact on the phenology of plant dormancy and freezing tolerance. To address these concerns, a thorough understanding of the genetics and mechanisms of plant freezing tolerance and dormancy is needed. With highbush blueberries, freezing injuries in winter and early spring are major concerns.

The CBF/DREB1 orthologues in blueberry
In woody fruit crops, constitutive expression of CBF1 and CBF4 or their orthologues has resulted in similar phenotypic changes to those observed in A. thaliana, suggesting that CBF/DREB1 mediated-freezing tolerance is conserved in plants [35,37,59,60]. In this study, our phylogenetic analysis of blueberry CBF/DREB1 proteins suggest the previous BB-CBF (an orthologue of CBF2) is more likely to be a DDF1 orthologue (VcDDF1) (Fig. 1). The orthologues of this VcDDF1 are in many deciduous woody plants but none of them was annotated as a CBF2 orthologue in GenBank (Fig. 1b). It is also interesting that we do not see CBF1 orthologues in our transcriptome reference. The low coverage of our transcriptome reference may have contributed to the lack of CBF1 orthologues but the lack of orthologues is probably due to the genome specificity of blueberry plants.
Regardless of distinction between CBF2 orthologue (BB-CBF) or DDF1 orthologue (VcDDF1), constitutive expression of this gene is anticipated for dwarfism and late flowering of transgenic plants if its designated function is conserved [45,46,49]. However, this is not the phenotypic change observed in transgenic Arabidopsis and blueberry plants, where the VcDDF1-OX enhanced plant freezing tolerance ( Fig. 2; Fig. 3) [50,51]. These results suggest the function of DREB1 orthologues in different plant species may vary from their functions designated in Arabidopsis (Fig. 6).

Effect of the VcDDF1-OX on plant flowering
In A. thaliana, the delayed flowering caused by overexpression of CBF1,2,3 was due to the increased expression of FLOWERING LOCUS C (FLC) and repressed expression of SOC1 [61]. In this study, VcDDF1-OX repressed the expression of VcSOC1, which is similar to the previous report [61]. However, none of the orthologues of FLC in blueberry, including the MADS-AF-FECTING FLOWERING 2-like gene (VcMAF2), MADS-AFFECTING FLOWERING 5-like gene (VcMAF5), and VERNALIZATION1-like gene (VcVRN1), showed differential expression. In addition, the expression of VcSVP, an orthologue of A. thaliana SVP that acts as a negative regulator of A. thaliana plant flowering [62], was repressed. These differences between the response of the flowering pathway genes to the VcDDF1-OX in blueberry and those of the overexpression of CBF1,2,3 in A. thaliana are responsible for the unchanged flowering phenotype caused by VcDDF1-OX.

Effect of the VcDDF1-OX on DELLA protein genes
In A. thaliana the delayed flowering and growth retardation caused by the overexpression of CBF1,2,3 are due to the changes of DELLA proteins [39,40]. The dwarf A. thaliana plants caused by overexpression of DDF1 was because of reducing bioactive gibberellin [46,48]. In this study, VcDDF1-OX induced little change to the expression of DELLA protein genes, providing further molecular evidence to support that normal growth and flowering of 'Legacy-VcDDF1-OX' plants,

Effect of the VcDDF1-OX on phytohormone genes
In A. thaliana, the cold-response of CBF1, CBF2, and CBF3 is ABA independent while the response of CBF4 is ABA dependent [41]. Additionally, ethylene signaling can affect expression of CBFs [63,64], and DELLA proteins responding to the overexpression of CBF/DREB1 genes are GA-related [46,48]. It seems that CBF/DREB1 genes interact with phytohormone genes to affect plant growth and development. In this study, VcDDF1-OX in tetraploid blueberry plants affected gene expression of the synthesis pathways of IAA, cytokinin, GA, and ethylene but not ABA ( Fig. 5a; Additional file 2: Table S2). The 49 DE transcript contigs of 36 gene contigs in the ethylene pathway could contribute to the increased freezing tolerance [63,64]. The 23 DE transcript contigs of 16 gene in the GA pathway did not show any changes over 4-folds and may have affected the minor changes in the two DE genes of DELLA proteins. Eleven DE transcript contigs of nine genes are orthologues of two IAA pathway genes of A. thaliana, including CYP83B1 and AUXIN TRANSPORTER PROTEIN 1(AUX1) (c90563_g2_i1); both genes are key regulators of root growth and development [65,66]. The orthologues of two cytokinin pathway genes A. thaliana B-TYPE RESPONSE REGULATOR18 (ARR18) and A. thaliana PSEUDO-RESPONSE REGULATOR 2 (APRR2) include five DE transcript contigs, the upregulated ARR18 orthologue could promote root elongation [67]. The DE transcripts involved in auxin and cytokinin pathways have likely altered plant tolerance to drought or high salinity pending on further investigations [45,46,48]. The analysis of DE transcripts of phytohormone genes in addition to DELLA protein genes provide molecular evidence to support that VcDDF1-OX was not associated with dwarf and delayed flowering in tetraploid blueberry plants.

Effect of the VcDDF1-OX on freezing tolerance
In terms of its role in freezing tolerance, VcDDF1 has the same function as DDF1 and other CBF/DREB1 genes in A. thaliana [48]. Based on both our previous electrolyte leakage assay of in vitro tissues [51] and freezing tolerance assay of intact plants in this study, we have demonstrated that VcDDF1-OX is able to enhance freezing tolerance in blueberry plants. In addition, the comparison of VcCOR genes in 'Legacy-VcDDF1-OX' plants with non-transgenic 'Legacy' has provided molecular evidence to support the role of overexpressed VcDDF1 in enhanced freezing tolerance (Fig. 2). Of the DE orthologues of CBF/DREB1 genes, VcDDF1-OX down-regulated VcCBF2, which did not alter plant growth and flowering.

Conclusion
In tetraploid blueberry plants, VcDDF1-OX resulted in enhanced freezing tolerance and normal plant growth and flowering compared to non-transgenic plants. The increased freezing tolerance is attributed to DE VcCOR genes, which are similar to the DDF1 and the other CBF/DREB1 genes (Fig. 6). In contrast to dwarf plant and delayed flowering associated with overexpression of DDF1 or other CBF/DREB1 [45,46,48], normal phenotypes with regards to plant growth and flowering was due to minimal effect of overexpressed VcDDF1 on the expression of DELLA proteins, flowering pathway genes, and other phytohormone genes related to plant growth (Fig. 6). The DE genes in phytohormone pathways of auxin and cytokinin imply that VcDDF1-OX might enhance plant tolerance to drought and high salinity.
This is the first known investigation of a DDF1 orthologue in any crop. More importantly, this is the first time the overexpression of a CBF/DREB1 orthologue was found to enhance plant freezing tolerance without altering plant growth and flowering time. This finding opens a new approach to increase freezing tolerance of deciduous woody plants by using overexpression of VcDDF1 or its orthologues.

Plant materials
A southern highbush blueberry cv. Legacy is tetraploid and needs over 800 chilling units (CU) for normal flowering. The 'Legacy' plants used in this study was original derived from the blueberry cultivar collections growing in a research field of the Horticulture Teaching and Research Center of Michigan State University. Transgenic 'Legacy' plants (herein 'Legacy-VcDDF1') contain a blueberry derived CBF gene (AVI45245.1), which was designated as BB-CBF [50,51] and renamed as VcDDF1 in this report. The 'Legacy-pCAMBIA' is a transgenic control for the VcDDF1. 'Legacy-VcDDF1-OX' is a representative transgenic 'Legacy-VcDDF1' that was used for RNA sequencing. The 'Legacy-VcDDF1-OX' (named as II7) contains a single copy of transgenes and showed high freezing tolerance [51]. Production of the 'Legacy- Fig. 6 Schematic diagram illustrating the effect of overexpressing a blueberry DDF1 orthologue on plant freezing tolerance. The diagram illustrating the effect of overexpression of CBF/DREB1 on freezing tolerance, plant growth and development in A. thaliana was derived from previous reports [40,41,46]. The overexpression of VcDDF1 showed little effect on DELLA protein genes but could affect plant tolerance to drought and high salinity through altered gene expression in auxin and cytokinin pathways VcDDF1' and Legacy-pCAMBIA' containing the binary vector pCAMBIA2301 was described in our previous report [51].
All non-transgenic and transgenic plants were obtained through micropropagation of in vitro cultured shoots. Plant age was determined based on the time after the shoot was rooted in soil. Rooting of in vitro cultured shoots and plant growth in the greenhouse were performed according the protocols established by Song [68]. All plants were grown normally and were fully vernalized unless otherwise mentioned. For full vernalization in winter, plants were potted and grown in a non-heated hoop house or in a secured courtyard under natural light conditions at Michigan State University, East Lansing, Michigan (latitude 42.701847, longitude −84.482170). The average low and high temperatures in January are −10.6°C and −1.8°C, respectively (http://www.usclimatedata.com/climate/east-lansing/michigan/united-states/usmi0248). RNA preparation, sequencing, and de novo transcriptome assembly Floral buds were collected in November 2013 before the plants were exposed to a non-heated greenhouse for chilling treatments. All tissues collected were frozen immediately in liquid nitrogen and stored at −80°C.

Plant growth and flowering
Total RNA isolation, RNA sequencing using the Illumina HiSeq2500 platform, de novo transcriptome assembly using the Trinity platform (trinity/20140413p1) [56] were described in our recent report [54].

Differential expression analysis and transcriptome annotation
RNA-seq reads of three biological replicates for each of 'Legacy' and 'Legacy-VcDDF1-OX' were analyzed. Two technical replicates were sequenced for each biological replicate and were combined together for analysis. The paired reads, two sets for each biological replicate, were aligned to the transcriptome reference developed for 'Legacy' [54] and the abundance of each read was estimated using the Trinity command "align_and_estimate_abundance.pl". The Trinity command "run_DE_analysis.pl -method edgeR" was used for differential expression analysis. The differentially expressed (DE) (relative to non-transgenic 'Legacy' unless other mentioned) genes or transcripts with false discovery rate (FDR) values below 0.05 were used for further analyses.

Phylogenetic analysis of VcDDF1
Representative nucleotide sequences of CBF/DREB1 of five A. thaliana CBF/DREB1 genes were retrieved using The A. thaliana Information Resource (TAIR) server (https://www.arabidopsis.org/tools/bulk/index.jsp). Orthologues of A. thaliana CBF/DREB1 genes in blueberry were identified from our annotated transcripts. The selected transcripts were converted to amino acid sequences based on BLAST results retrieved using the NCBI server (http://blast.ncbi.nlm.nih.gov/Blast. cgi). Selected nucleotide sequences from both A. thaliana and blueberry were aligned using Clustal Omega multiple sequence alignment program at EBI with default parameters (http://www.ebi.ac.uk/Tools/msa/clustalo/). Phylogenetic trees were generated using MEGA 6.06 software [71].
The VcDDF1 protein sequence (AVI45245.1) was used to search for VcDDF1 orthologues using the NCBI server. The selected protein sequences were aligned using Clustal Omega multiple sequence alignment program at EBI with default parameters.

Gene network construction
Annotated transcripts were imported to Cytoscape 3.4.0 under BiNGO's default parameters with selected ontology file 'GOSlim_Plants' and selected organism A. thaliana [72,73].

Identification of the selected pathway genes
Representative protein sequences of selected genes of A. thaliana were download from the TAIR server (https:// www.arabidopsis.org/tools/bulk/sequences/index.jsp). The retrieved sequences were used to search the transcriptome reference of blueberry (herein refTrinity) using the tblastn command of BLAST+. The resultant transcripts that show e-value lower than −20 were used to screen the DE transcript list of non-acclimated floral buds.
The 2637 cold-regulated genes (CORs) identified in wild-type A. thaliana plants and 172 CORs differentially expressed at a warm temperature (22°C) in transgenic A. thaliana plants overexpressing CBF1, CBF2 or CBF3 were obtained from Park et al. [42]. These CORs were used to identify their orthologues in blueberry (VcCORs), which was used to analyze the effect of VcDDF1-OX on VcCORs. The blueberry floral genes identified in our previous study [54] were used to analyze flowering pathway genes affected by VcDDF1-OX.
Quantitative RT-PCR (qRT-PCR) of DE transcripts Reliability of DE genes/transcripts identified through RNA-seq was evaluated through qRT-PCR analysis of six selected transcripts (Additional file 3: Table S3). These transcripts are from the representative DE genes in auxin, ethylene, cytokinin, and GA pathyways. They have high fold changes (>2) and sequence specificity (based on alignment result of different isoforms) for PCR amplification. Eukaryotic translation initiation factor 3 subunit H was the internal control (Additional file 3: Table S3).
The RNA samples used for RNA-sequencing, including samples of three biological replicates for each of 'Legacy' and 'Legacy-VcDDF1-OX' , were used for cDNA preparation. Reverse transcription of RNA to cDNA was performed using SuperScript II reverse transcriptase (Invitrogen, Carlsbad, CA, USA). The resulting cDNA of one micro gram of RNA was diluted (volume 1: 4) in water and 1 μl/sample (25 ng) was used for PCR reactions.
The primers were designed using the online tool provided by Integrated DNA Technologies, Inc. (https:// www.idtdna.com/Primerquest/Home/Index), where the primers were synthesized (Additional file 3: Table S3). qRT-PCR was performed in triplicate on an Agilent Technologies Stratagene Mx3005P (Agilent Technologies, Santa Clara, CA) using the SYBR Green system (Life Technologies, Carlsbad, CA). In each 25 μl reaction mixture, 25 ng cDNA, 200 nM primers and 12.5 μl of 2× SYBR Green master mix were included. The reaction conditions for all primer pairs were 95°C for 10 min, 40 cycles of 30 s at 95°C, 60 s at 60°C and 60 s at 72°C, followed by one cycle of 60 s at 95°C, 30 s at 55°C and 30 s at 95°C . The specificity of the amplification reaction for each primer pair was determined by the melting curve. Transcript levels within samples were normalized to the eukaryotic translation initiation factor 3 subunit H. Fold changes were calculated using 2 -ΔΔCt , where ΔΔCt = (Ct GOI -Ct nom ) Legacy-VcDDF1-OX -(Ct GOI -Ct nom ) Legacy for each transgenic 'Legacy-VcDDF1-OX' versus a non-transgenic 'Legacy' sample (n = 3) [79]. In addition, regular RT-PCR was also used for selected transcripts. The reaction conditions using 50 ng cDNA per reaction for all primer pairs were 94°C for 2 min, 35 cycles of 45 s at 94°C, 60 s at 60°C and 60 s at 72°C, with a final 10 min extension at 72°C. RT-PCR products were separated on 1.0% agarose gel containing ethidium bromide, visualized, and photographed under UV light.