Transcriptome-Based Identification of Candidate Flowering-Associated Genes of Blueberry in a Plant Factory with Artificial Lighting (PFAL) under Short-Day-Length Conditions

In blueberry (Vaccinium corymbosum L.), a perennial shrub, flower bud initiation is mediated by a short-day (SD) photoperiod and buds bloom once the chilling requirement is satisfied. A plant factory with artificial lighting (PFAL) is a planting system that can provide a stable and highly efficient growing environment for blueberry production. However, the characteristics of bud differentiation of blueberry plants inside PFAL systems are poorly understood. To better understand flower bud initiation and the flowering mechanism of blueberry in PFAL systems, the anatomical structure of apical buds under SD conditions in a PFAL system was observed using the southern highbush cultivar ‘Misty’ and a transcriptomic analysis was performed to identify the candidate flowering genes. The results indicated that the apical bud of ‘Misty’ differentiated gradually along with SD time course and swelled obviously when chilling was introduced. A total of 39.28 Gb clean data were generated, and about 20.31–24.11 Mb high-quality clean reads were assembled, yielding a total of 17370 differentially expressed genes (DEGs), of which 9637 were up-regulated and 7733 were down-regulated. Based on the functional annotation, 26 DEGs were identified including 20 flowering-related and 6 low-temperature DEGs, out of which the expressive level of four flowering-related DEGs (VcFT2, VcFPA, VcFMADS1, and VcCOP1) and two low-temperature-induced DEGs (VcTIL-1 and VcLTI 65-like) were confirmed by qRT-PCR with a good consistency with the pattern of transcriptome. Functional analysis indicated that VcFT2 was highly conserved with nuclear and cytoplasmic subcellular localization and was expressed mainly in blueberry leaf tissue. In Arabidopsis, ectopic overexpression of VcFT2 results in an early flowering phenotype, indicating that VcFT2 is a vital regulator of the SD-mediated flowering pathway in blueberry. These results contribute to the investigation of photoperiod-mediated flowering mechanisms of blueberry in PFAL systems.


Introduction
Blueberry (Ericaceae, Vaccinium spp.) is an economically important small woody perennial bush that includes several ecotypes, such as lowbush, semi-highbush, highbush, and rabbiteye blueberries, and interspecific hybrids, and the cultivation area of blueberry in China has statistically achieved about 66,400 ha with an annual yield of 0.35 million tons [1,2].Blueberry fruits contain exceptional amounts of anthocyanidins, flavonoids, organic acids, vitamin C, and antioxidants, and are famous for their high nutritional value and validated benefits for human health [3,4].Generally, the timing of flowering varies according to plant species, plant age, intrinsic genetic components, various environmental factors, and even the cultivar and the corresponding cultivation techniques used [5][6][7].For blueberry, flowering time is significantly associated with plant productivity and marketing time, and thus with the sale price of blueberry fruits [8,9].To pursue substantial economic benefits, an increasing number of blueberry plants have been cultivated by forcing culture using very early to early ripening cultivars, such as 'Misty', 'Emerald', 'O'Neal', etc., under protected conditions [10,11] or in greenhouses and controlled rooms [12].These conditions not only ensure safe overwintering but also earlier flowering and market time [13].Recently, blueberry cultivation using a plant factory with artificial lighting (PFAL) system was introduced to improve fruit yields and quality, generally, in China, and the average blueberry fruit yield in greenhouse or open-field farms was approximately 5.27 tons per hectare [14,15].The PFAL system is a modern and highly developed planting platform [16,17], inside which the growth parameters are always independent of seasonal changes and can be artificially manipulated according to plant growth dynamics [18,19].However, the annual growth cycle of blueberries and especially the flowering event inside the PFAL system is significantly different from that of traditional open-field farms or greenhouses.Knowledge about flower bud initiation, differentiation, development, and blooming, as well as the possible regulatory mechanism associated with the biological flowering process of blueberry in PFAL systems is still very limited.Therefore, research into blueberry flower bud development and elucidation of the underlying mechanism that controls the flowering event in closed PFAL environments is urgently required.
Previous studies have confirmed that blueberry flower bud differentiation and development is a phytochrome-mediated short-day (SD) response, which usually occurs under an 8 h photoperiod, but temperature also has a profound effect on the photoperiodic induction of flower buds [20,21].For instance, flower bud initiation in the lowbush blueberry (Vaccinium angustifolium Ait.) occurred under 11 or 13 h photoperiods with temperatures of either 10 or 21 • C, and flowers were produced under an 8 h photoperiod SD conditions for 2-4 weeks at constant 18 • C [22].In rabbiteye blueberry (Vaccinium ashei Reade) [23], 6 weeks of inductive short days are required for normal flower bud formation.In addition, blueberry is a chilling-mediated flowering species that requires a certain number of chilling hours before flowering; fully chilled blueberry plants flower normally, whereas non-chilled plants might not flower [24].Under natural conditions, most blueberry cultivars initiate flower buds in autumn and bloom the following spring [25].In Shanghai, China, the blueberry buds begin to differentiate in July as soon as the harvest is finished, under a SD length rhythm.In this period, the subtropical hot and humid climate forces the new shoots to stop expanding and quickly cap off, and the process of flower bud differentiation lasts until late November or early December.The buds then undergo a shallow dormancy, begin to swell rapidly, and bloom once the chilling requirement has been satisfied.To reveal the genetic and physiological mechanisms of blueberry flowering events, numerous studies have been conducted in recent decades, and a large number of flowering-related genes have been successfully identified, such as VcFT [26,27], SOC1-like MADS-box [28], VcRR2 [29], SPL [13], and VcSVP [30].FT is a major integrator of signaling that stimulates the transition of meristem tissue into flower buds [31], is expressed mainly in leaves, and the resultant FT protein travels through the phloem to the shoot apical meristem, where it induces flower determination [32].Constitutive expression of VcFT enables precocious flowering in transgenic blueberry plants, indicating that VcFT positively induces flower bud formation [26,27]; however, VcFT does not seem to be the only factor affecting chillingmediated blueberry flowering [24].Thus, the relationship between flowering and chilling levels needs further study to determine whether the SD length conditions could trigger bud differentiation in closed PFAL systems, and whether the differentiated buds could bloom successfully when chilling is introduced under SD conditions.
Transcriptome analysis is a powerful tool to identify differentially expressed genes (DEGs) and regulatory networks involved in complex biological processes.RNA-sequencing (RNA-Seq) technology has been successfully applied in blueberry to identify candidate genes or transcripts involving multiple specific traits, such as flowering [30,33], adventitious root formation [34], anthocyanin accumulation [35,36], flavonoid biosynthetic pathway [37], and fruit development [38,39].Therefore, to better understand the process of flower bud development, the effects of SD length and chilling on flower bud initiation of blueberry in a PFAL system were surveyed in this study, and transcriptomic analysis was carried out to identify the key candidate flowering genes.This is the first report to reveal the molecular regulatory basis of flowering genes in blueberry using the PFAL system, which will contribute to the understanding of blueberry plant growth and flower development in plant factories with artificial lighting.

