Transcriptomic Identification of Floral Transition and Development-Associated Genes in S tyrax japonicus

: Styrax japonicus ( S. japonicus ) is an important flowering tree species in temperate regions, and it is regarded as a nectariferous plant. However, there have been few studies to date analyzing floral development in this species. In order to understand gene expression dynamics during S. japonicus flower development, we; therefore, prepared cDNA libraries from three distinct stages of S. japonicus. Illumina sequencing generated 31,471 differentially expressed unigenes during flower development. We additionally conducted pathway enrichment analyses using the GO and KEGG database in order to assess the functions of genes differentially expressed during different stages of the floral development process, revealing these genes to be associated with pathways including phytohormone signaling, Transcription factor, protein kinase, and circadian rhythms. In total, 4828 TF genes, 8402 protein kinase genes, and 78 DEGs related to hormone pathways were identified in flower development stages. Six genes were selected for confirmation of expression levels using quantitative real-time PCR. The gene expression data presented herein represent the most comprehensive dataset available regarding the flowering of S. japonicus , thus offering a reference for future studies of the flowering of this and other Styracaceae species.


Introduction
Styrax japonicus (S. japonicus) is a deciduous species within the Styracaceae family that is distributed from the Qinling Mountains to the south of the Yellow River in China [1]. S. japonicus is an important flowering tree species in temperate regions, where it regarded as a nectariferous plant. Owing to its ornamental value and practical characteristics, S. japonicus is widely used in city greening efforts. However, there have been relatively few studies conducted to date regarding the basic biology and genetics of this species, with previous molecular biology studies on this plant only focused on the development of simple sequence repeat markers [2].
In higher plants, flower development is a key life cycle stage, especially in ornamental plants, and it is influenced by both endogenous and exogenous factors that modulate gene expression [3]. A great deal of research has been conducted examining the flowering mechanisms of model plants, but comparatively little research has focused on woody flowering plants [4]. Unlike annual plants, perennial woody plants bloom only under favorable conditions, with some meristematic tissues transforming into flower organs, while other meristematic tissues remain vegetative. In temperate regions, most perennial plants flower seasonally, including S. japonicus, which blooms in the spring. In Arabidopsis, there are six key pathways regulating floral development: The photoperiod, autonomous, vernalization, gibberellin (GA), ambient temperature, and age pathways [5]. Several key transcription factors (TFs) regulate floral development, including FLOWERING LOCUST (FT), CONSTANS (CO), SUPPRESSOR OF OVEREXPRESSION OF CO1 (SOC1), and FLOWERING LOCUS C (FLC). These floral integrators trigger floral meristem identity genes, including LEAFY (LFY), APETALA1 (AP1), and CAULIFLOWER (CAL) [6][7][8][9], which then determine whether a given tissue will undergo the transition from vegetative to reproductive growth.
Vernalization is the process whereby plants are exposed to cold for extended periods of time in order to complete their floral transformation [10]. Photoperiod regulation and vernalization are the primary means of preventing early flowering. GAs can promote cell growth, but studies have shown that GA can also induce flowering in many plants [11,12]. In the woody flowering plant Amygdalus persica, GA treatment can inhibit flower bud differentiation [13]. In sugar apples, GA, abscisic acid (ABA), and cytokinin (CK) have all been found to play key roles in flower development [14]. In Agapanthus praecox, GAs can promote early flowering, while IAA can delay this flowering [15]. Studies of aspen trees have found that CO mRNA accumulation peaks at the end of the day, and that CO is able to bind to the FT promoter and to thereby induce its expression, resulting in a circadian pattern of FT expression [16]. Diurnal CO and FT expression have significant effects on the growth and development of certain plants [17]. Indeed, many similar studies have been carried out in species including Cymbidium sinense, Citrus limon, and Vernicia fordii [18][19][20].
In the present study, we used RNA-seq technology to construct cDNA libraries corresponding to the three different flowering stages in S. japonicus, yielding a dataset that was sufficiently comprehensive to allow us to identify all flowering-related genes and major metabolic pathways encoded in the transcriptome of these plants. Overall our annotation of the S. japonicus transcriptome and analysis of differential gene expression profiles in these plants serve as a valuable genomic resource for future studies of the flowering mechanisms of S. japonicus and related species.

Materials and Methods
Samples of S. japonicus were obtained from the Rongcheng Donglin Seedling planting cooperative (37°25′12″ N, 122°17′40″ E). Samples obtained included small buds (1-2 mm) (C1), flower buds (7-9 mm) (C2), and open flowers (C3). Three groups of samples were taken for each of these three periods and from three plants. The samples were snap frozen prior to storage at −80 °C. Paraffin sections preparation and RNA extraction were performed at the same time ( Figure 1). Samples from the C1, C2, and C3 stages were collected on 12 March 2017, 10 April 2017, and 25 April 2017, respectively. Total RNA was extracted from these developing flower tissues using TRIzol (Promega, Beijing, China). Equal quantities of three RNA samples for each stage of floral development were mixed prior to treatment with DNase I (Takara, Dalian, China) for 30 min. Agarose gel electrophoresis was used to assess RNA quality, while a 2100 Bioanalyzer device was used for assessing RNA quality and concentrations. High-quality samples with 260/280 nm values between 1.8 and 2.0 were used for library preparation. Briefly, oligo (dT) was used to separate total RNA, which was then fragmented into −200 nt fragments in fragmentation buffer. These mRNA fragments were then used as templates for first-strand cDNA synthesis using random primers and reverse transcriptase. The second cDNA strand was synthesized from this first cDNA strand. AMPure XP beads were then used to purify double-stranded cDNA, and these cDNA fragments were then subjected to purification, end pairing, and a single adenine residue was added to the end of each fragment prior to Illumina adapter ligation. Following ligation, these fragments were separated via agarose gel electrophoresis, with appropriately sized fragments being collected for PCR amplification prior to sequencing on an Illumina HiSeq 2500 platform, with sequencing being performed by Gene De-novo Co. (Guangzhou, China). The transcriptomic sequencing dataset for S. japonicus produced in the present study is available from the NCBI Short Read Archive (SRA) with the accession number SRR7888585-7888593.
Raw high throughput sequencing data were filtered to yield high-quality clean data. Initial filtering focused on the removal of reads that contained adapters, were of low quality (>20% base mass value < 10), had >5% unknown bases, or were <60 nucleotides long. The remaining sequences were then assembled into contigs using the Trinity application [21], yielding a transcriptome reference database. Clean reads next underwent de novo splicing for transcript assembly, and the longest sequence of a given transcript was considered to be a unique gene ("Unigene"). The BLASTX alignment algorithm was then used to align these unigenes to proteins with an E-value ≤ 1 × 10 −5 [22]. Differentially expressed genes were analyzed based upon fragments per million bases (FPKM) values, with the Benjamini and Hochberg method used for p-value correction (false discovery rate, FDR) [23].
A cuffdiff analysis was used to analyze levels of differentially expressed genes (DEG) in these samples, with levels being used to map DEG tags, after which the transcriptomic data were filtered to remove tags that had only a single copy, tags that were of low quality, tags with unknown sequences, and empty tags corresponding only to adapter sequences. The number of mapped clean reads for each unigene was then determined and normalized based upon FPKM numbers in order to calculate unigene expression [24].
RT-qPCR was then used to identify the relative expression levels of the candidate genes during each floral developmental stage. RT-qPCR was conducted using the SYBR Green (TaKaRa, Beijing, China) fluorescent dye on an ABI 7500 Real-Time PCR System (Applied Biosystems, Foster City, CA, USA). The RT-qPCR primers were designed by Primer Premier 5.0 software. GAPDH was selected as the reference gene. Each sample was analyzed in three technical replicates. Relative expression levels were calculated using the 2 −ΔΔCt method [25].

