Skip to main content
  • Research article
  • Open access
  • Published:

Comprehensive transcriptome analysis reveals distinct regulatory programs during vernalization and floral bud development of orchardgrass (Dactylis glomerata L.)



Vernalization and the transition from vegetative to reproductive growth involve multiple pathways, vital for controlling floral organ formation and flowering time. However, little transcription information is available about the mechanisms behind environmental adaption and growth regulation. Here, we used high-throughput sequencing to analyze the comprehensive transcriptome of Dactylis glomerata L. during six different growth periods.


During vernalization, 4689 differentially expressed genes (DEGs) significantly increased in abundance, while 3841 decreased. Furthermore, 12,967 DEGs were identified during booting stage and flowering stage, including 7750 up-regulated and 5219 down-regulated DEGs. Pathway analysis indicated that transcripts related to circadian rhythm, photoperiod, photosynthesis, flavonoid biosynthesis, starch, and sucrose metabolism changed significantly at different stages. Coexpression and weighted correlation network analysis (WGCNA) analysis linked different stages to transcriptional changes and provided evidence of inner relation modules associated with signal transduction, stress responses, cell division, and hormonal transport.


We found enrichment in transcription factors (TFs) related to WRKY, NAC, AP2/EREBP, AUX/IAA, MADS-BOX, ABI3/VP1, bHLH, and the CCAAT family during vernalization and floral bud development. TFs expression patterns revealed intricate temporal variations, suggesting relatively separate regulatory programs of TF modules. Further study will unlock insights into the ability of the circadian rhythm and photoperiod to regulate vernalization and flowering time in perennial grass.


Flowering is a critical developmental stage of most higher plants, which makes the plants produce seeds, thus to pass the genetic information from one to the next generation. The timing of flowering is regulated by multiple genetic and environmental factors. Specifically, the plant genes controlling flowering are induced via synchronization of climatic and environmental conditions. Thus, flowering time varies extensively according to climate and latitudinal or altitudinal gradients. For crop production, it is important to coordinate the flowering time with changes of environment to avoiding adverse natural conditions during flower differentiation [1]. With the advance of molecular biology, genes controlling flowering time in annual crops such as Oryza sativa, Triticum aestivum, Hordeum vulgare, and especially in model plants Arabidopsis thaliana, have been identified, which mainly belong to four pathways interacting to each other including photoperiod pathway [2], GA pathway [3], vernalization pathway, and autonomous pathway [1, 4]. For winter cereals, low temperatures and long days following low temperatures are the major driving factors for plant development and flowering.

Vernalization or jarovization was first defined by Lysenko et al. in 1986 when it was observed that wheat varieties required cold for stem elongation and flowering [5]. Vernalization in winter cereals depends on a highly integrated and complex system of development, involving numerous structural and regulatory genes in an intricate network of signaling pathways. Flowering time of certain species strongly correlated with day time length and temperature during vernalization [6]. In wheat, chromosomes 5A and 5D play an important role in determining the responses to vernalization [7]. In barley, a series of quantitative trait loci (QTLs) in clusters associated with heading dates have been found on chromosome 7 [8]. Subsequent studies revealed other colinear regions of chromosomes 5A [9], 5B [10], and 5D in wheat; the genes were cloned in sequential studies and named the vernalization family (VRN) since they are extensively involved in the vernalization response. The latter demonstrated that VRN3 is completely linked to a gene similar to the FLOWERING LOCUS T (FT) in A.thaliana [11]. A model was further developed to elaborate the function of VRN1, VRN2, and the FT genes in response to vernalization. Prior to cold exposure, VRN2 represses the expression of FT. During cold exposure, VRN1 expression increases, resulting in the repression of VRN2, which in turn allows activation of FT during long days, thus inducing flowering [12]. In A. thaliana, the vernalization pathway converge on FLOWERING LOCUS C (FLC), a MAD-box transcription regulator; furthermore, its activator FRIGIDA represses flowering via increase of FLC mRNA abundance [13, 14]. The upstream gene WRKY34-induced and CULLIN3A (CUL3A)-dependent on FRIGIDA, modulate flowering in response to vernalization [15]. A latest study on transcriptome sequencing and chromatin immunoprecipitation sequencing (ChIP-seq) suggested that VRN1 binds the promoter of FLOWERING LOCUS T-like, and also targets VRN2 and ODDSOC2 [16]. This reasearch provide more information on flowering regulation of cereal crop.

With the known genetic background and natural variations, A. thaliana is an excellent tool in studying the complexity of vernalization regulation. Currently, the full flowering network can only be approached in A. thaliana. Studies on Oryza sativa, Triticum aestivum, and Hordeum vulgare have led to the identification of components within individual signaling pathways that affect flowering as well as their positioning within molecular hierarchies. Even more importantly, the current study concentrates on annuals. In contrast to annuals, perennials require a long vegetative phase to accumulate and achieve the transition to the reproductive stage [17]. In addition, perennials still retain meristems in a vegetative state once they become adult [18]. Obvious differences and more complex molecular mechanisms are involved in vernalization, due to the differential life-cycles.

Orchardgrass (Dactylis glomerata L.) is a winter perennial gramineae grass, which is native to northern Africa, western and central Europe, and temperate Asia [19]. As an economically important perennial forage crop, orchardgrass has been widely used for a long time due to its good adaptability, abundant nutrients, and high biomass production [20]. Vernalization is an important physiological characteristic of orchardgrass, which decide the timing of flowering. A clear understanding of the molecular mechanism of flowering can provide a basis for targeted forage breeding. However, the molecular mechanisms that regulate vernalization of perennial grass species remain unclear. In this study, we generated a comprehensive and continuous transcriptome of orchardgrass to elaborate dynamic changes of gene expression among different stages and constructed the co-expression network of genes responding to vernalization and flowering regulating. This research will facilitate the understanding of flowering mechanisms of perennial grasses in general.


Plant material and RNA extraction

A orchardgrass variety DONATA (Registered No. 398) was grown in a greenhouse of the Sichuan Agricultural University located in Chengdu (30°42′N, 103°51′E; Chengdu city, Sichuan, China) under natural light conditions, then was transferred to field in Chongzhou farm (30°37′N, 103°40′; Chongzhou city, Sichuan, China) prior to temperature reduction in the year of 2016. For sampling, mixed young leaf samples were collected from three biological replicates at six different stages and immediately frozen in liquid nitrogen and stored in −80 °C freezers before being used. The six sampling stages included before vernalization January 4th, vernalization February 2nd, after vernalization March 2nd, vegetative growth stage March 24th, before heading April 9th, and heading stage May 5th. We defined these six sampling stages as before vernalization stage (BV_DON), vernalization stage (V_DON), after vernalization stage (AV_DON), vegetative growth stage (VG_DON), before heading stage (BH_DON), and heading stage (H_DON) using a sampling timeline (Additional file 1: Figure S1). From BV_DON to V_DON, this was defined as phase 1 (P1), V_DON to AV_DON was defined as phase 2 (P2), AV_DON to VG_DON was defined as phase 3 (P3), VG_DON to BH_DON was defined as phase 4 (P4), and BH_DON to H_DON was defined as phase 5 (P5) in our data analysis. DON refers to the orchardgrass cultivated variety DONATA.

Total RNA was extracted from samples by using the RNAprep pure plant kit (Tiangen Biochemical Technology Company, Beijing, China). RNA samples were treated with DNaseIto remove DNA, then RNA samples were monitored on 1% agarose gels. RNA purity was verified using the NanoPhotometer® spectrophotometer (IMPLEN, CA,USA) and the concentration was measured via the Qubit® RNA Assay Kit in Qubit® 2.0 Fluorimeter (Life Technologies, CA, USA). RNA integrity was assessed using the RNA Nano 6000 Assay Kit of the Agilent Bioanalyzer 2100 system (Agilent Technologies, CA, USA). RNA samples with 260/280 ratio between 1.8 and 2.0, 260/230 ratio between 2.0 and 2.5, and RNA integrity number above 8.0, were chosen for further library construction and sequencing.

Transcriptome sequencing and data analysis

