Whole transcriptome profiling of the vernalization process in Lilium longiflorum (cultivar White Heaven) bulbs
© Villacorta-Martin et al. 2015
Received: 14 August 2014
Accepted: 1 June 2015
Published: 28 July 2015
Vernalization is an obligatory requirement of extended exposure to low temperatures to induce flowering in certain plants. It is the most important factor affecting flowering time and quality in Easter lily (Lilium longiflorum). Exposing the bulbs to 4 °C gradually decreases flowering time up to 50 % compared to non-vernalized plants. We aim to understand the molecular regulation of vernalization in Easter lily, for which we characterized the global expression in lily bulb meristems after 0, 2, 5, 7 and 9 weeks of incubation at 4 °C.
We assembled de-novo a transcriptome which, after filtering, yielded 121,572 transcripts and 42,430 genes which hold 15,414 annotated genes, with up to 3,657 GO terms. This extensive annotation was mapped to the more general GO slim plant with a total of 94 terms. The response to cold exposure was summarized in 6 expression clusters, providing useful patterns for dissecting the dynamics of vernalization in lily. The functional annotation (GO and GO slim plant) was used to group transcripts in gene sets. Analysis of these gene sets and profiles revealed that most of the enriched functions among genes up-regulated by cold exposure were related to epigenetic processes and chromatin remodeling. Candidate vernalization genes in lily were selected based on their sequence similarity to known regulators of flowering in other species.
We present a detailed analysis of gene expression dynamics during vernalization in Lilium, covering several time points and accounting for biological variation by the use of replicates. The resulting collection of transcripts and novel isoforms provides a useful resource for studying the changes occurring during vernalization at a fine level. The selected potential candidate genes can shed light on the regulation of this process.
Lilium longiflorum (Easter lily) is a leading bulbous crop worldwide and is produced as cut flower, potted plant, garden plant and as dry cell bulb . Like many other ornamental bulbs , L. longiflorum flowering requires cooling of the bulbs to meet the obligatory vernalization requirement of this plant species [3, 4]. L. longiflorum, (cultivar White Heaven) plants developing from non-vernalized bulbs grown at a constant temperature of 25 °C produced only leaves and did not flower over a period of more than 15 months (Ram et al., in preparation), confirming the obligatory cold requirement of this cultivar. Vernalization is also the main parameter involved in flowering time regulation in Easter lily and therefore has been the focus of a considerable amount of research related to physiological aspects of this species’ development, in order to reach flowering at specific dates [3, 5]. Typically, cold exposure of L. longiflorum bulbs at 2 to 10 °C quantitatively hastens flowering time while decreasing height, leaf and flower number, up to a saturation point of 6 weeks, after which additional cold exposure will not have a further effect on these parameters [3, 5–10]. In a previous study on L. longiflorum cultivar White Heaven , we found that bulb exposure to 4 °C for one week induced a decrease of about 20 % in the time from planting to floral transition and to flowering. Additional cold exposure led to a gradual decrease up to about 80 % and 55 % for floral transition and flowering, respectively, after nine weeks at 4 °C.
Despite the importance of vernalization in Lilium flowering, the molecular regulation of this mechanism is largely unknown in this species and other ornamental flowering bulbs. Most of the information available on molecular control of vernalization comes mainly from work performed on Arabidopsis, cereals and sugar beet and revealed that, while the general mechanism of vernalization is conserved among distant species, the sequence of the main regulatory genes is not [12–17].
In Arabidopsis, FLOWERING LOCUS C (FLC), a MADS-box gene encoding a potent repressor of flowering, is active in meristems in autumn. Flowering repression by FLC is mediated by its binding to major genes that promote flowering, such as FLOWERING LOCUS T and D (FT and FD, respectively) and SUPPRESSOR OF OVEREXPRESSION OF CONSTANS1 (SOC1) [18, 19]. While FLC represses genes that induce meristems to form flowers, it relies on FRIGIDA (FRI) to elevate its autumnal expression to a level that prevents flowering [19, 20]. During winter, vernalization causes the acquisition of meristem competence to flower by repressing FLC expression. Once it has been repressed by vernalization, FLC remains off for the rest of the plant’s life cycle after the return of warm conditions, i.e. the repression is epigenetic in the sense that it is mitotically stable in the absence of the inducing signal (cold exposure). The mechanism of epigenetic repression of FLC involves histone modifications that convert FLC into a heterochromatin-like state. A key player in the vernalization-mediated silencing of FLC is VERNALIZATION INSENSITIVE 3 (VIN3), which is required for all FLC chromatin modifications associated with vernalization-mediated silencing and as a measure of the cold period . Recently, it was shown that all members of the VIN3 family act together to repress FLC family members during vernalization . In addition, the non-coding (nc) antisense transcript COOLAIR and the intronic long ncRNA COLDAIR are upregulated at different points during cold exposure and are apparently playing a role in the epigenetic regulation of FLC [23–25]. Altogether, this measure of gradual cold acquisition ensures that only a prolonged cold exposure (the winter season) will lead to activation of the vernalization process.
In winter cereals, which require vernalization, a system similar to that in Arabidopsis exists. Specifically, a flowering repressor prevents flowering prior to cold exposure and the expression of this repressor is turned off by cold. In wheat, the repressor is a zinc-finger type protein VERNALIZATION 2 (VRN2). One of the genes repressed by VRN2 is VERNALIZATION1 (VRN1), which encodes a MADS-box protein that promotes flowering [26, 27]. When expressed in the leaves at high levels, VRN2 also represses FT , thereby playing a similar role as FLC in Arabidopsis. However, these two genes are unrelated and no FLC orthologues have been isolated from grasses . In sugar beet, three genes have been found to regulate the vernalization response, BvBTC1, BvFT1, and BvFT2 . Both BvFT1 and BvFT2 belong to the phosphatidylethanolamine-binding protein (PEBP) family and act in an antagonistic way: overexpression of the flower repressor BvFT1 leads to the repression of BvFT2 (the homologue of FT) and delays flowering time. Furthermore, like FLC in Arabidopsis, BvFT1 is downregulated by cold exposure. Different alleles of the BvBTC1 locus are associated with the expression level of BvFT1 and BvFT2 and with the vernalization response and flowering habit of various beet genotypes .
High throughput sequencing can produce a wealth of information on the genes involved in a certain process. Individual genes can be annotated based on the prediction of open reading frames and by comparison with expressed sequence tag (EST) collections. For example, such technology was successfully used to analyse the transcriptome of vernalization and gibberellin responses in sugar beet, revealing connections between gene expression patterns, genotypes and treatments, as well as potential new functions of the RAV1-like AP2/B3 domain protein in the vernalization response . In Lilium, several studies based on ESTs for genetic analyses and marker development are available [29–32]. A report about the transcriptome of Asiatic lily hybrid (cultivar Tiny Ghost) bulbs 27 days after exposure to 25 °C or 4 °C was recently published. The reported results revealed changes in the expression of a large number of genes belonging to several metabolic pathways and orthologs to vernalization genes in other plant species . In that study, the length of cold exposure of the bulbs was established according to the peak in sugar content in the scales, however, the report does not specify the effect of cold exposure on flowering time, therefore, a definite association between the results and the actual vernalization process is not explicit.
The search for genes regulating vernalization in lily, which has a huge genome and no vernalization mutants available, can be narrowed with a selection of genes differentially expressed during cold exposure. In this study, we use RNA sequencing to investigate gene expression patterns in the meristems of lily bulbs at several time points during cold exposure, aiming to characterize the gradual changes occurring in the bulb during vernalization, in view of its quantitative effect on flowering.
Design of the experiment, plant material and RNA collection
Previous experiments conducted in our lab  have shown that exposure of L. longiflorum cultivar White Heaven bulbs to 4 °C reduces the time from planting to flowering in a quantitative manner. Therefore, several time points between 0 and 9 weeks of cold exposure were selected for the transcriptome construction in order to cover a wide range of vernalization related transcription patterns useful for later inference.
In order to obtain sufficient RNA for each extraction, five meristems were pooled. The number of individuals to be pooled is a compromise between the technical requirements of the library preparation protocol, the need to reduce the impact of potential outliers and the soundness of downstream statistical analysis: pooling too many individuals could lead to an underestimation of the variability of expression between replicates and thus make many false positive rejections of the null hypothesis when comparing different treatments, inflating the number of differentially expressed genes.
Two biological replicates were made for each time point. This detail of the experimental design is important, since one of the central goals of our study is to compare the expression in different stages of the vernalization process, for which we need to infer the expression variance in each condition . Moreover, only using biological replication can we test whether the difference between two samples with different treatment is stronger than what we expect to see between two samples that are replicates, and therefore, whether that difference can be attributed to the cold treatment.
RNA-seq sample preparation and sequencing
RNA was isolated using the AurumTM Total RNA mini kit (Biorad # 732-6820) following manufacturer’s instructions. The procedure was followed by an additional DNase treatment with TURBO™ DNase (Life technologies # AM2238) to remove all genomic DNA still present after the isolation. The RNA integrity was confirmed using the 2100 Bioanalyzer (Agilent Technologies).
The samples were prepared for transcriptome sequencing using the Illumina kit (TruSeq RNA Sample Preparation Kit v2, # RS-122-2001) following the manufacturer’s recommendations. The obtained libraries were adjusted and pooled at a concentration of 20 nmol and sequencing was performed with 100-bp, paired-end reads on the HiSeq 2000 (Illumina).
Quality control and data pre-processing
After base calling via the Illumina pipeline (HiSeq Built-In Software), a total of 528 million of paired-reads (MPR) was obtained. A stringent quality control was then applied by removing the first three bases in the 5’ ends of each read and also the bases at the 3’ ends with a quality score smaller than 20. Reads containing primer/adaptor sequences were removed as well as the ones with long mono and di-nucleotide repeats.
De-novo transcriptome assembly
Four methods for the de-novo transcriptome assembly were tested: tIDBA , trans-Abyss , Oases  and Trinity . The resulting assemblies were compared using a reference set of L. longiflorum ESTs and the criteria proposed by Martin et al. .
Considering the resulting metrics of accuracy, completeness and contiguity, Trinity was chosen to perform the final assembly because it showed the best performance for this dataset . In order to obtain a transcriptome reference sequence that is common to all the data points, the reads were merged from all the different treatments before the final assembly. Since all samples were from the same lily cultivar, which is vegetatively propagated, merging different samples did not increase the complexity of the assembly. The cut-off for minimum contig length was set at 200 bp. Parameters for the assembly and further bioinformatic analyses are provided in Additional file 4.
Mapping, summarization, filtering and exploratory analysis
Each sample was mapped separately using a Bowtie  wrapper bundled with Trinity. After this, the counts of mapping reads per contig and per gene were summarized. Variation between replicates and distance between samples was explored using heatmaps and principal components analysis (PCA). For this purpose we used a subset of the top 500 genes with most variable expression after variance stabilizing transformation (which renders the expression values homoscedastic in relation to the fitted expression mean). This allowed us to rule out unexpected batch effects.
In order to reduce the amount of chimeric transcripts and noise, sequences whose total read count over all conditions was below the 80th percentile of read counts were filtered out. That resulted in 42,430 genes passing this filter, a number coherent with the amount of genes found in other plant species. Further filtering was applied to remove isoforms supported by less than 4 % of the reads mapping in any given gene. The purpose of independent filtering is to get rid of genes that, if tested, would have no chance of showing significant evidence of differential expression. By using a statistically independent filter (one that does not depend on the test statistic), a strong multiple-test correction (with the procedure described by Benjamini and Hochberg, ) was avoided and thus, statistical power was increased.
Differential expression analysis
Differential expression between each pair of conditions was analyzed following the methods as implemented in DESeq v1.10.1 . Robust estimates of the variability in the expression of each gene were calculated doing an average of its variability over each condition (weighted by the number of samples for each combination of factors, which is the same for all conditions). For a confident dispersion estimate, a conservative approach was used: estimates of high dispersion genes (those above the fitted regression line) were accepted but estimates of low dispersion genes were pulled up to the fitted line . We used a corrected  p-value of 0.05 as threshold to call differentially expressed genes.
After analyzing differential expression, we extracted genes beyond a threshold of a 4 fold-change in expression between any pair of time points and clustered them according to their patterns of differential expression across time. This was made using both the standard K-means clustering algorithm and pre-computed models for time series expression data in STEM . In particular, this last method first defines a set of distinct and representative models of gene expression, all of whom use a value of zero as a starting point. The number of real gene profiles assigned to each model is then computed. The number of genes expected to be assigned to a profile is estimated by randomly permuting the original time point values, renormalizing the gene’s expression values and then assigning genes to their most closely matching model profiles, repeating this for a large number of permutations. The average number of genes assigned to a model profile over all permutations is the expected number of genes for that model. This provides a framework to test the statistical significance of profile clusters. The biological significance of the model or cluster was investigated by means of GO  enrichment analysis, as explained in the gene set analysis section. For this clustering, those profiles that did not show a minimum absolute expression change of 4 log-fold (based on the difference between the minimum and maximum) were filtered out.
Function prediction was done in three steps: first inferring ORFs with Transdecoder (bundled with Trinity, ), second, calculating sequence similarity to a protein sequence database (UniProtKB) and third, to a protein profile database .
Computational approaches to annotation is not straightforward due to ambiguities in the model sequence-structure-function. Therefore, it is important to back inferences in different sources that should be integrated according to the semantics and strength of the evidence for each annotation source. For this integration purpose, we used Argot2 web-service . We provided as an input the results of sequence-based methods (Blast to a relevant subset of UniprotKB: plant sequences with a GO annotation) and functional domains assignments (HMMER3  to PfamA profiles ).
For those genes whose ORFs could not be inferred or could not be annotated, an additional similarity search (DNA to protein database) was done using Blastx and mapping the GO terms to the Blast results with Blast2GO . Once all the relevant GO annotations were obtained, these were mapped to GOslim Plants to achieve a reduced representation of GO terms.
Gene set analysis
The purpose of gene set analysis (GSA) [50, 51] is to shift the focus from single genes to sets of related genes. GSA was performed using two very different approaches: first, over gene profiles (clusters of genes with similar time-series expression patterns) and second, over pairwise comparisons of differential expression performed between different vernalization treatment time-points.
For the first approach, GSA was performed using time-series profiles with STEM . In this way, the cluster significance was determined first (checking whether the number of genes in each profile is higher than expected by chance alone), and afterwards their biological significance using an enrichment calculation based on their GO annotation and on the expected number of genes annotated with that specific GO given that particular cluster size. This calculation involves a comparison to a base set of genes (i.e. background or complement used as a reference) which in this case was comprised by all the genes in the dataset. Gene sets were based on the functions and pathways represented by the initial GO annotations (each of the 3657 ontology terms defines a gene set) and also on the transposons and other repetitive elements annotated with RepeatMasker .
In the second approach, GSA was performed with the package for R/Bioconductor GSVA 1.4.4 . This uses as input a gene expression matrix in the form of read counts and a database of gene sets, which were defined over the whole transcriptome, using GO slim plants (a summarized version of the 3657 GO annotations that results in 94 more general terms) and a subset of the transcriptome corresponding to genes annotated as transposons, using RepeatMasker  for that purpose (as in the previous method). This package implements a non-parametric unsupervised method for assessing gene set enrichment. It calculates gene expression level statistic for each sample and then it orders them by rank. Next, a Kolmogorov-Smirnov-like rank statistic is calculated which is subsequently used to obtain the resulting GSVA enrichment scores.
Orthologue vernalization gene finding
In order to find genes related to vernalization, several approaches were followed. First we used an “omics” approach, already described in the functional annotation sections. We selected genes whose inferred ORFs were annotated with terms related to vernalization. From these sequences we selected those whose domains matched the ones present in genes known to be involved in vernalization from other species (Arabidopsis , rice  and wheat ). Then, the protein sequence of each of the lily candidate genes was blasted  against the UniProtKB database. Top hits were then filtered based on the highest percentage of hit coverage combined the highest percentage of sequence similarity. Expression profile from the top hits was used to compare it against the expression of their putative orthologues and used it as another filtering criterion. Additionally we blasted protein sequences of known genes from other species, involved in vernalization, against the transcriptome data using BioEdit software  Top hits were then selected and filtered as described previously.
The initial Trinity assembly yielded a total of 329,599 transcripts belonging to 212,304 components (potential genes). By filtering out genes with lower count numbers and removing isoforms that represent just a small fraction of a gene’s total counts (smaller than 4 %) we reduced the total number to 121,572 transcripts and 42,430 genes. The total size of the unfiltered transcriptome was 263,039,169 bp. The N50 was 1,443 bp. The longest contig was 16,512 bp and the shortest, 201 bp. G + C content was 42.3 %.
Differential expression analysis
Final numbers of differentially expressed genes and isoforms differentially expressed (DE). Between pairs of vernalization treatments (0 to 9 weeks at 4 °C) using the most conservative parameters
Clustering and time-series analysis
Significant temporal expression profiles and the genes associated with these profiles were identified by clustering and visualization of the time series expression data. We integrated gene ontology data to perform enrichment analyses for sets of genes having the same temporal expression pattern.
Overrepresented gene sets in the upregulated and downregulated expression time-series clusters. Genes are categorized in sets according to their Gene Ontology. Column “Observed” is the number of genes from each category or set assigned to the cluster. ”Expected” is the number of genes that was calculated in the permutation simulation based on the total size of each category and the uniqueness of the profile. The results are order according to the Fold change. Only significant results are shown (Bonferroni corrected p-value < 0.05)
GO enrichment for Cluster 1 upregulated (Profiles 41 and 42)
DNA replication initiation
DNA packaging complex
DNA bending complex
protein-DNA complex assembly
protein-DNA complex subunit organization
chromatin assembly or disassembly
DNA-dependent DNA replication
DNA conformation change
regulation of cyclin-dependent protein serine/threonine kinase activity
regulation of protein serine/threonine kinase activity
protein heterodimerization activity
protein kinase binding
regulation of kinase activity
regulation of protein kinase activity
regulation of cell cycle
cellular macromolecular complex assembly
protein complex assembly
protein complex biogenesis
mitotic cell cycle process
macromolecular complex assembly
protein complex subunit organization
macromolecular complex subunit organization
cell cycle process
cellular component assembly
cellular response to DNA damage stimulus
cellular component biogenesis
DNA metabolic process
intracellular non-membrane-bounded organelle
cellular component organization
cellular component organization or biogenesis
GO enrichment for Cluster 0 downregulated (Profiles 0 and 4)
RNA-dependent DNA replication
RNA-directed DNA polymerase activity
While other profiles (like the #38, Fig. 4) showed patterns of expression matching our expectations of vernalization as a gradual cumulative process, the number of genes assigned to them in our clustering algorithm was not significant (a random pattern could have similar numbers of genes assigned just by chance) and was therefore not further explored at this stage.
In profile #41 (Fig. 4), most of the significantly enriched GO terms (corrected p-val < 0.05) (significant and upregulated, see Fig. 4) were related to epigenetic processes and chromatin remodeling (Table 2). From these profiles, the expression of several genes was validated by qPCR. In general, there were good correlations between the results obtained by RNA-seq and qPCR (Additional file 1: Figure S1).
Annotation statistics and functional enrichment
From the final 121,572 transcripts, a total of 82,219 ORFs were inferred (of which 67,658 ORFs had an annotation). These ORFs correspond to 15,414 genes, which were annotated with up to 3,657 GO terms.
An inherent problem of using GO annotations to define gene sets, is the strongly unbalanced sizes of the resulting groups (some GO terms apply to thousands of genes, while other, more specific functions, are represented by fewer than 10 genes). This disparity is inconvenient for sound statistical analysis (see gene set analysis).
In order to further facilitate the description of the functions in the transcriptome and decrease the effect of expression heterogeneity in the smaller sets, we translated the previous GO annotations into GO slim plant ontologies. GO slims are cut-down versions of the GO ontologies containing a subset of the terms in the whole GO. They give a broad overview of the ontology content without the detail of the specific fine grained terms. This turns out particularly useful as a summary of the results of GO annotation.
Instead of creating a custom GO slim vocabulary or a generic “all-species” one, we used an intermediate approach: a standard GO slim plant vocabulary developed by The Arabidopsis Information Resource (TAIR)  in order to facilitate comparisons with other plant species in future studies. With this new ontology, the number of functions was summarized in 94 annotations (Additional file 2: Figure S2), which could be further tested for overrepresentation using gene set analysis.
Gene set analysis
Gene set enrichment analysis is a way of abstracting the differential expression analysis from the level of gene expression profiles into a pathway or functional level. It provides a reduction dimensionality of the sample (testing a few hundred gene sets, rather than many thousand individual genes), as well as greater biological interpretability.
Specifically, with GSVA (Gene Set Variation Analysis)  we calculated sample-wise gene set enrichment scores as a function of genes inside and outside the gene set, analogously to a competitive gene set test (see Methods). Moreover, we estimated the variation of gene set enrichment over the samples independently of any class label. Conceptually, this is equivalent to a change in coordinate systems for gene expression data, from genes to gene sets.
Our first enrichment analysis (see Clustering and time-series analysis results above) was done for each of the time-course clusters using the complete GO annotation, and only 2 clusters (those with profiles 41/42 and with profiles 4/0) had functions significantly overrepresented on them and are therefore shown in detail in Table 2. Other clusters might simply not provide enough power for the hypergeometric enrichment test to detect any effect due to the smaller number of genes included or to the small number of genes annotated with each term.
This period between a non-vernalized and vernalized state, which is critical for the plant’s ability to flower in view of its obligatory requirement for cold, was characterized by major changes in gene expression (Fig. 3), which therefore might be related, among other processes, to the flowering mechanism.
Ortholog vernalization genes finding
Lily putative orthologues of vernalization-related genes and their expression in bulb meristems during cold exposure. Gene expression during cold exposure is schematically represented based on the log-fold change from week 0 to 9 at 4 °C
Putative ortholog in lily
Protein identity (%)
Expression profile in lily meristem during cold exposurea
FD / Arabidopsis
FT / Arabidopsis
SOC1 / Arabidopsis
VIP4 / Arabidopsis
VRN A1, VRN1 / Wheat, Barley
VRN B3, VRN H3 / Wheat, Barley
VRN2 / Arabidopsis
EMF2 / Arabidopsis
VIN3 / Arabidopsis
PHD-type Zinc finger
VIP6 (ELF8) / Arabidopsis
SUF4 / Arabidopsis
C2H2-type Zinc finger
SVP / Arabidopsis
Importance of the assembled transcriptome as reference for further research
In bulbous species, large genome sizes (36 Gbp for lily) make full genome sequencing expensive and error-prone . Transcriptome sequencing provides a proxy for high throughput comparisons of the exomes of crops with large genomes. In this way, ESTs of lily and tulip cultivars have been assembled using 454 pyro-sequencing, essentially for marker development . The transcriptome produced in this study provides additional molecular data for lily and a reliable quantitative profile of genes differentially expressed during cold exposure, enabling the molecular analysis of the vernalization response at a global scale.
While the physiological response of the plant to cold exposure has been well characterized in the literature [3, 5–10, 59] and by our own lab on this specific cultivar , the transcriptome analysis provides deeper insights into the molecular regulation of this vernalization response in a major flowering species. It will also enable the isolation of genes regulating flowering, which could be used to reduce vernalization requirements and to develop molecular markers for optimal bulb cold treatment, aiming at year-round production, high flower yield and quality, and reduced economic and environmental costs. In addition, it can serve to address fundamental questions regarding the conservation of the vernalization response among higher plants.
The emergence of systemic approaches for researching vernalization
From the methodological point of view, there is notable interest in importing the knowledge of vernalization acquired in model-species like Arabidopsis and wheat into other commercially relevant crops. This can be done using a singular, targeted approach -from genes to systems, or a more comprehensive, integrative approach -from systems to genes .
In our approach to tackle the vernalization process in Lilium, we have focused in the latter: the functional genomics of vernalization. To our knowledge, this is the first time in which vernalization of a bulb species is studied as a true dynamic process (including several time points), with a comprehensive level for comparative genomics (sequencing of more than 500 million read pairs yielding 121,572 isoforms) and a realistic model of expression variation (using biological replicates). With the assembly of a comprehensive de-novo transcriptome we provide the ground to develop and test hypotheses related to the dynamics of vernalization in a poorly studied, yet economically important group of plants.
During the annotation we identified novel isoforms for many potential orthologues of vernalization genes. A previous study measured expression of Asiatic lily bulbs before and after vernalization (27 days at 4 °C) using RNA-seq . Both studies show that vernalization induces many cellular, biological and metabolic changes and similar genes, such as the SOC1 homologue, are found upregulated after cold treatment in the bulb meristem hinting that this gene, known to be involved in flowering, might play a role in the vernalization pathway.
A recent report by Huang et al.,  attempted to understand the process of vernalization using what they refer to as “dynamic transcriptomes”. Their expression study involved a comparison between two samples taken before and after vernalization takes place but not during the vernalization process itself. Moreover, that comparison was done without considering biological variation in the expression. The approach followed in our research was to monitor the changes in expression throughout the whole process of vernalization. It is important to mention that the need of biological replication (sequencing different cDNA libraries extracted from different individuals) for assessing the effect of a given treatment in gene expression is critical and has been addressed before . The lack of this biological replication might be an explanation for the striking number of 68116 genes “differentially expressed”  found in the previously mentioned report by Huang et al., .
By using multiple time-points and biological replication in our research, we can say that the present work is the first dynamic study of the effect of vernalization in gene expression of lily bulbs.
The differential expression analysis for multiple time points during cold exposure allowed us to point at all the individual genes affected by vernalization at each stage. This information is critical for prioritization of downstream analysis and validation. On the other hand, gene set analysis allowed us to gain insight into the functions being influenced during vernalization and also the global effect that it has over the transcribed fraction of the transposable elements in the genome. Gene set analysis provides a scalable way to look at the changes that take place during vernalization at a high level.
The clustering performed for the time-series profiles also shows a general view on the number of genes that are up and downregulated, their functions and the patterns across time. None of these questions could be addressed adequately using an approach focused solely on individual genes. Additionally, unlike microarray studies, the use of high throughput sequencing technologies gave us margin to discover and measure new isoforms that would otherwise be unaccounted in the probes of array platforms.
Functional enrichment during cold exposure
The difference in the number of enriched functions that were upregulated versus downregulated by cold was strikingly high (see results of gene set analysis). In agreement with experiments in wheat and Arabidopsis, which had previously shown that vernalization is an epigenetic phenomenon , Table 2 illustrates empirically that the clusters found during bulb vernalization are consistent with theoretical expectations: the vast majority of enriched genes were annotated to functions related to chromatin and chromosome modifications and to functions related to cell division. Indeed the “memory of the cold” in the vernalization process is closely linked with cell division . Moreover, cell division is a pre-requisite for the integration and the maintenance of the vernalization signal [64, 65]. Therefore, upregulation of cell division-related genes during Lilium bulb vernalization is to be expected and is observed here.
Chromatin remodeling is an important regulatory mechanism of gene expression, well conserved among eukaryotes . The series of epigenetic events taking place during cold accumulation by the plant seem to be a common mechanism in flowering plants responding to vernalization . Indeed, chromatin modifications are tightly associated with the vernalization pathway in the context of the epigenetic repression of floral inhibitors. In the case of Arabidopsis, repression of the major floral repressor FLC by vernalization shows a pattern of epigenetic silencing comprising downregulation during cold exposure and continuation of silencing during the warm period following cold . This mechanism involves the action of chromatin-remodeling complexes, which are responsible for histone methylation at the FLC locus . In cereals, the vernalization response also involves epigenetic regulation. However, this regulation is not targeted to a floral repressor, but rather to the flowering enhancer VRN1, whose expression is silenced by histone methylations before cold exposure [68, 69]. In view of the conservation of cold-related epigenetic regulation of flowering among such distant species as Arabidopsis and cereals, it is probable that such a regulation exists in Lilium as well. The enrichment of chromatin related functions observed in the lily transcriptome during cold exposure could also be related to the silencing of transposable elements , whose expression decreased during cold exposure, as shown in the gene set analysis. Chromatin-related gene enrichment could also be attributed to the plant’s reponse to cold as an abiotic stress, as was observed for the ATP-dependent chromatin remodeling factors Snf2, in rice .
Ortholog and candidate vernalization genes
One of our goals was to identify candidate genes that play a role in lily vernalization. Based on their expression profile and/or sequence similarity to proteins involved in vernalization in other species, we identified 12 genes. These genes were classified then into enhancers or repressors (expected to be upregulated and downregulated respectively after cold induction) based on their reported function in other species. The expression profile of some of the candidate genes in lily did not match the expression reported in other species, suggesting that they may have a different role in the pathway. These results are similar to findings in beet, in which two paralogs of the FT gene (in Arabidopsis) have antagonistic functions. BvFT2 is functionally conserved with FT but BvFT1 represses flowering and its downregulation is essential for vernalization response in this species . In addition these data strongly suggests that, while the general mechanism of vernalization is relatively conserved among these distant species (a floral repressor is downregulated), the genes that mediate this process are not homologs and some of these components have resulted from convergent evolution [10–16]. For example, in wheat, although VRN2 is considered to have the same function as FLC as a floral repressor, it is evolutionarily unrelated, and so far there are not known orthologues of FLC in monocots . New functions for known genes are also hypothesized and will be subsequently validated.
In the differential expression analysis a relatively high number of genes were found to have significant changes during cold induction, which may only hint that these genes are candidates of regulators of vernalization in lily. Thus, a more thorough analysis of these data is needed. Additional experiments, including expression analysis of these and other genes in meristems, scales and leaves of lily plants grown under different environmental conditions are currently being performed, in order to further investigate the involvement of the genes in the vernalization pathway and in plant development. Furthermore, functional analysis of the selected candidate genes in lily and in Arabidopsis, also under way, will provide more information on how and when they are involved in the vernalization response. This will serve as a starting point to elucidate the vernalization pathway in this species and to have a detailed understanding of the molecular changes occurring in the bulb under cold induction.
In this study, we presented a detailed analysis of gene expression dynamics, covering a series of time points during bulb cold exposure in Lilium longiflorum, a leading ornamental crop. The resulting collection of transcripts and their novel isoforms provides a valuable data base for the exploration of fine molecular changes occurring during vernalization and the selected potential candidate genes can lead to the elucidation of the molecular regulation of vernalization in lily.
Overall, this study constitutes an important contribution to our current understanding of the molecular regulation and evolution of the vernalization response and can serve as a basis for similar research in other flowering bulbs.
Sequence data are available in the ArrayExpress database (www.ebi.ac.uk/arrayexpress) under accession number E-MTAB-2825. Parameters used for bioinformatics analyses are provided as Supplementary Information.
The authors thank Noam Rubin, Shani Sa’don-Shitrit and Iftach Mazor for performing the qPCR analyses and Sil'it Lazare and Rina Jager for meristem pictures. The study was funded by the Marie Curie Actions—Industry-Academia Partnerships and Pathways (IAPP) (Project number 286175).
- Grassotti A, Gimelli F. Bulb and Cut Flower Production in the Genus Lilium: Current Status and the Future. In: Grassotti A, Burchi G, editors. II International Symposium on the Genus Lilium, vol. 900. 2011. p. 21–35.Google Scholar
- Kamenetsky R, Zaccai M, Flaishman M. Florogenesis. In: Kamenetsky RaO, Hiroshi, editor. Ornamental Geophytes: From Basic Science to Sustainable Production. Boca Raton, FL: CRC Press Taylor and Francis Group; 2012. p. 197–232.View ArticleGoogle Scholar
- Miller WB. Lilium longiflorum. In: Nard DHaL, editor. The Physiology of Flower Bulbs. Amsterdam, London, New York, Tokyo: Elsevier; 1993. p. 391–422.Google Scholar
- Zlesak DC, Anderson NO. Inheritance of non-obligate vernalization requirement for flowering in Lilium formosanum Wallace. Isr J Plant Sci. 2009;57(4):315–27.View ArticleGoogle Scholar
- Roh SM, Willkins HF. Temperature and photoperiod effect on flower numbers in Lilium longiflorum Thunb. J Amer Soc Hort Sci. 1977;102(3):235–42.Google Scholar
- Dole JM, Wilkins HF. Interaction of bulb vernalization and shoot photoperiod on “Nellie White” Easter lily. Hortscience. 1994;29(3):143–5.Google Scholar
- Dole JM, Wilkins HF. Lilium, Easter. 1999.Google Scholar
- Holcomb EJ, Berghage R. Photoperiod, chilling, and light quality during daylight extension affect growth and flowering of tissue-cultured easter lily plants. Hortscience. 2001;36(1):53–5.Google Scholar
- Roh SM, Willkins HF. The effects of bulb vernalization and shoot photoperiod treatments on growth and flowering in Lilium longiflorum Thunb. Cv. Nellie White. J Amer Soc Hort Sci. 1977;102(3):229–35.Google Scholar
- Roh SM, Willkins HF. Comparison of continuous and alternating bulb temperature treatments on growth and flowering in Lilium longiflorum Thunb. J Amer Soc Hort Sci. 1977;102(3):242–7.Google Scholar
- Lugassi-Ben Hamo M, Martin CV, Zaccai M. Characterization of expressed sequence tags from Lilium longiflorum in vernalized and non-vernalized bulbs. J Plant Physiol. 2015;173:72–81.View ArticleGoogle Scholar
- Distelfeld A, Li C, Dubcovsky J. Regulation of flowering in temperate cereals. Curr Opin Plant Biol. 2009;12(2):178–84.PubMedView ArticleGoogle Scholar
- Kim DH, Doyle MR, Sung S, Amasino RM. Vernalization: winter and the timing of flowering in plants. In: Annual Review of Cell and Developmental Biology, vol. 25. Palo Alto: Annual Reviews; 2009. p. 277–99.Google Scholar
- Pin PA, Benlloch R, Bonnet D, Wremerth-Weich E, Kraft T, Gielen JJL, et al. An Antagonistic Pair of FT Homologs Mediates the Control of Flowering Time in Sugar Beet. Science. 2010;330(6009):1397–400.Google Scholar
- Pin PA, Zhang WY, Vogt SH, Dally N, Buttner B, Schulze-Buxloh G, et al. The role of a pseudo-response regulator gene in life cycle adaptation and domestication of beet. Curr Biol. 2012;22(12):1095–101.Google Scholar
- Ream TS, Woods DP, Amasino RM. The molecular basis of vernalization in different plant groups. Cold Spring Harb Symp Quant Biol. 2012;77:105–15.PubMedView ArticleGoogle Scholar
- Samach A. Control of flowering. Waltham, MA and London: Elsevier/Academic Press; 2012.View ArticleGoogle Scholar
- Helliwell CA, Wood CC, Robertson M, Peacock WJ, Dennis ES. The Arabidopsis FLC protein interacts directly in vivo with SOC1 and FT chromatin and is part of a high-molecular-weight protein complex. Plant J. 2006;46(2):183–92.PubMedView ArticleGoogle Scholar
- Sheldon CC, Finnegan EJ, Dennis ES, Peacock WJ. Quantitative effects of vernalization on FLC and SOC1 expression. Plant J. 2006;45(6):871–83.PubMedView ArticleGoogle Scholar
- Michaels SD, Amasino RM. FLOWERING LOCUS C encodes a novel MADS domain protein that acts as a repressor of flowering. Plant Cell. 1999;11(5):949–56.PubMed CentralPubMedView ArticleGoogle Scholar
- Sung SB, Amasino RM. Vernalization in Arabidopsis thaliana is mediated by the PHD finger protein VIN3. Nature. 2004;427(6970):159–64.PubMedView ArticleGoogle Scholar
- Kim DH, Sung S. Coordination of the vernalization response through a VIN3 and FLC gene family regulatory network in Arabidopsis. Plant Cell. 2013;25(2):454–69.PubMed CentralPubMedView ArticleGoogle Scholar
- Heo JB, Sung S. Encoding memory of winter by noncoding RNAs. Epigenetics. 2011;6(5):544–7.PubMedView ArticleGoogle Scholar
- Swiezewski S, Liu FQ, Magusin A, Dean C. Cold-induced silencing by long antisense transcripts of an Arabidopsis Polycomb target. Nature. 2009;462(7274):799–U122.PubMedView ArticleGoogle Scholar
- Turck F, Coupland G. When Vernalization Makes Sense. Science. 2011;331(6013):36–7.PubMedView ArticleGoogle Scholar
- Yan L, Fu D, Li C, Blechl A, Tranquilli G, Bonafede M, et al. The wheat and barley vernalization gene VRN3 is an orthologue of FT. Proc Natl Acad Sci U S A. 2006;103(51):19581–6.Google Scholar
- Yan LL, Loukoianov A, Blechl A, Tranquilli G, Ramakrishna W, SanMiguel P, et al. The wheat VRN2 gene is a flowering repressor down-regulated by vernalization. Science. 2004;303(5664):1640–4.Google Scholar
- Mutasa-Gottgens ES, Joshi A, Holmes HF, Hedden P, Gottgens B. A new RNASeq-based reference transcriptome for sugar beet and its application in transcriptome-scale analysis of vernalization and gibberellin responses. BMC Genomics. 2012;13:99.PubMed CentralPubMedView ArticleGoogle Scholar
- Du F, Wu Y, Zhang L, Li XW, Zhao XY, Wang WH, et al. De novo assembled transcriptome analysis and SSR marker development of a mixture of six tissues from Lilium Oriental Hybrid ‘Sorbonne’. Plant Mol Biol Rep 2015;33(2):281–293.Google Scholar
- Shahin A, Arens P, van Heusden AW, van der Linden G, van Kaauwen M, Khan N, et al. Genetic mapping in Lilium: mapping of major genes and quantitative trait loci for several ornamental traits and disease resistances. Plant Breed. 2011;130(3):372–82.Google Scholar
- Shahin A, van Kaauwen M, Esselink D, Bargsten JW, van Tuyl JM, Visser RGF, et al. Generation and analysis of expressed sequence tags in the extreme large genomes Lilium and Tulipa. BMC Genomics. 2012;13:640.Google Scholar
- Yuan SX, Ge L, Liu C, Ming J. The development of EST-SSR markers in Lilium regale and their cross-amplification in related species. Euphytica. 2013;189(3):393–419.View ArticleGoogle Scholar
- Huang J, Liu X, Wang J, Lü Y. Transcriptomic analysis of Asiatic lily in the process of vernalization via RNA-seq. Mol Biol Rep. 2014;41(6):3839–52.PubMedView ArticleGoogle Scholar
- Anders S, Huber W. Differential expression analysis for sequence count data. Genome Biol. 2010;11(10):R106.PubMed CentralPubMedView ArticleGoogle Scholar
- Peng Y, Leung HC, Yiu SM, Chin FY. Meta-IDBA: a de Novo assembler for metagenomic data. Bioinformatics. 2011;27(13):i94–101.PubMed CentralPubMedView ArticleGoogle Scholar
- Simpson JT, Wong K, Jackman SD, Schein JE, Jones SJ, Birol I. ABySS: a parallel assembler for short read sequence data. Genome Res. 2009;19(6):1117–23.PubMed CentralPubMedView ArticleGoogle Scholar
- Schulz MH, Zerbino DR, Vingron M, Birney E. Oases: robust de novo RNA-seq assembly across the dynamic range of expression levels. Bioinformatics. 2012;28(8):1086–92.PubMed CentralPubMedView ArticleGoogle Scholar
- Grabherr MG, Haas BJ, Yassour M, Levin JZ, Thompson DA, Amit I, et al. Full-length transcriptome assembly from RNA-Seq data without a reference genome. Nat Biotechnol. 2011;29(7):644–52.Google Scholar
- Martin JA, Wang Z. Next-generation transcriptome assembly. Nat Rev Genet. 2011;12(10):671–82.PubMedView ArticleGoogle Scholar
- Villacorta Martín C, Huijben K, Zaccai M, de Haan J, Passarinho P. RNA-seq priming bias and quality assessment of de-novo transcriptome assemblies. Poster session presented at: Benelux Bioinformatics Conference 2012 (BBC12). Nijmegen, Netherlands; 2012.Google Scholar
- Langmead B, Trapnell C, Pop M, Salzberg SL. Ultrafast and memory-efficient alignment of short DNA sequences to the human genome. Genome Biol. 2009;10(3):R25.PubMed CentralPubMedView ArticleGoogle Scholar
- Benjamini Y, Hochberg Y. Controlling the false discovery rate - A practical and powerful approach to multiple testing. J R Stat Soc B Methodological. 1995;57(1):289–300.Google Scholar
- Ernst J, Bar-Joseph Z. STEM: a tool for the analysis of short time series gene expression data. BMC Bioinformatics. 2006;7(1):191.PubMed CentralPubMedView ArticleGoogle Scholar
- Consortium TGO. Gene Ontology: tool for the unification of biology. Nat Genet. 2000;25:25–9.View ArticleGoogle Scholar
- Haas BJ, Papanicolaou A, Yassour M, Grabherr M, Blood PD, Bowden J, et al. De novo transcript sequence reconstruction from RNA-seq using the Trinity platform for reference generation and analysis. Nat Protoc. 2013;8(8):1494–512.Google Scholar
- Punta M, Coggill PC, Eberhardt RY, Mistry J, Tate J, Boursnell C, et al. The Pfam protein families database. Nucleic Acids Res. 2012;40(D1):D290–301.Google Scholar
- Falda M, Toppo S, Pescarolo A, Lavezzo E, Di Camillo B, Facchinetti A, et al. Argot2: a large scale function prediction tool relying on semantic similarity of weighted Gene Ontology terms. BMC Bioinformatics. 2012;13 Suppl 4:S14.Google Scholar
- Finn RD, Clements J, Eddy SR. HMMER web server: interactive sequence similarity searching. Nucleic Acids Res. 2011;39 suppl 2:W29–37.PubMed CentralPubMedView ArticleGoogle Scholar
- Conesa A, Götz S, García-Gómez JM, Terol J, Talón M, Robles M. Blast2GO: a universal tool for annotation, visualization and analysis in functional genomics research. Bioinformatics. 2005;21(18):3674–6.PubMedView ArticleGoogle Scholar
- Goeman JJ, Buhlmann P. Analyzing gene expression data in terms of gene sets: methodological issues. Bioinformatics. 2007;23(8):980–7.PubMedView ArticleGoogle Scholar
- Nam D, Kim SY. Gene-set approach for expression pattern analysis. Brief Bioinform. 2008;9(3):189–97.PubMedView ArticleGoogle Scholar
- Smith AFA, Hubley R, Green P. RepeatMasker Open-3.0. 1996-2010. http://www.repeatmasker.org.
- Hanzelmann S, Castelo R, Guinney J. GSVA: gene set variation analysis for microarray and RNA-seq data. BMC Bioinformatics. 2013;14:7.PubMed CentralPubMedView ArticleGoogle Scholar
- Lamesch P, Berardini TZ, Li D, Swarbreck D, Wilks C, Sasidharan R, et al. The Arabidopsis Information Resource (TAIR): improved gene annotation and new tools. Nucleic Acids Res. 2012;40(Database issue):D1202–1210.Google Scholar
- Kawahara Y, de la Bastide M, Hamilton JP, Kanamori H, McCombie WR, Ouyang S, et al. Improvement of the Oryza sativa Nipponbare reference genome using next generation sequence and optical map data. Rice. 2013;6(1):4.Google Scholar
- Brenchley R, Spannagl M, Pfeifer M, Barker GL, D’Amore R, Allen AM, et al. Analysis of the bread wheat genome using whole-genome shotgun sequencing. Nature. 2012;491(7426):705–10.Google Scholar
- Altschul SF, Gish W, Miller W, Myers EW, Lipman DJ. Basic local alignment search tool. J Mol Biol. 1990;215(3):403–10.58.Google Scholar
- Hall TA. BioEdit: a user-friendly biological sequence alignment editor and analysis program for Windows 95/98/NT. Nucl Acids Symp Ser. 1999;41:95–8.Google Scholar
- Amasino RM, Michaels SD. The timing of flowering. Plant Physiol. 2010;154(2):516–20.PubMed CentralPubMedView ArticleGoogle Scholar
- Leeggangers HA, Moreno-Pachon N, Gude H, Immink RG. Transfer of knowledge about flowering and vegetative propagation from model species to bulbous plants. Int J Dev Biol. 2013;57(6-8):611–20.PubMedView ArticleGoogle Scholar
- Robasky K, Lewis NE, Church GM. The role of replicates for error mitigation in next-generation sequencing. Nat Rev Genet. 2014;15(1):56–62.PubMed CentralPubMedView ArticleGoogle Scholar
- Song J, Angel A, Howard M, Dean C. Vernalization - a cold-induced epigenetic switch. J Cell Sci. 2012;125(16):3723–31.PubMedView ArticleGoogle Scholar
- Sung S, Amasino RM. Remembering winter: Toward a molecular understanding of vernalization. Annu Rev Plant Biol. 2005;56:491–508.PubMedView ArticleGoogle Scholar
- Burn JE, Bagnall DJ, Metzger JD, Dennis ES, Peacock WJ. DNA methylation, vernalization, and the initiation of flowering. Proc Natl Acad Sci. 1993;90(1):287–91.PubMed CentralPubMedView ArticleGoogle Scholar
- Wellensiek SJ. Dividing cells as the prerequisite for vernalization. Plant Physiol. 1964;39(832-835):832.PubMed CentralPubMedView ArticleGoogle Scholar
- Huff JT, Zilberman D. Regulation of biological accuracy, precision, and memory by plant chromatin organization. Curr Opin Genet Dev. 2012;22(2):132–8.PubMedView ArticleGoogle Scholar
- Zografos BR, Sung SB. Vernalization-mediated chromatin changes. J Exp Bot. 2012;63(12):4343–8.PubMedView ArticleGoogle Scholar
- Oliver SN, Finnegan EJ, Dennis ES, Peacock WJ, Trevaskis B. Vernalization-induced flowering in cereals is associated with changes in histone methylation at the VERNALIZATION1 gene. Proc Natl Acad Sci. 2009;106(20):8386–91.PubMed CentralPubMedView ArticleGoogle Scholar
- Trevaskis B. The central role of the VERNALIZATION1 gene in the vernalization response of cereals. Funct Plant Biol. 2010;37(6):479–87.View ArticleGoogle Scholar
- Hu Y, Zhu N, Wang X, Yi Q, Zhu D, Lai Y, et al. Analysis of rice Snf2 family proteins and their potential roles in epigenetic regulation. Plant Physiol Biochem. 2013;70:33–42.Google Scholar
- Abe M, Kobayashi Y, Yamamoto S, Daimon Y, Yamaguchi A, Ikeda Y, et al. FD, a bZIP protein mediating signals from the floral pathway integrator FT at the shoot apex. Science. 2005;309(5737):1052–6.Google Scholar
- Araki T, Kobayashi Y, Kaya H, Iwabuchi M. The flowering-time gene FT and regulation of flowering in Arabidopsis. J Plant Res. 1998;111(1102):277–81.View ArticleGoogle Scholar
- Lee H, Suh SS, Park E, Cho E, Ahn JH, Kim SG, et al. The AGAMOUS-LIKE 20 MADS domain protein integrates floral inductive pathways in Arabidopsis. Genes Dev. 2000;14(18):2366–76.Google Scholar
- Zhang H, van Nocker S. The VERNALIZATION INDEPENDENCE 4 gene encodes a novel regulator of FLOWERING LOCUS C. Plant J. 2002;31(5):663–73.PubMedView ArticleGoogle Scholar
- Dubcovsky J, Lijavetzky D, Appendino L, Tranquilli G. Comparative RFLP mapping of Triticum monococcum genes controlling vernalization requirement. Theor Appl Genet. 1998;97(5-6):968–75.View ArticleGoogle Scholar
- De Lucia F, Crevillen P, Jones AME, Greb T, Dean C. A PHD-Polycomb Repressive Complex 2 triggers the epigenetic silencing of FLC during vernalization. Proc Natl Acad Sci U S A. 2008;105(44):16831–6.PubMed CentralPubMedView ArticleGoogle Scholar
- Kim SY, Zhu T, Sung ZR. Epigenetic regulation of gene programs by EMF1 and EMF2 in Arabidopsis. Plant Physiol. 2010;152(2):516–28.PubMed CentralPubMedView ArticleGoogle Scholar
- He YH, Doyle MR, Amasino RM. PAF1-complex-mediated histone methylation of FLOWERING LOCUS C chromatin required for the vernalization-responsive, winter-annual habit in Arabidopsis. Genes Dev. 2004;18(22):2774–84.PubMed CentralPubMedView ArticleGoogle Scholar
- Kim SY, Michaels SD. SUPPRESSOR OFFRI 4 encodes a nuclear-localized protein that is required for delayed flowering in winter-annual Arabidopsis. Development. 2006;133(23):4699–707.PubMedView ArticleGoogle Scholar
- Hartmann U, Hohmann S, Nettesheim K, Wisman E, Saedler H, Huijser P. Molecular cloning of SVP: a negative regulator of the floral transition in Arabidopsis. Plant J. 2000;21(4):351–60.PubMedView ArticleGoogle Scholar
This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly credited. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated.