Identification and characterization of regulatory pathways involved in early flowering in the new leaves of alfalfa (Medicago sativa L.) by transcriptome analysis

Alfalfa (Medicago sativa L.) is a perennial legume extensively planted throughout the world as a high nutritive value livestock forage. Flowering time is an important agronomic trait that contributes to the production of alfalfa hay and seeds. However, the underlying molecular mechanisms of flowering time regulation in alfalfa are not well understood. In this study, an early-flowering alfalfa genotype 80 and a late-flowering alfalfa genotype 195 were characterized for the flowering phenotype. Our analysis revealed that the lower jasmonate (JA) content in new leaves and the downregulation of JA biosynthetic genes (i.e. lipoxygenase, the 12-oxophytodienoate reductase-like protein, and salicylic acid carboxyl methyltransferase) may play essential roles in the early-flowering phenotype of genotype 80. Further research indicated that genes encode pathogenesis-related proteins [e.g. leucine rich repeat (LRR) family proteins, receptor-like proteins, and toll-interleukin-like receptor (TIR)-nucleotide-binding site (NBS)-LRR class proteins] and members of the signaling receptor kinase family [LRR proteins, kinases domain of unknown function 26 (DUF26) and wheat leucine-rich repeat receptor-like kinase10 (LRK10)-like kinases] are related to early flowering in alfalfa. Additionally, those involved in secondary metabolism (2-oxoglutarate/Fe (II)-dependent dioxygenases and UDP-glycosyltransferase) and the proteasome degradation pathway [really interesting new gene (RING)/U-box superfamily proteins and F-box family proteins] are also related to early flowering in alfalfa. Integrated phenotypical, physiological, and transcriptomic analyses demonstrate that hormone biosynthesis and signaling pathways, pathogenesis-related genes, signaling receptor kinase family genes, secondary metabolism genes, and proteasome degradation pathway genes are responsible for the early flowering phenotype in alfalfa. This will provide new insights into future studies of flowering time in alfalfa and inform genetic improvement strategies for optimizing this important trait.


Background
Alfalfa (Medicago sativa L.) is a perennial leguminous forage extensively planted as livestock feed around the world due to the high crude protein and nutrition contents in leaves [1,2]. While flowering time is an important determinant of harvesting time, which contributes to forage quality, silage, yield, and production in alfalfa [3]. Early flowering can improve the number of cuttings, forage yield, and year-round production of alfalfa; it can also help to avoid the rainy season or drought stress for seed production, and enhance winter hardiness related to the temperate semi-arid continental climate in China [4]. Despite the importance, the underlying molecular mechanisms involved in flowering time control in alfalfa remain elusive.
Past endeavors of flowering time studies discovered six major regulatory pathways in model plants Arabidopsis, rice (Oryza.sativa L.), wheat (Triticum aestivum L.), and barley (Hordeum vulgare L.) [5,6]. These pathways can be categorized into two groups: the photoperiod, vernalization, and ambient pathways that respond to environmental conditions, and the autonomous, gibberellin, and age pathways that respond to endogenous cues associated with the up/downregulation of floral integrator genes, including FLOWERING LOCUS T (FT), FLOWERING LOCUS C (FLC), CONSTANS (CO) and SUPPRESSOR OF OVEREXPRESSION OF CONSTANS1 (SOC1) [7][8][9]. Furthermore, small non-coding RNAs and epigenetic pathways play important roles in regulating flowering time in plants [10]. For instance, in the plant adult phase, miR156 targets SQUAMOSA PROMOTER BINDING-LIKE (SPL) and induces flowering by upregulation expression of flower-promoting genes in the shoot apex [11], and FLOWERING LOCUS D, histone demethylase, downregulates FLC expression associated with mediating H3K4 demethylation and therefore promotes flowering in Arabidopsis [12]. An earlier study also indicated that jasmonate-activated (JA) transcription factors MYC2, MYC3, and MYC4 negatively regulate flowering time by interacting with FT, and exogenous JA could delay flowering time in Arabidopsis [13]. Applying exogenous abscisic acid (ABA) or inducing the overexpression of ABI5, an ABSCISIC ACID-INSENSITIVE MUTANT 5, through FLC expression can also trigger delayed flowering in Arabidopsis [14,15]. In addition, receptor-like kinase (RLK) family proteins play a key role in the perception of internal and external signals involved in flowering time regulation [16]. For instance, FLOR1, a leucine-rich repeat (LRR) protein expressed in the shoot meristem, physically interacts with AGA-MOUS (AG, a MADS-box transcription factor) and is induced in the early inflorescence meristem, and the Arabidopsis flor1 mutant exhibits a delayed flowering phenotype [17]. CURVY1, a novel receptor-like protein kinase, also responds to flowering time and seed production in Arabidopsis [18]. Furthermore, the ubiquitinproteasome system is a major catabolic pathway involved in plant growth, development, and physiological processes, including floral development. Tagging of the target proteins with ubiquitin is mainly performed by E1 (ubiquitin-activating enzyme), E2 (ubiquitin-conjugating enzyme), and E3 (ubiquitin-protein ligase enzyme), followed by degradation mediated by the 26S proteasome [19]. For example, PUB13, a U-box protein possessing E3 ubiquitin ligase, physically interacts with HFR1 (long hypocotyl in far-red light1) to regulate photomorphogenesis and flowering time in Arabidopsis [20]. However, whether their alfalfa homologs participate in flowering time regulation is not clear.
Previous studies using the shoot tip or leaf samples of alfalfa plants of different flowering stages have indicated that proteins involved in carbohydrate metabolism, phenylpropanoid biosynthesis and immunity are responsive to flowering time [4,21]. However, little is known about the underlying molecular mechanisms associated with flowering time control in new leaves, which represent the transition from the vegetative to the floral phase in alfalfa. It has been shown that leaf development and floral organ initiation are controlled by common homeotic factors [22]. For instance, knocking down the AG gene resulted in the reprogramming of leaf development and affected the expression of key regulatory genes (TRIPTYCHON, APETALA3 and PISTILLATA) in Arabidopsis leaves [22,23]. Therefore, we hypothesize that these metabolic genes may play essential roles in new leaves of alfalfa during the early flowering stage. Here, integrated phenotypical, physiological, and transcriptome analyses showed that hormone biosynthesis and signaling pathway genes, pathogenesis-related (PR) genes, signaling receptor kinase family genes, secondary metabolism genes, and proteasome degradation pathway genes might contribute to the early flowering phenotype in alfalfa.