Approximately, 3 μg RNA for each sample was used for the RNA sample preparations and a total of eighteen sequencing libraries were generated using NEBNext® Ultra™ Directional RNA Library Prep Kit for Illumina® (NEB, USA) following manufacturer’s recommendations. Library construction and transcriptome sequencing were conducted by the Novogene Bioinformatics Institute (Beijing, China) on Illumina Hiseq 4000 platform (Illumina, San Diego, CA, USA). Without reference genome, the clean reads were assembled as a reference genome via Trinity software, according to standard parameters for next analysis. Read counts per gene were expressed as the expected number of Fragments Per Kilobase of transcript sequence per Million base pairs sequenced (FPKM). To depict the global abundance of gene expression, the FPKM values of transcripts were accessed via box-plot. All the assembled unigenes were searched and annotated against the publicly available protein databases including Nr, Nt, Pfam, KOG/COG(EuKaryotic Orthologous Groups/ Clusters of Orthologous Groups, Swiss-prot, KEGG(Kyoto Encyclopedia of Genes and Genomes), and GO (Gene Ontology), using BLASTx analysis with an E-value cut-off of 1.0E-05.

DEGs were identified via pairwise sample comparisons in adjacent sampling points (V_DON-BV_DON, AV_DON-V_DON, VG_DON-AV_DON, BH_DON-VG_DON and H_DON-BH_DON) by using DESeq R package (1.18.0) [21]. The expression levels were identified as significant DEGs by applying the cutoff at a fold change greater than 2 and p-value below 0.05(FDR P value < 0.01) [22].

To further and systematically predict the complex biological functions of genes and to identify active biological pathways among developmental phases, the assembled unigenes were mapped against both GO and KEGG databases. Genes with an adjusted P-value below 0.05 found by DESeq were assigned as differentially expressed and employed for GO and KEGG analyses. GO enrichment analysis of DEGs was implemented with the GOseq R package, in which the gene length bias was corrected. GO terms with corrected P-values below 0.05 were considered significantly enriched by DEGs [23]. KEGG enrichment analysis of DEGs in KEGG database ( and KOBAS software was used to test the statistical enrichment of DEGs in KEGG pathways [24].

The WGCNA package was used for a weighted gene coexpression network analysis in R (v3.3.0) as described by Langfelder et al. [25]. A total of 25,071 genes with FPKM values above 0 in more than two sampling points were used. A gene expression adjacency matrix was constructed for analyzing network topology, the soft thresholding power = 7 in our analysis. The blockwiseModules was used to obtaining the modules by default setting. A topological overlap matrix was calculated by comparing the connectivity similarities of each pair of probes among all genes. An eigen-gene network was constructed to represent the relationships among the modules. The networks were visualized using Cytoscape v.3.0.0.

Principal component analysis (PCA) and hierarchical clustering were performed to assess transcriptome similarity among samples and to evaluate sampling between biological replicates. PCA was performed based on expressed genes in different samples by using the R program with the default parameter. Hierarchical clustering of samples via the complete linkage method showed the change of gene expression levels across different stages. Self-organizing map (SOM) neural network was used for clustering analysis in search of co-expression DEGs groups [26]. Both of these clustering analysis were performed using the R language (v3.3.0). The heat map was obtained, using the OmicShare tools, which is a free online platform (

Unigenes were examined from the iTAK database for families of transcription factors or regulatory motifs ( [27, 28].The local BLASTx similarity search was performed at E-value below 1.0E-05 against flowering-related genes sequences which were downloaded from NCBI database [29].

Verification by qRT-PCR

The validity of RNA-seq were verified by quantitative real-time PCR(qRT-PCR). Eighteen flowering related genes of interest were validated using qRT-PCR, including FRI (c149523_g1), LHP1 (c143664_g4), VRN1(c147469_g1), VIP1(c137956_g2), LHY (c146679_g3), COL1 (c151793_g1), WNK1 (c148831_g2), CDF2 (c127155_g1), GI (c151406_g1), COL16(c140653_g3), CRY1 (c137241_g1), FD (c128431_g1), GAI (c140748_g1), FVE (c137063_g1), SPL (c131707_g1), SPL3 (c134262_g2), SPL9 (c150483_g1) and FLC-like (c147268_g1). First-strand cDNAs were synthesized using Prime Script™ RT Master Mix Kit (RR036A) (Takara, Japan). The qRT-PCR reaction was performed using Bio-Rad CFX96 following the instructions for the SsoFast™ EvaGreen® Supermix Kit (SYBR Green) (#1725200AP) (Bio-Rad, USA). Three biological replicates were sampled were performed on each sample. The primers for unigenes were designed using online Primer BLAST program( Glyceraldehyde 3-phosphate dehydrogenase (GAPDH) was selected as reference gene [30]. Relative quantitative data were calculated based on the ΔΔCT method: normalization: (ΔCT = CT (sample) – CT (GAPDH); ΔΔCT = ΔCT (samples) - ΔCT (sample1); relative quantification = 2-ΔΔCT [31]. All primers were listed in Additional file 2.

Microscopic examination

Anatomic and microscopic examinations were conducted, using an Olymus CX41 microscope to distinguish the different stages and to accurately identify the floral initiation and development. The tissue was fixed with 75% (v/v) anhydrous ethanol and 25% (v/v) glacial acetic acid prior to microscopic examination. Pictures were obtained with an Olymus SZX12 stereo-microscope system (Olymus, Japan).


Identification of morphological characteristics and RNA-seq data statistics

We collected the samples weekly from December 2015 to May 2016 and determined the developmental stages via morphological and microscopic examination. The distinct stages were identified by changes of floral bud size (Fig. 1a, b, c, and d). Considering the development stage and climatic conditions, we identified six representative stages covering low temperature induction, booting stage, and flowering stage. In the first two phases (BV_DON and V_DON), orchardgrass plants have thickened lower internodes that form a crown (Additional file 3: Figure S2). The shoot apex and short stem are enclosed within the whorl of older leaf sheaths and near the ground level. After the flowering stimulus (AV_DON and VG_DON), the elongated stems (culms) of orchardgrass are divided into distinct nodes and hollows, but may be pithy or solid internodes, and these structures differentiate from the inflorescence (Fig. 1e, f, g, and h). The booting stage is critical for inflorescence formation. This stage was recognized by initiation and differentiation of floral organs from the flower primordium at the top of the reproductive stem.

Fig. 1
figure 1

The morphological photographs of young inflorescence and stem. a Photographs of flower primordium after vernalization. b, c and d, the young inflorescence after flower primordium initiation for (b) 5 days (c) 10 days and (d) 3 weeks. e and g, Photographs of (g) stem and (e) transection of stem before young inflorescence formation. f and h, Photographs of (h) internodes and (f) transection of internodes after young inflorescence formation

To obtain an overview of the Dactylis glomerata L. transcriptome information during vernalization and at the floral stage, we generated a set of consecutive and comprehensive transcriptome. Three biological replicates at each sampling point, a total of 18 samples, were collected for RNA-seq. A total of 1,108,545,336 raw reads were generated (Additional file 4: Table S1) and were deposited in the National Center for Biotechnology Information in Sequence Read Archive with accession number of SRR5341102. After trimming, high-quality clean reads were retrieved ranged from 48 M to 69 M for each sample. The clean reads were assembled into a total of 499,624 transcripts (> 200 bp) with a N50 of 1356 bp and a N90 of 339 bp. The assembled transcript dataset included 270,221 sequences with a N50 of 1000 bp and a N90 of 263 bp (Additional file 5: Figure S3A and B). In total, 133,371 unigenes were annotated in at least one of these databases, 13,701 unigenes in all databases (Additional file 6: Table S2). The expression of unigenes were showed in Additional file 7: Figure S4.

Identification of DEGs among stage comparisons

To investigate molecular differences among developmental stages, 21,508 DEGs were found in five pairwise stage comparisons, including 3247 in V_DON-BV_DON (V vs BV), 5292 in AV_DON-V_DON (AV vs V), 4410 in VG_DON-AV_DON (VG vs AV), 7859 in BH_DON-VG_DON (BH vs VG), and 700 in H_DON-BH_DON (H vs BH), respectively. Moreover, the number of up−/down-regulated genes were also presented (Fig. 2a). The number of DEGs in BH vs VG was remarkably higher than other pairwise comparisons, indicating the involvement of complex developmental procedures between stages VG_DON and BH_DON. The PCA analysis and hierarchical clustering revealed that the 18 samples clustered in six corresponding discrete groups. The sampling groups BV_DON and V_DON as well as AV_DON and VG_DON revealed a shorter distance than other stages, which was also showed in hierarchical clustering (Fig. 2b and c). The FPKM values of all 18 samples were accessed via box plots (Fig. 2d). These data constituted a consecutive and distinct set, highlighting the specialized nature of their transcriptomes related to their biological characteristics.

Fig. 2
figure 2

Transcriptional relationship between samples. a Number of up- and down-regulated genes in five pairwise sampling stages, including stage V_DON vs stage BV_DON (V vs BV), stage AV_DON vs stage V_DON (AV vs V), stage VG_DON vs stage AV_DON (VG vs AV), stage BH_DON vs stage VG_DON (BH vs VG), and stage H_DON vs stage BH_DON (H vs BH). b Principal component analysis based on all expressed genes, showing six distinct groups of samples: BV_DON_rep1, BV_DON_rep2, BV_DON_rep3 (light green circle); V_DON_rep1, V_DON_rep2, V_DON_rep3 (blue regular triangle); AV_DON_rep1, AV_DON_rep2, AV_DON_rep3 (brown cross); VG_DON_rep1, VG_DON_rep2, VG_DON_rep3 (purple cross); BH_DON_rep1, BH_DON_rep2, BH_DON_rep3 (khaki rhombus), H_DON_rep1, H_DON_rep2, H_DON_rep3 (green inverted triangle). c Cluster dendrogram, showing the global relationship between biological replicates. The y-axis shows the degree of variance. d Number of expressed unigenes and respective expression levels in each sample type, based on the FPKM of biological replicates. Sample labels are as follows: BV, before vernalization; V, vernalization; AV, after vernalization; VG, vegetative growth; BH, before heading; H, heading. DON refers to the orchardgrass cultivated varity DONATA (Registered No.398)

The GO and KEGG analysis indicated, the pathways of plant hormone signal transduction, circadian rhythm-plant, and flavonoid biosynthesis were significantly enriched in P1, relevant to hormone signal, cold response, light signaling, and vernalization induction in up/down-regulated DEGs. After more than a month of vernalization, the cold induction was disappeared as temperatures rose up during P2. The pathways of starch and sucrose metabolism were significantly enriched both in up or down-regulated DEGs. Furthermore, genes related to phenylpropanoid biosynthesis expressed abundantly in the group of up-regulated DEGs, and the result was consistent with morphological observations and gene expression profiling (Additional file 8: Figure S5A, B and Additional file 9: Figure S6A, B).

The transition of vegetative growth to reproductive growth occurred at P3. Large amounts of DEGs were identified involving in cellular biosynthetic process, intracellular component biosynthesis, macromolecule biosynthetic, and many other processes related to cell proliferation and differentiation (Additional file 8: Figure S5-C and Additional file 9: Figure S6-C). The floral primordia mainly generated and differentiated at P4; the DEGs within the pathways of starch and sucrose metabolism, ribosome and phenylpropanoid biosynthesis were presented at this stage. In contrast, the DEGs related to the pathways of circadian rhythm-plant, carbon metabolism, photosynthesis-antenna proteins, and biosynthesis of amino acids were down-regulated (Additional file 8: Figure S5-D and Additional file 9: Figure S6-D). The inflorescence protruded from the leaf sheath at P5, and the foremost biological process during this stage was pollen formation. The pathways of cell cycle, meiosis, DNA replication, mismatch repair, nucleotide excision repair, purine metabolism, and pyrimidine metabolism, which related to chromosome replication and meiosis were significantly enriched in up-regulated DEGs (Additional file 8: Figure S5-E and Additional file 9: Figure S6-E). These data suggest that specific regulators were required for different developmental phases, which converged to orchestrate the development of whole vernalization induction, floral bud formation, and heading.

Temporal gene expression dynamics

To identify the major transcriptional dynamics associated with low temperature induced vernalization and the transition from vegetative growth to reproductive growth, the SOM-clustering approach were used to group genes according to similar expression profiles. A total of 6 coexpression clusters with 30 subclusters were defined (Additional file 10: Data S1. SOM-clustering results). Most subclusters showed a distinct expression peak at these 6 sampling points; however, the subcluster_5_5 and subcluster_5_6 were exceptional with two peaks. With the environment factors and growth situation within the gene expression dynamics, these sub-clusters were chosen and were looked into further (Fig. 3).

Fig. 3
figure 3

Clusters expression associated with stages. The x-axis shows the sampling stages, y-axis shows the relative expression.Sample labels are as follows: BV, before vernalization; V, vernalization; AV, after vernalization; VG, vegetative growth; BH, before heading; H, heading

Supercluster 1 contained subcluster_1_1, subcluster_1_2, subcluster_1_3, subcluster_1_4, and subcluster_1_5, which comprised genes primarily expressed in phase V_DON with maximum expression. Subcluster_1_1 had overrepresentation of genes from processes associated to transport, regulation of transcription, metabolic processes, lipid and fatty biosynthetic processes, and signal transduction. The signal transduction categories contained many genes involved in protein phosphorylation, oxidation-reduction process response to cold, oxidative stress, and other abiotic stress. Genes response to light stimulus and photosynthesis were also identified in a variety of biosynthetic pathways; furthermore, several of the genes in circadian regulation of gene expression were found in this cluster, such as CCA1 [32].

Subcluster_1_2 had an overrepresentation function category of ATP/ADP binding, protein binding, oxidative phosphorylation, and transcription factor activity mostly involved in DNA-templated transcription with the energy-related metabolism. Large numbers of DEGs in subcluster_1_3 were significantly enriched in plant hormone signal transduction, photosynthesis-antenna proteins, and circadian rhythm-plant pathways. Unigenes for hormone like auxin and ethylene were present in this cluster for functional action, including AUXIN RESPONSE FACTOR 18 (ARF18) and AUXIN RESPONSE FACTOR 22 (ARF22), which encode transcription factors involved in auxin signaling. In addition, the auxin responsive protein INDOLE 3 ACETIC ACID 16 (IAA16), IAA18, IAA28, ethylene responsive transcription factor 061(ERF061) and ERF073 were also found in this cluster.. Compared to subcluster_1_2, the pathways of photosynthesis and plant hormone signal transduction enriched more significantly in this cluster, as a response to light stimulus and abiotic stress.

For subcluster_1_4, a large proportion of genes in the sucrose metabolic process, starch metabolic process, and lipid metabolic process were linked. Many pathways such as cytoskeleton organization, cell wall macromolecule catabolic process, and acyl-carrier-protein biosynthetic process were involved in cell replication. In this cluster, several genes associated with calcium ion binding were expressed, which constituted part of the calcium signaling system involved in both the photosystem and late embryogenesis. Subcluster_1_5 comprised genes that expressed in molecular function of ATP binding, zinc ion binding, and DNA binding were involved in transferase activity, protein kinase activity, oxidoreductase activity, and hydrolase activity. The minimum expression peak of subcluster_4_2, subcluster_4_4, subcluster_4_5, subcluster_5_3, and subcluster_5_4 were also appeared in phase V_DON. After this time point, the gene expression demonstrates an increasing tendency in the subsequent developmental phases. These clusters including unigenes category are similar to supercluster 1, but have a completely opposite expression tendency in phase V_DON.

Subcluster_2_2, subcluster_2_6, subcluster_3_1, subcluster_3_6, and subcluster_4_6 showed an expression peak in phase AV_DON, but relatively low expression during other developmental phases. Genes in these clusters encode proteins with similar functions in signal transduction, protein phosphorylation, photosynthesis, oxidation-reduction process, and defense responses. While investigating the profile in subcluster_2_2, we identified several genes, encoding photosystem II reaction centre N protein (psbN). In subcluster_2_6 and subcluster_3_6, the genes showed the minimum expression in biological pathways such as oxidation-reduction processes. The genes of oxidation-reduction process, transport, sucrose metabolic process, carbohydrate metabolic process, and starch metabolic process were greatly changed relative to phase V_DON. Genes involved in cell morphogenesis and embryonic morphogenesis had minimal expression levels in subcluster_4_6 and genes response to auxin, reduction process, and nitrogen compound metabolic process had a maximum fold change compared to phase V_DON.

Prior to the heading phase, we focused on subcluster_1_6, subcluster_5_1, subcluster_5_2, subcluster_5_3, and subcluster_5_4, which all showed drastic expression changes. 290 genes in subcluster_5_2 followed an ascending expression trend overall, but with three expression peaks in both phases V_DON and BH_DON. Interestingly, the genes in V_DON showed an apparent polarization. For example, the expression transcription factor bHLH35 in protein dimerization activity, and the AP2 domain for transcription factor activity sharply decreased at this stage. Several genes started to express and some other genes had maximum expression. We found many TFs, such as bHLH35, bHLH92, BIM2, and WRKY transcription factor 40 (WRKY40) and WRKY70. Ethylene-responsive transcription factor 2(ERF2) and the MADS-box transcription factor 34 were also highly expressed in stage BH_DON.

In subcluster_1_6, the genes relevant to stress responses, cell proliferation, morphological development, and iron uptake were detected,such as sigma-70(σ70) and sigma-54(σ54). Photomorphogenesis and photosystem regulated genes from the cytochrome family were also found in abundance in this cluster, including cytochrome P450, cytochrome oxidase C, cytochrome B6-F, and cytochrome b562.. In addition, subcluster_5_1 also expressed several genes of the cytochrome family, encoding cytochrome b and cytochrome b245.

Subcluster_5_1 and subcluster_5_2 showed peak expression in phase BH_DON with an increasing tendency from phase BV_DON. Functional categories of unique genes associated with genes encoding heat shock protein (Hsp) responding to stress were Hsp 70, Hsp 71, Hsp 83, Hsp 90, and Hsp 98 which were also presented in this cluster, as well as the identified salt stress-induced protein. In subcluster_1_6, several genes in abiotic stress categories contained multiple genes encoding cold shock proteins, dehydrin, and G-protein, which may be induced by various stresses during development. Furthermore, a large proportion of elongation factors (EF) such as EF1-alpha, EF1-beta, EF1-delta, EF1-gamma, EF2, elongation factor G, and elongation factor-Tu were contained in this cluster.These clusters comprised several genes encoding lipases, containing the GDSL motif and lipid-transfer protein, which is involved in the lipid metabolic process. Many genes in the cell organization category were related to the cytoskeleton, including several tubulin alpha-1/−2, tubulin beta-1, and kinesin genes. In addition, genes encoding cell division cycle protein 48, cell division control protein 6, and cell division control protein 45 were identified in this cluster. Furthermore, a number of genes associated with cellulose synthase, cell wall, and cell cycle were expressed in this phase, possibly indicating a major shift in the regulation of cell morphogenesis.

Subcluster_5_3 and subcluster_5_4 showed low expression peaks in phase V_DON, followed by rapid increases from phase VG_DON to phase BH_DON. To evaluate the transcriptional changes associated with the morphologic changes after vernalization, we examined a number of genes in these clusters via enrichment in GO and KEGG. Genes in these clusters were predominantly enriched in phenylpropanoid and starch and sucrose metabolism in KEGG; also, the pathways of plant hormone signal transduction were involved. The functional category protein binding, transferase activity, and catalytic activity were revealed via GO enrichment. An overrepresentation of transcription factors were enriched, including AUXIN(AUX)/INDOLE-3-ACETIC ACID (IAA) family genes IAA6, IAA13, IAA21, and IAA31, the AUXIN RESPONSE FACTOR1(ARF1), and WRKY transcription factor family genes WRKY 27, WRKY 30, WRKY 40, WRKY 46, and WRKY 50. Otherwise, the primordia in embryos and flower development related transcription factors, no apical meristem (NAM) also involved in this cluster.

Subcluster_5_5 and subcluster_5_6 showed unique expression patterns, which have expression peaks in phases V_DON and VG_DON by negative regulation, and were also strongly expressed in phase BH_DON via positive regulation. These clusters contained a large number of genes involved in the protein process in the endoplasmic reticulum and phenylpropanoid biosynthesis, which played a role in cell morphogenesis and structural molecule activities. Subcluster_5_5 contained an overrepresentation of transcription factor induced by various abiotic factors, including ethylene-responsive transcription factor11 (ERF11), ethylene-responsive transcription factor34 (ERF34), and WAX INDUCER1 (WIN1).

Coexpression network analysis

To further investigate the coexpression network of candidate genes, a weighted correlation network analysis (WGCNA) was adopted. Based on pairwise correlations and gene expression trends in all samples, coexpression networks were constructed using the normalized microarray expression data of these 25,071 probes from all 18 samples (three biological replicates for six time points) via R library. Different color represents a specific module, containing a cluster of highly correlated genes. This analysis resulted in 14 distinct modules (Fig. 4a). Remarkably, 6 out of 14 coexpression modules comprise genes that are highly expressed in a single stage (Fig. 4b). For example, the brown module, included 14,121 genes that were significantly enriched in stage H_DON. Alignment of these genes in GO classes covered terms associated with signal transduction, transcription regulation, carbohydrate metabolic, and photosynthesis. Specific terms relevant to floral bud development and heading were abundant, such as cell wall synthesis and remodeling, including cell wall pectin biosynthetic processes, cell wall macromolecule catabolic processes, and cell wall modification. Moreover, the terms of cell morphogenesis, cell proliferation, cell growth, and cell division were also abundant in this module. These terms were identified and can be further used to predict key regulatory genes involved in floral bud development and heading. The red module contained the genes mainly expressed during stage BH_DON and GO terms in this module were similar to the brown module, which may demonstrate that these two modules shared a common gene network and thus respond to similar physiological processes. In the same way, we can identify the genes that significantly impact vernalization inducted by low temperatures. The genes in saddlebrown module and sienna 3 module, expressed lower than other modules, however, the skyblue module contained genes with higher expression in stage V_DON and may provide information for the key regulators in vernalization.

Fig. 4
figure 4

WGCNA of genes at six stages. a Hierarchical cluster tree shows coexpression modules, identified via WGCNA. Each leaf in the tree represents one gene. Major tree branches constitute 14 modules labeled by different colors. b This heatmap shows the gene relative expression of different modules in six stages, the y-axis represents the relative expression of modules and the x-axis shows the different stages. c Visualization of the eigengene network represents the relationships among the modules and the clinical trait weight. The hierarchical clustering dendrogram of the eigengenes shows the relationships among the modules. Heatmap shows the correlation of different modules, and the deeper red color represents the higher correlation. Sample labels are as follows: BV, before vernalization; V, vernalization; AV, after vernalization; VG, vegetative growth; BH, before heading; H, heading

In addition, the correlation in different modules was also considered, and 7 broad clades were identified in 14 modules (Fig. 4c). A heat map showed that there was a high degree of correlation among the darkgreen module, the royalblue module, and the skyblue module. The genes in these modules were massively expressed in stage BV_DON and V_DON. GO analysis in these modules identified series terms, involved in response to stress, signal transduction, and protein phosphorylation, indicating that these clusters may participate in the process of response stimulation and may guide vernalization. Some genes were coexpressed in the grey60 module, darkorange module, and orange module. In addition, the heatmap showed that the brown, red module, and darkmagenta modules coexpressed in stages BH_DON and H_DON. A host of transcription factor families were identified in these modules, such as the zinc finger transcription factor family (MIZ-type, C2H2-type, and C3HC4-type), the WRKY transcription factor family, and the NRT1/PTR family. In particular, we identified a series of genes that encode the YABBY protein, which plays a role in reproductive development [33].

The expression pattern of transcription factors and regulatory motifs at different stages

A total of 3079 transcription factors were identified and classified in 74 transcription factor families. We are concerned about the following nine families: WRKY, NAC, AP2/EREBP, Alfin-like, AUX/IAA, MADS-BOX, bHLH, ABI3/VP1 and CCAAT. In our data, we noticed that some transcription factor families expressed in specific period. For example, the Alfin-like transcription factor family peaked in phase H_DON and were only minor expressed in other phases (Fig. 5d). The basic Helix-Loop-Helix (bHLH) family, comprising transcription factors that are known to be active during flower development in A. thaliana, peaked in the phase of H_DON [34] (Fig. 5g). The nuclear transcription factor Y contains the three subunits NF-YA, NF-YB, and NF-YC that bind with highly specific CCAAT motifs in a variety of genes. We found these transcription factors overrepresented in phase H_DON (Fig. 5i). In addition to the above three transcription factor families, the AUX(AUXIN)/IAA(INDOLE-3-ACETIC ACID) transcription factors also have low expression level in other phases and with high expression level in phase H_DON (Fig. 5e).

Fig. 5
figure 5

Transcription factor profiling according to RNA-seq data. The heat map shows the expression of transcription factor families that are overrepresented among coexpression clusters. The y-axis represents transcription factor families: WRKY transcription factor families (a); NAC transcription factor families (b); AP2/EREBP transcription factor families (c); Alfin-like transcription factor families (d); AUX/IAA transcription factor families (e); MADS-BOX transcription factor families (f); bHLH transcription factor families (g); ABI3/VP1 transcription factor families (h); CCAAT transcription factor families (i), and the x-axis shows the different stages. Sample labels are as follows: BV, before vernalization; V, vernalization; AV, after vernalization; VG, vegetative growth; BH, before heading; H, heading

Otherwise, some transcription factor families have two expression peak among developmental stages. ABI3/VP1 transcription factor family encoded the B3 domain-containing protein, were the key components of auxin signaling include the AUX/IAA coreceptors and ARF transcription factors. Our data showed this type of transcription factors was strongly expressed in phase V_DON and phase H_DON (Fig. 5h), as well as the MADS-box transcription factors (Fig. 5f). The WRKY transcription factor family was widely distributed in the regulation of the transcriptional response to stress. In our data, we found the WRKY transcription factors expressed during different phases, especially in phase BH_DON and H_DON (Fig. 5a), which have the similar expression pattern of NAC transcription factor family (Fig. 5b). Unlike other transcription factor families, AP2/EREBP had expression peaks in phase V_DON (cold stress), AV_DON, and VG_DON (heat stress), as well as in BH_DON and H_DON (floral development) (Fig. 5c).

Identification of expression patterns of key flowering regulators in Orchardgrass

Due to the importance of flowering for plant reproduction, many studies report the gene function and genetic network of flowering regulation in A. thaliana [35]. And more than 200 flowering associated genes have been identified and characterized [36, 37]. In this study, we captured these candidate genes in multiple different flowering-related pathways to understand the dynamic change during the developmental switch. A total of 77 reported flowering-related genes were identified in our transcriptome data and the great majority of these were located in the pathways of vernalization and photoperiod (Fig. 6, Additional file 11: Table S3, Additional file 12). According to previous research, 24 potential genes including VRN1, VIN3, FRI, SVP, VIP1, VIP2, and SUF4 were implicated in vernalization. Furthermore, 22 candidate genes were involved in the photoperiod pathway, including MAF1, CDF2, NF-YB1, NFYB2, TIC, COL, and FD. In addition, the pathways of circadian clock, autonomous, GA, and aging contained 10, 8, 6, and 4 candidate genes, respectively. Most vernalization pathway related genes with high expression levels during phase 2 and phase 5 include VRN1 and FRI. VRN1 is a key regulator in the vernalization pathway. The high expression level of VRN1 in phase 2 matches the observation that low temperatures induce VRN1 [38]. In A.thaliana, FRI is an upstream regulator, repressing the expression of integrators FT via activating FLC, thus causing late flowering. Our data shows the FRI with high expression during phase 2, followed by a drastic decrease during phase 3. This result is consistent with research, showing that cold stress induces the degradation of FRI [15]. Within the photoperiod pathway, genes are enriched in short-day conditions (phase 2) and long-day conditions (phase 5), suggesting that genes in this pathway may respond to different length of hours of sunlight.

Fig. 6
figure 6

Putative schematic network of bolting and flowering regulation in orchardgrass. Arrows indicate positive regulation. The heat map shows the relative expression of candidate genes in different stages, the y-axis represents identified candidate genes, and the x-axis shows the different stages. These identified candidate genes are involved in various flowering pathways according to previous reports. Sample labels are as follows: BV, before vernalization; V, vernalization; AV, after vernalization; VG, vegetative growth; BH, before heading; H, heading

To confirm the RNA-seq results, 18 unigenes were selected for qRT-PCR. In general, the relative expression based on qRT-PCR was consistent with the sequencing data (Fig. 7). But there were also some inconsistent trends between qRT-PCR and sequencing data at specific stages. For example, the qRT-PCR results showed that the expression of COL1 and SPL3 were down-regulated in stage BH_DON comparing to stage H_DON, while the sequencing data showed up-regulated (Fig. 7F and P). In addition, the expression of CRY1 was up regulated from stage VG_DON to BH_DON base on the qRT-PCR results, while the sequencing data showed the down regulation (Fig. 7K). On the contrary, the expression of GAI was down-regulated in stage VG_DON relative to BH_DON, while the sequencing data showed up regulation(Fig. 7M).

Fig. 7
figure 7

qRT-PCR validation. The abscissa represents the different sampling time, ordinate represents the expression. The bar with oblique stripes represents the relative expression base on qRT-PCR results and the bar with faillette base on RNA-seq results. The bottom title represents the different genes. Sample labels are as follows: BV, before vernalization; V, vernalization; AV, after vernalization; VG, vegetative growth; BH, before heading; H, heading


Flowering time is one of the major factors affecting forage quality in orchardgrass. A comprehensive transcriptome profile of the critical periods, comprising vernalization, young spike differentiation, and flower development was obtained via RNA-seq. DEGs were identified and functional characteristics also elaborated. The analysis of expression patterns during six key developmental stages revealing the abundant diversity in gene expression. Stage-specific profiles indicated a dramatic shift in gene expression during different developmental stages, and while many genes were expressed constitutively, very few of these were expressed evenly throughout all stages. The analysis captured the distinct molecular characteristics during vernalization and flower formation in orchardgrass, which could be further explored to help understand flowering regulation in response to seasonal cues in grass.

Genome-wide expression profiling uncovered developmental modules in vernalization and flower formation

The number of DEGs between the adjacent sampling points indicates that there were much more DEGs in phase 4 (P4) than any other phases, and more remarkably, phase 5 only had the least DEGs. This demonstrated that P4 involves far molecular processes and most of these processes have had a regulatory effect before heading. This result agrees with the markedly morphological alterations during P4 (inflorescence formation) and P5 (heading). The plant perception of vernalization and outgrowth of floral buds, involved numerous processes, both temporal and spatial. GO and KEGG analysis of DEGs throughout all five phases showed that the processes focused differently, which also supported this notion. The protein related processes, such as protein phosphorylation and kinases, and protein modification, were highlighted during P1, while the functional hormone related categories were also overrepresented. Protein phosphorylation plays a central role in regulating many cellular processes in eukaryotes. In particular, protein phosphorylation is a major currency of signal transduction pathways [39]. This result indicates that the plant responded to multiple environmental stimuli during this phase. Compared to P1, the functional categories of carbohydrate metabolic processes and sucrose metabolic processes were significantly enriched and previous research suggests that these two processes may associate with gene expression throughout dormancy in crown buds [40]. Cold temperature was involved in breaking dormancy, inhibiting shoot growth, and inducing flowering competence [41]. These are critical phases of plant floral bud formation and development during P4 and P5. The functional categories including cell wall, cytoskeletal, polysacharide metabolic, and catalytic activity, collectively coordinate organogenesis. These biomolecular processes potentially relate to subtle morphological changes. KEGG analysis of DEGs highlighted the otherness of temporal and spatial.

We found various pathways that focus differently throughout developmental phases. In P1 and P2, plant hormone signal transduction, NOD-like receptor signaling, and MAPK signaling were emphasized. Plant hormones exhibit strong functions in response to endogenous signals and environmental cues. Furthermore, the ubiquitin biosynthesis pathways were also identified in this phase, and ubiquitin-mediated proteolysis has been shown to regulate the different steps of plant hormone signal transduction [42]. Ethylene is an important plant stress and defense hormone and MAPK modules act downstream of the receptors, thus mediating ethylene signaling [43]. Recently, researchers suggested that the starch and sucrose metabolism, carbon metabolism, and plant hormone signal transduction all mediate the timing of transition from vegetative to reproductive phase [44]. This report is consistent with our KEGG enrichment results during P3 and P4. The functional category of circadian rhythm-plant, flavonoid biosynthesis, and starch and sucrose metabolism were enriched in most phases, suggesting that biological processes during specific stages and non-specific stages play equally important roles.

Stage-specific analysis outlines the core molecular processes related to vernalization and flower formation

Coexpressed gene cluster analysis illustrated the temporal distribution of core molecular events during vernalization and floral bud formation. Our results suggested that transcriptome specialization is established at a specific stage for several biological processes. At vernalization stage, genes related to signal transduction and stress response were overrepresented. When temperature dropped, several mechanisms have been described to enhance the freeze tolerance, including changes in lipid composition, increases in active-oxygen-scavenging enzymes, anthocyanin accumulation, and carbohydrate metabolism, which are consistent with our transcriptome sequencing results [45]. Signal transduction includes protein phosphorylation, calcium ion signaling, plant hormone signaling pathways, and many metabolic processes. These various pathways reflect the biological processes of plants at this stage to receive and cope with the external multiple stimuli. For perennial grasses, vernalization induced by cold stimulation is critical for plant flowering. We can capture potential genes that mediate vernalization in these clusters. Our results also indicated that photosynthesis-related genes are show strong stage dependency and with high expression levels during V_DON stage, which may indicate these kind of genes are tightly co-expressed. Under normal conditions, low temperature causes photoinhibition of photosynthesis [46]; however, the elevated expression of these genes is consistent with previous studies. This suggests that the photosynthetic contribution (reproductive assimilation) of reproductive structures propagates production and influences flower initiation [47]. Previous research demonstrated that CO2 assimilation may play a regulatory role during early reproductive stages [48].

To identify genes associated with flower formation, we examined genes with expression peak in BH_DON and H_DON. Several clusters highlighted a number of specialized features of transcriptome data. Subclusters_1_6, subcluster_5_1, and subcluster_5_3 revealed the importance of cell wall remodeling by the active cellulose metabolism in the BH_DON stage, which agrees with reports of cellulose as a cell wall polymer as found in eukaryotes [49]. Plant cell wall-loosening proteins promote cell growth and are essential for floral bud initiation and differentiation. During the plant cell expansion process, microfibers are usually functional in the direction of cell expansion, which were laying down secondary cell walls [50]. Furthermore, genes related to cellular aromatic compound metabolic processes, cell wall modification, cell wall macromolecule catabolic processes, and cell wall biogenesis also play pivotal roles in cell growth [51]. In addition, genes within cell morphogenesis, cell division, cell cycle, and cell differentiation categories with higher expression levels, are indicative of multiple processes that actively coordinate tissue development and organ morphogenesis during floral bud formation [52].

This is exemplified by subclusters_5_5 and subcluster_5_6, which indicate a functional partition of different developmental stages by the expression pattern. These subclusters had higher expression during stages V_DON and VG_DON compared to other stages. The subclusters were characterized by genes related to signal transduction and stress responses, likely reflecting the fact that the external environment significantly influenced plant growth during these stages. A distinct gene set associated with cell proliferation and differentiation during BH_DON and H_DON provides evidence that genes with specific function express during a specific period. In other cases, expression preference exists even though gene clusters were expressed throughout the growth stage.

WGCNA analysis enabled the identification of specific modules during critical developmental stages. The skyblue module and darkmagenta module connect with vernalization induction in stages V_DON and AV_DON. Furthermore, the grey60 module, red module, and brown module, may contain numerous genes that are associated with growth transduction, floral organ development, and flowering control, respectively. Furthermore, WGCNA also provided evidence in the interaction of different modules during certain stages. This clearly revealed a network of genes instead of individual genes, which may help to uncover the molecular mechanisms underlying vernalization and flowering. Furthermore, the conjoint analysis of WGCNA and clustering analysis supplied a more efficient method to search core regulators in diverse and abundant data.

Transcriptional regulation among different stages

An analysis of genes that are associated with transcription regulation showed that TFs play an important role during vernalization and floral bud formation. Expression profiles during different developmental stages suggest temporal specialization of TFs and the data set can be used to identify key TFs relevant for flowering time control. For example, WRKY, MADS-BOX, and ABI3/VP1 transcription factor families were found to preferentially accumulate during vernalization, floral organ formation, and heading, which could be part of the transcriptional regulatory complex regulating stress response and flower development. Several WRKY transcription factors expressed in phase BH_DON and H_DON, which is consistent with their reported role in Miscanthus, including WRKY 12, which participates in cell wall formation and promotes flowering [53]. WRKY20, WRKY23, WRKY30, WRKY40, WRKY50, WRKY53, WRKY57, and WRKY70 were involved in ABA and jasmonic acid signaling pathways in moderating flowering [54]. The ABI3/VP1 transcription factor family encoded the B3 domain-containing protein and play a central role in embryogenesis and dormancy by regulating the gene expression during the plant embryo maturation via ABA responsive elements or conserved RY repeats [55]. In a previous study, the B3 domain-containing protein demonstrated to participate in the early development of flowers [56], including the reproductive meristem gene (REM) [57], which was preferentially expressed in reproductive meristems.

The MADS-box family has been defined on the basis of primary sequence similarity amongst numerous proteins [58]. In a previous study, the MADS-box was considered to play an important role during the origin and evolution of flower development [59]. MADS-box transcription factor suppressor of overexpression of CO1 (SOC1) signals from the GA-dependent pathway, which influences the flowering time during short days [60]. The presence of MAD-box transcription factors is consistent with the overrepresentation of the regulation of transcription functional categories.

In addition, Alfine-like and CCAAT transcription factor family were overrepresented in stages BH_DON and H_DON. Alfin-like transcription factors encoded several PHD finger proteins, including vernalization insensitive 3 (VIN3) and vernalization 5 (VRN5) [61, 62], which mediated both vernalization and photoperiod. The nuclear transcription factor Y contains the three subunits NF-YA, NF-YB, and NF-YC that bind with highly specific CCAAT motifs in a variety of genes. The NF-YA regulate transcriptionally and post transcriptionally to promote drought resistance [63, 64]. NF-YB plays a role in the regulation of the flowering time in A. thaliana [65]. Furthermore, members of NF-YC can physically interact with constans (CO), and are genetically required for CO-mediated floral promotion.

Previous research demonstrated basic helix-loop-helix (bHLH) transcription factor family to participate in controlling cell proliferation and development of specific cell lineages [66]. Our study found bHLH TFs significant enriched during booting stage and heading time, suggesting this type of transcription factor family may play a similar role in the regulation of flower formation and heading. NAC transcription factors (NAM, ATAF, and CUC2) were described during recently years [67], only a proportion of NAC proteins have been studied. NAC genes are tightly associated with embryonic, floral, and vegetative development [68]. We found NAC genes concentrated expression in phase BH_DON and H_DON in our data, meaning that this type of TFs are produced to join the pathways that lead to flower organ formation and flowering, which was consistent with the previous reports.

The ethylen-responsive transcription factor (ERF) is an important member of AP2/EREBP and encoded a type of AP2 containing protein. Many ERFs have been reported to be involved in the regulation of floral development and stress response [69]. This transcription factor family consists of a large quantity of genes that encode IAA, which have an important function in floral organ development and plant growth regulation [70]. Hormone-related TF families AUX/IAA and AP2/EREBP were found in abundance at stage V_DON, VG_DON, BH_DON and H_DON, this suggested that ethylene biosynthesis and signaling have great function in stress response and are particularly associated with flower development [71]. The auxin response factor (ARF) family comprises transcription factors that are known to act during different phases in floral organ formation and are operated by a complex transcriptional network [72]. Furthermore, ARF transcription factors contribute to the response of abiotic stresses, containing ARF7 and ARF19. ARF is an important regulator of auxin activity, and these genes have maximal expression during stage H_DON, suggesting common centers for auxin biosynthesis and transduction during heading time. The diverse expression profiles of gene related hormone signaling may constitute the basis for this stage-specific response, also supplying evidence for the importance of hormones in vernalization and flower development. Timing of CAB expression 1(TOC1) is a widespread study circadian clock gene belongs to AP2/EREBP family, which has been suggested to be a component of the central oscillator of controlling flowering time in A. thaliana. TOC1 is designated as Arabidopsis pseudo-response regulators (APRR1/TOC1), a circadian-associated transcription factor family. Other than APRR1/TOC1, most APRR family members have been implicated in the mechanisms underlying the circadian rhythm; in particular APRR5, APRR7, ARR4, ARR3, and ARR9 in A. thaliana [73], OsPRR73 and OsPRR95 in rice [74], and SbPRR37 in sorghum [75].

Most MYB proteins function in a variety of plant-specific processes especial in controlling plant development, responses to stresses [76]. The MYB-related gene Circadian 1 (CIR1) was identified in the transcriptome data and has been reported to affect several circadian-regulated processes [77]. In addition, MYB proteins mainly participated in signal transduction such as phytochrome A signaling, jasmonate signaling, ABA signaling, and GA signaling, which is involved in a variety of processes [78, 79]. Otherwise, the photoperiod and the circadian clock pathway gene STO (BBX24) that belongs to B-Box Family can affect the key flowering time genes FLC and FT/SOC1, thus regulating the flowering time [80]. The zinc finger proteins constans-like (COL6) and the SBP-box gene family squamosa promoter binding protein-like3 (SPL3) [81], which are involved in flowering regulation were also identified in transcription factors examination.

Key flowering regulators in orchardgrass

Studies on the flowering network of A. thaliana and other species uncovered several critical regulators for the integration of multiple flowering pathways, including FT and SOC [82]. These regulators are mediated by two dominant upstream genes, which implicate FLC in pathways of vernalization, autonomous, and aging. The previous study indicated the upstream gene FLC have negative regulation on FT. We found the FLC have low expression in vernalization stage and then increase gradually until heading stage. On the contrary, the FT have high expression in vernalization.while lowly expressed in other stages.This results were consistent with the existing conclusion that low temperature inhibit the expression of FLC and release the suppression for FT, which cause flowering transition [83]. Otherwise, the other upstream gene CO positively regulate the FT [84]. Our results showed that the CO and FT highly expressed in vernalization stage, which may provide support that the CO promote the expression of FT in flowering regulation. Identified these key regulators in RNA-seq data of D.glomerata demonstrated that different species may share a highly conservative flowering genetic network and several homologous critical candidate genes likely have a similar function. Furthermore, we present these flowering pathway related genes via expression dynamics, providing an intuitive display of when the genes are active during developmental phases. Furthermore, our results provide more essential information for functional analysis of flowering regulatory pathways in perennial grasses.


In conclusion, an RNA-seq approach was utilized to investigate the patterns of gene expression during six key flowering developmental stages involved in vernalization and flower development, revealing novel networks and key regulators. Stage-specific profiling provided biological information of molecular events. This evidence included the process of signal transduction, stimulation of the vernalization response, of hormone control, cell proliferation, and differentiation in floral organ formation. Furthermore, this study added insight into the vital function of transcriptional factors in plant growth as well as valuable information for plant biology in the area of flower development. The WGCNA approach revealed a tightly co-expressed gene clusters and highly ordered gene expression networks that control plant growth. Our work highlighted the effectiveness of the clustering analysis intersected with the WGCNA analysis tool in multi-sample and high-volume data analysis.



Abscisic acid


Abscisic acid insensitive 3 or maize ortholog viviparous 1


Apetala2 and ethylene-responsive element binding proteins


Arabidopsis pseudo-response regulators


Auxin response factor


Auxin/indole-3-acetic acid


After vernalization


Before heading


Basic helix-loop-helix


Before vernalization


Chromatin immunoprecipitation sequencing


Circadian 1




Cluster of orthologous group




Different expressed genes


The orchardgrass varity DONATA


Elongation factors


Ethylene-responsive transcription factor


False discovery rate


Flowering locus C


Fragments Per Kilobase of transcript sequence per Million base pairs sequenced


Flowering locus T




Germins and germin-like protein


Gene Ontology




Heat shock protein


Indole-3-acetic acid


Inducer of CBF expression 1


Identification and classification of plant transcription factors and protein kinases database


Jasmonic acid


Kyoto Encyclopedia of Genes and Genomes


Kegg Orthology Based Annotation System


Eukaryotic orthologous groups


Long after far-red light 1




Myeloblastosis, a transcription factor family




No apical meristem




Non-redundant protein sequence


Oryza sativa pseudo-response regulators


Protein family


Photosystem II reaction centre N protein


Quantitative trait locus


Reproductive meristem gene


Ribonucleic acid sequencing


Salicylic acid


Suppressor of overexpression of CO1


SBP-box gene family squamosa promoter binding protein-like3


Annotated protein sequence database


Transcription factors


Timing of CAB expression 1




Vegetative growth


Vernalization insensitive 3




Vernalization 5


APETALA1(AP1)-like MADS box gene in wheat


Weighted correlation network analysis


Wax inducer 1


  1. Amasino R. Seasonal and developmental timing of flowering. Plant J. 2010;61:1001–3.

    Article  CAS  PubMed  Google Scholar 

  2. Mouradov A, Cremer F, Coupland G. Control of flowering time: interacting pathways as a basis for diversity. Plant Cell. 2002;14(1):S111–30.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  3. Metzger JD. Thermoinductive regulation of gibberellin metabolism in Thlaspi Arvense L. Plant Physiol. 1990;94(1):151–6.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  4. Levy YY, Mesnage S, Mylne JS, Gendall AR, Dean C. Multiple roles of Arabidopsis VRN1 in vernalization and flowering time control. Science. 2002;297:243–6.

    Article  CAS  PubMed  Google Scholar 

  5. Flood RG, Halloran GM. Genetics and physiology of Vernalization response in wheat. Adv Agron. 1986;39:87–125.

    Article  Google Scholar 

  6. Huang HR, Yan PC, Lascoux M, Ge XJ. Flowering time and transcriptome variation in Capsella Bursa-Pastoris (Brassicaceae). New Phytol. 2012;194:676–89.

    Article  CAS  PubMed  Google Scholar 

  7. Law CN, Worland AJ, Giorgi B. The genetic control of ear-emergence time by chromosomes 5A and 5D of wheat. Heredity. 1976;36:49–58.

    Article  Google Scholar 

  8. Hayes PM, Blake T, Chen THH, Tragoonrung S, Chen F, Pan A, Liu B. Quantitative trait loci on barley (Hordeum vulgareL.) chromosome 7 associated with components of winterhardiness. Genome. 1993;36:66–71.

    Article  CAS  PubMed  Google Scholar 

  9. Galiba G, Quarrie SA, Sutka J, Morgounov A, Snape JW. RFLP mapping of the vernalization ( Vrn1 ) and frost resistance ( Fr1 ) genes on chromosome 5A of wheat. Theor Appl Genet. 1995;90:1174–9.

    Article  CAS  PubMed  Google Scholar 

  10. Iwaki K, Nishida J, Yanagisawa T, Yoshida H, Kato K. Genetic analysis of Vrn-B1 for vernalization requirement by using linked dCAPS markers in bread wheat ( Triticum Aestivum L.). Theorapplgenet. 2002;104:571–6.

    CAS  Google Scholar 

  11. Yan L, Fu D, Li C, Blechl A, Tranquilli G, Bonafede M, Sanchez A, Valarik M, Yasuda S, Dubcovsky J. The wheat and barley vernalization gene VRN3 is an orthologue of FT. Proc Natl Acad Sci U S A. 2006;103:19581–6.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  12. Woods D, McKeown M, Dong Y, Preston JC, Amasino RM. Evolution of VRN2/GhD7-like genes in vernalization-mediated repression of grass flowering. Plant Physiol. 2016;170:2124.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  13. Choi K, Kim J, Hwang H-J, Kim S, Park C, Kim SY, Lee I. The FRIGIDA complex activates transcription of FLC, a strong flowering repressor in Arabidopsis, by recruiting chromatin modification factors. Plant Cell. 2011;23:289–303.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  14. Johanson U, West J, Lister C, Michaels S, Amasino R, Dean C. Molecular analysis of FRIGIDA, a major determinant of natural variation in Arabidopsis flowering time. Science. 2000;290:344–7.

    Article  CAS  PubMed  Google Scholar 

  15. Hu X, Kong X, Wang C, Ma L, Zhao J, Wei J, Zhang X, Loake GJ, Zhang T, Huang J, et al. Proteasome-mediated degradation of FRIGIDA modulates flowering time in Arabidopsis during vernalization. Plant Cell. 2014;26:4763–81.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  16. Deng W, Casao MC, Wang P, Sato K, Hayes PM, Finnegan EJ, Trevaskis B. Direct links between the vernalization response and other key traits of cereal crops. Nat Commun. 2015;6:5882.

    Article  PubMed  Google Scholar 

  17. Hsu CY, Yuceer C. FLOWERING LOCUS T duplication coordinates reproductive and vegetative growth in perennial poplar. Proc Natl Acad Sci U S A. 2011;108:10756–61.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  18. Huang YJ, Liu LL, Huang JQ, Wang ZJ, Chen FF, Zhang QX, Zheng BS, Chen M. Use of transcriptome sequencing to understand the pistillate flowering in hickory (Carya Cathayensis Sarg.). BMC Genomics. 2013;14(1):299.

    Article  Google Scholar 

  19. Huang L, Yan H, Jiang X, Zhang Y, Zhang X, Ji Y, Zeng B, Xu B, Yin G, Lee S, et al. Reference gene selection for quantitative real-time reverse-transcriptase PCR in orchardgrass subjected to various abiotic stresses. Gene. 2014;553:158–65.

    Article  CAS  PubMed  Google Scholar 

  20. Casler MD, Fales SL, Mcelroy AR, Hall MH, Hoffman LD, Leath KT. Genetic progress from 40 years of orchardgrass breeding in North America measured under hay management. Can J Plant Sci. 2000;40:1019–25.

    Google Scholar 

  21. Wang L, Feng Z, Wang X, Wang X, Zhang X. DEGseq: an R package for identifying differentially expressed genes from RNA-seq data. Bioinformatics. 2010;26:136–8.

    Article  PubMed  Google Scholar 

  22. Anders S, Huber W. Differential expression analysis for sequence count data. Genome Biol. 2010;11:1–12.

    Article  Google Scholar 

  23. Young MD, Wakefield MJ, Smyth GK, Oshlack A. Gene ontology analysis for RNA-seq: accounting for selection bias. Genome Biol. 2010;11:R14.

    Article  PubMed  PubMed Central  Google Scholar 

  24. Mao X, Cai T, Olyarchuk JG, Wei L. Automated genome annotation and pathway identification using the KEGG Orthology (KO) as a controlled vocabulary. Bioinformatics. 2005;21:3787–93.

    Article  CAS  PubMed  Google Scholar 

  25. Langfelder P, Horvath S. WGCNA: an R package for weighted correlation network analysis. BMC Bioinf. 2008;9:1–13.

    Article  Google Scholar 

  26. Mangiameli P, Chen SK, West D. A comparison of SOM neural network and hierarchical clustering methods. Eur J Oper Res. 1996;93:402–17.

    Article  Google Scholar 

  27. Pérez-Rodríguez P, Riaño-Pachón DM, Corrêa LG, Rensing SA, Kersten B, Mueller-Roeber B. PInTFDB: updated content and new features of the plant transcription factor database. Nucleic Acids Res. 2010;38(Database issue):822–7.

    Article  Google Scholar 

  28. Jin J, Zhang H, Kong L, Gao G, Luo J. PlantTFDB 3.0: a portal for the functional and evolutionary study of plant transcription factors. Nucleic Acids Res. 2014;42(Database issue):1182–7.

    Article  Google Scholar 

  29. Nie S, Li C, Xu L, Wang Y, Huang D, Muleke EM, Sun X, Xie Y, Liu L. De novo transcriptome analysis in radish (Raphanus Sativus L.) and identification of critical genes involved in bolting and flowering. BMC Genomics. 2016;17:389.

    Article  PubMed  PubMed Central  Google Scholar 

  30. Reid KE, Niclas O, James S, et al. An optimized grapevine RNA isolation procedure and statistical determination of reference genes for real-time RT-PCR during berry development. BMC Plant Biol. 2006;6(1):27.

    Article  PubMed  PubMed Central  Google Scholar 

  31. Schefe JH. Quantitative real-time RT-PCR data analysis: current concepts and the novel "gene expression's CT difference" formula. J Mol Med. 2006;84:901–10.

    Article  CAS  PubMed  Google Scholar 

  32. Alabadí D, Oyama T, Yanovsky MJ, Harmon FG, Más P, Kay SA. Reciprocal regulation between TOC1 and LHY/CCA1 within the Arabidopsis circadian clock. Science. 2001;293:880–3.

    Article  PubMed  Google Scholar 

  33. Meister RJ, Oldenhof H, Bowman JL, Gasser CS. Multiple protein regions contribute to differential activities of YABBY proteins in reproductive development. Plant Physiol. 2005;137:651–62.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  34. Zhou HJ, Yu-Cui WU, Jin XX, Rao LH, Wang ZZ. Cloning and expression pattern analysis of transcription factor SmbHLH93 from salvia miltiorrhiza. Chinese Trad Herb Drugs. 2014;45:3449–55.

    CAS  Google Scholar 

  35. Amasino RM, Michaels SD. The timing of flowering. Plant Physiology. 2010;154:516.

  36. Srikanth A, Schmid M. Regulation of flowering time: all roads lead to Rome. CellMol Life Sci. 2011;68:2013–37.

    Article  CAS  Google Scholar 

  37. Fornara F, De MA, Coupland G. SnapShot: control of flowering in Arabidopsis. Cell. 2010;141:1–2.

    Article  Google Scholar 

  38. Campoli C, Korff MV. Chapter five – genetic control of reproductive development in temperate cereals. Adv Bot Res. 2014;72:131–58.

    Article  CAS  Google Scholar 

  39. Hunter T. Protein kinases and phosphatases: the yin and yang of protein phosphorylation and signaling. Cell. 1995;80:225–36.

    Article  CAS  PubMed  Google Scholar 

  40. Anderson JV, Gesch RW, Jia Y, Chao WS, Horvath DP. Seasonal shifts in dormancy status, carbohydrate metabolism, and related gene expression in crown buds of leafy spurge. Plant Cell Environ. 2005;28:1567–78.

  41. Moe R, Wickstrøm A. The effect of storage temperature on shoot growth, flowering, and carbohydrate metabolism in tulip bulbs. Physiol Plant. 1973;28:81–7.

    Article  CAS  Google Scholar 

  42. Giovanna F, Nam-Hai C. Ubiquitin-mediated proteolysis in plant hormone signal transduction. Trends Cell Biol. 2002;12:308–11.

    Article  Google Scholar 

  43. Yoo SD, Sheen J. MAPK signaling in plant hormone ethylene signal transduction. Plant Signal Behav. 2008;3:848–9.

    Article  PubMed  PubMed Central  Google Scholar 

  44. Xing L, Zhang D, Li Y, Zhao C, Zhang S, Shen Y, An N, Han M. Genome-wide identification of vegetative phase transition-associated microRNAs and target predictions using degradome sequencing in Malus Hupehensis. BMC Genomics. 2014;15:1–22.

    Article  Google Scholar 

  45. Allen DJ, Ort DR. Impacts of chilling temperatures on photosynthesis in warm-climate plants. Trends Plant Sci. 2001;6:36–42.

    Article  CAS  PubMed  Google Scholar 

  46. Gombos HW Z, Murata N. The recovery of photosynthesis from low-temperature photoinhibition is accelerated by the unsaturation of membrane lipids: a mechanism of chilling tolerance. Proc Natl Acad Sci U S A. 1994;91:8787–91.

    Article  Google Scholar 

  47. Karlsson H, Erwin B, Carlson B. Temperature and photosynthetic photon flux influence chrysanthemum shoot development and flower initiation under short-day conditions. J Am Soc Hortic Sci. 1989;114:158–63.

    Google Scholar 

  48. Vu JC, Yelenosky G, Bausher MG. Photosynthetic activity in the flower buds of ;Valencia' Orange (Citrus Sinensis [L.] Osbeck). Plant Physiol. 1985;78:420–3.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  49. Delmer DP. Cellulose biosynthesis: exciting times for a difficult field of study. Plant Biol. 1999;50:245–76.

    CAS  Google Scholar 

  50. Lipchinsky A. How do expansins control plant growth? A model for cell wall loosening via defect migration in cellulose microfibrils. Acta Physiol Plant. 2013;35:3277–84.

    Article  CAS  Google Scholar 

  51. Belisle JT, Vissa VD, Sievert T, Takayama K, Brennan PJ, Besra GS. Role of the major antigen of mycobacterium tuberculosis in cell wall biogenesis. Science. 1997;276:1420–2.

    Article  CAS  PubMed  Google Scholar 

  52. Fu Y, Gu Y, Zheng Z, Wasteneys G, Yang Z. Arabidopsis Interdigitating cell growth requires two antagonistic pathways with opposing action on cell morphogenesis. Cell. 2005;120:687–700.

    Article  CAS  PubMed  Google Scholar 

  53. Yu Y, Hu R, Wang H, Cao Y, Guo H, Fu C, Zhou G. MlWRKY12, a novel Miscanthus transcription factor, participates in pith secondary cell wall formation and promotes flowering. Plant Sci. 2013;212:1–9.

    Article  CAS  PubMed  Google Scholar 

  54. Jiang Y, Liang G, Yang S, Yu D. Arabidopsis WRKY57 functions as a node of convergence for jasmonic acid- and auxin-mediated signaling in jasmonic acid-induced leaf senescence. Plant Cell. 2014;26:230–45.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  55. Sakata Y, Kuwano T, Takaoka S, Takashima N, Taji T, Takenaga H, Tanaka S. Possible involvement of a bZIP protein in the repression of ABI3/VP1-mediated Chymotrypsin inhibitor gene expression at the late-seed maturation in winged bean. J Agric Sci Tokyo Nogyo Daigaku. 2006;51:113–21.

    CAS  Google Scholar 

  56. Wang Y, Deng D, Zhang R, Wang S, Bian Y. Systematic analysis of plant-specific B3 domain-containing proteins based on the genome resources of 11 sequenced species. Mol Biol Rep. 2012;39:6267–82.

    Article  CAS  PubMed  Google Scholar 

  57. Francozorrilla JM, Cubas P, Jarillo JA, Fernándezcalvín B, Salinas J, Martínezzapater JM. AtREM1, a member of a new family of B3 domain-containing genes, is preferentially expressed in reproductive meristems. Plant Physiol. 2002;128:418–27.

    Article  CAS  Google Scholar 

  58. Shore P, Sharrocks AD. The MADS-box family of transcription factors. Eur J Biochem. 1995;229:1–13.

    Article  CAS  PubMed  Google Scholar 

  59. Theissen G, Becker A, Di RA, Kanno A, Kim JT, Münster T, Winter KU, Saedler H. A short history of MADS-box genes in plants. Plant Mol Biol. 2000;42:115–49.

    Article  CAS  PubMed  Google Scholar 

  60. Zhao S, Luo Y, Zhang Z, Xu M, Wang W, Zhao Y, Zhang L, Fan Y, Wang L. ZmSOC1, a MADS-box transcription factor from Zea Mays, promotes flowering in Arabidopsis. Int J Mol Sci. 2014;15:19987–20003.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  61. Sung S, Amasino RM. Vernalization in Arabidopsis Thaliana is mediated by the PHD finger protein VIN3. Nature. 2004;427:159–64.

    Article  CAS  PubMed  Google Scholar 

  62. Sung S, Schmitz RJ, Amasino RM. A PHD finger protein involved in both the vernalization and photoperiod pathways in Arabidopsis. Genes Dev. 2006;20:3244–8.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  63. Yan-Jie Li YF, Ya-Ru F, Huang J-G, Wu C-A, Zheng C-C. NFYA1 is involved in regulation of Postgermination growth arrest under salt stress in Arabidopsis. PLoS One. 2013;8:e61289.

    Article  PubMed  PubMed Central  Google Scholar 

  64. Ni Z, Zheng H, Jiang Q, Hui Z. GmNFYA3, a target gene of miR169, is a positive regulator of plant tolerance to drought stress. Plant Mol Biol. 2013;82:113–29.

    Article  CAS  PubMed  Google Scholar 

  65. Kumimoto RW, Adam L, Hymus GJ, Repetti PP, Reuber TL, Marion CM, Hempel FD, Ratcliffe OJ. The nuclear factor Y subunits NF-YB2 and NF-YB3 play additive roles in the promotion of flowering by inductive long-day photoperiods in Arabidopsis. Planta. 2008;228:709–23.

    Article  CAS  PubMed  Google Scholar 

  66. Miao L, Zhang L, Raboanatahiry N, Lu G, Zhang X, Xiang J, Gan J, Fu C, Li M. Transcriptome analysis of stem and globally comparison with other tissues in Brassica Napus. Front Plant Sci. 2016;7(R106):1403.

    PubMed  PubMed Central  Google Scholar 

  67. Nuruzzaman M, Manimekalai R, Sharoni AM, Satoh K, Kondoh H, Ooka H, Kikuchi S. Genome-wide analysis of NAC transcription factor family in rice. Gene. 2010;465:30–44.

    Article  CAS  PubMed  Google Scholar 

  68. Olsen AN, Ernst HA, Leggio LL, Skriver K. NAC transcription factors: structurally distinct, functionally diverse. Trends plant Sci. Trends Plant Sci. 2005;10:79–87.

    Article  CAS  PubMed  Google Scholar 

  69. Licausi F, Ohmetakagi M, Perata P. APETALA2/ethylene responsive factor (AP2/ERF) transcription factors: mediators of stress responses and developmental programs. New Phytol. 2013;199:639–49.

    Article  CAS  PubMed  Google Scholar 

  70. Wang J, Yan DW, Yuan TT, Gao X, Lu YT. A gain-of-function mutation in IAA8 alters Arabidopsis floral organ development by change of jasmonic acid level. Plant Mol Biol. 2013;82:71–83.

    Article  CAS  PubMed  Google Scholar 

  71. Riechmann JL, Meyerowitz EM. The AP2/EREBP family of plant transcription factors. Biol Chem. 1998;379:633–46.

    CAS  PubMed  Google Scholar 

  72. Liu N, Wu S, Van HJ, Wang Y, Ding B, Fei Z, Clarke TH, Reed JW. Van dKE. Down-regulation of AUXIN RESPONSE FACTORS 6 and 8 by microRNA 167 leads to floral development defects and female sterility in tomato. J Exp Bot. 2014;65:2507–20.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  73. Ishida K, Yamashino T, Mizuno T. Expression of the cytokinin-induced type-a response regulator gene ARR9 is regulated by the circadian clock in Arabidopsis Thaliana. Biosci Biotechnol Biochem. 2008;72:3025–9.

    Article  CAS  PubMed  Google Scholar 

  74. Murakami M, Tago Y, Yamashino T, Mizuno T. Characterization of the Rice circadian clock-associated pseudo-response regulators in Arabidopsis Thaliana. Biosci Biotechnol Biochem. 2007;71:1107–10.

    Article  CAS  PubMed  Google Scholar 

  75. Murphy RL, Mullet JE. Coincident light and clock regulation of pseudoresponse regulator protein 37 (PRR37) controls photoperiodic flowering in sorghum. Proc Natl Acad Sci U S A. 2011;108:16469–74.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  76. Dubos C, Stracke R, Grotewold E, Weisshaar B, Martin C, Lepiniec L. MYB transcription factors in Arabidopsis. Trends Plant Sci. 2010;15:573–81.

    Article  CAS  PubMed  Google Scholar 

  77. Zhang X, Chen Y, Wang ZY, Chen Z, Gu H, Qu LJ. Constitutive expression of CIR1 (RVE2) affects several circadian-regulated processes and seed germination in Arabidopsis. Plant J. 2007;51:512–25.

    Article  CAS  PubMed  Google Scholar 

  78. Song S, Qi T, Huang H, Ren Q, Wu D, Chang C, Peng W, Liu Y, Peng J, Xie D. The Jasmonate-ZIM domain proteins interact with the R2R3-MYB transcription factors MYB21 and MYB24 to affect Jasmonate-regulated stamen development in Arabidopsis. Plant Cell. 2011;23:1000–13.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  79. Jaradat MR, Feurtado JA, Huang D, Lu Y, Cutler AJ. Multiple roles of the transcription factor AtMYBR1/AtMYB44 in ABA signaling, stress responses, and leaf senescence. BMC Plant Biol. 2013;13:1–19.

    Article  Google Scholar 

  80. Li F, Sun J, Wang D, Bai S, Clarke AK, Holm M. The B-box family gene STO (BBX24) in Arabidopsis Thaliana regulates flowering time in different pathways. PLoS One. 2014;9:e87544.

    Article  PubMed  PubMed Central  Google Scholar 

  81. Kim JJ, Lee JH, Kim W, Jung HS, Huijser P, Ji HA. The microRNA156-SQUAMOSA PROMOTER BINDING PROTEIN-LIKE3 module regulates ambient temperature-responsive flowering via FLOWERING LOCUS T in Arabidopsis. Plant Physiol. 2012;159:461–78.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  82. Nie S, Li C, Xu L, Wang Y, Huang D, Muleke EM, Sun X, Xie Y, Liu L. De novo transcriptome analysis in radish ( Raphanus Sativus L.) and identification of critical genes involved in bolting and flowering. BMC Genomics. 2016;17:1–16.

    Article  Google Scholar 

  83. Seo E, Lee H, Jeon J, Park H, Kim J, Noh YS, Lee I. Crosstalk between cold response and flowering in Arabidopsis is mediated through the flowering-time gene SOC1 and its upstream negative regulator FLC. Plant Cell. 2009;21:3185–97.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  84. Seonghoe J, Virginie M, Panigrahi KCS, Stephan W, Wim S, Deng XW, Federico V, George C. Arabidopsis COP1 shapes the temporal pattern of CO accumulation conferring a photoperiodic flowering response. EMBO J. 2008;27:1277–88.

    Article  Google Scholar 

Download references


We are grateful to B. Shaun Bushman from the USDA-FRRL for technical support in determined the sampling time point.


This research work was funded by the National Basic Research Program (973 program) in China (No. 2014CB138705), the National Natural Science Foundation of China (NSFC 31372363) and the earmarked fund for Modern Agro-industry Technology Research System (No. CARS-35-05). All these funding play roles in the design of the study and collection, analysis, and in writing the manuscript.

Availability of data and materials

All the data pertaining to the present study has been included in table and/or figure form in the manuscript and authors are pleased to share analyzed/raw data and plant materials upon reasonable request; The raw RNA-seq reads have been deposited in the NCBI database (accession SRR5341102, The plant materials were provided by Department of Grassland Science, College of Animal Science and Technology, Sichuan Agricultural University, Chengdu, China.

Author information

Authors and Affiliations



GF, LH and XZhang conceived and designed the experiments. XW and TH prepared samples for small RNA sequencing and analysed the data. XZhao performed the qRT-PCR. GF and LX contributed in dissection and microexamination, and drafted the manuscript. JW, JL and LP revised the manuscript. All authors have read and approved the final manuscript.

Corresponding author

Correspondence to Xinquan Zhang.

Ethics declarations

Ethics approval and consent to participate

Not applicable

Consent for publication

Not applicable

Competing interests

The authors declare that they have no competing interests.

Publisher’s Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Additional files

Additional file 1: Figure S1.

The photo of orchardgrass in different stages. Including stage before vernalization (BV_DON); vernalization (V_DON); after vernalization (AV_DON); vegetative growth (VG_DON); before heading (BH_DON); heading (H_DON). (TIFF 8381 kb)

Additional file 2:

The primers information for qRT-PCR. (XLSX 16 kb)

Additional file 3: Figure S2.

The crown at the bottom of stem. A, showed the overall view. B, showed the anatomical structure of the plant. C, showed the magnified structure. The white arrow point the crown. (TIFF 7662 kb)

Additional file 4: Table S1.

RNA-seq statistics. (DOCX 17 kb)

Additional file 5: Figure S3.

Statistics of de novo assembly of transcriptome. A, Transcript length distribution. B, Unigene length distribution. (TIFF 8933 kb)

Additional file 6: Table S2.

Statistics of annotation analysis of unigenes. (DOCX 16 kb)

Additional file 7: Figure S4.

The box-plot describing the FPKM distribution of expressed transcripts after filtering in different samples. Sample labels are as follows: BV, before vernalization; V, vernalization; AV, after vernalization; VG, vegetative growth; BH, before heading; H, heading. DON refers to the orchardgrass cultivated varity DONATA (Registered No.398). (TIFF 9497 kb)

Additional file 8: Figure S5.

GO functional classification of DEGs in five pairwise sampling stages. Including stage V_DON vs stage BV_DON(A), stage AV_DON vs stage V_DON(B), stage VG_DON vs stage AV_DON(C), stage BH_DON vs stage VG_DON(D) and stage H_DON vs stage BH_DON(E). (TIFF 10209 kb)

Additional file 9: Figure S6.

KEGG functional classification of DEGs in five pairwise sampling stages. Including stage V_DON vs stage BV_DON(A), stage AV_DON vs stage V_DON(B), stage VG_DON vs stage AV_DON(C), stage BH_DON vs stage VG_DON(D) and stage H_DON vs stage BH_DON(E). (TIFF 9769 kb)

Additional file 10:

Data S1 SOM-clustering results. (RAR 14661 kb)

Additional file 11: Table S3.

Identified flowering-related gene in orchardgrass. (DOCX 19 kb)

Additional file 12:

The annotation information for flowering-related genes identified in orchardgrass base on RNA-seq data. (XLSX 39 kb)

Rights and permissions

Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (, which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The Creative Commons Public Domain Dedication waiver ( applies to the data made available in this article, unless otherwise stated.

Reprints and permissions

About this article

Check for updates. Verify currency and authenticity via CrossMark

Cite this article

Feng, G., Huang, L., Li, J. et al. Comprehensive transcriptome analysis reveals distinct regulatory programs during vernalization and floral bud development of orchardgrass (Dactylis glomerata L.). BMC Plant Biol 17, 216 (2017).

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: