Comparative Transcriptomic Analysis of Differentially Expressed Transcripts Associated with Flowering Time of Loquat ( Eriobotya japonica Lindl.)

: Flowering is an important phenophase of plant species, however, knowledge about the regulatory mechanism controlling ﬂowering cues in loquat is limited. To identify candidate genes regulating ﬂowering time in loquat, we used RNA-Seq technology to conduct a comparative transcriptome analysis of differentiating apical buds collected from the early-ﬂowering variety ‘Baiyu’ and the late-ﬂowering variety ‘Huoju’. A total of 28,842 differentially expressed transcripts (DETs) were identiﬁed. Of these, 42 DETs controlled ﬂowering time while 17 other DETs were associated with the ABA signaling pathway. Compared with those in ‘Huoju’, EjFT , EjFY , EjFLK , and EjCAL1-like were signiﬁcantly upregulated in ‘Baiyu’. Moreover, transcripts of the ABA 8 (cid:48) -hydroxylases ( EjABH2 , EjABH4 , and EjABH4-like2 ), the ABA receptors ( EjPYL4/8 ), and the bZIP transcription factor EjABI5-like were upregulated in ‘Baiyu’ compared with ‘Huoju’. Hence, they might regulate loquat ﬂowering time. There was no signiﬁcant difference between ‘Baiyu’ and ‘Huoju’ in terms of IAA content. However, the ABA content was about ten-fold higher in the apical buds of ‘Baiyu’ than in those of ‘Huoju’. The ABA:IAA ratio sharply rose and attained a peak during bud differentiation. Thus, ABA is vital in regulating ﬂoral bud formation in loquat. The results of the present study help clarify gene transcription during loquat ﬂowering.


Introduction
Flowering is the onset of plant reproduction and an important phenophase. Understanding the mechanisms that control flowering events is essential for efficient fruit reproduction [1]. Flowering has evolved to maximize outcrossing and mitigate the damaging effects of frost. It involves elegant mechanisms that integrate complex signals including photoperiod and temperature [2,3]. Flower bud differentiation and, by extension, flowering itself, are influenced by sophisticated regulatory networks [4][5][6][7]. Several major pathways in Arabidopsis include floral integrator genes that govern flowering time. Photoperiod as well as circadian, ambient temperature, vernalization, autonomous, gibberellin, and aging pathways are implicated [8]. FLOWERING LOCUS T (FT) is a crucial mobile signal [9]. FT and FT-like play pivotal roles in the flowering of fruit trees in the Rosaceae family such as apple [10,11], pear [12], peach [13], and loquat [14][15][16]. Flowering time is affected by environmental conditions and phytohormones. Previous studies have demonstrated phytohormone crosstalk in flower induction [17,18]. ABA is generally considered a flowering suppressor [19]. In loquat, however, ABA plays a positive role in flower bud

Plant Materials
Two loquat (Eriobotrya japonica Lindl.) varieties with distinct flowering traits were used in this study. 'Baiyu' is an early-flowering cultivar while 'Huoju' is a late-flowering variety. Both were planted at the Jinshan Fruit Tree Experimental Station of the Shanghai Academy of Agricultural Sciences, Shanghai, China (30 • 47 27 ; 121 • 8 6 ). During flower bud differentiation (early July to late August in Shanghai), about 30 apical buds were randomly picked at 7-d intervals between 13 July 2017 and 24 August 2017. All samples were immediately frozen in liquid nitrogen and stored at −80 • C before being analyzed.

