Transcriptomic Insight into Underground Floral Differentiation in Erythronium japonicum

Erythronium japonicum Decne (Liliaceae) flowers in early spring after overwintering. Its sexual reproduction process includes an underground development process of floral organs, but the underlying molecular mechanisms are obscure. The present study is aimed at exploring the transcriptional changes and key genes involved at underground floral developmental stages, including flower primordium differentiation, perianth differentiation, stamen differentiation, and pistil differentiation in E. japonicum. Multistage high-quality transcriptomic data resulted in identifying putative candidate genes for underground floral differentiation in E. japonicum. A total of 174,408 unigenes were identified, 28,508 of which were differentially expressed genes (DEGs) at different floral developmental stages, while only 44 genes were identified with conserved regulation between different stages. Further annotation of DEGs resulted in the identification of 270 DEGs specific to floral differentiation. In addition, ELF3, PHD, cullin 1, SE14, ZSWIM3, GIGNATEA, and SERPIN B were identified as potential candidate genes involved in the regulation of floral differentiation. Besides, we explored transcription factors with differential regulation at different developmental stages and identified bHLH, FAR1, mTERF, MYB-related, NAC, Tify, and WRKY TFs for their potential involvement in the underground floral differentiation process. Together, these results laid the foundation for future molecular works to improve our understanding of the underground floral differentiation process and its genetic regulation in E. japonicum.


