De Novo Transcriptome Analysis Reveals Flowering-Related Genes That Potentially Contribute to Flowering-Time Control in the Japanese Cultivated Gentian Gentiana triflora

Japanese cultivated gentians are perennial plants that flower in early summer to late autumn in Japan, depending on the cultivar. Several flowering-related genes, including GtFT1 and GtTFL1, are known to be involved in regulating flowering time, but many such genes remain unidentified. In this study, we obtained transcriptome profiling data using the Gentiana triflora cultivar ‘Maciry’, which typically flowers in late July. We conducted deep RNA sequencing analysis using gentian plants grown under natural field conditions for three months before flowering. To investigate diurnal changes, the plants were sampled at 4 h intervals over 24 h. Using these transcriptome data, we determined the expression profiles of leaves based on homology searches against the Flowering-Interactive Database of Arabidopsis. In particular, we focused on transcription factor genes, belonging to the BBX and MADS-box families, and analyzed their developmental and diurnal variation. The expression levels of representative BBX genes were also analyzed under long- and short-day conditions using in-vitro-grown seedlings, and the expression patterns of some BBX genes differed. Clustering analysis revealed that the transcription factor genes were coexpressed with GtFT1. Overall, these expression profiles will facilitate further analysis of the molecular mechanisms underlying the control of flowering time in gentians.


Introduction
Japanese cultivated gentians (Gentiana triflora, Gentiana scabra, and their hybrids) are the flowers most frequently used to decorate graves on special occasions, such as the Bon Festival, 'Obon', in mid-August and Equinox week, 'Ohigan,' in late September [1]. Therefore, the flowering time of gentians is an important factor for farmers seeking to ship the flowers to market in a timely manner. Except for a minority of potted plants grown in greenhouses, gentians are usually cultivated in the field; thus, controlling their growth and flowering time is challenging. Several pre-and post-harvest studies have been conducted to assess the quality and longevity of cut gentian flowers [2][3][4], and the results depend on cultivar and species; therefore, a universal method for extending gentian flower life has not been established. In contrast to chrysanthemums, which are short-day plants for which flowering time can be controlled using photoperiodic lighting (e.g., by day-length extension or night breaks) [5], Japanese cultivated gentians are considered day-neutral plants [6], like tomato and cucumber [7,8]. Therefore, managing gentian flowering time using light is relatively difficult, and most farmers instead cultivate several cultivars with different flowering times, thereby accounting for the year-by-year fluctuations in flowering time in the field. However, this strategy is wasteful, and gentian flower farmers would genes involved in the biosynthesis of active ingredients [29,30]. We also applied RNAseq analysis to Japanese gentians in relation to flower color and flower opening [31][32][33]. However, RNA-seq analysis has yet to be targeted at flowering in Gentiana species.
Thus, the aim of the present study was to produce gene catalogs and expression profiles of the developmental and diurnal changes in flowering-related genes in a gentian species. To this end, field-grown G. triflora plants were the subjects of RNA-seq analysis for three months prior to flowering. To evaluate the usefulness of the RNA-seq-based gene catalog, we first analyzed B-BOX (BBX) family genes, including the CO/COL family, and MADS-box genes, as this gene family is known to regulate flowering in other plant species. Next, using in-vitro-grown seedlings, we tested the photoperiodic response in selected BBX genes, confirming that some of these genes have photoperiodic responses under different conditions. Subsequently, we identified the transcription factor genes coexpressed with GtFT1 via clustering analysis, finding genes that were likely related to flowering and plant growth. Overall, this study revealed candidate genes that may be involved in the regulation of flowering time in gentians. These data will facilitate further research on flowering in gentians at the molecular level and improve our understanding of the molecular mechanisms underlying the control of flowering time in these plants, which will enable the development of molecular markers and breeding of suitable cultivars for flowering control.