Differential Gene Expression during Floral Development
The basic transcriptome information of S. japonicus have been reported in a previous study [2]. By comparing unigene transcript levels in samples from these three stages of floral development, we were able to explore the transcriptomic dynamics of this flowering process. We used a false discovery rate (FDR) ≤ 0.05 and an absolute log2 ratio value ≥ 1 as criteria to identify DEGs. Using this approach, a total of 14,763 unigenes were differentially expressed between the C1 and C2 stages (7879 up-and 6884 down-regulated, respectively), while 21,617 unigenes were differentially expressed between the C1 and C3 stages (9967 up-and 11,650 down-regulated, respectively) and 15,401 genes were differentially expressed between the C2 and C3 stages (6569 up-and 8830 down-regulated, respectively) ( Figure 2). There were; thus, the fewest differentially expressed genes between stages C1 and C2, suggesting that of the three analyzed stages of floral development, stages C1 and C2 are the most similar. Gene expression differences in different flowering stages were next compared using a Venn diagram ( Figure 3). Among the 31,471 total differentially expressed unigenes, 17,707 (56.26%) were shared across at least two comparisons, and 2606 (8.3%) were differentially expressed across all three comparisons. The overlap between the C1 vs. C2 and C1 vs. C3 comparisons consisted of 9119 (28.9%) unigenes, while there were 5021 (15.9%) unigenes overlapping between C1 vs. C2 and C2 vs. C3, and 8779 (27.8%) unigenes overlapping between C1 vs. C3 and C2 vs. C3. During these three stages of floral development, 957 unigenes were progressively up-regulated, while 598 were progressively down-regulated. This suggests that these genes may play a key role in regulating the growth and reproduction of this species.

Differential Expression of Transcription Factors and Protein Kinases
Transcription factors (TFs) and protein kinases are known to play major roles in regulating flower development. Therefore, we studied the expression dynamics of TF genes in S. japonicus. In total, 4828 TF genes and 8402 protein kinase genes were identified in flower development stages. The 10 most predominant TF families identified during floral development included the basic helix-loophelix, NAC, B3, MYB-related, bZIP, WRKY, ERF, FAR1, C2H2, and MYB TF families, with 626, 479, 45, 10, 58, 357, 347, 428, 25, and 422 members, respectively, identified via de novo transcriptome assembly ( Figure 5A). Among these, some TF families were significantly up-regulated. For example, the bHLH, MYB, and bZIP families were significantly up-regulated after the C1 stage, while the NAC and C2H2 families were up-regulated after the C2 stage. In contrast, other TF families, such as WRKY, ERF, and FAR, were significantly down-regulated at these same time points. A heat map depicting the trends in TF family gene expression during flower development was constructed using the MeV software (Dana-Farber Cancer Institute and Harvard School of Public Health, Boston, USA). These results suggested that complex transcriptional reprogramming is involved in S. japonicus flower development ( Figure S1).
The most abundant group of protein kinases in these sequenced samples were members of the serine/threonine-protein kinase family (3387, 40%). Other common protein kinases in these plant samples included receptor-like protein kinases (1673, 20%), LRR receptor-like serine/threonineprotein kinases (1093, 13%), and CBL-interacting protein kinases (306, 4%). We also found that CBLinteracting serine/threonine-protein kinases (305, 4%) and Wall-associated receptor kinases (227, 3%), which have been reported to involved in flower development and maturation [27,28], were also highly abundant in our transcriptomic dataset ( Figure 5B). These results suggested that complex transcriptional reprogramming is involved in S. japonicus floral growth and development.

Transcriptome Analysis of Phytohormone Pathway Genes during Flower Development
Several phytohormones, including auxin (AUX), ABA, GA, and cytokinin (CTK), are closely associated with flower transition and development. In order to investigate the relationship between phytohormones and flower development in S. japonicus, changes in the expression of unigenes involved in phytohormone signaling pathways were next analyzed ( Figure 6). Hormones including AUX, ABA, CTK, and GA serve as internal cues that govern petal expansion and flower development. In our transcriptomic dataset, we identified 78 DEGs related to these hormone pathways which were arranged into a heat map ( Figure 6). A total of 24 floral unigenes were annotated as being associated with the ABA signaling pathway, of which 15 were expressed at higher levels in stages C2 and C3. In addition, we identified 22 DEGs that were involved in the AUX signaling pathway, 20 that were involved in the CTK signaling pathway, and 12 that were involved in the GA signaling pathway. These genes exhibited strong and specific expression patterns during flower development. MIF1 is involved in the ABA, GA, CTK, and AUX hormone pathways. Four genes in the ABA signaling pathway were also involved in the AUX signaling pathway (MIF1, DRB1, ETR1, ABA3). There were additionally two overlapping genes between the AUX and GA signaling pathways (MIF1, ETR1), and two overlapping genes between the ABA and CTK signaling pathways (DRB1, MIF1).