Background
Erythronium japonicum Decne (Liliaceae) is a spring ephemeral plant commonly known as the Asian fawn lily [1]. The known geographic origins of E. japonicum are Northeast China, Japan, and Korea [1,2]. E. japonicum produces an eye-appealing florescence with reddish-purple flowers [3]. Its vernal characteristics and florescence in early spring make it a perfect ornamental plant. In ephemeral spring plants, flower buds are usually initiated before dormancy induction and continue during the dormancy period [4]. Many studies have been conducted to understand the life cycle, growth habits, reproduction, morphological distinctions, and environmental dynamics in E. japonicum [2,3,[5][6][7][8]. However, there is an apparent lack of studies concerning the molecular mechanisms underlying the underground floral differentiation in E. japonicum.
The floral structures are originated in the floral primordium; however, the specific differentiation of stamens and pistils governs the further floral development [9]. Therefore, it is important to understand the regulatory pathways underlying floral differentiation. Floral differentiation has been widely studied in many plant species viz., Jatropha curcas [9], Brassica napus [10], Camellia sinensis [11], Populus [12] Dianthus caryophyllus [13], Litsea cubeba [14], Rosa chinensis [15], Lilium [16], and Juglans regia [17]. In Arabidopsis, multiple pathways have been identified responsible for floral differentiation including, gibberellic acid (GA), vernalization pathways, aging pathway, and sugar signaling pathway [18][19][20]. FLOWERING LOCUS T (FT) is the integral component in flowering regulation, and most of the flowering-related pathways converge to FT regulation [11]. Regulation of flowering is a complex mechanism and is generally triggered by environmental variables, i.e., temperature and humidity [21]. Furthermore, an overlap between pathways governing flowering and dormancy has been reported [21]. For instance, FLOWERING LOCUS C (FLC) and FRI-GIDA (FRI) have been reported with reduced expression during vernalization [22,23]. The considerable overlap between flowering and dormancy needs to be explored further to exploit their regulation.
This study investigated the transcriptional changes during floral differentiation in E. japonicum at four developmental stages viz., flower primordium differentiation, perianth differentiation, stamen differentiation, and the pistil differentiation period. Our analysis of underground the floral differentiation in E. japonicum provides an overview of differentially expressed genes and their roles in developing flower organs after overwintering.     (from June 9th to June 17th), and stamen differentiation (from June 15th to June 20th) and the pistil differentiation period (from June 17th to June 30th), according to a previously published report [45]. The sample collection stages have been elaborated in Figure 1. The samples of E. japonicum were collected in three biological replicates from different plants. During the differentiation period, samples were collected from the under-forest plot every ten days. After collecting the floral organ meristem samples, samples were wrapped in aluminum foil and frozen in liquid nitrogen immediately. Later, the samples were stored in a refrigerator at -80°C until further use. During sampling, the phenological phase of each plant was recorded. A total of 12 samples were used for transcriptome sequencing analysis. All samples were obtained from the wild, and no permissions are necessary to collect such samples. The formal identification of the samples was conducted by Prof Rengui Zhao, and novoucher specimens have been deposited.

Methods
2.2. RNA Extraction, Library Preparation, and Sequencing. Transcriptome sequencing was performed by constructing four libraries corresponding randomly collected flower bud samples, each with three replicates. Total RNA was extracted using TRIzol reagent (TaKaRa, China). To access, the quality of extracted RNA contamination and RNA integrity number was checked using 1% agarose gel and Agilent 2100 Bioanalyzer system (Agilent Technologies, CA, USA), respectively. Pair end sequencing libraries were constructed using 3 μg RNA for each sample. Further, libraries were generated using NEBNext® UltraTM RNA Library Prep Kit for Illu-mina® (NEB, USA) following the manufacturer's instructions. Illumina HiSeq platform was utilized for RNS sequencing and was performed by company Novogene (https://en.novogene.com/). Following, the libraries were sequenced by paired-end sequencing on Illumina Hiseq.
Low-quality reads and short sequence reads (<50 bp) were removed using FastQC and in-house Perl scripts. Finally, clean reads were de novo assembled using Trinity v.2.6.6 [46].

Gene Expression Quantification and Differential
Expression Analysis. The mapped reads numbers were calculated using featureCounts v1.5.0-p3 [47]. Then, calculating the expected number of fragments per kilobase of exon model per million reads mapped (FPKM) of each gene based on the length of each gene and reads count mapped to the gene. Differentially expressed genes (DEGs) between samples from different floral developmental stages were identified using the DESeq R package (v1.18.0) [48]. The false discovery rate (FDR) method was used to estimate the p value threshold in   [49]. The DEGs were classified using GO [50] and KEGG enrichment analysis [51]. The annotated DEGs were further screened for their functions related to floral bud differentiation, and their corresponding expression levels at different stages were compared. Furthermore, transcription factors were identified using iTAK by integrating PlnTFDB and PlantTFDB [52]. The principle is to identify TF by hmmscan comparison by using the TF family information.

Gene Expression Validation
Using qRT-PCR. Quantitative real-time PCR (qRT-PCR) was performed for selected genes to verify the transcriptomic data and their corresponding gene expressions and different floral differentiation stages. Tiangen RNAprep Pure Plant kit (Tiangen biotech., Beijing, China) was used to isolate total RNA from samples. Eighteen genes related to floral differentiation and flowering-related pathways were selected, and corresponding primers were designed for qRT-PCR using the Oligo-7 software (Table S1). The primers were synthesized by Sangon Biotech (Shanghai, China). The cDNA was extracted from RNA and used as a template to make the reaction for qRT-PCR by using Takara qPCR kit SYBR Premix Ex TaqTM II (Tli RNaseH Plus). Three biological repeats were used for each qRT-PCR reaction and analysis was performed using 2 −ΔΔCt method [53].

Results
Based on morphological distinction, we divided the floral bud differentiation of E. japonicum into four stages, viz.,  Figure 4: Major GO terms, viz., biological processes, molecular functions, and cellular components associated with the DEGs. The x-axis represents the major GO terms, and y-axis represents the count number for each GO term.       BioMed Research International flower primordium differentiation, perianth differentiation, stamen differentiation, and pistil differentiation. Tissue samples from each stage, in three replicates, were subjected to RNA-sequencing. A separate transcriptome from each stage was subsequently analyzed to identify the molecular regulation of floral differentiation in E. japonicum. Approximately 548 million raw reads were filtered for clean reads (520 million). After filtering for unqualified reads, 78.04 Gb of clean bases were obtained where Q20% was above 97.38%, and Q30% was above 92.72%. The GC contents ranged from 48.35% to 50.11% (Table S2). After de novo assembly, 263,291 transcripts and 174,408 unigenes were identified with a mean length of 561 and 706, respectively. The spliced transcripts were sorted lengthwise, and N50 distribution was estimated to be 727 (Transcripts) and 880 (Unigenes). The estimated sequence length distribution of all unigene and transcripts has been presented in Figure 2 Figure S1). Before proceeding to the comparative analysis of transcriptomes, the individual transcriptome data were analyzed for quantitation. Correlation and principal component analysis (PCA) were performed. All samples showed highly significant correlations (Figure 2(b)), while PCA based on FPKM values depicted uniform distribution of replicates for each sample. However, all four samples were distributed as separate groups, suggesting different transcriptional regulations at each stage (Figure 2(c)).