Phenotypes and endogenous hormone contents
For the flowering time assays under controlled growth chamber conditions, plants were trimmed to maintain the same height. Genotype 80 first flowered at 25 days after planting with a mean plant height of 30.95 cm, when the late-flowering genotype 195 had an average plant height of only 19.73 cm (Fig. 1a). Genotype 195 then flowered at 50 days after planting. Endogenous hormone contents in mature leaves (ML), new leaves (NL) including apical meristem, and flower buds (FB) were determined in both genotypes at 25 days after planting when genotype 80 flowered ( Fig. 1b and c; Fig. 2a-d).
We observed significantly higher auxin (IAA), ABA, salicylic acid (SA), and JA contents in NL compared with ML in both genotypes, while, the differences in hormone levels between NL and ML was smaller in genotype 80 than in 195 ( Fig. 2a-d). In addition, except for IAA, hormone contents in FB were significantly higher than those in ML or NL in genotype 80 ( Fig. 2a-d). Together, these results suggest a higher translocation efficiency of hormones from new leaves to flower buds in genotype 80 than genotype 195.

Transcriptome profiling of leaves and flower buds
To further investigate the molecular basis of the earlyflowering phenotype in genotype 80, we performed transcriptome profiling in the ML, NL, and FB of both genotypes (Table 1). A total of 101.13 Gb of clean reads were obtained, the Q30 percentage (the percentage of bases with a quality score of 30 or higher) and GC (guanine and cytosine) percentage were 93.69 and 42.28%, respectively. A total of 53,897 unigenes were yielded from the assembly (N50 of 1770 bp, mean length of 1106 bp) and were annotated using seven functional databases, including the NCBI non-redundant protein database (NR), Nucleotide sequence database (NT), Clusters of Orthologous Groups of proteins (COG), Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG), Swissprot and Interpro (E-value < 10 − 5 ). Then, 31,962 coding sequence transcripts (N50 of 1866 bp, mean length of 1336 bp) were used for further functional analysis (Table 1). We identified 2148 (1454 upregulated/694 downregulated) and 1972 (1142 upregulated/ 830 downregulated) differentially expressed genes (DEGs) between NL and ML in 195 and genotype 80 respectively, among which 942 and 766 DEGs were exclusively expressed in genotype 195 and 80, respectively ( Fig. 3a and b; Additional file 1: Table S1). In addition, 4423 (2300 upregulated/2123 downregulated) DEGs were identified between FB and NL in genotype 80, and 6288 (3615 upregulated/2673 downregulated) DEGs were identified between FB and ML in genotype 80, among which 2802 common DEGs were found to respond to early flowering in genotype 80 ( Fig. 3a and b; Additional file 1: Table S1).