Endogenous Phytohormone Determination by LC-MS
About 0.5 g of fresh apical bud was ground in liquid nitrogen and homogenized in 5 mL of ethyl acetate. The suspension was vortexed for 1 min and stored in the dark at 4 • C for 12 h. The mixture was centrifuged at 10,000× g and 4 • C for 10 min. The supernatant was collected and the residue was extracted twice with 3 mL ethyl acetate at 10,000× g and 4 • C. The solution was vacuum-concentrated to dryness at 30 • C and the residue was dissolved in 250 mL chromatographic methanol. The solution was then passed through a 0.22-µm membrane filter and a 10 µL aliquot was injected into a liquid chromatograph-mass spectrometer (LC-MS). Changes in the levels of endogenous phytohormones including abscisic acid (ABA) and indoleacetic acid (IAA) at each sampling period were analyzed according to the methods of Niu et al. (2014) [28]. The mobile phase consisted of methanol and 0.5% (v/v) acetic acid dissolved in redistilled water. The standards were IAA (I2886) and ABA (A1094) (Sigma-Aldrich Corp., St. Louis, MO, USA). Calibration curves were plotted using 0 ng/mL, 2.5 ng/mL, 5.0 ng/mL, 10 ng/mL, 12.5 ng/mL, 25 ng/mL, 50 ng/mL, and 100 ng/mL of each standard. Linear regressions of the calibration curves are shown in Table S1.

Total RNA Extraction and RNA-Seq
Total RNA was extracted with a mirVana miRNA Isolation Kit (Thermo Fisher Scientific, Waltham, MA, USA) according to the manufacturer's protocol. RNA integrity was validated with an Agilent 2100 Bioanalyzer (Agilent Technologies, Santa Clara, CA, USA). RNA-Seq libraries were constructed with a TruSeq Stranded mRNA LTSample Prep Kit (Illumina, San Diego, CA, USA) and applied to an Illumina HiSeq X Ten sequencing platform for RNA-Seq analysis by Shanghai OE Biotech. Co. Ltd. (Shanghai, China).

Quality Control and de Novo Assembly
Raw data (reads) were processed with Trimmomatic [29]. Low-quality reads and those containing poly-N were removed to obtain clean reads. Adaptor and low-quality sequences were removed and the clean reads were assembled into expressed sequence tag clusters (contigs). They were then assembled de novo into transcripts by the paired-end method using Trinity v. trinityrnaseq_r20131110 [30]. The longest transcript was selected as a unigene based on sequence similarity and length, and used in the subsequent analysis.

Functional Annotation
Unigene function was annotated by aligning the unigenes with the databases NCBI non-redundant (NR), SwissProt, and Clusters of Orthologous Groups for Eukaryotic Complete Genomes (KOG). For this purpose, Blastx was used at a threshold E-value of 10 −5 [31]. Proteins with the highest numbers of hits to the unigenes were used to assign functional annotations to them. Gene ontology (GO) classification was performed by the mapping relationship between SwissProt and the GO terms. The unigenes were mapped to the Kyoto Encyclopedia of Genes and Genomes (KEGG) database to annotate their potential metabolic pathways [32].

Identification of Differentially Expressed Transcripts (DETs) Involved in Flowering Time
Fragments per kilobase per million mapped reads (FPKM) and read count values for each unigene were calculated with bowtie2 and eXpress [33][34][35]. Differentially expressed transcripts (DETs) were identified with the DESeq (2012) functions SizeFactors and nbinomTest [36]. p < 0.05 and FoldChange > 2 were set as the thresholds for significant differential expression. A hierarchical cluster analysis of the DETs was performed to explore transcript expression patterns. GO and KEGG pathway enrichment analyses of the DETs were performed on the hypergeometric distribution.

Gene Expression Analysis by qRT-PCR
Total RNA was extracted from the loquat apical buds collected at each sampling date. First-strand cDNA was synthesized from 1 g total RNA using the PrimerScript RT Reagent Kit with gDNA Eraser (RR047; TaKaRa Bio Inc., Shiga, Japan), diluted tenfold with ddH 2 O, and used as templates for qRT-PCR. The qRT-PCR mixture (10 µL) comprised 5 µL SYBR Premix ExTaq (RR820; TaKaRa), 0.5 µL of each primer (10 µM), 1 µL cDNA, and 2.5 µL RNase-free water. The qRT-PCR was performed on a LightCycler480 instrument (Roche Diagnostics, Basel, Switzerland) according to the manufacturer's instructions. The two-step qRT-PCR program was as follows: 95 • C for 30 s, followed by 40 cycles at 95 • C for 5 s and 60 • C for 20 s [24]. Template-free controls were included in each run for each primer pair. β-actin was the internal reference. Degenerate primer sequences for the qRT-PCR are shown in Table S2.