Flower Bud Differentiation of 'Misty' in the PFAL System
To clarify the mechanism of blueberry flower bud differentiation in the PFAL system, we followed the bud developmental process of 'Misty' plants histologically.At the B1 stage, when the 'Misty' plants had been just switched from long-day (LD) length to SD length conditions, the apical meristem was small and closed (Figure 1a).At B2, 10 d after SD treatment, the bud scales began to loosen, and the apical meristems swelled and thickened to a hemispherical shape (Figure 1b).At B3, 20 d after SD treatment, the apical meristems further enlarged and became nearly elliptical (Figure 1c).At B4, 30 d after SD treatment, the flower bud primordium formed (Figure 1d).At B5, 50 d after SD treatment, when chilling was introduced for 20 d (SD + chilling treatment), the pistil primordium emerged (Figure 1e).At B6, after 60 d of SD treatment + 30 d of chilling, the differentiated flower buds, including apical and lateral flower buds, were observed (Figure 1f).The results suggest that blueberry flower buds can be initiated and differentiated in the PFAL system under SD length conditions combined with chilling treatment.
Int. J. Mol.Sci.2024, 25, x FOR PEER REVIEW 3 of 20 pathway [37], and fruit development [38,39].Therefore, to better understand the process of flower bud development, the effects of SD length and chilling on flower bud initiation of blueberry in a PFAL system were surveyed in this study, and transcriptomic analysis was carried out to identify the key candidate flowering genes.This is the first report to reveal the molecular regulatory basis of flowering genes in blueberry using the PFAL system, which will contribute to the understanding of blueberry plant growth and flower development in plant factories with artificial lighting.

Flower Bud Differentiation of 'Misty' in the PFAL System
To clarify the mechanism of blueberry flower bud differentiation in the PFAL system, we followed the bud developmental process of 'Misty' plants histologically.At the B1 stage, when the 'Misty' plants had been just switched from long-day (LD) length to SD length conditions, the apical meristem was small and closed (Figure 1a).At B2, 10 d after SD treatment, the bud scales began to loosen, and the apical meristems swelled and thickened to a hemispherical shape (Figure 1b).At B3, 20 d after SD treatment, the apical meristems further enlarged and became nearly elliptical (Figure 1c).At B4, 30 d after SD treatment, the flower bud primordium formed (Figure 1d).At B5, 50 d after SD treatment, when chilling was introduced for 20 d (SD + chilling treatment), the pistil primordium emerged (Figure 1e).At B6, after 60 d of SD treatment + 30 d of chilling, the differentiated flower buds, including apical and lateral flower buds, were observed (Figure 1f).The results suggest that blueberry flower buds can be initiated and differentiated in the PFAL system under SD length conditions combined with chilling treatment.

RNA-Sequencing and Identification of Differentially Expressed Genes
To understand the molecular basis of flowering bud differentiation of 'Misty' in PFAL under SD length conditions, six samples (B1, B2, B3, B4, B5, and B6) were subjected to total RNA extraction and transcriptome sequencing.After removing the adaptor and low-quality sequence, a total of 39.28 Gb of clean data were generated, containing 25. 81

Sequencing Annotation
All unique sequences were further annotated based on BLASTx searches against eight databases, i.e., clusters of orthologous genes (COG), GO, KEGG, KOG, NR, Pfam, Swiss-Prot, and evolutionary genealogy of genes: Non-supervised Orthologous Groups (eggNOG) databases.Using a cut-off E-value of 10 −5 , a high number of differentially expressed genes (DEGs) was identified, the number of annotated DEGs were shown detail in Table 2. Notes: Library (i.e., B1-vs-B2, B2-vs-B3, etc.) indicate pairwise between each two groups; the second column to the last represents the number of differentially expressed genes (DEGs) annotated by the functional databases separately.
GO analysis was performed to classify the functions of the 21,287 predicted unigenes into three main categories: biological process (3468), cellular component (5480), and molecular function (12,339) (Figure S1).In the biological process category, unigenes involved in 'metabolic process', 'cellular process', and 'single-organism process' were dominant.In the cellular component category, a great number of unigenes were mainly categorized as 'cell', 'cell part', 'membrane', 'organelle', and 'membrane part'.In the molecular function category, most of the unigenes were involved in 'catalytic activity' and 'binding' (Figure S1).

Flowering-Related DEGs of Blueberry in PFAL
Short-day length is one of the key environmental factors that induce blueberry apical buds to complete flower bud differentiation.In addition, differentiated flower buds

Flowering-Related DEGs of Blueberry in PFAL
Short-day length is one of the key environmental factors that induce blueberry apical buds to complete flower bud differentiation.In addition, differentiated flower buds require certain chilling or cooling hours conditions to bloom.Therefore, to identify the candidate DEGs involved in flowering and the chilling pathway, the distribution and functional annotation of unigenes was calculated based on the FPKM values and GO and KEGG enrichment with p = 0.05, and a total of 26 DEGs were identified, including 20 floweringrelated and six low-temperature DEGs, and their gene expression patterns are illustrated in a heatmap (Figure 3).require certain chilling or cooling hours conditions to bloom.Therefore, to identify the candidate DEGs involved in flowering and the chilling pathway, the distribution and functional annotation of unigenes was calculated based on the FPKM values and GO and KEGG enrichment with p = 0.05, and a total of 26 DEGs were identified, including 20 flowering-related and six low-temperature DEGs, and their gene expression patterns are illustrated in a heatmap (Figure 3).4).The observation indicated that the expression patterns of VcFT2, VcFMADS1, and VcFPA were similar, and that the three genes were up-regulated gradually with SD treatment (Figure 4a-c).VcCOP1 expression exhibited a trend opposite to those of VcFT2, VcFMADS1, and VcFPA.VcCOP1 was expressed highly during the B1-B4 treatments, but down-regulated markedly in the B5 and B6 treatments (Figure 4d), suggesting that VcCOP1 might play a negative regulatory role during blueberry flower bud differentiation in the PFAL system.The expression levels of VcTIL-1 and VcLTI 65-like were relatively low; however, they were both up-regulated significantly during the chilling phase (B6 treatment) (Figure 4e,f), indicating that the two genes might be important for the chilling requirement stage (i.e., low-temperature-induced vernalization pathway).qRT-PCR analysis was performed to confirm the expression of the six DEGs, indicating a good reproducibility with the RNA-Seq data (Figure 4).4).The observation indicated that the expression patterns of VcFT2, VcFMADS1, and VcFPA were similar, and that the three genes were up-regulated gradually with SD treatment (Figure 4a-c).VcCOP1 expression exhibited a trend opposite to those of VcFT2, VcFMADS1, and VcFPA.VcCOP1 was expressed highly during the B1-B4 treatments, but down-regulated markedly in the B5 and B6 treatments (Figure 4d), suggesting that VcCOP1 might play a negative regulatory role during blueberry flower bud differentiation in the PFAL system.The expression levels of VcTIL-1 and VcLTI 65-like were relatively low; however, they were both up-regulated significantly during the chilling phase (B6 treatment) (Figure 4e,f), indicating that the two genes might be important for the chilling requirement stage (i.e., low-temperature-induced vernalization pathway).qRT-PCR analysis was performed to confirm the expression of the six DEGs, indicating a good reproducibility with the RNA-Seq data (Figure 4).

Gene Expression Analysis
Based on the reported chilling-mediated flowering pathway of blueberry reported in the literature [24] and the candidate DEGs identified in present study, a potential regulatory network that controls flower bud differentiation and thus flowering under SD conditions in PFAL system is derived.It is speculated that the SD length conditions would induce flower bud differentiation and trigger the expression of VcFT2, VcFMADS1, VcFPA, VcCOP1, and other flowering-related genes.When the temperature-associated genes or transcript factors, such as VcTIL-1 and VcLTI 65-like genes, were stimulated for high expression following introduction of chilling, the blueberry plants eventually flowered (Figure 5).Based on the reported chilling-mediated flowering pathway of blueberry reported in the literature [24] and the candidate DEGs identified in present study, a potential regulatory network that controls flower bud differentiation and thus flowering under SD conditions in PFAL system is derived.It is speculated that the SD length conditions would induce flower bud differentiation and trigger the expression of VcFT2, VcFMADS1, VcFPA, VcCOP1, and other flowering-related genes.When the temperature-associated genes or transcript factors, such as VcTIL-1 and VcLTI 65-like genes, were stimulated for high expression following introduction of chilling, the blueberry plants eventually flowered (Figure 5).

Isolation and Characteristics Analysis of the VcFT2 Gene
The ORF of the VcFT2 gene was 525 bp long and encoded a protein of 174 amino acids.For the subcellular localization assay, the coding region of VcFT2 was fused to GFP to form the fusion gene VcFT2::GFP driven by 35S promoter in Vector 1300.Both the fusion gene and GFP control plasmid were expressed transiently in Nicotiana benthamiana leaves.The fluorescence of VcFT2::GFP was observed simultaneously in the nucleus and cytoplasm (Figure 6a).The putative VcFT protein showed a high similarity (Figure 6b, 99.81%) to the known Vaccinium corymbosum Flowering locus T (MN708233.1),and when the nucleotide sequence was aligned against the NCBI database using BLASTn, only a single nucleotide difference was observed between the two VcFTs; therefore, the VcFT gene ob-

Isolation and Characteristics Analysis of the VcFT2 Gene
The ORF of the VcFT2 gene was 525 bp long and encoded a protein of 174 amino acids.For the subcellular localization assay, the coding region of VcFT2 was fused to GFP to form the fusion gene VcFT2::GFP driven by 35S promoter in Vector 1300.Both the fusion gene and GFP control plasmid were expressed transiently in Nicotiana benthamiana leaves.The fluorescence of VcFT2::GFP was observed simultaneously in the nucleus and cytoplasm (Figure 6a).The putative VcFT protein showed a high similarity (Figure 6b, 99.81%) to the known Vaccinium corymbosum Flowering locus T (MN708233.1),and when the nucleotide sequence was aligned against the NCBI database using BLASTn, only a single nucleotide difference was observed between the two VcFTs; therefore, the VcFT gene obtained in the present study was named VcFT2 and deposited in GenBank with accession number OR660085.A comparison of the deduced amino acid sequences of VcFT2 in GenBank indicated that VcFT2 had high similarity with FT-like homologous of Camellia sinensis (AB741571.1,85.90%) [40], Actinidia chinensis (KX611595.1,84.38%) [41], Fagus crenata (AB775532.1,82.29%) [42], Hydrangea macrophylla (MF374627.1,82.29%), and Betula luminifera (JQ951966.1,81.00%) [43], suggesting that the VcFT2 obtained in the present study belonged to the FT-like superfamily.Subsequently, seven known FT genes, including Vaccinium corymbosum (VcFT, QNM38066) [27], Arabidopsis thaliana (AtFT, AAF03936) [44], Solanum lycopersicum (SlSP3D, AAO31792) [45], Populus deltoids (PdFT1, AAS00056) [46], Poplar tremula (PtFT1, ABD52003) [47], Malus domestica (MdFT3, ACV92037) [48], and Pyrus pyrifolia (PpFT2atw, BAK74837), were used to analyze the phylogenetic relationships among VcFT2 and FT-like proteins, revealing that VcFT2 had the closest kinship relationship with blueberry VcFT (Figure 6b).Multiple sequence alignments based on putative amino acid sequences indicated that VcFT2 possessed typical features of the FT-like protein subfamily (Figure 6c).It not only contained the conserved amino acid residues Tyr (Y) and Gln (Q), but also a highly conserved 14-amino acid stretch (P-loop domain) and LYN motif (Figure 6c), suggesting that VcFT2 is a putative ortholog of FT-like genes in blueberry.To determine the expression patterns of VcFT2 in blueberry, the tissue-specific ex- To determine the expression patterns of VcFT2 in blueberry, the tissue-specific expression of VcFT2 was analyzed using quantitative real-time reverse transcription (qRT)-PCR.This indicated that VcFT2 could be detected in all investigated tissues, with large differences in expression levels (Figure 7).VcFT2 was expressed primarily in leaves, followed by higher expression in mature fruits, young fruits, opened flowers, and unopened flowers, with lower expression in root and stem tissues (Figure 7).