Transcription Factors Associated with Floral
Differentiation. Understanding the developmental regulatory networks is essential to comprehend the specific developmental process. Transcription factors (TF) play a pivotal role in the developmental process. Therefore, we explored the DEGs and identified 213, 181, 397, 303, 191, and 316 TFs in Az vs. Bz, Az vs. Cz, Az vs. Dz, Bz vs. Cz, Bz vs. Dz, and Cz vs. Dz, respectively (Table S10). AP2/ERF-ERF, bHLH, FAR1, mTERF, MYB-related, NAC, Tify, and WRKY were the most prominent TFs differentially expressed at the different floral differentiation stages. Further annotation of these TFs families identified TFs associated with floral differentiation. Based on corresponding annotation results, 41 TFs were identified (Table S11). The differential expression of these TFs has been presented in Table S10. Twelve TFs,  FPKM values of identified 28,508 were subjected to Kmean clustering analysis to identify coexpressed TF genes with DEGs at the four floral development stages ( Figure 6). We identified subclasses 1, 3, 6, and 8 containing the most number of structural DEGs related to flower development. These subclasses contains several TFs, including AP2/ERF-ERF, bHLH, FAR1, mTERF, MYB-related, NAC, Tify, and WRKY. Further molecular characterization of identified coexpressed TFs with structural genes can potentially narrow down the DEGs involved in floral differentiation in E. japonicum.

qRT-PCR-Based Verification of Expression Pattern of
Identified Genes. Based on the transcriptome analysis and further bioinformatics analyses, we identify 18 genes potentially associated with floral differentiation in E. japonicum. To validate the transcriptome data and corresponding expression of selected genes at different floral development stages, we performed qRT-PCR-based validation. As a result, the expression profile of selected genes confirmed the transcriptome's reliability and demonstrated the differential expression pattern at the four floral developmental stages (Figure 7).

Discussion
Floral organ development in ornamental plants is a key process shaping their commercial value [47,56]. To meet the everincreasing demand in floriculture industry, many wild flowers have been domesticated for their commercial use [57]. Erythronium japonicum Decne (Liliaceae) is one such example that is native to Asia [2]. E. japonicum is an early spring ephemeral, and the initial flower development phase starts underground without photoperiod induction and vernalization [58]. Therefore, it is important to understand and explore the floral differentiation process in E. japonicum to better utilize its commercial value. The present study is aimed at exploring diverse transcriptomic landscape pertaining to different developmental stages involved in floral morphogenesis, viz., including primordium differentiation, perianth differentiation, stamen differentiation, and pistil differentiation stage.
Floral induction is a series of developmental processes, with each stage providing substantial inputs to govern the overall process [59]. Based on a previous report [45], we selected four stages of underground flower development in E. japonicum, and samples from each stage were subjected to transcriptomic profiling. A similar approach for exploring transcriptomic profile for floral differentiation has been adopted in different species such as Ranunculus glacialis [60], Chrysanthemum morifolium [61], Chrysanthemum lavandulifolium [62], Staphisagria Ranunculaceae [63], rice [64], and Delphinieae [63]. The obtained results in this study suggested significant variation in the expression profiles at the different stages.

BioMed Research International
We identified putative genes responsible for floral differentiation in E. japonicum and observed upregulated expression of ELF3, PHD, cullin1, SE14, ZSWIM3, GIDNATEA, and SERPIN B from the initial primordium differentiation stage to the perianth differentiation stage. ELF3 plays a crucial role in regulating the circadian clock and is responsible for many downstream regulatory pathways [65]. Furthermore, ELF3 interacts with ELF4, LUX, and other proteins to regulate hypocotyl extension, thermo-morphogenesis, and flowering time [66,67]. In addition, ELF3 gene has been reported to suppress cell elongation under increasing temperature [65]. Thus, further functional analysis of ELF3 can provide valuable insights into its role in regulating floral differentiation in E. japonicum. Similarly, other DEGs identified with differential expression at early stages of floral differentiation, including PHD [68][69][70], cullin 1 [71], SE14 [72], ZSWIM3 [73], GIGNATEA [74], and SERPIN B [75], have been characterized for their positive role in the regulation of flowering in different plant species with their potential involvement in circadian pathways.
Comparative transcriptomic profile suggested that expression of ELF3 and FT, cullin 1, GLP1, and CONSTANS significantly increased at later stages, suggesting an enhanced role in later flower differentiation stages, viz., stamen differentiation, and pistil differentiation. Similarly, stage-specific genes were identified for each transition stage during floral organ development. Similar results have been reported previously, suggesting stage-specific regulatory genes [76,77]. MADS domain protein APETALA1 (AP1) and LEAFY (LFY) are generally considered as master regulators of flowering in plants [78]; however, activation of these genes is dependent on MADS domain proteins, including SUPPRESSOR OF OVEREXPRESSION OF CONSTANS 1 (SOC1), FRUITFULL (FUL), and AGAMOUS-LIKE 24 (AGL24) [79,80].
TFs and their roles in the developmental process have been extensively studied over the past decades [81][82][83]. Several TFs such as AP2/ERF, MYB, bHLH, MADS-box, and NAC have been previously characterized for their active role in development and flower initiation in plants. Our study identified TFs such as AP2/ERF [84], bHLH [85], FAR1 [86], mTERF [87], MYB-related [88], NAC [89], Tify [90], and WRKY [91] as major regulators involved in floral differentiation in E. japonicum. Ethylene-Responsive Factor (ERF) gene family is known for its diverse role in plant developmental process, including germination, flowering, maturation, and senescence [92,93]. Ethylene has been reported with its regulatory role in the transition to flowering phase [94,95], and ethylene regulation in flowering plants is controlled by ERF gene family [96]. The bHLH TF family regulates CONSTANS in Arabidopsis, which is crucial for photoperiodic flowering. Similarly, FAR1, an important regulator in the photo-sensitive circadian clock, regulates ELF4 by directly binding to FBS cis-elements and promotes flowering [97]. However, the flowering phase in E. japonicum starts underground in the absence of light. Therefore, further characterization of bHLH and FAR1 and their relationship with CONSTANS in E. japonicum may yield a potential breakthrough in activating flower organs of underground bulbs in ephemeral plants. Myb-related protein positively regulates flowering by activating FLOWERING LOCUS T and FLOWERING LOCUS T INTERACTING PROTEIN 1 [88]. Furthermore, based on GO terms associations, we identified twelve upregulated TFs at the four stages of flower initiation. Therefore, we speculated that these TFs play a crucial role in underground floral differentiation in E. japonicum.