Statistical Analysis
The data were subjected to one-way analysis of variance (ANOVA) in SPSS v. 18.0 (SPSS Inc., Chicago, IL, USA). Differences between treatment means were considered significant at p ≤ 0.05. Graphics were created with GraphPad Prism v. 6.0 (GraphPad Software Inc., La Jolla, CA, USA).

Loquat Apical Bud Phytohormone Analyses
The IAA levels were low in the 'Baiyu' and 'Huoju' apical buds throughout bud differentiation. They did not significantly differ from each other (Figure 1a). The ABA level in 'Baiyu' was about tenfold higher than that in 'Huoju' and reached maxima at the sampling dates of 3 August and 10 August (Figure 1b). The ABA:IAA ratios varied widely between cultivars (Figure 1c). In 'Baiyu', ABA:IAA exhibited a 'down-up-down' trend. It was relatively lower at the earlier stage (13-27 July, i.e., T0-T2), but sharply rose at sampling time (3 August, T3). Hence, this period might have been the key flower bud differentiation stage (Figure 1c). ABA:IAA did not fluctuate in 'Huoju' throughout the experimental period (Figure 1c). The foregoing results indicated that elevated ABA content and ABA:IAA balance might be required for loquat flower bud differentiation.

Statistical Analysis
The data were subjected to one-way analysis of variance (ANOVA) in SPSS v. 18.0 (SPSS Inc., Chicago, IL, USA). Differences between treatment means were considered significant at p ≤ 0.05. Graphics were created with GraphPad Prism v. 6.0 (GraphPad Software Inc., La Jolla, CA, USA).

Loquat Apical Bud Phytohormone Analyses
The IAA levels were low in the 'Baiyu' and 'Huoju' apical buds throughout bud differentiation. They did not significantly differ from each other (Figure 1a). The ABA level in 'Baiyu' was about tenfold higher than that in 'Huoju' and reached maxima at the sampling dates of 3 August and 10 August (Figure 1b). The ABA:IAA ratios varied widely between cultivars (Figure 1c). In 'Baiyu', ABA:IAA exhibited a 'down-up-down' trend. It was relatively lower at the earlier stage (13-27 July, i.e., T0-T2), but sharply rose at sampling time (3 August, T3). Hence, this period might have been the key flower bud differentiation stage (Figure 1c). ABA:IAA did not fluctuate in 'Huoju' throughout the experimental period (Figure 1c). The foregoing results indicated that elevated ABA content and ABA:IAA balance might be required for loquat flower bud differentiation.

Sequence Analysis and Assembly
To obtain a transcriptome landscape for flower bud differentiation, cDNA samples were extracted from the 'Baiyu' and 'Huoju' apical buds at each sampling date and sequenced on the Illumina HiSeq 4000 platform. Table 1  Next-generation short-read sequences were assembled de novo with Trinity (v. trinityr-naseq_r20131110) into 162,080 unigenes with mean size = 1109.34 bp and N50 length = 1625 bp ( Table 2). Unigene length was in the range of 301-2000 nt, and the total length was 87,192,678 nt (Table 2; Figure S1). Distribution of the unigene lengths are shown in Figure S1.