Ectopic Expression of VcFT2 Promotes Flowering in Arabidopsis
To further investigate the function of VcFT2, transgenic Arabidopsis plants overexpressing VcFT2 combined with 35S promoter were generated, and three independent T3 lines, i.e., VcFT2-OE3, VcFT2-OE4, and VcFT2-OE5, were selected for flowering time analysis.The transgenic VcFT2-OE plants flowered at least 8 d earlier than the non-transgenic wild-type plants (Figure 8).Furthermore, both the number of flower buds and flower bud rate of the three lines VcFT2-OE3, VcFT2-OE4, and VcFT2-OE5 were significantly higher than those of the wild-type plants (Figure 9a,b), suggesting that VcFT2 in Arabidopsis might act as a flowering activator, thus leading to an earlier flowering phenotype.As for the transgenic lines, no significant differences in flower bud number and rate were observed among VcFT2-OE3, VcFT2-OE4, and VcFT2-OE5 on day 0 to day 2, and day 10; however, the flower bud number and rate of VcFT2-OE4 on day 4 to day 8 was significantly higher than those of VcFT2-OE3 and VcFT2-OE5 (Figure 9a,b).Vegetative growth of all transgenic lines of VcFT2-OE Arabidopsis was suppressed after approximately 63 d (Figure 10a).The plant height, stem diameter, rosette leaf number, and leaf area of VcFT2-OE3, VcFT2-OE4, and VcFT2-OE5 were significantly lower than those of the non-transgenic wild-type plants (Figure 10b-e), indicating slower vegetative growth or dwarf traits of plants after VcFT2 overexpression.As for the transgenic lines, the stem height of VcFT2-OE3 was significantly higher than that of VcFT2-OE4, and the differences in stem height between VcFT2-OE4 and VcFT2-OE5 as well as that between VcFT2-OE3 and VcFT2-OE5 were not significant (Figure 10b).The stem diameter, rosette leaf number, and leaf area of VcFT2-OE3 were significantly higher than those of VcFT2-OE4 and VcFT2-OE5; however, there were no significant differences in stem diameter, rosette leaf number, and leaf area between VcFT2-

Ectopic Expression of VcFT2 Promotes Flowering in Arabidopsis
To further investigate the function of VcFT2, transgenic Arabidopsis plants overexpressing VcFT2 combined with 35S promoter were generated, and three independent T3 lines, i.e., VcFT2-OE3, VcFT2-OE4, and VcFT2-OE5, were selected for flowering time analysis.The transgenic VcFT2-OE plants flowered at least 8 d earlier than the non-transgenic wild-type plants (Figure 8).Furthermore, both the number of flower buds and flower bud rate of the three lines VcFT2-OE3, VcFT2-OE4, and VcFT2-OE5 were significantly higher than those of the wild-type plants (Figure 9a,b), suggesting that VcFT2 in Arabidopsis might act as a flowering activator, thus leading to an earlier flowering phenotype.As for the transgenic lines, no significant differences in flower bud number and rate were observed among VcFT2-OE3, VcFT2-OE4, and VcFT2-OE5 on day 0 to day 2, and day 10; however, the flower bud number and rate of VcFT2-OE4 on day 4 to day 8 was significantly higher than those of VcFT2-OE3 and VcFT2-OE5 (Figure 9a,b).Vegetative growth of all transgenic lines of VcFT2-OE Arabidopsis was suppressed after approximately 63 d (Figure 10a).The plant height, stem diameter, rosette leaf number, and leaf area of VcFT2-OE3, VcFT2-OE4, and VcFT2-OE5 were significantly lower than those of the non-transgenic wild-type plants (Figure 10b-e), indicating slower vegetative growth or dwarf traits of plants after VcFT2 overexpression.As for the transgenic lines, the stem height of VcFT2-OE3 was significantly higher than that of VcFT2-OE4, and the differences in stem height between VcFT2-OE4 and VcFT2-OE5 as well as that between VcFT2-OE3 and VcFT2-OE5 were not significant (Figure 10b).The stem diameter, rosette leaf number, and leaf area of VcFT2-OE3 were significantly higher than those of VcFT2-OE4 and VcFT2-OE5; however, there were no significant differences in stem diameter, rosette leaf number, and leaf area between VcFT2-OE4 and VcFT2-OE5 (Figure 10c-e).

Discussion
Bud differentiation in plant species such as blueberry is a complex physiological and biochemical process influenced by multiple factors including temperature, photoperiod, and day length.Blueberry is a photoperiod-sensitive berry shrub, and flower bud initiation appears to be promoted by exposure to short days, typically with a critical day length