Conclusions
This study investigated the transcriptional profiles of underground floral differentiation in E. japonicum at four developmental stages. Through a comparative transcriptome analysis, we identified several putative candidate genes, including ELF3, PHD, cullin 1, SE14, ZSWIM3, GIGNATEA, SERPIN B, bHLH, FAR1, mTERF, MYB-related, NAC, Tify, and WRKY. Further functional characterization of these putative candidate genes can provide a better understanding of the process of underground floral organ differentiation in E. japonicum. Furthermore, the excavated information can be used as a base study for further characterization of floral differentiation in spring ephemeral plants.

DEG:
Differentially expressed gene FPKM: Fragments per kilobase of exon model per million reads mapped FRD: False discovery rate GO: Gene ontology KEGG: Kyoto Encyclopedia of Genes and Genomes PCA: Principal component analysis qRT-PCR: Quantitative real-time PCR TF: Transcription factor.

Data Availability
The raw RNA-seq data has been submitted to NCBI SRA under the project number: PRJNA730644.

Conflicts of Interest
The authors declare no conflict of interest.

Authors' Contributions
Conceptualization was done by H.W., R.Z., and J.Z.; methodology was done by H.W.; software was done by H.W. and L.Z.; validation was done by H.W., P.S., and X.L.; formal analysis was done by H.W.; investigation was done by H.W., L.Z., P.S., and X.L.; resources was done by H.W., L.Z., P.S., and X.L.; data curation was done by H.W.; writing-original draft preparation was done by X.X.; writing-review and editing was done by R.Z. and J.Z.; supervision was done by R.Z. and J.Z.; project administration was done by R.Z. and J.Z. All authors have read and agreed to the published version of the manuscript.

BioMed Research International
Supplementary Materials Figure S1: annotation of unigenes from different databases; Table S1: list of primers corresponding to selected DEGs used for qRT-PCR; Table S2: statistics of transcriptome sequencing; Table S3: list of differentially expressed genes corresponding to different floral developmental stages; Table  S4: differentially expressed genes related to floral differentiation identified in comparison of Az (flower primordium differentiation) and Bz (perianth differentiation); Table S5: differentially expressed genes identified related to floral differentiation in comparison of Az (flower primordium differentiation) and Cz (stamen differentiation); Table S6: differentially expressed genes identified related to floral differentiation in comparison of Az (flower primordium differentiation) and Dz (pistil differentiation); Table S7: differentially expressed genes related to floral differentiation identified in comparison of Bz (perianth differentiation) and Cz (stamen differentiation); Table S8: differentially expressed genes related to floral differentiation identified in comparison of Bz (perianth differentiation) and Dz (pistil differentiation); Table S9: differentially expressed genes related to floral differentiation identified in comparison of Cz (stamen differentiation) and Dz (pistil differentiation); Table S10: differentially expressed transcription factors between different groups;