Sequence Annotation
Several complementary approaches were utilized to annotate the assembled sequences. The sequences were aligned against the NR, Swissport, KEGG, KOG, eggNOG, GO, and Pfam databases using the Diamond program. The threshold E-value was 10 −5 , and 45,107 NR, 31,795 Swissport, 15,890 KEGG, 23,886 KOG, 38,848 eggNOG, 28,705 GO, and 57 Pfam sequences were identified (Table 3). All unigenes were annotated and matched against loquat genome database (Table S3), and distribution of the top ten species was determined by the Pollution_test based on BLASTx and the NT database. The cutoff E-value was 1 × e −10 , and coverage was >80%. About 47.96%, 46.53%, and 2.21% of the unigenes matched sequences from Malus domestica, Pyrus × bretschneideri, and Eriobotrya japonica, respectively ( Figure S2). GO and KEGG analysis categorized the functions of the predicted loquat genes (Tables S4 and S5). The unigenes (28,705) were classified into three main categories, including 'biological process' (24,099), 'cellular component' (25,113), and 'molecular function' (24,690) ( Figure 2). In the 'biological process' category, the unigenes implicated in 'cellular process' and 'metabolic process' predominated. There were also high percentages of unigenes under the categories 'response to stimulus', 'biological regulation', and 'regulation of biological process'. For the 'cellular component' category, many unigenes were categorized under 'cell', 'cell part', and 'organelle'. Under the 'molecular function' category, most of the unigenes were involved in 'binding', 'catalytic activity', and 'transporter activity'. The GO enrichment of whole uingenes in each library is shown in Figures S9-S15. GO and KEGG analysis categorized the functions of the predicted loquat genes (Tables S4 and S5). The unigenes (28,705) were classified into three main categories, including 'biological process' (24,099), 'cellular component' (25,113), and 'molecular function' (24,690) ( Figure 2). In the 'biological process' category, the unigenes implicated in 'cellular process' and 'metabolic process' predominated. There were also high percentages of unigenes under the categories 'response to stimulus', 'biological regulation', and 'regulation of biological process'. For the 'cellular component' category, many unigenes were categorized under 'cell', 'cell part', and 'organelle'. Under the 'molecular function' category, most of the unigenes were involved in 'binding', 'catalytic activity', and 'transporter activity'. The GO enrichment of whole uingenes in each library is shown in Figures S9-S15.  (25,113), and 'molecular function' (24,690). Note: The left x-axis represents the % of a specific class of genes in the main category. The right y-axis represents the number of genes in each category.

Identification of DETs Involved in Flowering Time
We calculated the distributions of the unigenes in the various pathways to identify the DETs controlling flowering time in 'Baiyu' and 'Huoju'. Based on the FPKM values,  Table S7).
Flowering locus T (EjFT, CL3783Contig2), flowering time control protein FY-like (EjFY, CL21741Contig1), flowering locus K homology domain-like (EjFLK, CL9491Contig1), and truncated transcription factor CAULIFLOWER A-like (EjCAL1, Comp12801_c0_seq1_1) were selected to analyze the relative differences in their expression during 'Baiyu' and 'Huoju' flower bud differentiation. A qRT-PCR evaluated the expression profiles obtained from RNA-Seq data ( Figure 5). EjFT, EjFY, EjFLK, and EjCAL1-like were significantly upregulated in the apical buds of 'Baiyu' compared with those of 'Huoju', especially in the early stage of 13-27 July (T0-T2) which was the key period of loquat flower bud differentiation. Hence, these genes play vital roles in flowering time regulation ( Figure 5).
Transcriptome data showed that the EjFY and EjFLK expression levels were higher in 'Baiyu' than 'Huoju' throughout apical bud differentiation ( Figure S4b,c).  Table S7). Flowering locus T (EjFT, CL3783Contig2), flowering time control protein FY-like (EjFY, CL21741Contig1), flowering locus K homology domain-like (EjFLK, CL9491Contig1), and truncated transcription factor CAULIFLOWER A-like (EjCAL1, Comp12801_c0_seq1_1) were selected to analyze the relative differences in their expression during 'Baiyu' and 'Huoju' flower bud differentiation. A qRT-PCR evaluated the expression profiles obtained from RNA-Seq data ( Figure 5). EjFT, EjFY, EjFLK, and EjCAL1-like were significantly upregulated in the apical buds of 'Baiyu' compared with those of 'Huoju', especially in the early stage of 13-27 July (T0-T2) which was the key period of loquat flower bud differentiation. Hence, these genes play vital roles in flowering time regulation ( Figure 5). Transcriptome data showed that the EjFY and EjFLK expression levels were higher in 'Baiyu' than 'Huoju' throughout apical bud differentiation ( Figure S4b,c).

Identification of DETs Involved in ABA Signalling
ABA plays important roles in regulating loquat apical bud differentiation and flowering time. Based on the FPKM values and the GO or KEGG enrichment analyses, we identified DETs involved in ABA signaling. A total of 17 ABA signaling pathways were

Discussion
Early flowering is a useful trait in high-value crops as it may enable early fruit harvest and timely market distribution. However, early-flowering loquat cultivars have low yield because they are exposed to low wintertime temperatures. Extension or delay of flowering can reduce the risk of cold stress injury and increase or stabilize production. Nevertheless, flowering time is a complex biological process influenced by numerous factors. Identifying the environmental and genetic mechanisms regulating loquat flowering time could expedite breeding and enhance fruit production (Lin et al. 1999). In the present study, we performed a comparative transcriptome analysis of differentiating loquat flower buds. To the best of our knowledge, this study is the first to use RNA-Seq to identify candidate genes regulating flowering time in loquat.
Certain flowering genes such as EjFT1/2 [15], EjAP1 [14], EjSOC1/2 [44], EjGI, and EjCO [45] were recently identified in loquat and cloned from it. However, little is known about the mechanisms regulating the transformation from vegetative to flower bud in loquat. To identify flowering time-related genes and clarify their transcription during bud differentiation in loquat, we conducted a comparative transcriptome analysis on early-flowering 'Baiyu' and late-flowering 'Huoju'. We generated 40.85-46.91 M and 46.26-47.56 M pairs of 150-bp clean reads from 'Baiyu' and 'Huoju', respectively (Tables 1 and 2). After functional annotation, 28,842 DETs were identified (Table 3; Figures 2 and 3; Figures S1 and S2). These data revealed the DET transcription characteristics and biological processes involved in loquat flower bud differentiation.
Several flowering loci have been identified in Arabidopsis. Overexpression of FLOW-ERING T (FT) leads to early flowering and loss-of-function in the flowering time repressor FLOWERING LOCUS C (FLC). In this manner, the late-flowering phenotype is eliminated [46]. EjLFY-1 in strawberry is homologous to LEAFY (LFY) in loquat. Liu et al. reported that EjLFY-1 overexpression accelerates strawberry flowering, shortens the time required for flower induction, and maintains the early-flowering trait in the asexual progeny [47]. Our transcriptome data identified certain known flowering genes such as flowering locus T-like (EjFT-like) [15] and flowering time control protein FY-like (EjFY) [47]. They also revealed other genes possibly regulating the early flowering trait such as flowering locus K (EjFLK), EARLY FLOWERING-LIKE protein (EjELF), CAULIFLOWER A-like (EjCAL1), and others ( Figure 4). Our transcriptome data revealed that the expression pattern of unigenes varied significantly between 'Baiyu' and 'Huoju' ( Figure S7) and the qRT-PCR indicated that that the flowering-related genes EjFT, EjCAL1-like, EjFY, and EjFLK were upregulated in 'Baiyu' compared with 'Huoju', especially in the early stage of bud differentiation ( Figure 5 and Figure S4). Thus, early flowering might have been partially mediated by upregulation of these genes in the buds. However, the functional characteristics and molecular mechanisms of these DETs in the control of loquat flowering time remain to be elucidated.
Flowering is a complex biological process integrating multiple endogenous and environmental signals that ensure timely flowering. This process is regulated by phytohormones [48,49]. Several studies demonstrated the importance of interactions among various phytohormones in flower induction [50,51]. In Arabidopsis, ABA is a floral repressor, and flowering is delayed by exogenous ABA application and overexpression of ABSCISIC ACID INSENSITIVE MUTANT 5 (ABI5). The latter is mediated by upregulation of FLC and the ABI5 promoter encoding bZIP-type transcription factors [52,53]. It was previously reported that ABA plays the opposite role in loquat, and its level sharply rises during flower bud differentiation [20]. Flowering time was advanced when 200 ppm ABA was sprayed onto loquat apical buds (unpublished data). In the present study, IAA content was low in 'Baiyu' and 'Huoju' apical buds during floral bud development (Figure 1a). Nevertheless, both the ABA content and the ABA:IAA ratio were significantly higher in 'Baiyu' than 'Huoju' (Figure 1b,c). Hence, the early flowering trait observed in 'Baiyu' may be partially explained by massive ABA accumulation in its buds. The results obtained in the present study supported the theory that in loquat, an increase in the ABA content in the bud is conductive to flower bud formation [20]. We speculated that ABA accumulation during loquat bud differentiation could be promoted by certain environmental conditions such as elevated temperatures and moderate drought. Drought and high temperatures cause water loss, induce stomatal closure, and trigger the drought escape response which, in turn, favors early flowering [54,55]. Flowering induction in the short-day plant Lemna aequinoctialis was associated with high ABA content [56]. Leng et al. reported that ABA and its signaling play important roles in regulating floral development and fruit set in sweet cherry. ABA accumulation increased and peaked with floral development. ABA synthetase PaNCED1 expression was consistent with changes in ABA accumulation [57]. In the present study, most ABA signaling genes such as EjABHs, EjPYLs, and EjABI5like showed higher transcription levels in the buds of 'Baiyu' than in those of 'Huoju' (Figures 6-8, Figures S5 and S6). Thus, ABA signaling has a potential regulatory role in loquat flowering time. However, it is unclear whether this positive effect on the expression of genes regulating flowering time can be directly or indirectly attributed to ABA. Furthermore, it is unknown how ABA signaling genes interact with flowering genes and their regulatory networks in loquat. Other stress-related signaling pathways and factors may determine whether flowering will be early or late [58]. Drought stress promoted flowering under long-day conditions in Arabidopsis [59,60] and Citrus latifolia [61]. Salt stress delayed flowering in Arabidopsis [62][63][64]. Recent research has focused on the potential mechanisms by which stress regulates flowering. Substances produced under stress conditions were considered transmissible flowering stimuli. Salicylic acid was regarded as the most likely compound involved in stress-induced flowering [65]. Nevertheless, the physiological responses and molecular regulatory networks implicated in stress-induced loquat flowering remain obscure. In addition, flowering is influenced by multiple factors such as tree age and nutrition, phytohormones, sunshine duration, ambient temperature, and other environmental conditions. A great number of transcription factors (TFs) have been identified, of which AP2, ARF, ARR-B, C2H2, ERF, MYB family, WORK, bHLH, and bZIP were the predominated subgroups ( Figure S8, Tables S9-S15), these TFs may participate in the loquat flowering pathway. However, functional roles of these DETs and TFs associated with these pathways have yet to be identified.

Conclusions
To the best of our knowledge, this report is the first to provide comprehensive transcriptome and DET profiling data associated with loquat floral bud differentiation. Approximately 28,842 DETs involved in multiple regulatory pathways were significantly differentially regulated between early-flowering and late-flowering loquat cultivars. Fortytwo candidate DETs were identified that might play important roles in flowering time control. Seventeen candidate DETs involved in the ABA signaling pathway were identified. The results of the present study could help clarify the gene transcription and regulatory networks involved in loquat flowering time. They also provide an empirical basis for further study on the progression of bud development and the determination of flowering time in loquat.