Discussion
Bud differentiation in plant species such as blueberry is a complex physiological and biochemical process influenced by multiple factors including temperature, photoperiod, and day length.Blueberry is a photoperiod-sensitive berry shrub, and flower bud initiation appears to be promoted by exposure to short days, typically with a critical day length of less than 12 h [20,21,49].Under controlled conditions, 8 weeks of 8, 10, or 12 h photoperiods will trigger and drive flower bud development in blueberries [19,22,23], and the flower bud number is positively associated with the length of exposure to short days [50].Although highly advanced plant factory systems can provide a stable and high-efficiency growing environment for producing high-yield and high-quality blueberry fruits [14,15,51], the characteristics of the growth cycle of blueberry plants inside the PFAL system are poorly understood.In addition, it is uncertain whether the buds of blueberry plants can transition autonomously from the vegetative to the reproductive phase under controlled PFAL conditions.Sophisticated studies on temperature, day length, and photoperiod settings are needed to reveal the underlying flowering mechanisms in closed PFAL environments.In this study, an 8 h SD length (8/16 h day/night) was artificially set to induce flower bud differentiation using the low-chilled southern highbush 'Misty'.Our observations of histoanatomical structures indicated that the apical bud of 'Misty' differentiated gradually along with the SD time course (Figure 1a-f), and when chilling (8/16 h day/night, 10 • C about 10-14 d) was introduced to satisfy the cold requirement levels of 'Misty', the differentiated buds visibly swelled.Moreover, when the indoor PFAL parameters were returned to 16/8 h (day/night) and 22 • C (the temperature was increased gradually), the differentiated buds burst and bloomed normally, a high fruit setting rate was recorded, and these blueberry fruits matured gradually (Figure S2).This demonstrates the possibility of inducing bud differentiation artificially in a PFAL system under certain SD photoperiod conditions.
Transcriptome analysis is a useful tool for identifying potential flowering pathway genes [52][53][54].Using this approach, a large number of DEGs related to the flowering regulatory network have been identified in blueberry, such as the conserved flowering pathway genes VcFT, VcFD, VcTFL1, VcLFY, VcAPP6, and MADS-box genes; the photoperiod-related genes VcCOL2 and VcCOL5; and some downstream genes of VcFT including VcSOC1, VcAP1, VcCAL1, and VcFUL [24,55].These candidate DEGs have provided vital evidence for studying the mechanism of flowering in blueberry [56].However, transcript information on blueberry flower bud development under SD length conditions in a PFAL system has not been reported.In the present study, comparative transcriptome analyses performed against different stages of floral bud in a PFAL system generated 20.31-24.11Mb clean reads (Table 1) and numerous unigenes participating in blueberry floral bud development were screened out (Table 2, Figures 2 and 3).After functional annotation, six key candidate DEGs involved in the flowering pathway were identified, namely VcFT2, VcFMADS1, VcFPA, VcCOP1, VcTIL-1, and VcLTI 65-like genes.These DEGs were significantly up-regulated or down-regulated along with the time course of SD length (Figures 4 and 5), indicating their potential roles in the SD-mediated blueberry flowering pathways.These results were partly consistent with those reported by Song and Chen [24].However, the functional characteristics and molecular mechanisms of the flowering DEGs identified in this study remain to be elucidated.
FLOWERING LOCUS T (FT) is a major flowering pathway integrator and a universal promoter of flowering, and functional analyses of FT and FT-like genes in numerous plant species have been reported [57,58].FT protein has a highly conserved P-loop domain with a 14 amino acid (LGRQTVYAPGWRQN) motif and conserved YLH and LYN triads [59,60].To verify FT-like genes in blueberry, in this study, a 525 bp VcFT2 homolog was cloned from the flower bud of 'Misty' based on RNA-Seq data (Figure 6).Our results showed subcellular localization of VcFT2 in the nucleus and cytoplasm (Figure 6a).The phylogenetic tree suggested that VcFT2 clustered closely with FT or FT-like genes identified in the blueberry cultivar 'Bluecrop' and other species, and the sequence of VcFT2 was highly conserved (Figure 6b,c).Alignment-based nucleotide sequences suggested that the CDs of VcFT and VcFT2 only had a single nucleotide difference, possibly because the blueberry cultivar used to isolate FT was different, or because of the heterozygosity of the blueberry chromosome.Day length or photoperiod is perceived by leaves and induces a special signal called a florigen that moves through the phloem to the shoot apex [61][62][63].Many FT or FT-like genes have been found to be highly expressed in leaves and are regulated by photoperiod, such as FaFT1 in strawberries [64] and GmFT2a in soybeans [65].Although FT is a mobile signal originating from leaves [66,67], a large number of FT or FT-like genes have been found mainly in reproductive tissues, such as JcFT, which is highly expressed in the flower buds of Jatropha curcas L. [68] and GhFT1 expressed primarily in the stamens, sepals, petals, and carpels of Gossypium hirsumtum L. [69].In the present study, VcFT2 was highly expressed in leaves and secondarily in fruits and flowers (Figure 7).This is partly consistent with the findings in other species described above, indicating not only divergent expression patterns of VcFT homologs in blueberries, but also suggesting that FT genes might play multifaceted roles in plant development in addition to flowering time control [58,70].In addition, ectopic overexpression of VcFT2 in Arabidopsis thaliana L. resulted in an early flowering phenotype (Figures 8-10), indicating that the VcFT2 gene isolated in this study might act as a floral stimulator, that is, promote the transition from the vegetative to the reproductive phase.The results obtained in this study were consistent with those previously published by the Song research group [27,28,55].However, the related functional genes or transcription factors (TFs) upstream and downstream of VcFT2 remain unclear, and further study of VcFT2 and its related flowering regulatory genes or TFs is required.Future studies should focus on mining and verifying the crosstalk between VcFT2 and its upstream and downstream interacting genes.Furthermore, how VcFT2 participates in the response of blueberry plants to photoperiod (SD) and temperature (chilling) conditions, and thus, the regulation of flowering, is unknown.