RT-qPCR Validation of Selected Flowering-Related Genes
To confirm the reliability of our RNA-Seq data, we additionally performed RT-qPCR assays in order to analyze gene expression in independent samples collected from these S. japonicus flowers during different developmental stages (C1, C2, and C3). Six flowering-related genes were selected for these RT-qPCR analyses. We found that the expression levels of these selected genes were strongly correlated with our RNA-seq results, verifying the reliability of our RNA-Seq data (Figure 7).

Transcription Factors Associated with Floral Development
TFs are key regulators of flower development; therefore, we assessed their expression dynamics during floral development in S. japonicus. MYB genes contain DNA binding domains that can control plant growth and development [29]. In Arabidopsis, stamen maturation is regulated by several MYB genes, with MYB108 and MYB24 acting similarly downstream of MYB21 to facilitate transcriptional events leading to stamen and pollen maturation [30]. Recently, the apple R2R3-MYB transcription factor MdMYB1 was found to be a master regulator of anthocyanin biosynthesis and fruit coloration [31]. C2H2 zinc finger (C2H2-ZF) proteins are a large gene family in plants that participate in various aspects of plant growth and development. C2H2 genes can be expressed in carpel primordium of flower buds and regulate stamen development [32]. Members of the WRKY family modulate many biological processes, including responses to biotic and abiotic stress, dormancy, and development [33]. With respect to the regulation of flower development, WRKY71 accelerates flowering in A. thaliana via its activation of FT and the floral meristem identity gene LEAFY [34]. During S. japonicus floral development, we found certain TFs to be substantially up-regulated. For example, the bHLH, MYB, and bZIP TF families were significantly up-regulated between stages C2 and C3, whereas the NAC and C2H2 family members were expressed reached to the highest in stage C3. In contrast, the WRKY, ERF, and ARF family members were down-regulated.

Time-Associated Gene Expression during S. japonicus Flowering
In woody flowering plants, photoperiod regulation is the most important pathway governing the flowering process, and it consists of three modules: Photoreceptors, a circadian clock, and a flowering-specific signaling pathway that is linked to the activity of this circadian clock [30]. Light signals are received by phytochromes (PHYA, PHYB, PHYC, PHYD, PHYE) and cryptochromes (CRYl, CRY2), creating an internal circadian clock that coordinates plant activity [35,36]. In Dimocarpus longan, the circadian clock can measure changes in day length so as to regulate the expression of FKF1 and GI [37], which work together to drive the up-regulation of the TF GO that in turn drives FT and SOC1 expression. ESD4 further regulates SUMO conjugation, thereby controlling the time of flowering [38]. The synthesis of the protein kinase WNK1 is also regulated by circadian rhythms [39]. The FVE, FPA, and FY genes all affect plant circadian rhythms and function in a manner similar to the FCA genes [40,41]. LD serves primarily to repress FLC expression, while FVE, FPA, and FY can all repress FLC and serve as positive regulators of the integration pathway [42]. These genes are redundant in the context of voluntary floral development. As experiments have shown that FLC insertion mutations do not change flowering time [43], FLC has been determined not to be the master regulator of this process. In the present study, we identified genes homologous to LD, FRI, FPA, and FLC in our S. japonicus transcriptome database. We also detected many unigenes associated with the circadian rhythm process, which is a key regulator of flowering. During the C1 stage, we observed the highest expression of the CK2α, CK2β, COP1, CHE, and FT, genes, whereas the genes encoding PHYA, CRY, PRR3, EARLY FLOWERING 3, and TOC1, were most highly expressed during stage C2, and those encoding PHYB, LHY, PRR9, and GI were most highly expressed during stage C3.

Identification of Hormone Unigenes Related to Flowering Stage
Through KEGG pathway enrichment analyses, we determined that genes differentially expressed during S. japonicus flowering were significantly enriched in hormone metabolism and signal transduction pathways. Auxin plays an important role in flower development, especially in regulating the formation of flower primordium and the transformation of flower organs [44]. During the flowering of S. japonicus, auxin-related unigene levels changed significantly, with several genes linked to ARFs and Aux/IAAs being differentially expressed at different time points. IAA10 is an early mediator of auxin regulation [45]. In S. japonicus, peak IAA10 expression was observed during the C3 stage of flower development. This suggests that IAA10 may play a role in flower maturation. Many ARFs are known to be involved in floral development. For example, in Arabidopsis, ARF3 integrates AGAMOUS and APETALA2 functions during floral meristem determinacy [46]. Two other members of the ARF family, ARF6 and ARF8, have been shown to play some role in controlling vegetative and floral organ growth and development. We found that in S. japonicus, ARF1 was highly expressed primarily during the C2 stage, while ARF19 was highly expressed during the C1 stage. This suggests that ARF19 and ARF6 play different roles in flower development.
Studies have shown that certain auxin transporter genes are important mediators of floral development [47,48]. The regulation of auxin flux by ATPIN1, for example, is known to impact male and female gamete development in Arabidopsis [49]. ATPIN1 expression is controlled by ATMLO4, which encodes a heptahelical protein found on the plasma membrane [50]. In this study, we identified the PIN7 and MLO4 auxin transporter genes. Interestingly, the expression levels of PIN7 increased during flower development in S. japonicus, while MLO4 expression exhibited a trend of first increasing and then decreasing.
GA, ABA, and cytokinin are also key mediators of floral development [51]. GA acts primarily through the DELLA proteins GIBBERELLIC ACID INSENSITIVE (GAI), REPRESSOR of ga1-3 (RGA), RGA-LIKE1 (RGL1), RGL2, and RGL3 [52]. ABA can promote flower bud differentiation in some woody plants [53]. Cytokinin drives the flowering of Arabidopsis via regulating the expression of TSF, which is paralogous to FT [54]. GA2OX family genes, encoding the 2oxoglutarate-dependent dioxygenases that catalyze later steps in the biosynthetic pathway of GA, have been identified in different plant species [55]. Interestingly, the expression of the GA biosynthesis gene GA3OX1 was down-regulated in this study, while that of GA2OX6 was initially down-and then up-regulated during the flowering process in this study, reaching peak expression at stage C3. These differential expression dynamics of these GA2OX family members suggest that diverse GA biosynthesis mechanisms are engaged during floral development in S. japonicus. As further research regarding the mechanisms of flowering is required, trends in hormone gene expression and the specific functions of these genes in S. japonicus flowering remain to be further clarified. Furthermore, the MIF1 and ETR1 regulators are involved in the auxin and ABA signaling pathways. ETR1 is regulated by both GA and auxin signaling, while MIF1 is controlled by ABA and cytokinin signaling.

Conclusions
In conclusion, by generating and sequencing transcriptomic libraries for three different stages of S. japonicus flowering, we were able to identify patterns of differential gene expression during the flowering process. Based on the results of GO and KEGG pathway analyses, we were able to highlight particular patterns of gene expression that were linked to this floral development process. We further focused on the identification of particular TFs, circadian rhythm-associated genes, genes linked with flowering time, and hormone-related genes that were engaged during the flowering process. Some of these key genes were further selected for RT-qPCR-mediated confirmation of their expression levels. Together these identified genes will help to support future efforts to more thoroughly understand the molecular basis for the flowering of S. japonicus.
Supplementary Materials: The following are available online at www.mdpi.com/xxx/s1. Figure S1: Identification and analysis of floral transition and flower development-associated transcription factor genes, Table S1: GO annotation results of differentially expressed genes in three stages, Table S2: KEGG annotation results of differentially expressed genes in three stages.
Author Contributions: K.W. and C.Z. conceived and designed the experiments; W.L. performed the experiments, conducted the data analysis. and Z.X. wrote the manuscript; X.J. revised the manuscript and contributed the plant material. All authors have read and agreed to the published version of the manuscript.