An invaluable transgenic blueberry for studying chilling-induced flowering in woody plants

Background Many deciduous woody crops require a minimum level of chilling to break dormancy and allow the seasonal growth of vegetative and floral buds. In this study, we report the discovery of an invaluable transgenic event of the blueberry cultivar ‘Legacy’ (hereafter, Mu-Legacy) for studying chilling-induced flowering in woody plants. Mu-legacy and its progeny provide a unique material to study the unknown mechanism of chilling-mediated flowering in woody plants. Results Unlike nontransgenic ‘Legacy’ and plants of 48 other transgenic events, Mu-Legacy plants were able to flower under nonchilling conditions and had early flower bud formation, reduced plant size, and reduced chilling requirement for normal flowering. These characteristics were heritable and also observed in self-pollinated, transgenic T1 progenies of Mu-Legacy. A 47-Kbp genomic sequence surrounding the transgene insertion position was identified. RNA-sequencing data showed increased expression of a RESPONSE REGULATOR 2-like gene (VcRR2), located adjacent to the insertion position in Mu-Legacy and likely driven by the CaMV 35S promoter of the transgene. The Mu-Legacy showed 209 differentially expressed genes (DEGs) in nonchilled flower buds (compared to nontransgenic ‘Legacy’), of which only four DEGs were in the flowering pathway. This suggests altered expression of these few genes, VcRR2 and four flowering DEGs, is sufficient to significantly change flowering behavior in Mu-Legacy. Conclusions The significance of VcRR2 in Mu-Legacy suggests that the VcRR2-involved cytokinin pathway likely contributes to the major differences in chilling-mediated flowering between woody and herbaceous plants. More importantly, Mu-Legacy shows increased yield potential, a decreased chilling requirement, and better winter hardiness than many low-chilling cultivars growing in southern warm winter conditions. Electronic supplementary material The online version of this article (10.1186/s12870-018-1494-z) contains supplementary material, which is available to authorized users.