Plant Samples and Experimental Design
Three-year-old southern highbush blueberry cultivar (Vaccinium corymbosum 'Misty') plants were used in the present study (Figure 11).In total, 24 'Misty' plants were grown in room 1 of the PFAL system.All the plants were potted in the medium containing peatmoss/perlite/vermiculite (1:1:1, v/v/v) [14] and were irrigated daily using an automatic AEtrium-4 system (AEssence Grows, Santa Clara, CA, USA) with modified Hoagland nutrient solution (N:P:K = 1:1:1, pH = 5.0 ± 0.5, EC-value was 1500 ± 200 µS/cm).The plant height, stem diameter, shoot number, and shoot length of 'Misty' plants used in the study were 45.88 ± 12.13 cm, 10.91 ± 0.57 mm, 16.75 ± 3.89, and 24.86 ± 6.54 cm, respectively.Before treatment, the PFAL environmental parameters were set as follows: room temperature of 25 ± 2 • C/22 ± 2 • C (day/night), photoperiod of 16/8 h (day/night), CO 2 concentration of 400 µmol mol −1 , relative air humidity of 60-70%, and light intensity over blueberry canopy of 200 µmol m −2 s −1 [14].Generally, in greenhouse or open-field farms, blueberry cultivars begin to differentiate flower buds under an SD length rhythm; however, the differentiated flower buds require a fixed quantity of chilling hours during winter endo-dormancy for vernalization [71], and over 80% of flower buds cannot bloom normally without a chilling period [24].However, the characteristics flower bud formation and flowering events in PFAL systems are not clear.Therefore, in the present study, SD in combination with chilling treatments were designed to induce blueberry flower bud differentiation.The present study was carried out from March to May in 2020, when new shoots of 'Misty' stopped expanding.The photoperiod was changed to an SD length, i.e., 8/16 h (marked as 0 d, B1, collected on 3 March 2020), to induce flower bud differentiation, and the other growth parameters remained unchanged.Apical buds from the shoots were sampled every 10 d for one month; i.e., the buds were sampled at 10 d (B2, collected on 13 March 2020), 20 d (B3, collected on 23 March 2020), and 30 d (B4, collected on 2 April 2020).Subsequently, all SD-treated 'Misty' plants were transferred into room 3 of the PFAL for chilling at 10 ± 2 • C (day/night) under the same photoperiod, i.e., 8/16 h (day/night).Apical buds were sampled at 20 d (B5, collected on 22 April 2020) and 30 d (B6, collected on 2 May 2020) after the chilling treatment.

Observation of Blueberry Bud Differentiation in the PFAL System
To observe the process of blueberry flower bud formation and differentiation in the PFAL system, ten apical buds were randomly collected at B1, B2, B3, B4, B5, and B6, as described above.The samples were fixed using formaldehyde-acetic acid-alcohol (FAA) fixation solution, which was consisted by formaldehyde, acetic acid and 70% alcohol with the ratio of 5:5:90 (v/v/v).Before observation, all the sampled buds were softened for about 14 d within 4% ethylenediamine solution following the method of An et al. [34].Subsequently, dehydration was carried out using a graded ethanol series.Vitrification was performed using dimethylbenzene before embedding the samples in paraffin.Finally, 10 µmthick sections were cut using a rotary microtome (LEICA, RM2265), and photographs were captured under a light microscope (NIKON ECLIPSE E200).

Total RNA Extraction, RNA-Sequencing, and De Novo Transcriptome Assembly
Ten apical buds were collected at each sampling date (B1, B2, B3, B4, B5, and B6), total RNA was extracted using the Polysaccharide-and Polyphenolic-rich RNAprep Pure Plant Kit (TIANGEN, Beijing, China) according to the manufacturer's protocol.The RNA integrity was validated using an Agilent 2100 Bioanalyzer (Agilent, Santa Clara, CA, USA).RNA-Seq libraries were constructed using a TruSeq Stranded mRNA LT Sample Prep Kit (Illumina, San Diego, CA, USA) and applied to an Illumina NovaSeq 6000 sequencing system guide for RNA-Seq analysis by Biomarker Technologies Corporation (Beijing, China).Due to limited space in the blueberry plant factory, a small amount of apical bud samples (10 buds at each sampling period) were collected; therefore, the RNA-Seq analysis was performed only once with no biological replicates.

Sequence Annotation and Identification of DEGs
Gene functions were annotated by aligning the genes with the nucleotide (NT) sequences of the National Center for Biotechnology Information (NCBI), non-redundant (NR) [72], Swiss-Prot [73], Pfam [74], Kyoto Encyclopedia of Genes and Genomes (KEGG) [75], Clusters of Orthologous Groups of proteins (COGs) [76], eggNOG [77], and Gene Ontology (GO) databases [78].Gene expression levels were calculated using the fragments per kilobase per million reads (FPKM) method, which is the mostly used RNA-Seq analysis method because both sequencing depth and unigene length are considered by the method.Differential expression analysis was performed using the edgeR package [79].Genes with fold change ≥ 2 and FDR ≤ 0.05 were identified as DEGs.

Observation of Blueberry Bud Differentiation in the PFAL System
To observe the process of blueberry flower bud formation and differentiation in the PFAL system, ten apical buds were randomly collected at B1, B2, B3, B4, B5, and B6, as described above.The samples were fixed using formaldehyde-acetic acid-alcohol (FAA) fixation solution, which was consisted by formaldehyde, acetic acid and 70% alcohol with the ratio of 5:5:90 (v/v/v).Before observation, all the sampled buds were softened for about 14 d within 4% ethylenediamine solution following the method of An et al. [34].Subsequently, dehydration was carried out using a graded ethanol series.Vitrification was performed using dimethylbenzene before embedding the samples in paraffin.Finally, 10 µm-thick sections were cut using a rotary microtome (LEICA, RM2265), and photographs were captured under a light microscope (NIKON ECLIPSE E200).

Total RNA Extraction, RNA-Sequencing, and De Novo Transcriptome Assembly
Ten apical buds were collected at each sampling date (B1, B2, B3, B4, B5, and B6), total RNA was extracted using the Polysaccharide-and Polyphenolic-rich RNAprep Pure Plant Kit (TIANGEN, Beijing, China) according to the manufacturer's protocol.The RNA integrity was validated using an Agilent 2100 Bioanalyzer (Agilent, Santa Clara, CA, USA).RNA-Seq libraries were constructed using a TruSeq Stranded mRNA LT Sample Prep Kit (Illumina, San Diego, CA, USA) and applied to an Illumina NovaSeq 6000 sequencing system guide for RNA-Seq analysis by Biomarker Technologies Corporation (Beijing, China).Due to limited space in the blueberry plant factory, a small amount of apical bud samples (10 buds at each sampling period) were collected; therefore, the RNA-Seq analysis was performed only once with no biological replicates.