Response of hormone-related genes to early flowering
The significant change in hormone levels prompted us to further examine the transcript levels of hormonerelated genes. The analysis detected 27 and 34 DEGs (between NL and ML) involved in the IAA, ABA, SA and JA metabolism and signaling pathways in genotype 195 (NL195/ML195) and genotype 80 (NL80/ML80), respectively (Figs. 4a, b and 5; Additional file 2: Table S2). Among these DEGs, 11 of the 16 IAA metabolism and  Table S2). Seven of the 13 IAA-related genes were downregulated in NL compared with ML in genotype 195, such as those encoding SAUR family proteins, o-fucosyltransferase, auxin efflux carrier component 5, and thromboxane-A synthase ( Fig. 5; Additional file 2: Table S2). Furthermore, ten DEGs in genotype 195 (five downregulated/five upregulated) and eight DEGs in genotype 80 (three downregulated/five upregulated) were enriched in NL compared with ML ( Fig. 5; Additional file 2: Table S2). These were primarily involved in ABA metabolism and signaling pathways, including the ABA-responsive protein, 9-cisepoxycarotenoid dioxygenase, and HVA22-like protein encoding-genes ( Fig. 5; Additional file 2: Table S2). Two SA-related genes that encode SA carboxyl methyltransferase (SAMT) and lipoxygenase were significantly downregulated in NL compared with ML in genotype 80 but were not enriched in genotype 195 ( Fig. 5; Additional file 2: Table S2). Seven of the nine JA biosynthetic genes were significantly downregulated in NL compared with ML in genotype 80; while two of the four JArelated genes were downregulated in NL compared with ML in genotype 195, including those encoding 12oxophytodienoate reductase-like protein, allene oxide synthase, and lipoxygenase ( Fig. 5; Additional file 2: Table S2). These results further indicate that the downregulation of these hormone-related genes may have a negative impact on early flowering in genotype 80. In  addition, 16 and 27 DEGs related to hormone metabolism and signaling pathways were identified between FB and NL and between FB and ML, respectively, in genotype 80; these include genes involved in ABA signal transduction and genes encoding auxin-responsive proteins, SAMT, allene oxide synthase, and lipoxygenase 3 ( Fig. 5; Additional file 2: Table S2).
Response of genes-related to pathogenesis and metabolism to early flowering   between FB and NL and 78 (50 upregulated and 28 downregulated) DEGs between FB and ML were identified in genotype 80 (Additional file 2: Table S2). Altogether, these findings indicate that compared with genotype 195, more PR genes were downregulated, and more metabolismrelated genes were upregulated in genotype 80 upon flowering.
Responses of receptor-like kinase family genes to early flowering  Table S2). Taken together, these results reveal that more receptor-like kinase genes were downregulated in response to early flowering in genotype 80 than in genotype 195.

Response of proteasome-related genes to early flowering
Thirty-two proteasome-related DEGs between NL and ML (21 up-regulated and 11 down-regulated) were identified in genotype 195, and 25 such DEGs (nine upregulated and 16 downregulated) DEGs were identified genotype 80. These genes mainly encode F-box family proteins, really interesting new gene (RING) /U-box superfamily proteins, and S-ribonuclease binding proteins (Fig. 4a, b and 8; Additional file 2: Table S2). Meanwhile, ten DEGs between FB and NL (four upregulated and six down-regulated) and 32 DEGs between FB and ML (16 upregulated and 16 downregulated) were  Table S2). Together, these results reveal that more proteasome related genes were downregulated in response to early flowering in genotype 80.

qRT-PCR validation of DEGs identified by RNA-seq
We randomly selected 15 DEGs between NL and ML identified by transcriptome profiling from genotypes 195 for qRT-PCT validation. Six of the 15 DEGs showed good agreement between the qRT-PCR analysis and the RNAseq data. These include genes involved in ABA signal transduction, and those encoding beta-amyrin synthase R, a tyrosine kinase family protein, an LRR receptor-like kinase, and a RING/U-box superfamily protein. For genotype 80, the expression of all DEGs between NL and ML exhibited a high correlation between the RNAseq data and qRT-PCR analysis; similarly, six DEGs between FB and NL and 11 DEGs also showed good agreement between the RNA-seq data and qRT-PCR analysis ( Fig. 10; Additional file 3: Table S3).

Discussion
Early flowering of alfalfa is related to hormone metabolism Hormones play distinctive roles in controlling plant growth and development, especially in flowering processes [24,25]. For instance, disruption of auxin biosynthesis, signaling, or polar auxin transport can extensively inhibit the initiation of inflorescence [15,26] whereas exogenous ABA or JA applications can significantly delay floral transition [13,14]. Exogenous SA is known to induce flowering by suppressing FLC expression and increasing FT, CO, and SOC1 expression in tobacco (Nicotiana tabacum L.) [27], Lemna paucicostata [28] and Arabidopsis [29]. In the present study, the ratio of hormone contents (IAA, ABA, SA, and JA) in new leaves compared with mature leaves was lower in the early-flowering genotype 80 compared with the late-flowering genotype 195. Consistent with this observation, 11 auxin-related genes were downregulated in NL compared with ML in phenotype 80, these include those encoding SAUR family proteins, o-fucosyltransferase, auxin efflux carrier component 5, and thromboxane-A synthase. Seven JA-related genes were downregulated in NL compared with ML in genotype 80, including those encoding the lipoxygenase, the 12-oxophytodienoate reductase-like protein, and SAMT. Meanwhile, we observed the upregulation of genes involved in JA biosynthesis in the FB of genotype 80, such as genes encoding SAMT and lipoxygenase3/5. In addition, we detected higher ABA and SA contents as well as higher transcript levels of the HVA22 family protein and SAMT encoding genes in the FB of genotype 80 compared with genotype 195. ABA is known to induce AtHAV22 expression mainly in flower buds and inflorescence tissues, whereas SAMT catalyzes the production of methyl-salicylate from SA, which confers pathogen resistance in Arabidopsis and white tea (Camellia sinensis L.) [30][31][32]. Therefore, the early-flowering phenotype of genotype 80 might be negatively regulated by JA biosynthetic genes, and the interactive networks of these hormones may also contribute to early flowering in alfalfa.
The early flowering of alfalfa is associated with pathogenesis-related and signaling receptor kinase genes Plants often respond to biotic and abiotic stress by diverting more resources away from growth and development to promote early flowering [9]. For instance, Arabidopsis plants infected by fungal and bacterial pathogens showed an accelerated flowering time and increased transcript levels of floral integrator genes FT and GIGANTEA. Further investigation indicated that 19 genes encoding receptor-like proteins were upregulated in the Arabidopsis-Fusarium oxysporum-resistant mutants [9,33] and ERECTA lead to early flowering and increased slique size in Arabidopsis and tomato, respectively, and overexpression of the tomato LRR receptor-like kinase enhances drought tolerance [34,35]. In addition, overexpression of AtLRK10L1.2, the homolog of wheat LRK10-like kinase, confers positive response toward ABA signaling and drought tolerance, features associated with early flowering in Arabidopsis [16]. In line with these published results, we found that genotype 80 suffered from pathogen attacks and abiotic stress while exhibiting early flowering during our field experiments. Consistently, pathogenesisrelated and signaling receptor kinase genes were significantly downregulated in genotype 80 compared with the later-flowering genotype 195.
The early flowering of alfalfa is related to secondary metabolism and protein degradation 2-oxoglutarate/Fe (II)-dependent dioxygenases play versatile roles in secondary metabolite synthesis and catalyze epimerization and demethylation during floral development [36,37]. For example, PKDM78 and JMJ14, which 2OG oxygenases, suppress the expression of floral integrator genes FT, AP1, SOC1, and LFY by mediating H3K4 demethylation of chromatin, thereby inhibiting floral transition in Arabidopsis [38,39]. In this study, eight of ten 2OG oxygenases were significantly downregulated in the new leaves of genotype 80, indicating that 2OG oxygenases might negatively affect early-flowering. Furthermore, a previous study reported that glycosyltransferases catalyze the glycosylation of plant secondary compounds to maintain cell homeostasis and regulate plant growth and development [40]. Ectopic overexpression of glycosyltransferase genes UGT87A2 and PtGT1 were found to inhibit FLC expression and increase floral integrator gene expression, thus resulting in an earlyflowering phenotype in Arabidopsis and tobacco, respectively [40,41]. Here, seven UDP-glycosyltransferases were significantly upregulated in the new leaves of genotype 80 compared with those of genotype 195; these genes were also significantly upregulated in the flower buds compared with mature or new leaves in genotype 80. Together, these results indicate the UGT genes might play critical roles in regulating early flowering in genotype 80. The widely existed ubiquitin/proteasome system affects plant development by participating in signaling transduction cascades, pathogen defense, and biotic and abiotic stress response [19]. For instance, HOS1 (high expression of osmotically responsive genes1), a RING finger-containing E3 ubiquitin ligase, negatively regulates CONSTANS abundance and delays flowering time in Arabidopsis [42]. By contrast, the FKF1 F-box protein physically interacts with CDF1 (CYCLING DOF FACTOR 1) to induce its degradation, thereby increasing the expression of CO, which promotes earlyflowering in Arabidopsis [43]. In addition, the F-box protein UFO (UNUSUAL FLORAL ORGANS) associates with the AP3 promoter to enhance LFY transcription [44]. Consistent with these published data, we found that nine of the 11 RING/U-box superfamily proteins were significantly upregulated in new leaves in genotype 195 compared with genotype 80, indicating that these ubiquitin/proteasome-related genes may be negative regulators of flowering time. Six of seven F-box family proteins were significantly downregulated in the new leaves of genotype 195 compared with genotype 195, suggesting that these F-box proteins may promote early flowering in alfalfa. These findings suggest that secondary metabolism genes and protein degradation genes may be essential for flowering time regulation in alfalfa.