Background
Deciduous fruit-bearing crops that grow in temperate climates come mainly from nine families, including Rosaceae (apples, pears, quinces, almonds, apricots, plums, cherries, peaches, raspberries, blackberries, loquats, and strawberries) [1], Fagaceae (chestnuts), Betulaceae (filberts), Juglandaceae (pecans and walnuts), Ebenaceae (persimmons), Moraceae (figs and mulberries), Vitaceae (grapes), Ericaceae (blueberries and cranberries), and Grossulariaceae (currants and gooseberries). These crops often require a certain period of cold exposure (chilling) to stimulate dormancy release and induce floral buds to blossom in their life cycles. Healthy flowering under the changing climate is of critical importance for sustainable production of temperate woody crops. Generally, the chilling requirement and plant hardiness, both of which vary among plant species, are the main consideration in regard to low temperatures when temperate woody crops are adopted.
Vaccinium is a genus of terrestrial shrubs in the family Ericaceae (Syn. Heath) containing approximately 450 species [2]. Highbush blueberries (Vaccinium corymbosum L.), including northern and southern ecotypes, are the major cultivated blueberries. Most of the commercial highbush cultivars require chilling to ensure normal flowering. The northern highbush cultivars require more than 800 chill units (CU) to break dormancy and flower in the spring and have generally better winter hardiness. In contrast, the southern highbush blueberry often needs 150 to 800 CU, but has greater high-temperature tolerance in the summer. Blueberry floral bud initiation often starts before endo-dormancy. Enough chilling accumulation during endo-dormancy is critical for floral bud formation and bud-break in the Spring for decidous woody fruit and nut crops. Insufficient chilling hours prevent bud-break and often lead to reduced fruit production. Manipulation of chilling requirement are considered to be long-term solution secure deciduous fruit production under the changing climate [3].
To date, extensive studies have elucidated the molecular basis of vernalization and flowering in Arabidopsis [4][5][6][7][8][9][10][11][12], the cereals [8,9,13] and beets [14]. In contrast, far less progress has been made in unveiling the molecular pathways of chilling-induced dormancy release in woody perennials, due mainly to the high complexity of their genomes. It has been suggested that overlap may exist between components of the herbaceous vernalization pathways and the dormancy release pathways of woody plants [15]. However, neither the FRIGIDA (FRI) [16] or FLOWER-ING LOCUS C (FLC)-determined [4,5,17] vernalization models in Arabidopsis nor the VERNALIZATION gene (VRNs)-mediated vernalization models in cereals and their relatives have proven sufficient to form a model of cold-dependant flowering in woody plants [8,9,13]. In general, the chilling requirement for dormancy-break is species-and genotype-dependent [14,18].
Functional analysis of a blueberry FLOWERING LOCUS T (FT) gene (VcFT) has been conducted to study flowering mechanisms in highbush blueberry (Vaccinium corymbosum L.) [19][20][21]. Overexpression of VcFT (approx. 2900-fold increase in leaf tissues) caused continuous and precocious flowering in in vitro shoots and in one-year-old ' Aurora' plants [19]. However, two-and three-year-old VcFT-overexpressing plants did not flower normally, and the majority of the flower buds did not open under greenhouse conditions without chilling. It is interesting that overexpression of VcFT was not sufficient to completely release all blueberry floral buds from dormancy, suggesting that the molecular pathways for floral transition and breaking of seasonal dormancy only partially overlap in this species [20,21].
Blueberry floral initiation occours before bud dormancy release. Transcriptome analysis using RNA sequencing data from nonchilled, chilled, and late pink buds of southern highbush blueberry 'Legacy' has revealed genes associated with chilling accumulation and bud break [3]. It is interesting that VcFT expression did not show a differential expression in chilled flower buds (compared to nonchilled flower bud) but were up-regulated in late-pink buds (compared to chilled flower bud. DORMANCY ASSOCIATED MADS-box (DAM) genes which are a cluster of six MADS-box transcription factors. The loss of all or part of the DAMs resulted in the non-vernalizaed peach evergrowing mutant [22,23]. The DAM genes are considered alternatives to FLC in regulating vernalization-mediated chilling requirement and flowering [22,24]. However, the DAM genes show high similarities to A. thaliana AGAMOUS-LIKE 24 (AGL24) and SHORT VEGETATIVE PHASE (SVP) genes [24,25]. Additionally, functional analysis of DAMs to reveal their roles in chilling-mediated flowering through reverse genetics has not been reported in peach. VcSVP showed differential expression in chilled and late-pink buds. In chilled blueberry flower buds, both VcFLC and VcSVP homologues decreased but increased in late-pink buds [3]. The functional orthologues of FLC and AGL24 were not detected in blueberry, suggesting that the vernalization/chilling-mediated flowering pathway of blueberry is different from A. thaliana.
Expression of blueberry DWARF AND DELAYED FLOWERING gene (VcDDF1) is not cold-inducible; however, overexpression of the VcDDF1 increases freezing tolerance in transgenic leaves and flower buds [26][27][28]. One transgenic event of blueberry 'Legacy' (herein Mu-Legacy) transformed with VcDDF1 was identified from 49 transgenic events because it can flower under nonchilling conditions whereas all the nontransgenic and transgenic plants of 48 other events could not. This Mu-Legacy is a desirable material for revealing the mechanism of chilling-induced flowering in woody plants.

Phenotypes of mu-legacy
Sixteen nontransgenic and 202 transgenic 'Legacy' plants were grown in a heated greenhouse in Michigan in 2009. These plants were developed from in vitro shoot cultures and were grown in the greenhouse for approximately 14 months after they were rooted. The transgenic plants, three to five plants per transgenic event, were from 49 transgenic events, each containing an overexpression vector for VcDDF1 [27,28]. All three plants of one VcDDF transgenic event, which was designated Mu-Legacy, started flowering in the heated greenhouse in mid-October (Fig. 1a, b), while all 16 nontransgenic 'Legacy' and 199 transgenic plants of the other 48 VcDDF1 transgenic Legacy events could not flower without chilling. The Mu-Legacy had an obviously shortened juvenility phase, and could flower without the typical dormancy period and subsequent release following chilling. These results suggest that the overexpressed VcDDF1 is probably not the only factor causing the changes in Mu-Legacy plants. As of this report, 200 plants for each of the Mu-Legacy and nontransgenic  (Fig. 1d). The difference in flowering between long-day and short-day grown nonchilled Mu-Legacy plants is that the long-day grown plants showed flowered on newly formed soft-wood shoots, which was not observed in plants grown under short day (Fig. 1g). These results suggest that compared to nontransgenic 'Legacy' the Mu-Legacy plants have a mutation that contributes to altered chilling requirment (Fig. 1). However, the Mu-Legacy still requires some accumulation of chilling for normal flowering.
When looking at one year-old plants with over 1000-CU of cold exposure, none of the nontransgenic 'Legacy' plants flowered, wheras 84.9% of Mu-Legacy plants flowered. This indicates early floral bud formation in the Mu-Legacy plants with a reduced juvenility phase. Chilling accumulation influences the bloom period of both Mu-Legacy and nontransgenic 'Legacy' plants. For Mu-Legacy, lower levels of chilling accumulation (133 or 300 CU) were associated with longer bloom periods when compared with greater chilling accumulation (500 or 850 CU) (Fig. 2). Flowering of nontransgenic 'Legacy' plants was not observed until chilling accumulation reached 500 CU (Fig. 2c). With 850-CU chilling, Mu-Legacy and nontransgenic 'Legacy' plants showed a similar, very condensed bloom period, but Mu-Legacy flowered two days earlier than nontransgenic 'Legacy' plants ( Fig. 2d).
Mu-Legacy had a smaller plant size, greater flower bud formation, and a reduced chilling requirement when compared to either nontransgenic 'Legacy' (Fig. 1e

Phenotypes of the T 1 progeny of mu-legacy
The chilling-independent flowering trait in Mu-Legacy was heritable and segregated in the T 1 generation. Of 36 self-pollinated T 1 seedlings of Mu-Legacy grown in a heated greenhouse for two years, 22 plants were able to flower prior to their exposure to any chilling while 14 plants did not flower (Additional file 1: Fig. S1). All of these 22 T 1 plants were PCR-positive for both the neomycin phosphotransferase II (nptII) and VcDDF1 transgenes; in contrast, none of the 14 non-flowering plants were PCR-positive. The transgene segration rate is 3:1 (chi square test, p < 0.01). Other phenotype variations, such as plant size, architecture, and leaf shape and size, were also observed in the 36 T 1 plants. One extremely dwarf transgenic seedling (herein Mu-Legacy-T 1 ) was identified (Fig. 4). These results suggest that the tranfer DNA (T-DNA) insertion is responsible for the phenotypic changes in the Mu-Legacy.

Genetic control of the mutation in mu-legacy
An electrolyte leakage assay suggested that overexpressing VcDDF1 resulted in a significant increases in freeze tolerance in leaf tissues of several transgenic events, i.e., the representative transgenic event VcDDF1-Legacy (previous code: II7) but not in Mu-Legacy (previous code: II3) [28]. In compared to nontransgneic 'Legacy' , VcDDF1 overexpression (VcDDF1_OX) in Mu-Legacy resulted in increases in VcDDF1 expression (c32575_g1_i1) at 108fold in nonchilled flower buds (Additional file 2: Table S1) and 40-fold in leaves (Additional file 3: Table S2), respectively. In compared to VcDDF1_Legacy, the VcDDF1_OX in Mu-Legacy did show differential expression in nonchilled flower buds (Additional file 4: Table S3) but showed a significant decrease in leaves (Additional file 5: (See figure on previous page.) Fig. 1 Flowering of Mu-Legacy and nontransgenic 'Legacy' plants. a, b, Flowering of two-year old, nonchilled Mu-Legacy (a) and 'Legacy' (b) plants under a short-day photoperiod (nine hours). Red arrows indicate flowers, fruits, or flower buds. c, d, Flowering of fully chilled two-year old Mu-Legacy (c) and 'Legacy' (d) plants under a long-day photoperiod (14 h). e, f, Pattern of flowering when three-year old Mu-Legacy and 'Legacy' plants (n = 6) were grown in greenhouses under nonchilling conditions with a short-day photoperiod (nine hours) (e) and a long-day photoperiod (16 h) (f). The primary y-axis is for the line chart and the secondary y-axis is for the column chart. The bars are showing the mean number of flower buds per plant and the number of flower buds that bloomed. Each data point is an average of six plants plus standard deviation bars. g, In addition to the flowering pattern observed under a short-day photoperiod (nine hours) (a), flowering was shown on the new shoots of Mu-Legacy under a longday photoperiod (16 h). h, A shorter Mu-legacy plant is showing more fruit production than a nontransgenic 'Legacy' plant. The five-year old plants were grown in a secured courtyard under natural environmental conditions. Significance (compared to nontransgenic 'Legacy') determined using a Student's t-test is denoted. One asterisk (*) indicates p < 0.05 and two asterisks (**) indicate p < 0.01 Table S4). Thus, little evidence supports that the overexpressed VcDDF1 is the only explanation for the mutations in Mu-Legacy.
Southern blot analysis indicated that Mu-Legacy contains one copy of the T-DNA (Fig. 5a). The T-DNA insertion position was initially found in a 1279 base pairs (bp) fragment. The left border of the inserted T-DNA has an 18-bp deletion, and the right border has a 68-bp deletion. The insertion broke a 354-bp reading frame, which shows the highest similarity at protein level to a retrotransposon protein of rice (Oryza sativa). Sequence analysis of RT-PCR products showed the presence of the 354-bp fragment in the leaf samples of both Mu-Legacy and nontransgenic 'Legacy' , suggesting that the 354-bp sequence has at least two identical alleles in nontransgenic 'Legacy'. Using the 1279 bp were grown under normal chilling conditions (> 1000 chilling hours) before they reached three-years old. The primary y-axis is for the column chart and the secondary y-axis is for the line chart. Each data point is a mean of five plants. Significant changes determined using a Student's t-test are denoted. Asterisk (**) indicates p < 0.01 sequence of the T-DNA insertion position of Mu-Legacy to search blueberry ESTs, we found a 781 bp EST contig CV091265.1 that has 199 bp overlap with the 1279 bp sequence. The 199 bp region (herein VcRR2) of this EST shows a high similarity to B-type RESPONSE REGULATOR 2 (ARR2).
Based on this EST sequence, the sequence at the insertion position could be extended to 3053 bp (Additional file 6 Supplemental information 1); however, the sequence information alone is insufficient to explain why the insertion is responsible for the changes in Mu-Legacy. Four DEGs of the flowering pathway in mu-legacy RNA sequencing data were generated for comparative analysis of nontransgenic'Legacy' and Mu-Legacy. The overrepresented gene ontology (GO) term "reproduction"  Table S1). These results suggest that the four DEGs of the flowering pathway in Mu-Legacy play significant roles in chilling-mediated flowering.

DEGs of major phytohormone-related genes and COR genes in mu-legacy
Due to the potential effect of phytohormones in chillingmediated blueberry flowering [3], DE orthologues of five major phytohormones of A. thaliana [i.e., abscisic acid (ABA), gibberellic acid (GA), auxin, cytokinin, and ethylene] were examined in Mu-Legacy tissues (Additional file 2: Tables S1, Additional file 3: Tables S2). GA is involved in the flowering pathway through its interaction with the SUPPRESSOR OF OVEREXPRESSION OF CONSTANS 1 (SOC1) and TFL1 [29,30]. In nonchilled flower buds, orthologues of ENT-KAURENOIC ACID HYDROXYLASE (KAO2) and GA REQUIRING 3 (GA3) in the GA biosynthesis pathway were both upregulated (|Log 2 FC = Log 2 fold changes | < 2) without being associated with a differential expression of blueberry SOC1 (VcSOC1). DE orthologues of the CYP83B1 PROTEIN (SUR2) in the auxin pathway were also upregulated. DE orthologues of two A. thaliana genes in the cytokinin pathway, one in the ethylene pathway, and one in the auxin pathway showed low fold changes (|Log 2 FC| < 2). Compared to the number in nonchilled flower buds (eight in total), fewer DE orthologues of these phytohormone genes were present in leaves (two) (Additional file 3: Table S2). These DE phytohormone genes in different tissues could play a role in the early flowering and reduced size of Mu-Legacy plants.
DE orthologues of A. thaliana cold-regulated genes (CORs) were present in Mu-Legacy plants (Additional file 2: Tables S1, Additional file 3: Tables S2). No significant increase in freezing tolerance was observed in leaf and bud tissues in the electrolyte leakage assay (previous code: II3) [28], suggesting that these DE CORs were insufficient to drive a significant increase in freezing tolerance.

DEGs caused by insertion position
Unlike Mu-Legacy plants, VcDDF1-Legacy plants, flowered similar to nontransgenic 'Legacy' [27], suggesting that the overexpressed VcDDF1 might not be the major factor for the altered flowering in the Mu-Legacy. To further explore the potential candidate genes that could be resposible for the mutation in Mu-Legacy, transcriptomic analyses of Mu-Legacy (compared to VcDDF1-Legacy) were also conducted in leaves and nonchilled flower buds. The comparisons resulted in 1108 and 1119 DE unique genes in leaves and nonchilled flower buds, respectively  Table S3, Additional file 5: Table S4). Further analysis of the four groups of DEGs identified in the comparions of Mu-Legacy vs. nontransgenic 'Legacy' (77 for leaves and 209 for flower buds) and Mu-Legacy vs. VcDDF1-Legacy (1108 for leaves and 1119 for flower buds) resulted in 18 shared DEGs and 31 unannotated transcripts by Trinotate, which were the major DEGs and DE transcripts caused by the transgene insertion position in Mu-Legacy ( Fig. 8; Table 2, Additional file 2: Table S1, Additional file 3: Table S2, Additional file 4: Table S3, Additional file 5: Table S4). The 18 DEGs and 31 DE transcripts were further annotated by searching the Arabidopsis protein database ( Table 2). Of the 18 DEGs, the upregulated c96767_g2_i12 is a VcRR2 gene adjacent to the insertion position in Mu-Legacy (Fig. 5, Tables 1 and 2). The repressed c49456_g1_i1 and c49456_g2_i2 (annotated to ACTS_RAT or actin-11) showed a high similarity to the Arabidopsis ARP6 gene, which regulates FLC independent of vernalization in Arabidopsis [31] Thus, the upregulated VcRR2 and the repressed ARP6 could be the key DEGs causing the changes (compared to both nontransgenic 'Legacy' and VcDDF1-Legacy) in chillingmediated flowering in the Mu-Legacy. It appears from this data that a small numbers of DEGs were responsible for the phenotypic changes in the Mu-Legacy ( Fig. 1-3), suggesting that the Mu-Legacy is an invaluable material to study the blueberry flowering mechanism.

Blueberry B-type RESPONSE REGULATOR (VcRR2) is likely responsible for the mutation of mu-legacy
In the 47-Kbp transgene insertion region, transgenes showed overexpression; more importantly, the VcRR2 was the only DE gene and it was upregulated in both leaves and buds (Table 1; Additional file 2: Table S1, Additional file 3: Table S2, Additional file 4: Table S3, Additional file 5: Table S4). Similar results were    observed in leaves of Mu-Legacy-T 1 (Fig. 9). Since the VcRR2 did not show differential expression in leaves and flower buds of VcDDF1-Legacy [27], the upregulated expression of VcRR2 could be responsible for the mutation of Mu-Legacy. It is likely that the CaMV 35S promoter of the gusA gene drove through the NOS terminator and then promoted expression of the adjacent VcRR2 gene.

T-DNA insertion upregulated expression of adjacent host genes
Regardless of transgene, a T-DNA integration is often associated with a position mutation by a random gene insertion [32]. In Mu-Legacy, the insertion broke a 354-bp retrotransposon (Fig. 5); a homologue of the retrotransposon was also identified. The role of the 354-bp retrotransposon in regulating expression of the VcRR2 is not known in this study. In the literature, there is almost no evidence showing that a retrotransposon itself contributes to regulating plant flowering. Since the retrotransposon is adjacent to VcRR2, it is a genetic marker of the VcRR2; however, whether it has some regulatory roles in VcRR2 expression or evolution is still to be determined.
Random T-DNA insertion is often associated with nucleotide deletion at the T-DNA borders [32], which was demonstrated in this study (Fig. 5). More interestingly, we found the constitutive 35S promoter driving through the NOS polyA terminator where a truncation occurred right after the stop code (Fig. 5), which enabled overexpression of a downstream VcRR2 gene in Mu-Legacy (Table 1, Fig. 9). Consequently, the Mu-Legacy and its transgenic progeny showed early flowering, altered chilling requirement, and reduced plant height (Figs. 1-4).

Roles of VcRR2 in plant flowering
In the dicot Arabidopsis, ARR2 is involved in cytokinin responses that affect a wide range of developmental processes in response to phytohormones [33][34][35][36][37][38]. In comparison to wild type Arabidopsis, loss-of-function mutants of ARR2 (arr2) show retarded growth and development, and early flowering [34,39]; the gain-of-function of ARR2 promotes in vitro shoot production and leaf differentiation, and delays leaf senescence [35]. In the monocots (e.g., rice and corn), two-component systems play a central role in cytokinin signaling, but the roles of ARR2 orthologues are not clear [40][41][42]. B-type RR orthologues are identified in woody plants, for example, in black cottonwood (Populus trichocarpa) and peach (Prunus persica) [43].
Here, for the first time in woody plants (blueberry), it was found that an enhanced expression of VcRR2 in Mu-Legacy plants, in addition to the overexpressed VcDDF1, resulted in reduced plant size, early flower bud initiation and flowering without chilling, and enhanced flower bud formation (Figs. 1-4, Table 2). Although the roles of the overexpressed VcDDF1 can not be excluded in the Mu-Legacy plants, our recent studies have demonstrated that the overexpressed VcDDF1 alone were  Table 1 not able to change chilling requirment in the other 48 VcDDF1-OX lines [26][27][28]. Functional analyses of VcRR2 through overexpression, gene knock-out, and phytohormone profiling are still needed to reveal the roles of VcRR2 in blueberry flowering.
Flowering pathway genes and chilling-mediated flowering of blueberry plants Comparative transcriptome analysis provides a powerful tool to study the differentially expressed genes associated with phenotypic changes in transgenic blueberry plants [20,21,26,27]. In terms of flowering time regulation, overexpression of a blueberry FLOWERING LOCUS T (VcFT) (an increase of approximately 2900-fold) was unable to completely eliminate the need for chilling in blueberry for normal flowering [20]. More recently, overexpression of of the VcSOC1K increases blueberry productivity by promoting early and more flower bud formation through the other DE MADS-box genes [44]. In this study, flowering of nonchilled Mu-Legacy plants has demonstrated non-VcFT-mediated plant flowering, where early flowering and reduced shilling requirement was associated with decreased expression of VcFD, VcTFL1, VcARP6, and VcDOF5.3 without the involvement of other major flowering pathway genes (Additional file 2: Table S1, Additional file 3: Table S2). Further transcriptomic comparsions of these transgenic plants will reveal the DE genes that are involved in chilling-mediated flowering in blueberries.
The transgene insertion position and VcDDF1-OX in Mu-Legacy caused 18 DEGs and 31 unannotated DE transcripts ( Fig. 8; Table 2). Most of the 18 DEGs have not been well-studied in plants. The obvious phenotypic changes associated with a small number of DE genes in the leaf and nonchilled bud tissues of Mu-Legacy make it invaluable for studying the chilling-mediated flowering mechanism in woody plants (Fig. 6, Additional file 1: Figure S1). Further investigations on the effect of chilling on expression of flowering pathway genes in the progenies of Mu-Legacy will allow the unravelling of the mechanism of chilling-mediated flowering in woody plants.

Conclusions
Mu-Legacy was identified from 49 VcDDF1 transgenic events. The most obvious phenotypic change is that Mu-Legacy was able to flower under nonchilling conditions, whereas nontransgenic 'Legacy' and the other transgenic events could not. In addition, transgenic progenies derived from the self-pollinated seeds were also able to flower under nonchilling conditions, but none of the nontransgenic segregants in progenies derived from Mu-legacy, or transgenic progenies from another VcDDF1 transgenic event, flowered. Since the mechanism of chilling-mediated flowering remains unknown in woody plants, Mu-Legacy and its progeny provide a unique material to study woody plant flowering. The significance of VcRR2 in Mu-Legacy suggests that the VcRR2-involved cytokinin pathway likely contributes to the major differences in chilling-mediated flowering between woody and herbaceous plants. More importantly, Mu-Legacy shows increased yield potential, a decreased chilling requirement, and better winter hardiness than many low-chilling cultivars growing in southern warm winter conditions.

Plant materials
A southern highbush blueberry cv. 'Legacy' is tetraploid and needs more than 800 CU for normal flowering. Forty-nine transgenic 'Legacy' plants contain an overexpressed VcDDF1 [28]. Both transgenic and nontransgenic plants were developed from in vitro cultured shoots; they were grown in the greenhouse (heated for winter) or the courtyard between two greenhouses under natural light conditions and a regular schedule of irrigation and fertilization using 0.2 g/L fertilizer (Nitrogen: Phosphorus: Potassium = 21: 7: 7).
Plant chilling was conducted in an unheated hoop house in winter under natural light conditions or in growth chambers at 4°C with a 12-h photoperiod. The conversion of selected temperatures to chill units for highbush blueberry is based on the equation: total chill units = 0.5 × number of hours with temperatures below 2.4°C and 9.2-12.4°C + 1 × number of hours with temperatures 2.5-9.1°C -0.5 × number of hours with temperatures 16-18°C -1 × number of hours with temperatures above 18°C [45].
One hundred transgenic plants of a representative VcDDF1 transgenic event (herein VcDDF1-Legacy) and two hundred plants at different ages, for each of the nontransgenic 'Legacy' and Mu-Legacy (a flowering mutant of transgenic legacy), were grown for phenotyping experiments under various chilling conditions, including chill units of zero, 133, 300, 500, and 850 under controlled conditions or above 850 CU in the open-air conditions in Michigan. For each treatment, three plants for VcDDF1-Legacy and at least five replicated plants for each of the nontransgenic 'Legacy' and Mu-Legacy were used.
Self-pollinated T 1 seeds of Mu-Legacy were stored in a refrigerator for six months prior to their germination either in soil or on half strength Murashige and Skoog (1962) medium (MS) [46]. The seedlings germinated on half strength MS were micropropagated prior to being grown in the greenhouse. All T 1 progenies were grown in the greenhouse and were not exposed to any chilling. Plant size and flowering time were documented.

Mapping of the T-DNA insert in mu-legacy
Identification of the T-DNA insert in Mu-Legacy was conducted using the PCR method [32]. DNA samples of nontransgenic 'Legacy' and Mu-Legacy were obtained using a cetyltrimethylammonium bromide (CTAB) method. Both HindIII-digested and EcoRI-digested DNA samples were ligated to adapters and then used for PCR using adapter primer AP1 and T-DNA border primers according to O'Malley et al. (2007) [32]. Eight primers to cover an approximately 200 bp region from each end of the T-DNA were designed based on the sequence of the T-DNA. Primer and adapter sequences used in this study are included in Additional file 7: Table S5. Target PCR products were recovered from gel, purified and ligated to a pCR 2.1-TOPO vector (Invitrogen, Carlsbad, CA, USA) for sequencing.
The identified 1279 bp sequence of the blueberry genome at the T-DNA insertion position was used to search the blueberry EST database (http://www.vaccinium.org). A 781 bp expressed sequence tag (EST) (CV091265) from nonacclimated floral buds, which has a 199 bp overlap with the 1279 bp sequence, was used as a reference to design the primers (i.e., E1F & E2R, E2F & E2R, and H2F & H2R) (Additional file 7: Table S5) for further extension of the 1279 bp sequence. PCR products were ligated to a pCR 2.1-TOPO vector (Invitrogen) for Sanger sequencing.

DNA sequencing and genome assembly
Total DNA was isolated from 200 mg young leaf tissues for each of the Mu-Legacy and Mu-Legacy-T 1 (a selected T 1 progeny of Mu-Legacy) sample, using a CTAB method [47]. The samples were purified using DNeasy Mini Kit (Qiagen, Valencia, CA, USA). The integrity of the DNA samples was assessed using electrophoresis. All samples were sequenced for 150-bp pair end reads with about 40-fold blueberry genome coverage using the Illumina HiSeq4000 platform at the Research Technology Support Facility at Michigan State University (East Lansing, MI, USA).
For each of the Mu-Legacy and Mu-Legacy-T 1 , SOAP-denovo2 and ABySS/1.9.0 (k = 64) were used to assemble genome sequences using the resources at the High Performance Computing Center at Michigan State University.

RNA preparation and sequencing
Nontransgenic 'Legacy' , Mu-Legacy, VcDDF1-Legacy (named II7 in our previous report [27,28])were used for RNA sequencing (RNA-seq). Three three-year old plants of each genotype were used for collecting leaf and flower bud tissues. Nonchilled flower buds (30 to 50 buds) were collected in November before the plants were exposed to an unheated hoop house for chilling. Fully chilled flower buds, 30-50, were harvested at the end of January from courtyard-chilled plants. Young leaf tissues were collected in February. All the collected tissues were frozen in liquid nitrogen immediately and stored at − 80°C.
Total RNA was isolated from 0.5 g tissues for each sample, using a CTAB method [47]. The samples were purified using the RNeasy Mini Kit (Qiagen, Valencia, CA, USA). The integrity of the RNA samples was assessed using electrophoresis. All samples had an RNA quality score above 8.0 prior to the submission for 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). The FastQC program (www.bioinformatics.babraham.ac.uk/projects/fastqc/) was used to assess the quality of sequencing reads for the per base quality scores ranging from 30 to 40.

Differential expression analysis
The RNA-seq reads of three biological replicates for each of the nontransgenic 'Legacy' Mu-Legacy, and VcDDF1-Legacy were analyzed. Each biological replicate was sequenced twice. The paired reads, two sets for each biological replicate, were aligned to the transcriptome reference Reftrinity [20] and the abundance of 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 differential expression analysis [48]. The differentially expressed (DE) genes or isoforms with the false discovery rate (FDR) values below 0.05 (p-value < 0.001) were used for further analyses of the flowering genes of blueberry. Most of the analyses were performed using the resources at the High Performance Computing Center at Michigan State University. Sequence alignment and phylogenetic tree analyses were conducted using CLC Sequence Viewer 7.

Quantitative reverse transcriptase polymerase chain reaction (qRT-PCR) of DE transcripts
The RNA samples used for RNA-sequencing, including samples of three biological replicates for each of nontransgenic 'Legacy' and Mu-Legacy, 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 microgram of RNA was diluted (volume 1: 4) in water and 1 μl/sample (25 ng) was used for PCR reactions. PCR primers that cover the VcRR2 region (E1F & E1R) were used. Eukaryotic translation initiation factor 3 subunit H was the internal control (1). QRT-PCR was performed in triplicate on an Agilent Technologies Stratagene Mx3005P (Agilent Technologies, Santa Clara, CA) using the SYBR Green system (Life Technologies,