Sequence Annotation and Identification of DEGs
Gene functions were annotated by aligning the genes with the nucleotide (NT) sequences of the National Center for Biotechnology Information (NCBI), non-redundant (NR) [72], Swiss-Prot [73], Pfam [74], Kyoto Encyclopedia of Genes and Genomes (KEGG) [75], Clusters of Orthologous Groups of proteins (COGs) [76], eggNOG [77], and Gene Ontology (GO) databases [78].Gene expression levels were calculated using the fragments per kilobase per million reads (FPKM) method, which is the mostly used RNA-Seq analysis method because both sequencing depth and unigene length are considered by the method.Differential expression analysis was performed using the edgeR package [79].Genes with fold change ≥ 2 and FDR ≤ 0.05 were identified as DEGs.

Cloning and Sequence Analysis of VcFT2
The open reading frame (ORF) region of VcFT2 was amplified from the cDNA of 'Misty' flower bud using a conserved degenerate primer pair of VcFT2-F: 5 ′ -ATGCCACGGGATAG-GGATCC-3 ′ ; VcFT2-R 5 ′ -GCGTCGTCGTCCTCCTGAG-3 ′ .The VcFT2-F/R primer was designed based on the reference FT sequence (unigene ID: CUFF.18244) annotated using transcriptome data.The obtained sequence of VcFT2 was analyzed using DNAman 6.0 (LynnonBiosoft, Quebec, QC, Canada) and aligned against the NCBI database using BLASTn (www.ncbi.nlm.nih.gov/blast.cgi,accessed on 10 June 2021).A neighbor-joining (NJ) phylogenetic tree of VcFT2 and known FTs from other species was constructed using MEGA 6.0 software [80].Bootstrap values were calculated from 1000 replications.

Subcellular Location
The CDs region of VcFT2 without stop codon was amplified, then were fused with Green Fluorescent Protein (GFP) and inserted into pCAMBIA1300 vector.The transformation was carried out accordingly to the method described by Zuo et al. [81].Agrobacterium (GV3101) transformed with the target vectors was suspended in the infiltration buffer with the value of OD 600 was 0.6.Next, Nicotiana benthamiana leaves were injected.Infected plants were incubated for 48 h at 30 • C before observation under a NIKON C2-ER laser scanning confocal microscope (Nikon C2-ER, Nikon, Tokyo, Japan).

Arabidopsis Transformation
The CDS sequence of VcFT2 was inserted after the CaMV-35S promoter in the pCAMBIA-1300-GFP vector to form the fusion vector 35S-VcFT2-GFP.To further test the potential molecular biological function of VcFT2, the Arabidopsis Col-0 plants were transformed using A. tumefaciens containing the vector 35S-VcFT2-GFP via the floral dip method [82].Three T3 transgenic Arabidopsis lines (VcFT2-OE3, VcFT2-OE4, and VcFT2-OE5) possessed relatively high transcriptional levels of VcFT2 and these transgenic lines were selected for subsequent experiments.The flowering time, flower bud number, and flower bud rate of VcFT2 over-expressing plants and wild-type (WT) plants were recorded over a time course of 0-10 d after the Arabidopsis plants with 3-4 cotyledons were transplanted from the medium to soil composed by vermiculite and peat (2:1, v/v).For phenotype analyses, the stem height (cm), stem diameter (mm), rosette leaf number, and leaf area (mm 2 ) of the 63-day VcFT2 over-expressing plants and WT plants were measured.Each transgenic line or wild-type plant has three biological repeats and each biological repeat contains ten 63-day seedlings (i.e., the 63rd days after seeding).

Gene Expression Analysis by qRT-PCR
qRT-PCR analysis was performed according to the method of An et al. [34].Firststrand cDNA was synthesized using the PrimerScript TM RT reagent Kit with gDNA Eraser (RR047, Takara, Japan), the qRT-PCR reaction was performed on a LightCycler 480 Real-Time PCR System (Roche, Basal, Switzerland), and the program was as follows: 95 • C for 30 s, then repeated for 40 cycles at 95 • C for 5 s and 60 • C for 20 s.VcGAPDH was used as the internal reference gene.Each qRT-PCR reaction was performed in triplicate and the expression level of each gene was calculated by the values of 2 −∆∆Ct .The primer sequences for qRT-PCR analysis are listed in Supplementary Table S1.

Statistical Analysis
Data analysis was performed with one-way Analysis of Variance (ANOVA) using IBM SPSS Statistics 18 (IBM Corp., Armonk, NY, USA).Significant differences were compared by Duncan's multiple range tests and p ≤ 0.05 was considered significant.Graphics were created using GraphPad Prism 6.0 (GraphPad Software Inc., San Diego, CA, USA).

Figure 1 .
Figure 1.Observation of blueberry flower bud developmental process under short-day-length conditions in the PFAL system.Notes: subfigures (a)-(f) indicate the corresponding anatomical structure at sampled date of B1-B6, respectively.B1 indicates the date when the photoperiod was changed to an SD length, i.e., 8/16 h (marked as 0 d); B2 indicates 10 days after SD treatment; B3 indicates 20 days after SD treatment; B4 indicates 30 days after SD treatment; B5 indicates 50 days after SD treatment + 20 days after chilling (10 ± 2 °C, day/night) treatment; B6 indicates 60 days after SD treatment + 30 days after chilling treatment.BS, bud scale; AM, apical meristem; FBP, flower bud primordium; PI, pistil; AFB, apical flower bud; LFB, lateral flower bud; the bar indicates the scale.

Figure 1 .
Figure 1.Observation of blueberry flower bud developmental process under short-day-length conditions in the PFAL system.Notes: subfigures (a-f) indicate the corresponding anatomical structure at sampled date of B1-B6, respectively.B1 indicates the date when the photoperiod was changed to an SD length, i.e., 8/16 h (marked as 0 d); B2 indicates 10 days after SD treatment; B3 indicates 20 days after SD treatment; B4 indicates 30 days after SD treatment; B5 indicates 50 days after SD treatment + 20 days after chilling (10 ± 2 • C, day/night) treatment; B6 indicates 60 days after SD treatment + 30 days after chilling treatment.BS, bud scale; AM, apical meristem; FBP, flower bud primordium; PI, pistil; AFB, apical flower bud; LFB, lateral flower bud; the bar indicates the scale.

Figure 2 .
Figure 2. Statistics of differently expressed genes (DEGs).Note: (a) significantly up-or down-regulated unigenes using the −log10(p-value) and |log2FoldChange| ≥ 1 thresholds in B4-vs-B5; (b) graphical representation of overall DEGs in each library; (c,d) Venn diagram showing the number of up-regulated and down-regulated DEGs in each library.