Conclusions
Based on the phenotypical, physiological and transcriptomic analyses, reduced JA content in new leaves and the downregulation of JA biosynthetic genes, including those encoding lipoxygenases, 12-oxophytodienoate reductase-like proteins, and SAMT may play essential roles in controlling the early-flowering phenotypes in alfalfa. Further analyses reveal that PR genes (including those encoding LRR family proteins, receptor-like proteins, and TIR-NBS-LRR class proteins), genes encoding signaling receptor kinase family members (including DUF26 and wheat LRK10-like kinases), secondary metabolism genes that encode 2OG oxygenases and UGTs, and proteasome degradation pathway genes (such as those encoding RING/U-box superfamily proteins and F-box family proteins) also contribute to the establishment of early-flowering in alfalfa (Fig. 11). In summary, our results provide insights into the regulatory mechanisms underlying early flowering in alfalfa and offer new target genes for future functional characterization and the genetic improvement of alfalfa.

Plant materials and growth conditions
Alfalfa genotype 80 (Accession No. PI183404, collected in India) and genotype 195 (Accession No. PI584512, collected in Canada) were provided by the USDA National Plant Germplasm System as previously described [45]. For quality control, raw reads of the FASTQ format were processed through in-house Perl scripts, after which reads containing adapter sequences were removed and 100.13 Gb of clean reads for fifteen alfalfa samples were obtained. Next, the Q30, GC-content, and sequence duplication level of the clean reads were calculated. All downstream analyses were performed on the clean reads with high quality. The de novo assembly of non-redundant unigenes was performed by using the trinity program (v20131110) [46].
Endogenous ABA, IAA, JA, and SA analyses in alfalfa Endogenous hormone contents were determined in the leaves (the same leaves as those used for the RNA-seq analysis) of genotypes 80 and 195. ABA, IAA, JA, and SA were extracted and purified using previously reported methods [45,48]. Briefly,~50 mg of frozen tissue was ground in liquid nitrogen and mixed with 500 μL of indolepropionic acid extraction buffer (2:1 dilution) and 50 μL of internal standards ( 2 H 6 ABA, 2 H 2 IAA, 2 H 5 JA, and 2 H 4 SA). The mixture was centrifuged at 15000 g for 5 min, the supernatant was mixed with 1 mL chloroform, and dissolved in 100 μL methanol and used for liquid chromatography coupled with tandem mass spectrometry (LC-MS) according to Cao et al. [48] . The content of each hormone was calculated according to the standard curve of the internal standards and normalized to fresh weight (FW) associated with three biological replicates.

qRT-PCR validation of sixteen DEGs
Total RNA was isolated from the same samples used for the RNA-seq analysis using an RNAprep Pure Plant Kit (Tiangen, Beijing, China) following the manufacturer's instructions. First-stand cDNA was synthesized using the Takara MLV-Reverse transcriptase (Takara Bio, Inc., Otsu, Japan) used for qRT-PCR analysis on the Bio-Rad CFX96 real-time PCR detection system with four biological replicates and using the SYBR Premix Ex Taq (Takara Bio, Inc., Otsu, Japan). Sixteen unigenes were randomly selected for the qRT-PCR analysis using gene-specific qRT-PCR primers (Additional file 3: Table S3). The 2 -ΔΔCT method [49] was used to calculate the relative expression level of each gene to the endogenous control Actin.

Statistical analysis
All data were subjected to the analysis of variance (ANOVA) in SAS 9.0 (SAS Institute, Cary, NC) using the general linear model. Fisher's least significant difference test (LSD, P < 0.05 was used to determine significant differences among treatments according to Ma et al. [50].
Additional file 1: Table S1. Annotation of DEGs identified in genotypes 195 and 80.
Additional file 2: Table S2. Mapman display of the functional categories of DEGs in genotypes 195 and 80.