Identification of Reference Genes Expressed in Gentian Leaves, and Gene Annotation via BLAST Analysis and Protein Motif Analysis
The leaves of field-grown gentian plants were sampled as shown in Figure 1A. A list of the bulks used for RNA-seq is shown in Supplementary Table S1. The filtered RNAseq reads were assembled de novo using Trinity [34], with 521,292 transcripts detected. The coding regions of the assembled transcripts were predicted using TransDecoder [35]. Subsequently, Cluster Database at High Identity with Tolerance (CD-HIT) [36] was used to reduce transcript redundancy and obtain unique genes. Finally, 37,919 transcripts were obtained. BUSCO analysis was performed to evaluate the completeness of the assembly against a dataset set of 1614 core genes in Embryophyta [37]; the completeness of the transcript set was 87.3%, which is comparable to 89.0% in our previous RNA-seq analysis of corolla [33]. RNA-seq results are summarized in Table 1. The transcripts were annotated using a local BLASTX search for Arabidopsis thaliana (downloaded from ftp.ensemblgenomes. org/pub/plants/release-42/fasta/arabidopsis_thaliana/pep/ (accessed on 1 September 2022)) and the UniProt (downloaded from ftp.ebi.ac.uk/pub/databases/uniprot/current_ release/knowledgebase/complete/ (accessed on 1 September 2022)) protein database. Gene expression levels were calculated from the transcripts per kilobase million (TPM) values using the RSEM software within the Trinity package. The expression levels of all contigs and the annotation results are shown in Supplementary Table S2.   We searched contigs known to be flowering-related genes in Arabidop Flowering-Interactive Database (FLOR-ID; http://www.phytosystems.ulg  (3-year-old) plants were used. The third to fourth fully opened leaves were sampled at seven time points and subjected to RNA-seq analysis. (B) Aseptic seedlings cultured in plant boxes under long-and short-day conditions were used. The third and fourth leaves from the top were sampled at seven time points and subjected to qRT-PCR analysis. White and black boxes indicate light and dark periods, respectively.

Expression Profiling of Flowering-Related Genes in Leaves over a Three-Month Period Prior to Flowering
We searched contigs known to be flowering-related genes in Arabidopsis using the Flowering-Interactive Database (FLOR-ID; http://www.phytosystems.ulg.ac.be/florid/ (accessed on 10 March 2021)) [38]. Among the genes present in FLOR-ID, our RNA-seq contigs had 212 hits (65.2%; Tables 2 and S3). This gene catalog covered 44.2-74.1% of FLOR-ID genes depending on the pathway. First, we analyzed the gene expression patterns of BBX and MADS-box family genes. A list of these genes is provided in Table S4. Seventeen BBX and thirteen MADS-box genes were included, and phylogenetic analysis using Arabidopsis genes indicated that gentian genes were distributed in most clades (Figure 2A,B). The conserved motifs in the BBX genes are shown in Figure S1; the assembled contigs were full length, except for GtCOL12. The expression profiles of the gentian BBX and MADS-box genes are shown in Figures 3 and 4, respectively. We performed real-time reverse transcription-PCR (qRT-PCR) to validate the expression levels (TPM) obtained via RNA-seq analysis. For this purpose, we searched for internal reference genes for the standardization of gene expression because the expression of the ubiquitin gene used previously fluctuated in samples collected on 14 July ( Figure S2). To select these reference genes, we used the following criteria: TPM values at all points were >100, and the value of each point was 0.67-1.50 times the average of all points. Using these criteria, three genes were found to be stably expressed in all samples: MALE GAMETOPHYTE DEFECTIVE 1 (GtMGP1), 20S PROTEASOME ALPHA SUBUNIT PAD1 (GtPAD1), and GLYCERALDEHYDE-3-PHOSPHATE DEHYDROGENASE C SUBUNIT 1 (GtGAPC1) ( Figure S3). Because the variation in GtMGP1 gene expression was the lowest among these genes, we selected GtMGP1 as an internal standard. We subjected GtFT1 and four BBX genes near to Arabidopsis CO, i.e., GtCO, GtCOL, GtCOL4, and GtCOL5, to qRT-PCR, finding that the expression profile patterns were similar according to RNA-seq ( Figure 3) and qRT-PCR ( Figure S2) analyses. Therefore, the TPM values obtained via RNA-seq were considered reliable and suitable for detailed expression profiling analysis.
GtFT1 expression increased from May to June, but clear variation in expression within one day was not observed. In contrast, GtFT1 expression levels increased markedly during the day in July (Figures 3 and S2). GtFT2, involved in the phase transition of overwintering buds, was not detected in our RNA-seq analysis. Seventeen BBX genes showed genedependent developmental and diurnal changes, but no genes showed similar expression patterns to those of GtFT1 ( Figure 3). The expression patterns of 13 MADS-box genes are shown in Figure 4. The expression profiles were also gene-dependent, but the expression patterns of SEPALLATA1 (GtSEP1) and APETALA1 (GtAP1) were similar to that of GtFT1, whereas the expression patterns of GtSVP-L1 and GtFLC-L1 differed from those of GtFT1. The expression levels of GtFUL increased slightly during development, whereas those of GtSOC1b decreased.

Effect of Photoperiod on the Expression of BBX Genes in In Vitro-Grown Seedlings
We further investigated gentian BBX genes by analyzing the effects of different day lengths on their expression levels. It was difficult to control the day-length conditions in the field-grown gentians; therefore, we used three-month-old in-vitro-grown gentian seedlings ( Figure 1B). Flowering was not observed in these in-vitro-grown gentians ( Figure S4). As expected, qRT-PCR analysis revealed that GtFT1 expression was not induced under either day-length condition, i.e., either 16 h light/8 h dark or 8 h light/16 h dark ( Figure 5). The expression profiles of representative BBX genes, including GtCO, GtCOL, GtCOL4, and GtCOL5, over 24 h are shown in Figure 5. GtCO, GtCOL, and GtCOL5 showed different expression patterns (repressed in the light) under long-and short-day conditions, whereas the expression patterns of GtCOL4 barely changed over 24 h. Thus, some BBX genes exhibited altered expression levels in response to day length, even when GtFT1 expression was not induced.     GtFT1 expression increased from May to June, but clear variation in expression within one day was not observed. In contrast, GtFT1 expression levels increased markedly during the day in July (Figures 3 and S2). GtFT2, involved in the phase transition of overwintering buds, was not detected in our RNA-seq analysis. Seventeen BBX genes showed gene-dependent developmental and diurnal changes, but no genes showed similar expression patterns to those of GtFT1 ( Figure 3). The expression patterns of 13 MADS-box genes are shown in Figure 4. The expression profiles were also gene-dependent, but the expression patterns of SEPALLATA1 (GtSEP1) and APETALA1 (GtAP1) were similar to that of GtFT1, whereas the expression patterns of GtSVP-L1 and GtFLC-L1 differed from

Coexpression Cluster Analysis
Because GtFT1 expression levels in July showed diurnal changes compared with those in May and June (Figures 3 and S3), we searched for the genes coexpressed with GtFT1 in July using coexpression analysis. To this end, transcription factor genes were selected among all contigs, and only genes with specific expression levels (TPM values of ≥1 at any point) were used. The list of genes used is shown in Table S5, and a heatmap including all 1095 contigs used is shown in Figure S5. The genes belonging to the same clade as GtFT1 are shown in Table 3 and in closeup in Figure 6A. The list includes MADS-box genes, GtSEP1, GtAP3a, and GtAP3b but does not include any BBX genes. The clade included a CRYPTOCHROME-INTERACTING basic helix-loop-helix 1 (CIB1) homolog, known to be a positive regulator of FT expression. Notably, several hormone signaling-related genes were also included in the gene list. The expression profiles of phytohormone-related genes found in Table 2 are also shown in Figure S7.
We further investigated gentian BBX genes by analyzing the effects of different day lengths on their expression levels. It was difficult to control the day-length conditions in the field-grown gentians; therefore, we used three-month-old in-vitro-grown gentian seedlings ( Figure 1B). Flowering was not observed in these in-vitro-grown gentians (Figure S4). As expected, qRT-PCR analysis revealed that GtFT1 expression was not induced under either day-length condition, i.e., either 16 h light/8 h dark or 8 h light/16 h dark ( Figure 5). The expression profiles of representative BBX genes, including GtCO, GtCOL, GtCOL4, and GtCOL5, over 24 h are shown in Figure 5. GtCO, GtCOL, and GtCOL5 showed different expression patterns (repressed in the light) under long-and short-day conditions, whereas the expression patterns of GtCOL4 barely changed over 24 h. Thus, some BBX genes exhibited altered expression levels in response to day length, even when GtFT1 expression was not induced.  Figure 1B. The Y-axis shows relative expression levels, which are shown with the short day, ZT0, value as 1. The X-axis shows time points. White and black boxes indicate light and dark periods, respectively.

Coexpression Cluster Analysis
Because GtFT1 expression levels in July showed diurnal changes compared with those in May and June (Figures 3 and S3), we searched for the genes coexpressed with GtFT1 in July using coexpression analysis. To this end, transcription factor genes were selected among all contigs, and only genes with specific expression levels (TPM values of ≥1 at any point) were used. The list of genes used is shown in Table S5, and a heatmap including all 1095 contigs used is shown in Figure S5. The genes belonging to the same clade as GtFT1 are shown in Table 3 and in closeup in Figure 6A. The list includes MADSbox genes, GtSEP1, GtAP3a, and GtAP3b but does not include any BBX genes. The clade included a CRYPTOCHROME-INTERACTING basic helix-loop-helix 1 (CIB1) homolog, known to be a positive regulator of FT expression. Notably, several hormone signaling-  Figure 1B. The Y-axis shows relative expression levels, which are shown with the short day, ZT0, value as 1. The X-axis shows time points. White and black boxes indicate light and dark periods, respectively.

Discussion
We performed transcriptome analysis of field-grown Japanese cultivated gentians using the typical cultivar 'Maciry,' with data collected over three months prior to flowering and over 24 h at 4 h intervals (Table S1). De novo assembly resulted in 37,919 contigs which included 212 FLOR-ID homologous genes (Table 2). Gene expression levels prior to floral initiation and with diurnal changes were obtained (Table S2), and the reliability of the results was validated by qRT-PCR analysis ( Figure S5). The obtained data improve our understanding of the molecular mechanisms underlying flowering in gentians; in particular, we obtained analysis results for the BBX and MADS-box gene families, which are well-known for their roles in the flower development process. To identify the genes coexpressed with GtFT1 during development over three months, clustering analysis of transcription factor genes was performed using data from the threemonth period. A heatmap of all genes is shown in Figure S6, and a list of all identified genes is shown in Table S6. Genes belonging to the GtFT1 clade are shown in Table 4 and in closeup in Figure 6B. The clade included GtBBX22, a homolog of AtBBX22 involved in light and phytohormone pathways, as well as a MADS-box gene, GtSEP1, which was also observed in the coexpression analysis of the July data.

Discussion
We performed transcriptome analysis of field-grown Japanese cultivated gentians using the typical cultivar 'Maciry,' with data collected over three months prior to flowering and over 24 h at 4 h intervals (Table S1). De novo assembly resulted in 37,919 contigs which included 212 FLOR-ID homologous genes (Table 2). Gene expression levels prior to floral initiation and with diurnal changes were obtained (Table S2), and the reliability of the results was validated by qRT-PCR analysis ( Figure S5). The obtained data improve our understanding of the molecular mechanisms underlying flowering in gentians; in particular, we obtained analysis results for the BBX and MADS-box gene families, which are well-known for their roles in the flower development process.
We identified 17 BBX genes belonging to group I-V (Figures 2A and S5). In other plant species, ca. 30 BBX genes are usually identified as B-box protein family genes (e.g., 32 in Arabidopsis, 29 in tomato, and 30 in rice) [39]. Therefore, more genes are likely to be present in gentians. The unidentified genes might be expressed in other organs or at other times, such as under stress-induced conditions, e.g., temperature stress, water stress, and pathogen attack. However, the CO/COL genes expressed in leaves are generally important for photoperiodic flowering control [40], and expression profiling data are useful for identifying the candidate genes involved in flowering control. The expression of some BBX genes was induced on 14 July (about 2 weeks prior to flowering) with the expression of GtFT1, but the expression patterns of the BBX genes did not match that of GtFT1 (Figure 3). Notably, GtCO and GtCOL, homologs of Arabidopsis CO and COL, showed similar expression patterns over three months, i.e., the expression levels decreased under light and increased under dark in the three-month period. Although the expression patterns of GtCOL5 did not change over three months, the expression of this gene was induced in the morning when the induction of GtFT1 expression began (07:30), suggesting that GtCOL5 rather than GtCO and GtCOL contributes to the induction of GtFT1 expression. It is possible that repressor genes were expressed in both May and June to suppress GtFT1 induction or that other cooperative genes are lacking, even when GtCOL5 is expressed. It may also be possible that the protein stability of GtCOL5 is involved in the activation of GtFT1 transcription as reported in Arabidopsis CO [41].
Because it was difficult to subject field-grown gentians to day-length treatments, we analyzed the photoperiod responses of in-vitro-grown seedlings ( Figure 1B). We found that three BBX genes, namely GtCO, GtCOL, and GtCOL5, were expressed at different levels in response to day length ( Figure 5). As expected, GtFT1 expression was not fully induced at the young seedling stage as the plants were not ready for flowering; nevertheless, the photoperiodic response of the BBX genes suggested their involvement in flowering control in response to different day lengths. However, the gene expression levels of these genes in field-grown gentians did not change over three months, reflecting the similar day lengths on 15 May (14 h 20 m), 16 June (14 h 56 m), and 14 July (14 h 41 m). Therefore, another study is necessary to reveal the responsiveness of these genes to photoperiod in mature fieldgrown gentians, which could be achieved using artificial shading or light supplementation. In a previous study, SlCOL, SlCOL4a, and SlCOL4b, among the 13 members in the CO/COL family, were identified as potential candidate activators through their interactions with SFT in day-neutral tomato plants [24]. These tomato genes had different diurnal rhythms under long-and short-day conditions, and they belonged to group I, like Arabidopsis CO, COL, and COL5, which function as flowering inducers. Whether gentian CO/COL family genes also regulate flowering in response to photoperiodic conditions remains unknown but warrants further study.
Thirteen MADS-box genes were identified via RNA-seq ( Figure 3B and Table S4). In previous studies focusing on the identification of genes that regulate floral organ identity in G. scabra [12,13], 14 MADS-box genes belonging to A-E classes were identified, some of which were expressed at low levels in leaves. In our RNA-seq, 7 of these 14 genes were identified, indicating that these genes were expressed to some extent in leaves. The expression patterns of GtAP1, GtAP3a, and GtSEP1 (closest to AtSEP4) were found to be similar to that of GtFT1 (Figure 4). Arabidopsis AP1, AP3, and SEP genes expressed in the floral meristem are known as identity genes that specify floral organs and meristem identity [42,43]. GtAP1, GtAP3a, and GtSEP1 are also expressed in floral organs as A, B, and E class genes, respectively [12]. Thus, it is likely that these homolog genes play the same roles in gentians. However, their expression was not limited to flowers, and the TPM value of GtSEP1 was >100 on 14 July (maximum value: 146.73 at 15:30), suggesting the existence of any functional roles in leaves. Arabidopsis SEP1-3 genes are expressed specifically in inflorescence tissues, whereas SEP4 (also known as AGL3) is also expressed in leaves and floral stems [44], although functional analysis of SEP4 in relation to flowering time was not conducted. In Isatis indigotica, the silencing of IiSEP4, a homolog of Arabidopsis SEP4, delayed flowering time, whereas the overexpression of IiSEP4 in Arabidopsis led to early flowering in addition to changes in floral organs [45,46]. Thus, further detailed analysis is necessary to gain insights into the involvement of GtSEP1 in gentian flowering control. Takahashi et al. [15] identified GtAP1, GtFUL, and GtFLC-L in G. triflora, and these genes were among our RNA-seq contigs (Table S4) However, some genes, such as TFL1 and LFY homologs, were not among our RNA-seq contigs (Table S4), indicating that genes expressed predominantly in the shoot apical meristem (SAM) were missing. Thus, to complete the transcriptome analysis in gentians, RNA-seq should be performed for other organs, including flowers, SAMs, and overwintering buds. We are preparing these samples with the aim of obtaining comprehensive transcriptome data for gentians.
According to coexpression cluster analysis of our RNA-seq data, 23 genes were coexpressed with GtFT1 on 14 July over 24 h (Table 3), including several flowering-related genes. CIB1, which is involved in CRY2-dependent regulation of flowering, was one of these genes. CIB1 and CO act together to regulate FT transcription and flowering [47]; thus, the CIB1 homolog in gentians likely has a similar function, although the interacting BBX genes remain to be identified. Floral homeotic genes, e.g., GtAP3a, GtAP3b, and GtSEP1, were also coexpressed, and these may regulate flowering time as well as floral organ identity (as discussed above). Ten genes were detected in our coexpression cluster analysis of data from the three-month period (Table 4), and these genes are considered candidates for the regulation of flowering time as well as plant development in gentians. The gene lists shown in Tables 3 and 4 also included several plant hormone-related genes involved in auxin, abscisic acid, and ethylene signaling. As shown in Figure S7, gene expression levels of gibberellin metabolism enzyme were not remarkably altered during development over three months, although some showed diurnal changes. Instead, the expression of one GA signaling pathway-related gene, GtPIF4 (PHYTOCHROME-INTERACTING FACTOR 4), increased and showed reverse diurnal expression in July compared with May and June. Arabidopsis PIF4 was revealed to bind the FT promoter and activate FT expression in cooperation with CO under warm ambient temperature [48]; therefore, GtPIF4 may also be involved in gentian flowering induction in July when the temperature rises. In contrast, GtSPL3 expression was clearly reduced in July, but Arabidopsis SPL3 was reported to act as the photoperiodic activator signal in the FT-FD module in flowering [49]. It is likely that GtSPL3 may have the opposite function in gentian flowering, although further analyses are necessary. Some other phytohormone-related genes showed developmental and diurnal changes; therefore, these genes are also candidates along with the BBX and MADS-box genes mentioned above.
Efficient genome-editing and viral vectors have been established in Japanese gentians [50,51]; hence, in time, it will be possible to analyze these candidate genes. In addition, association analysis of gentian cultivars/lines with various flowering times would improve our understanding of the involvement of these genes in flowering-time control. Therefore, we are analyzing the DNA polymorphisms and expression levels of the candidate genes to determine their contribution to flowering time in Japanese cultivated gentians.
The first reference gentian genome, that of Gentiana dahurica, was constructed using three technologies (Nanopore, Illumina paired-end, and Hi-C) and has a total size of around 1.5 Gb [52]. In contrast, the Japanese cultivated gentians G. triflora and G. scabra have about 3-fold larger genome sizes, i.e., ca. 5 Gb, calculated according to their nuclear DNA contents [53]. Thus, the genome assembly of Japanese cultivated gentians will be more difficult, even using the latest technologies. Nevertheless, we are currently constructing a draft genome sequence of Japanese gentians because this sequence will be required to achieve a complete understanding of the various gene functions, especially those involved in flowering, in these plants.

Plant Materials
Gentiana triflora 'Maciry,' bred in Iwate prefecture, was used in this study. This cultivar is an F 1 hybrid between two G. triflora breeding lines, and its flowering time is late July to early August depending on the year. The flowering time in 2020 was 30 July. Three-year-old plants were sampled in 2020. On 15 May, 16 June, and 14 July 2020, the fourth to fifth fully unfolded leaves from the top were collected from 10 independent plants for 24 100 mg) were excised, and the mid ribs were cut using scissors and placed into 2 mL screw cup tubes with beads. Five to ten samples were collected at each time point.

RNA-Seq Analysis
Total RNAs were isolated from the leaves using an RBC Total RNA Extraction Kit Mini (Plant) (SciTrove, Tokyo, Japan) according to the manufacturer's instructions. DNase solution (Deoxyribonuclease, Nippon Gene, Tokyo, Japan; Buffer RDD, QIAGEN, Hilden, Germany) was used to conduct on-column digestion of genomic DNA. Construction of the cDNA libraries, sequencing, assembly, and functional annotation were performed as described previously [33]. The RNA-seq data was deposited in the DNA Data Bank of Japan (DDBJ; accession no. DRA014495). A summary of the generated RNA-seq data and the accession numbers deposited in DDBJ is shown in Table S1. During the quality control step, we filtered and discarded reads shorter than 50 bases and those with an average read quality of <20 and trimmed the poly A and adapter sequences using FaQCs version 2.08 [54].

Validation via qRT-PCR Analysis
The total RNAs isolated as described above were used for cDNA synthesis. qRT-PCR was performed using a QuantStudio 3 system (Thermo Fisher Scientific, Waltham, MA, USA) with Luna Universal qPCR Maser Mix (New England Biolabs, Ipswich, MA, USA). Reference genes for field-grown gentians were selected from the RNA-seq data, and the gentian MGP1 gene was used as an internal control. The relative expression of genes was calculated using the 2 −∆∆Ct method. Primers are listed in Table S7.

Phylogenetic Analysis of the BBX and MADS-Box Family Proteins
BBX and MADS-box family proteins from gentians and Arabidopsis were used for phylogenic tree analysis. Protein sequences from Arabidopsis were obtained from the database of the National Center for Biotechnology Information (http://www.ncbi.nlm.nih. gov (accessed on 10 June 2022)). Maximum-likelihood phylogenic trees were constructed using MEGA X software (http://www.megasoftware.net/home (accessed on 1 September 2022)) with the JTT + G + I model and 1000 bootstrap replications.

Analysis of BBX Gene Expression in Response to Photoperiod in In-Vitro-Grown Seedlings
Seeds of 'Maciry' were aseptically sown on a plastic plate containing half-strength Murashige Skoog medium supplemented with 3% sucrose and 0.24% gellan gum at 22 • C for 29 days under short-day (8 h light/16 h dark) conditions (35 µmol m −2 s −1 ). Germinated seedlings were transferred to a plant box (5 shoots per box) and further cultured in a growth chamber for 40 days under the same short-day conditions. The seedlings were then cultured under short-day or long-day (16 h light/8 h dark) conditions (35 µmol m −2 s −1 ). After 32 days of incubation, the third to fourth leaves from the top were sampled at 4 h intervals over 24 h and stored at −80 • C until their use. All cultures were performed in a plant growth chamber (CLE-305, TOMY, Tokyo, Japan). The isolation of total RNAs and qRT-PCR analysis were performed as described above. Primers are listed in Table S7.

Gene Coexpression Analysis in Field-Grown Gentian Plants
To identify the genes coexpressed with GtFT1 during developmental changes and diurnal changes, we performed clustering analysis. Specifically, we targeted transcription factor genes using the Plant Transcription Factor Database ver. 5.0 (http://planttfdb.gaolab.org (accessed on 1 September 2022) [55][56][57][58]. Contigs with TPM values of >1.0 at any point were used. To select genes during the development stage, we used datasets from three time points (03:30, 11:30, and 19:30) on 15 May, 16 June, and 14 July (nine time points in total). For the analysis of diurnal change, seven time points (24 h) on 14 July were used when GtFT1 was highly expressed. Hierarchical cluster analysis was performed based on the TPM value of each gene using Pearson correlation with an average-linkage clustering method via online tools (http://www.heatmapper.ca/expression/ (accessed on 28 June 2022)) [59].

Conclusions
To reveal the molecular mechanisms underlying flowering control, it is a prerequisite to identify the genes expressed during the developmental stages up to flowering as well as their responses to a variety of environmental conditions. In this study, we produced gene catalogs and performed expression profiling of field-grown gentians using RNA-seq. We successfully demonstrated the usefulness of our gene catalog by analyzing the expression profiles of BBX and MADS-box genes, which involved the screening of candidate genes involved in flowering. Coexpression cluster analysis also revealed several transcription factor genes potentially involved in flowering and growth regulation. Overall, this study represents the first step toward unraveling the molecular mechanisms underlying floweringtime and photoperiodic control in gentians. Flowering control is a highly complex process; therefore, the molecular and functional characterization of each candidate gene is warranted in future studies.  Data Availability Statement: All data supporting the findings of this study are available within the article and within its supplementary data published online. The raw RNA-seq data in this paper were deposited in DDBJ Sequence Read Archive (DRA) under the accession number DRA014495.