Figure 2 .
Figure 2. Statistics of differently expressed genes (DEGs).Note: (a) significantly up-or downregulated unigenes using the −log10(p-value) and |log2FoldChange| ≥ 1 thresholds in B4-vs-B5; (b) graphical representation of overall DEGs in each library; (c,d) Venn diagram showing the number of up-regulated and down-regulated DEGs in each library.

Figure 3 .
Figure 3. Heatmap of the expression of flowering-related DEGs identified by transcriptomic analysis of 'Misty' under short-day length in the PFAL system.Note: heatmap was created normalized FPKM value data using the Log2(1 + FPKM) formula.

Figure 3 .
Figure 3. Heatmap of the expression of flowering-related DEGs identified by transcriptomic analysis of 'Misty' under short-day length in the PFAL system.Note: heatmap was created normalized FPKM value data using the Log 2 (1 + FPKM) formula.

Figure 5 .
Figure 5. Assumed gene networks that regulate flower bud differentiation of blueberry under shortday length and chilling conditions in PFAL system.

Figure 5 .
Figure 5. Assumed gene networks that regulate flower bud differentiation of blueberry under short-day length and chilling conditions in PFAL system.

Figure 6 .
Figure 6.Subcellular localization of VcFT2 in Nicotiana benthamiana leaves (a), phylogenetic analysis (b), and alignment of amino acid sequences of VcFT2 and FT-like proteins from other plant species (c), including blueberry (Vaccinium corymbosum: VcFT, QNM38066), Arabidopsis (Arabidopsis thaliana: AtFT, AAF03936), tomato (Solanum lycopersicum: SlSP3D, AAO31792), poplar (Populus deltoids: PdFT1, AAS00056; Poplar tremula: PtFT1, ABD52003), apple (Malus domestica: MdFT3, ACV92037) and pear (Pyrus pyrifolia: PpFT2atw, BAK74837).Note: The bar in subfigure a was 20 µm.The derived protein of VcFT2 was initially aligned to other known FT proteins using Clustal W. MEGA 6.0 was used to construct a phylogenetic tree using the neighbor-joining (NJ) method with 1000 bootstrap replications.Bootstrap percentages are shown at dendrogram branch points.The red circle in in subfigure b indicates the accession number of VcFT2.The unit for the scale bar in subfigure b displays branch lengths (0.02).Accession numbers are indicated at the right of the gene name.The asterisks indicate the conserved residues Tyr and Gln.The 14-amino-acid stretch (P-loop domain) and the LYN triad are marked with boxes.

20 Figure 7 .
Figure 7. Tissue-specific expression of VcFT2 in blueberry.Note: Roots, stems, and leaves were sampled at the same day when the new shoots stopped expanding.Flowers including unopened and opened flowers were collected when the blueberry plants began flowering.Young fruit was collected at the green fruit stage, and blue fruit indicated mature fruit.The bar indicates the standard error, * indicates significance at p ≤ 0.05.

Figure 7 .
Figure 7. Tissue-specific expression of VcFT2 in blueberry.Note: Roots, stems, and leaves were sampled at the same day when the new shoots stopped expanding.Flowers including unopened and opened flowers were collected when the blueberry plants began flowering.Young fruit was collected at the green fruit stage, and blue fruit indicated mature fruit.The bar indicates the standard error, * indicates significance at p ≤ 0.05.

20 Figure 8 .
Figure 8.Effect of ectopic expression of VcFT2 on phenotypic change in T3 transgenic Arabidopsis thaliana.Note: red arrow indicates the flower.

Figure 9 .
Figure 9. Statistics of flower bud number (a) and flower bud rate (b) for VcFT2-OE and WT plants.Notes: Different lowercase letters in column indicate significance at level of p ≤ 0.05.

Figure 8 . 20 Figure 8 .
Figure 8.Effect of ectopic expression of VcFT2 on phenotypic change in T3 transgenic Arabidopsis thaliana.Note: red arrow indicates the flower.

Figure 9 .
Figure 9. Statistics of flower bud number (a) and flower bud rate (b) for VcFT2-OE and WT plants.Notes: Different lowercase letters in column indicate significance at level of p ≤ 0.05.

Figure 9 .
Figure 9. Statistics of flower bud number (a) and flower bud rate (b) for VcFT2-OE and WT plants.Notes: Different lowercase letters in column indicate significance at level of p ≤ 0.05.

Figure 10 .
Figure 10.Phenotype observations of ectopically expressed VcFT2-OE and non-transgenic wild-type (WT) Arabidopsis thaliana 63 d after seed germination.Note: (a) growth phenotype transgenic VcFT2-OE and non-transgenic wild-type (WT) Arabidopsis thaliana; (b) stem height of VcFT2-OE and WT plants.(c) stem diameter of VcFT2-OE and WT plants; (d) rosette leaf number of VcFT2-OE and WT plants; (e) leaf area of VcFT2-OE and WT plants; Different lowercase letters in column indicate the significance at p ≤ 0.05.

Figure 10 .
Figure 10.Phenotype observations of ectopically expressed VcFT2-OE and non-transgenic wild-type (WT) Arabidopsis thaliana 63 d after seed germination.Note: (a) growth phenotype transgenic VcFT2-OE and non-transgenic wild-type (WT) Arabidopsis thaliana; (b) stem height of VcFT2-OE and WT plants.(c) stem diameter of VcFT2-OE and WT plants; (d) rosette leaf number of VcFT2-OE and WT plants; (e) leaf area of VcFT2-OE and WT plants; Different lowercase letters in column indicate the significance at p ≤ 0.05.

Figure 11 .
Figure 11.Blueberry plants grown under short-day-length conditions in the PFAL system.

Figure 11 .
Figure 11.Blueberry plants grown under short-day-length conditions in the PFAL system.

Table 1 .
). Summary of transcriptome of apical bud collected from blueberry cultivar Vaccinium corymbosum 'Misty' during flower bud differentiation in the PFAL system under short-day-length conditions.

Table 2 .
Summary statistics of unigenes annotated by different databases.