Digital gene expression analysis of male and female bud transition in Metasequoia reveals high activity of MADS-box transcription factors and hormone-mediated sugar pathways

Metasequoia glyptostroboides is a famous redwood tree of ecological and economic importance, and requires more than 20 years of juvenile-to-adult transition before producing female and male cones. Previously, we induced reproductive buds using a hormone solution in juvenile Metasequoia trees as young as 5-to-7 years old. In the current study, hormone-treated shoots found in female and male buds were used to identify candidate genes involved in reproductive bud transition in Metasequoia. Samples from hormone-treated cone reproductive shoots and naturally occurring non-cone setting shoots were analyzed using 24 digital gene expression (DGE) tag profiles using Illumina, generating a total of 69,520 putative transcripts. Next, 32 differentially and specifically expressed transcripts were determined using quantitative real-time polymerase chain reaction, including the upregulation of MADS-box transcription factors involved in male bud transition and flowering time control proteins involved in female bud transition. These differentially expressed transcripts were associated with 243 KEGG pathways. Among the significantly changed pathways, sugar pathways were mediated by hormone signals during the vegetative-to-reproductive phase transition, including glycolysis/gluconeogenesis and sucrose and starch metabolism pathways. Key enzymes were identified in these pathways, including alcohol dehydrogenase (NAD) and glutathione dehydrogenase for the glycolysis/gluconeogenesis pathway, and glucanphosphorylase for sucrose and starch metabolism pathways. Our results increase our understanding of the reproductive bud transition in gymnosperms. In addition, these studies on hormone-mediated sugar pathways increase our understanding of the relationship between sugar and hormone signaling during female and male bud initiation in Metasequoia.

Introduction

Metasequoia glyptostroboides Hu et Cheng, a famous living fossil (also known as the dawn redwood) and the sole extant species of Metasequoia, which has been merged with Cupressaceae as the subfamily Sequoia (Brunsfeld et al., 1994), is found in the Hubei and Hunan provinces in south-central China, and has been widely distributed in the Northern Hemisphere (Nelson, 1998; Kunzmann and Mai, 2011).

Although of great ecologic and economic importance, like most gymnosperm trees, Metasequoia has a long juvenile phase (20–30 years) before female bud appearance, making their seeds scarce. It is a monoecious species with male cones and female cones on different twigs of adult trees. In the past few decades, there have been studies on morphology and anatomy (Sterling, 1949; Dörken, 2012); however, the majority of these have focused on female cones. Few reports have examined male cones (e.g., Jin et al., 2012) because few trees produce them. In addition, despite its large numbers and worldwide distribution, most adult trees have low seed viability as a result of inbreeding depression (Li et al., 2005). Given the fact that male bud initiation is functionally adapted for successful pollination, exploring the gene expression involved in male reproductive development is important for understanding the molecular regulation of flowering and cone development, and for determining how and why there are few male cones on adult trees. At this time, few studies are available on M. glyptostroboides due to limited genomic resources. The only report on Metasequoia genomic resources is on an expressed sequence tag (EST) dataset built from both vegetative buds and reproductive buds during the late stages of bud development (Zhao et al., 2013). Thus, Metasequoia behind lagged behind other gymnosperm species such as Norway spruce, White spruce, and loblolly pine in genomic resources (Birol et al., 2013; Nystedt et al., 2013; Neale et al., 2014; Zimin et al., 2014).

Physiological and genetic studies have demonstrated that vegetative-to-reproductive transitions are not only regulated by environmental signals such as day length, light intensity, and ambient temperature (Trevaskis et al., 2006; Winfield et al., 2009), but are also regulated by endogenous signals such as plant hormones and age (Wang et al., 2011; Poethig, 2013). We previously found that spraying hormones (GA, IAA, ABA) on 5 year old trees could induce female and male buds and shorten this phase (Zhao et al., 2013). Recent studies have revealed the molecular mechanism of the phase between vegetative and reproductive transitions in annual species of Arabidopsis and maize, as well as some perennial wood species. In Arabidopsis, GA3 (Gibberellin) treatment can activate SUPPRESSOR OF CONSTANS1 (SOC1), which belongs to the MADS-box family of transcription factors, to induce flowering (Moon et al., 2003). Thus, MADS-box transcription factors are essential for the flowering program. Previous studies showed that many interactions between MADS-domain proteins are conserved among Arabidopsis thaliana, rice (Oryza sativa), petunia (Petunia hybrida), and Antirrhinum majus (Ferrario et al., 2004). Based on advances in our understanding of the genetic network controlling flowering in Arabidopsis, various biotechnological approaches have been used to shorten the juvenile phase in angiosperm trees (Bohlenius et al., 2006). However, recent studies have suggested that regulation of the vegetative-to-reproductive switch may differ in gymnosperms compared to angiosperms (Klintenas et al., 2012). Studies on the MADS-box gene family in gymnosperms suggest that some members in this family play an important role in shortening the juvenile phase and specifying reproductive identity in the male and female buds (Carlsbecker et al., 2004).

Sugar can provide energy for the plant, and act as a signal in coordination with hormone signaling pathways to control various plant physiological processes (Rolland et al., 2006). However, no sugar analysis was carried out in Metasequoia. Distinct glycolysis/gluconeogenesis, starch and sucrose metabolism pathways, and fructose signaling pathways have been examined (Cho and Yoo, 2011; Li et al., 2011). These pathways are regulated by various sucrose enzymes, including trehalose-6-phospath (T6p), UDPGluc, and glucose 6-phosphate (G6p) (Lunn et al., 2006; Huang et al., 2012). Sugar signaling plays an important role in controlling flowering time, and affects juvenile-to-adult phase transition in plants (Yu et al., 2013). Recently, the relationship between sugar metabolism/signaling and floral transition has been explored (Turnbull, 2011; Moghaddam and Van-den-Ende, 2013). However, the molecular mechanism controlling the transition to flowering through the hormone-mediated sugar pathway remains unclear.

In the present study, we performed digital gene expression (DGE) tag profiling using the Illumina next-generation sequencing (NGS) platform, taking advantage of its sensitivity for transcript analysis and high efficiency to sequence expressed genes compared to cDNA microarrays (Wang et al., 2009; Matsumura et al., 2010; Robinson et al., 2010). We mapped the Illumina reads to the previous 454 transcriptomic data, and obtained expression data for 69,520 putative unique transcripts (PUTs) and 243 associated Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways of Metasequoia. Furthermore, we compared the gene expression profiles of vegetative and reproductive buds (male and female) in different stages of bud formation, and identified numerous differentially and specifically expressed transcripts at different developmental stages of bud formation. Our study provides new information on the transcriptome analysis of Metasequoia buds, demonstrating the use of DGEs to examine transcript abundance of male and female buds in the early stage of bud initiation, and providing a valuable resource for studying gene function in gymnosperms.

Materials and Methods

Hormone Treatment and Sample Collection

Hormone (indole butyric acid, zeatin ribonucleoside, 0.5% Tween-20) treatments were conducted on 38 years old M. glyptostroboides trees located on the campus of Beijing Forestry University (latitude 40.0 N, longitude 116.2 E, 57 m above sea level) using the spray method, as previously described (Zhao et al., 2013). The treatments were divided into two groups in each tree, with three branches for each group. For group one, the hormone treatments were conducted at 5:00 p.m. on May 15 and May 25 in 2013. For group two, hormone treatments were conducted also at 5:00 p.m. on June 10, June 25, and July 10 in 2013. The treated shoots (TR) were sprayed until liquid drops started to drip down from the leaves, as previously described (Zhao et al., 2013). Untreated shoots (CK) in a similar location in each tree were used as the control.

Buds with the current year's stems (approximately 1 mm in length on TR and CK from the same tree) for three phenotypes were collected during bud differentiation when hormone treatment was efficacious for male bud induction on May 26 (samples TR1-1, TR1-2, TR1-3, CK1-1, CK1-2, CK1-3) and on June 11 (samples TR2-1, TR2-2, TR2-3, CK2-1, CK2-2, CK2-3), and was effective for female bud initiation on June 26 (samples TR3-1, TR3-1, TR3-3, CK3-1, CK3-2, CK3-3) and on July 11 (samples TR4-1, TR4-2, TR4-3, CK4-1, CK4-2, CK4-3). All of the samples were harvested from three individual trees (biological repeats) between 5:00 p.m. and 7:00 p.m., and were immediately frozen in liquid nitrogen. Samples were collected from all three trees for each time point, and were stored at −80°C for RNA extraction.

For the two group, some treated twigs remain for male bud and female percent calculation in winter and next spring, when it is easy for reproductive bud morphologically distinguishable with bare eyes. The estimate was done by the following formula: the percent of male (female) bud = the total number of male (female) buds/ the number of vegetative and male (female) buds.

RNA Extraction and Sequencing

Total RNA from all the samples was extracted separately using Trizol reagent (Invitrogen, CA, USA) according to the manufacturer's instructions. Total RNA quantity and purity were examined by Bioanalyzer (Agilent). Approximately 10 μg total RNA representing a specific type was subjected to mRNA isolation with poly-T oligo-attached magnetic beads (Invitrogen). Following purification, the mRNA was fragmented into small pieces using divalent cations under elevated temperatures. The cleaved RNA fragments were reverse-transcribed to create the final cDNA library in accordance with the protocol for the mRNA-Seq sample preparation kit (Illumina, San Diego, USA). 24 DGE libraries were constructed with Illumina's DGE Tag Profiling Kit according to the manufacturer's protocol (Version 2.1B). The eight (average of the three replicates) cDNA libraries in our study represented vegetative bud initiation (CK1 and CK2), male reproductive buds in early development (TR1 and TR2), vegetative buds in early development (CK3 and CK4), and female reproductive buds in early development (TR3 and TR4). Then, the 24 tag libraries underwent Illumina proprietary sequencing for cluster generation through in situ amplification. Single-end sequencing (36 bp) was performed on an Illumina Hiseq2500 at the LC-BIO (Hangzhou, China). Raw sequence data were available in GenBank Short Read Archive (Accession SRP051608).

DGE-Tag Profiling

The generated image files were processed to produce digital-quality sequence data. The raw data were filtered to obtain clean reads. The types of clean tags were represented as distinct clean tags. Subsequently, we classified the clean tags and distinct clean tags according to their copy number in the library, calculated their percentage in total clean and distinct tags, and analyzed library saturation. For annotation, all of the clean tags were mapped to unigenes of the Metasequoia database that we previously generated (http://www.genome.clemson.edu/node/273) with Bowtie2 (Langmead and Salzberg, 2012), because there is no Metasequoia genome sequence available.

Assembly of Different cDNA Libraries and De novo Transcript

All the mapped clean tags were assembled again according to the 96,565 PUTs. 454 original data without adaptors (Zhao et al., 2013) were De novo assembled again with Trinity (Haas et al., 2013).

Identification of DE Genes

The expression level of each gene was normalized to RPKM based on the number of clean tags. Genes were deemed significantly DE with a P < 0.05 (FDR < 0.1), and a relative change threshold of two-fold in the sequence counts across libraries. To obtain reliable and precise expression levels, the average of three replicates of RPKMs was used to determine the final fold change. All DE genes were divided into 16 comparisons (CK1 vs. CK2, CK1 vs. CK3, CK1 vs. CK4, CK2 vs. CK3, CK2 vs. CK4, CK3 vs. CK4; TR1 vs. TR2, TR1 vs. TR3, TR1 vs. TR4, TR2 vs. TR3, TR3 vs. TR4; CK1 vs. TR1, CK2 vs. TR2, CK3 vs. TR3, CK4 vs. TR4). All of the downregulated or upregulated genes in the following description were regulated in the second comparison component. Functional classification of DE genes was performed according to the functional categories of Gene Ontology (GO).

Biological Pathway Analysis

GO-defined pathways using the global test (available from Bioconductor: www.bioconductor.org) (Goeman et al., 2004) were significantly deregulated for comparisons. After summarization of the tags for each Entrez Gene entry, a global test was performed on the scaled and square root transformed data. P-values were calculated using the asymptotic method. Additional filtering of pathways was performed using the median of z-scores for each gene in the pathway to identify significant pathways.

Phylogenetic Analysis

Annotated MADS box genes with high similarity to Metasequoia PUTs putatively for MADS-protain from selected angiosperms and gymnosperms were retrieved from GenBank using BLAST searches. We obtained amino acid sequences and aligned them using the Clustal W program in Mega 5.0. A neighbor-joining (NJ) tree was obtained and local bootstrap probability of each branch was estimated with Mega 5.0 (Tamura et al., 2011).

Real-Time PCR Analysis

To verify the DGE results, a total of 32 Metasequoia PUTs were selected for quantitative RT-PCR (qRT-PCR) analysis using the same plant materials used for sequencing (three distinct biological samples). These transcripts were selected because of their high abundance in the cDNA libraries according to RNA sequencing. cDNA synthesis was performed using the first Strand cDNA Synthesis Kit (TaKaRa) with Oligo (dT) primers. The qRT–PCR experiments were performed on a StepOnePlus (Applied Biosystems) with SYBR Green supermix (Applied Biosystems). Three replicates were included for each sample. The reactions contained 45 cycles of 15 s at 95°C and 30 s at 60°C, followed by the thermal dissociation protocol for SYBR green detection. PCR Miner was used to analyze the raw data (Zhao and Fernald, 2005), and to calculate the cycle threshold (Ct-value) and PCR efficiency. The 2-ΔΔCT method (Schmittgen and Livak, 2008) was used to calculate the relative expression level of the gene of interest. 18S rRNA was used as the internal control. The primer sequences used are shown in Supplementary data 1.

Measurements of Non-Structural Carbohydrates

A total of 1.5 g of hormone-treated and untreated Metasequoia stems with buds (1 mm stem in length) were weighed and ground with mortar and pestle, separately. Five milliliters of 80% methanol (v/v) (0.1 M with pH 7) with 1.5 g samples was added into glass tube, incubated at 70°C for 15 min, and then centrifuged at 10,000 rpm for 10 min. This extraction was repeated twice. 0.5 ml of an internal standard (methyl-^a-D-glucopyranoside, 10 g/250 ml) was added to the supernatant, and then centrifuged at 10,000 rpm for 5 min. 80% methanol was then added until a volume of 25 ml. About 5 ml supernatant was concentrated to an anhydrous state after atmospheric distillation at 45°C for about 20 min. Then 1 ml derivatization reagent was added for solving for 15 min, and 0.5 ml was used for GC analysis (Bureau et al., 2013).

Results

Hormone-Treated Shoots Show Male Buds and Female Cones

According to our observation, in 2014, approximately 70% of female cones (Figures 1E,F) were found in June and July hormone treated shoots and the other half were vegetative buds. The shoots sprayed with hormone in May produced about 50% of male buds (Figures 1C,D), and the others were vegetative buds (Figures 1G,H). None of reproductive buds were observed on the control branches and on the shoots before hormone treatment (Figures 1A,B). The rectangle and elliptical female buds are 2–4 mm longer than the obovate male buds. The male cones were initiated in the axils of opposite leaves, with two male cones per leaf pair. Longitudinal sections of male buds shows pollen while female buds meristem were shown on the micrographs of female buds (Supplementary Figure S1).

FIGURE 1

Figure 1. Male buds and female cones on hormone-treated shoots. (A) Shoots before hormone treatment collected on May 24th, 2013. (B) Shoots before hormone treatment collected on June 11th, 2013. (C) Male buds on treated shoots collected on Sept 15th, 2013. (D) Male cones on treated shoots collected on Feb 5th, 2014. (E) Female buds on treated shoots collected on Oct16th, 2013. (F) Female buds on treated shoots collected on Mar 5th, 2014. (G) Untreated shoots (control) with vegetative buds collected on Oct 16th, 2013. (H) Untreated shoot (control) with vegetative buds collected on Mar 5th, 2014.

Sequence Analysis and Assembly of Different cDNA Libraries

Since there is no Metasequoia genome sequence available, the transcripts from different samples were reconstructed ab initio using a comprehensive set of 97,565 Metasequoia PUTs as a reference (Zhao et al., 2013). Based on 97,565 PUTs, we detected expression data for 69,520 genes, covering 71.3% of previous databases. After assembly, 32,186 PUTs at least 200 bp with an average of 740 bp (N50 = 1132) were generated. The original 454 data (Zhao et al., 2013) without adaptors were de novo assembled again. As a result, we generated 80,514 PUTs at least 200 bp in length with an average length of 807 bp (N50 = 1046). The size (bp) distribution of PUTs is represented in Figure 2. After assembly, 31,405 PUTs were mapped to the reference EST dataset of Metasequoia. The aligned sequences could be found in Supplementary data 2.

FIGURE 2

Figure 2. The size (bp) distribution of PUTs.

Our sequencing on an Illumina Genome Analyzer II generated a range of 25.7–50.2 million clean reads per library. When mapped to the 97,565 Metasequoia PUTs (http://www.genome.clemson.edu/node/273, Zhao et al., 2013), the number of alignable reads ranged from 14.6 million (CK1 on May 26, 56.7%) to 38.9 million (TR2 on June, 77.4%) based on clean reads. Among the alignable reads, 7.8–15.9 million were uniquely aligned reads, with CK1 making up the highest percentage (53.7%) (Table 1).

TABLE 1

Table 1. Summary of RNA-seq and mapping results of Metasequoia reads.

Go Annotation and DE Genes

To identify DE genes between reproductive and vegetative buds, we analyzed changes in expression levels among 69,520 transcripts. Statistical analyses with an adjusted P-value cutoff of 0.05 revealed PUTs with differential expression in all 16 comparisons (CK vs. CK, CK vs. TR, and TR vs. TR) (Figure 3). Comparisons with vegetative buds showed the most DE PUTs (6166). For CK and TR comparisons group, the highest number of DE PUTs was found in CK3 vs. TR3 (June 26) (1678), and CK1 vs. TR1 (May 26) showed the least (978) DE PUTs. For all of the comparisons, the most enriched categories in Biological Process, Cellular Component, and Molecular Functions were DNA integration, Nucleus, and ATP binding, respectively. Genes with expression responses covered a variety of regulatory and metabolic processes. All DE genes in each comparison were classified into several categories based on GO terms (http://www.geneontology.org/). In CK4 vs. TR4 comparison, genes related to DNA integration (5.87%) and RNA-dependent DNA replication (5.03%) were most abundant in DE genes. It should be noted that approximately 0.7% of DE genes were categorized as flower development (Figure 4).

FIGURE 3

Figure 3. Differentially expressed genes in different comparisons.

FIGURE 4

Figure 4. Gene classification based on gene ontology (GO) for differentially expressed genes in Metasequoia in CK4 vs. TR4 comparison in biological processes. The frequency of GO terms was analyzed using GO slim assignment. The x-axis and y-axis indicate the names of clusters and the percentage of each cluster, respectively.

DGE Analysis Shows Upregulation of MADS-Box Genes and Downregulation of a Flowering Time Control Protein (FPA) Involved in Early Male Bud Initiation

By comparing transcript abundance in samples from shoots undergoing male bud initiation with samples from CK, 978 (750 upregulated and 228 downregulated, CK1 vs. TR1) and 1675 (118 upregulated and 1557 downregulated, CK2 vs. TR2) significantly DE transcripts were identified in late May and early June, respectively (Figure 5A). When combining the two data sets, 35 transcripts were significantly upregulated or downregulated in TR for male bud initiation. Seven of these genes were also DE (CK1 vs. CK2) in control and (TR1 vs. TR2) treated plants. The gene ID and annotation of the 35 DE genes are listed in Supplementary data 3. Among the 978 DE genes in CK1 vs. TR1 comparison, we focused on two putative transcripts; isotig33092, and isotig32790, both of which were upregulated in TR. By examining these genes, we found a decrease in the transcript isotig33092 putatively encoding Agamous-like MADS-box protein (AGL21) in CK4 when compared to CK2. Another transcript isotig32790 putatively coding for MADS-box transcription factor 50 showed downregulation of TR3 compared to TR1, downregulation of TR4 compared to TR1, and downregulation of CK4 compared to CK2. Among the 1675 DE genes in CK2 and TR2, two transcripts: isotig12610 and isotig12611, both putatively encoding FPA, were downregulated in TR2. When compared among controls, these two transcripts were downregulated in CK4 compared to CK2, but were upregulated in CK2 compared to CK1. Based on these results, these six genes showed enhanced male bud initiation in Metasequoia (Table 2).

FIGURE 5

Figure 5. Venn diagrams with numbers of shared and unique differentially expressed genes. (A) DE genes between CK1 vs. TR1, CK2 vs. TR2, CK1 vs. CK2, and TR1 vs. TR2. (B) DE genes between CK3 vs. TR3, CK4 vs. TR4, CK3 vs. CK4, and TR3 vs. TR4.

TABLE 2

Table 2. Six differentially expressed genes identified in hormone-treated shoots undergoing male bud initiation.

DGE Analysis Indicates Upregulation of Zinc Finger CCCH Domain-Containing Protein Family Factors and FCA during the Female Bud Initiation Stage

By comparing transcript abundance in samples from shoots for female buds with samples from CK, 1678 (324 upregulated and 1354 downregulated, CK3 vs. TR3) and 1498 (1205 upregulated and 293 downregulated, CK4 vs. TR4) significantly DE transcripts were identified, respectively (Figure 5B). When combining the dataset of non-female bud formation comparisons with female bud shoots comparisons, 44 transcripts were common in CK3 vs. TR3, CK4 vs. TR4, and CK3 vs. CK4. We further characterized a transcript isotig14067, which was downregulated in TR3 compared to CK3, upregulated in TR4 compared to CK4, and downregulated in the vegetative bud formation stage (CK4 compared to CK3). This transcript putatively encoded zinc finger CCCH domain-containing protein 28, which plays a significant role in developing floral tissue. Excluding this gene, another six transcripts putatively coding for zinc finger proteins were DE in two stages for female bud formation comparisons (CK3 vs. TR3 and CK4 vs. TR4). Other than these transcripts, two similar transcripts isotig14235 and isotig14236 putative encoding FCA were upregulated in TR4 compared to CK4, but were downregulated in CK4 compared to CK3. Two transcripts (isotig13796 and isotig23456) putatively coding for EARLY FLOWERING 3 (ELF3) and FPA were downregulated in both comparisons (CK3 vs. TR3 and CK3 vs. CK4). A transcriptional co-repressor LEUNIG protein (LUG) was upregulated in TR3 compared to CK3. We believe that these 12 genes constitute a reasonable set of candidate genes involved in the early female bud initiation in Metasequoia (Table 3).

TABLE 3

Table 3. Twelve differentially expressed genes identified in hormone-treated shoots undergoing female bud initiation.

Upregulation of Genes Involved in the Vegetative-To-Reproductive Phase Transition

Based on GO annotation, we selected 7 DE genes with higher expression levels in the TR for male bud initiation (CK1 vs. TR1), and 10 DE genes were upregulated in TR for female bud initiation (CK3 vs. TR3 and CK4 vs. TR4). Based on our results, these genes were activated during the transition phase. Among the genes isotig14814 putatively encodes a LRR receptor-like protein involved in the vegetative-to-reproductive phase transition of meristem. Isotig29943 putatively codes RNA pseudourine synthase 2, which plays a role in spermatogenesis (Table 4).

TABLE 4

Table 4. Upregulated genes in hormone-treated shoots with annotations in the vegetative to reproductive phase transition and flowering development.

MADS-Box Genes Show High Expression Levels in Male Buds or Female Cone Setting Shoots, and Play Important Roles in the Transition from Vegetative to Inflorescence Meristem Identity

Several studies on various plant species have shown that MADS-box genes control various aspects of development and reproductive processes including flower formation. Among the DE genes in our comparisons, genes involved in flower development were identified, including 22 MADS-box family genes with 10 putatively coding for AGL MADS-proteins, five MADS-box proteins GGM13, and seven MADS-box transcription factors. To examine the expression patterns of the 22 MADS-box candidate genes in all samples, we clustered the genes in a heat-map. We observed that the upregulated genes were found in the male and female buds compared to non-cone setting samples, such as isotig16705 and isotig16706 putatively coding for AGL19 (which are upregulated in female buds), and isotig07445 putatively coding for AGL24 (which is highly expressed in male buds) (Figure 6). This analysis provides expression conditions for proteins involved in the transition of vegetative growth to reproductive phase and floral induction. Phylogenetic relationship of MADS box genes from angiosperms and gymnosperms show homolog to the putative transcripts from Metasequoia PUTs (isotig33092, isotig32790, isotig16705, isotig16706, and isotig07445). These five genes belong to a subclass of a MADS box family called MIKC type (Supplementary Figure S2).

FIGURE 6

Figure 6. Heat map showing the relative expression level for Metasequoia MADS-box genes in different stages during vegetative to reproductive phase transition based on DGE data analysis. The color scale (−1.5 to 1.5) represents the Z-score calculated by comparing Reads Per Kilobase of exon model per Million mapped reads (RPKM). Hierarchical clustering of genes and samples was shown in the dendrogram on the top and side of the heatmap using the complete linkage approach.

Identification of DE Genes Related to Hormone Induction

A total of 34 genes were involved in hormone regulation (Supplementary data 4), including Auxin response factor 6, ILR1-like 1, IAA-amino acid hydrolase ILR1, ABC transporter B and G family member, Auxin efflux carrier component 1, and Protein TORNADO 2. Two homeotic transcripts (isotig10488 and 10489) putatively coding for Auxin-responsive protein IAA2 were also upregulated in TR (CK2 vs. TR2). During the female bud initiation stage, two homologous genes (isotig15995 and isotig15996), both putatively coding for IAA-amino acid hydrolase ILR1-like 1, showed increased expression levels (by 1.72-fold) in female cone setting shoots in CK3 vs. TR3 (Table 5).

TABLE 5

Table 5. DE genes related to hormone signals in the vegetative to reproductive phase transition.

DE Genes Involved in Starch and Sucrose Metabolism and the Glycolysis/Gluconeogenesis Pathway

Annotated transcripts involved in important pathways were mapped to the KEGG (Ogata et al., 1999; Hashimoto et al., 2006) pathways using Blast2GO software (Conesa et al., 2005; Gotz et al., 2008). A total of 243 KEGG pathways (Supplementary data 5) were covered, among which 289 transcripts were found in the pathway of starch and sucrose metabolism. A total of 39 transcripts were significantly DE in TR compared to CK. One transcript was upregulated and 13 transcripts were downregulated in the TR for male bud initiation (CK1 vs. TR1 and CK2 vs. TR2). A total of 25 transcripts related to starch and sucrose metabolism were identified; among them, the transcript isotig37134, potentially coding for UDP-glycosyltransferase 85A1, was upregulated during the reproductive bud initiation stage (TR2 vs. TR3 and TR2 vs. TR4). In the glycolysis/gluconeogenesis pathway, 23 DE genes were identified in comparisons of TR with CK (CK1 vs. TR1, CK2 vs. TR2, CK3 vs. TR3, and CK4 vs. TR4). Among them, there was a common gene (isotig17002) putatively encoding bifunctional dihydroflavonol 4-reductase/flavanone 4-reductase in CK1 vs. TR1, CK2 vs. TR2, and CK3 vs. TR3. This gene showed upregulation in comparisons of CK1 vs. TR1 and CK3 vs. TR3, but downregulation in CK2 vs. TR2. This suggested that the glycolysis response for this gene was immediately induced after hormone treatment, but decreased after the transition phase (Supplementary data 6).

Validation of Illumina Sequencing Data by qRT-PCR

To validate the Illumina sequencing data, 32 DE genes identified by Illumina sequencing were tested by RT-PCR. The Ct-values of 18S rRNA of all samples ranged from 25.3 to 27.9 (Supplementary data 7). Selected genes included 17 upregulated genes and 15 downregulated genes. Among them, 29 genes (14 upregulated and 15 downregulated) were consistent between the two expression analysis platforms, with 92.3% showing the same trend of up- or down-regulation (Table 6). The correlation of expression measurements between these 29 consistent genes was 0.81 (R² = 0.81) (Figure 8). Overall, the qRT-PCR and RNA-Seq analyses were consistent.

TABLE 6

Table 6. Selected DE genes tested with RT-PCR.

Levels of Non-Structural Carbohydrates

In vegetative buds, fructose level did not fluctuate during May and June, with a decrease in July (Supplementary Figure S3). In contrast, the glucose level steadily increased during May to July. Similar to glucose, sucrose level increase in June followed by a decrease in July. During male bud initiation, the glucose level in male buds were three times higher than in vegetative buds 11 days after hormone treatment and then declined to the same level as vegetative buds. During female bud initiation, the glucose level in female buds were lower than in vegetative buds. Compared to vegetative buds, reproductive buds also had lower sucrose and fructose contents during May to July.

Discussion

To expand our understanding of the gene expression network involved in the transition from vegetative to reproductive phases in Metasequoia, tag sequencing was used as a tool to mine genes (Marioni et al., 2008), because it performs better in gene discovery and has a higher dynamic range (Robinson et al., 2010). Based on our previous report of the first Metasequoia EST dataset, we developed these expression profiles using DGE Illumina's sequencing to analyze DE genes. Our results demonstrate that deep sequencing of nucleotide DGE tags can be used to resolve expression profiles and to detect differences in transcript abundance over a broad dynamic range. In our test case, which compared hormone-treated shoots to control shoots in four stages for male bud and female bud transition, we identified 69,520 (71.3%) transcripts in a previous database of 24 DGE libraries. Results from this expression analysis provide testable hypotheses to increase our understanding of the functions of DE genes involved in the vegetative-to-reproductive phase transition (male and female bud transition stage) via hormonal induction. In addition, glycolysis/gluconeogenesis was shown to play a role during floral development.

We identified a large number of putative transcripts that were differentially expressed only in comparisons of TR vs. CK. In the male bud initiation stage, some transcripts putatively coding MADS-box transcription factors were actively expressed during early flower development. For example, isotig32790, putative coding for MADS-box transcription factor 50 was one of the candidates, and was downregulated in TR3 compared to TR1, TR4 compared to TR1, and CK4 compared to CK2, but was upregulated in TR1 compared to CK1. These results suggest that this gene may act as a vegetative growth repressor and an early flowering inducer, especially as the male bud inducer. Another candidate gene that could have been annotated as being upregulated in male bud setting shoots is isotig33092, putative coding for AGL21, which shows a similar expression pattern to isotig32790, and is downregulated in CK4 compared to CK2, but is upregulated in TR1 compared to CK1 (Table 2). Several angiosperm genes in MADS gene family have been shown to act as floral integrators during the transition from vegetative-to-reproductive phase (Dorca-Fornell et al., 2011), but other genes such as AGL14 and AGL19 play significant roles in root development (Alvarez-Buylla et al., 2000). In gymnosperms, DAL 19, which is in the same clade as AGL42, is thought to have a similar function and to be involved in the transition from vegetative-to-reproductive growth (Uddenberg et al., 2013). Interestingly, in Arabidopsis, AGL21 (belonging to Type I MADS domains) is expressed in developing embryos and roots after germination (Burgeff et al., 2002), and is expressed in the male bud initiation before the phase change from vegetative-to-reproductive growth. These results suggest that the expansion of isotig33092 putatively coding for AGL21 contributed to transition from vegetative-to-reproductive growth, especially during male floral initiation.

The regulation of flowering time is a complex process with many environmental and internal signals, and a switch from vegetative-to-generative development at the appropriate time. During this process, MADS-box genes have been shown to play important roles in most aspects of angiosperm development, especially in fruit development (de Folter et al., 2005; Gramzow and Theissen, 2010; Lovisetto et al., 2012). In a recent report on the genome of the Norway spruce, it was shown that 278 MADS-box homologs are involved in phase changes in gymnosperms (Nystedt et al., 2013). A total of 22 PUTs with clear homology to MADS boxes were identified in our Metasequoia assembly (Figure 6). Isotig07445 putatively coding for the MADS-box protein AGL24 showed the highest expression during the early stage of male bud formation. Previous studies have revealed that this gene is a positive regulator of flowering in Arabidopsis. Overexpression of AGL24 results in early flowering (de Folter et al., 2005). In the present study, isotig16705 and isotig16706 putatively coding for AGL19 was highly expressed in the female bud formation stage (TR3 and TR4). Previous reports have shown that this gene is involved in root development. Thus, it is important to identity these genes in gymnosperms. We propose that isotig07445, isotig16705, and isotig16706 may play an important role in reproductive bud formation during the transition phase of Metasequoia. However, further analyses about the transcript function are required to increase our understanding of the transition phase from vegetative to reproductive phases in gymnosperms.

Recent studies suggest that both FCA and FPA function in the independent (autonomous) pathway to control floral inhibitor FLOWERING LOCUS C (FLC) mRNA accumulation and flowering time (Schomburg et al., 2001; Feng et al., 2011). In our expression profiles, we did not identify Metasequoia FLC homologs or they did not show any DE. The Arabidopsis gene FCA encodes an RNA-binding protein that promotes floral transition and root development (Macknight et al., 2002). Interestingly, previous studies have shown that FCA is an ABA receptor (Hirayama and Shinozaki, 2007). In our experiments, isotig14235 and isotig14236 putatively coding FCA showed increased expression in hormone-treated shoots during female bud formation, while was downregulated during vegetative development (Table 3). Therefore, our data suggest that FCA may be activated in the female bud transition stage and is involved in the female bud formation development in gymnosperms. In Arabidopsis, FPA has been identified as a regulator of flower development and enables flowering by limiting expression of the floral repressor FLC (Michaels and Amasino, 2001; Schomburg et al., 2001). However, the biological role of FPA is not restricted to flowering, but also influences other aspects of development (Veley and Michaels, 2008). Our data suggest that isotig12610 and isotig12611 putatively encoding FPA in Metasequoia repress male bud initiation, but promote vegetative bud initiation. LEUNIG is a key regulator to repress the Arabidopsis floral homeotic gene AGAMOUS, which regulates the activity of both floral meristem identity genes as well as floral organ identity genes (Conner and Liu, 2000). LEUNIG_HOMOLOG and LEUNIG also play a role during Arabidopsis embryo and floral development by forming a putative corepressor complex with another protein, SEUSS (SEU) (Sitaraman et al., 2008). However, during the female transition stage, upregulation of isotig78340 putatively coding LUG in hormone-treated shoots compared to the control may suggest that it is a floral transition inducer. However, further studies are required to characterize the gymnosperm-specific genes involved in the female transition phase. Recent studies have shown that zinc finger CCCH domain-containing proteins regulate flowering time and play significant roles in developing floral tissue (Colasanti et al., 2006; Chao et al., 2014). In addition, the CONSTANS family of transcription factors, some of which play important roles in photoperiod-induced floral induction, have a unique combination of zinc fingers. For example, zinc finger protein CONSTANS-LIKE genes play a significant role in regulating flowering in Arabidopsis and rice (Griffiths et al., 2003; Wenkel et al., 2006; Tsuji et al., 2011). In our TR2 vs. TR3, TR1 vs. TR2, and TR1 vs. TR3 comparisons, all of the transcripts putatively coding the zinc finger CCCH domain-containing proteins were upregulated in the samples collected at earlier time points, while in TR3 vs. TR4, most of the zinc finger CCCH were also upregulated, suggesting that these transcripts may play a positive role in bud initiation.

In the starch and sucrose metabolism pathways (Figure 7, Supplementary data 6), we observed upregulation of isotig37134 putatively coding for UDP-glycosyltransferase 85A1 in TR2 vs. TR3 and TR2 vs. TR4, and upregulation of the transcripts putatively coding for ADP, ATP carrier protein in TR2 vs. TR4, indicating that these key transcripts are required for bud formation and flowering development. Glycosyltransferases (GTs) are enzymes involved in the modification of secondary metabolites, phytohormones, and xenobiotics (Kumar et al., 2011) that transport sugar molecules to the specific acceptors and form glycosidic bonds. The glycosylation reaction is often the final step in plant natural product synthesis, which stabilizes and provides a safe storage form for potentially toxic aglycones. The UGTs transfer the glycosyl group from an activated nucleoside diphosphate sugar donor (UDP-sugar) to a wide range of sugar or non-sugar acceptors (Sharma et al., 2014). Recently, the UDP-glycosyltransferase 85B1 model in Sorghum has been used to identify cations of amino acid residues important for the catalysis and binding of sugar donors and acceptors in the active site (Thorsoe et al., 2005). However, it remains unclear how the sugar donor specificity of UDP-glycosyltransferase 85A1 has changed. Identification of Metasequoia transcripts putatively coding UGTs can be used to develop plants that react efficiently to stress conditions and may be used for functional analysis. The 20 downregulated genes that potentially play a role in pectinesterase activity (GO:0030599) in TR3 vs. TR4, CK2 vs. TR2, CK3 vs. TR3, and TR2 vs. TR4 comparisons may block carbon flow to pectate (Zhang et al., 2014). During glycolysis/gluconeogenesis, downregulation of transcripts putatively coding alcohol dehydrogenase (NAD) and glutathione dehydrogenase in CK2 vs. TR2 may slow the entry of carbon into the citrate cycle, pyruvate metabolism, and propanoate metabolism, leading to less α-D-glucose-6P being converted into glycerate-3P and α-Dglucose-1P. As a result, the majority of α-D-glucose-6P is transported into the amyloplast for starch and sucrose metabolism. These results suggest that these enzymes are rate-limiting in two important pathways during vegetative bud formation. These findings may facilitate the breeding of Metasequoia to enhance starch accumulation during vegetative to reproductive bud transition and increase our understanding of biological processes in bud formation.

FIGURE 7

Figure 7. Multiple interactions among the genes involved in vegetative to reproductive phase transition and pathways of starch and sucrose metabolism and glycolysis/gluconeogenesis in response to hormone signaling. Colored rectangles correspond with the Matesequoia genes detected using DGE.

In Figure 7, we show that the effects of the sugar-hormone interplay may be mediated by hormones that enable tissues to respond to sugars, resulting in vegetative-to-reproductive phase transition. At the same time, hormone and sugar signaling may independently affect vegetative-to-reproductive phase transition (Cho and Yoo, 2011; Matsoukas, 2014). The actions of all pathways finally converge to control the expression of some floral pathway integrators. In our expression tags, genes putatively coding for auxin response protein (IAA) and ARFs respond to hormone regulation during male and female bud transition. In a previous study of endogenous hormones levels during initiation stage (early June to early August) (Liang and Yin, 1994), we found that while Metasequoia reproductive buds shared similar patterns of IAA, ABA, and GA₁₊₃ dynamics with vegetative buds, the former showing higher hormone contents and greater fluctuations with female buds having the highest contents of the analyzed hormones followed by male buds. The current study substantiated our previous findings with the identified DE transcripts involved in hormone synthesis pathways. For instance, most, if not all, of the DE transcripts with a GO function of gibberellin, IAA, or ABA biosynthetic process, were upregulated in reproductive buds when compared to vegetative buds of the same stage. pathway integrators including transcripts putatively coding AGL24, MADS50, FCA, and LUG were found to regulate male and female phase transition, and flowering locus T (MFT), SOC1, and LEAFY (LFY) identified in our previous report play roles in the late stages (Lee et al., 2008; Zhao et al., 2013). Here, we highlight the potential roles of sugar-hormone interactions in regulating the vegetative to reproductive phase transition; however, further studies are required to understand the antagonistic and agonistic interactions between the sugar and hormone signals in a spatiotemporal manner at the molecular level. Their association with other transcriptional regulatory networks remains unclear, so functional studies of these genes in Metasequoia remain to be done.

Biochemical or physiological changes, such as sugar, can be recognized during the transition of a vegetative apex into a floral apex. The non-structural carbohydrates tested in the present study (sucrose, glucose, and fructose) had a lower level in reproductive buds than in vegetative buds in Metasequoia, with the exception in late May and mid June for glucose. This is consistent with the finding of a decrease of sucrose level in strawberry shoot tips at the beginning of floral induction (Eshghi and Tafazol, 2006). We also found that the glucose level during early male bud initiation was thrice higher than in vegetative buds. These observations are supported by the differential expression of related genes. Isotig26509 involved in gluconeogenesis was found induced in comparisons of CK1 vs. TR1, CK1 vs. CK3, and CK1 vs. CK4. Isotig18149 with an annotation of galactinol-sucrose galactosyltransferase activity was also upregulated in female buds when compared to vegetative buds. Collectively, our data suggest that soluble carbohydrate contents may play a role in Metasequoia floral induction: lower levels of sucrose and fructose facilitate initiation of both types of reproductive buds, while elevated glucose content is favorable in male bud initiation and the opposite in female buds.

In a previous study, we identified 17,086 PUTs in late female bud initiation stage that were not found in vegetative buds during late August and early September (Zhao et al., 2013). In early spring (March) when buds rapidly enlarged, 25,643 PUTs were found uniquely in female cones. Our current study identified a total of 6166 DE PUTs between reproductive and vegetative bud initiation (late May–early July), which was the most DE PUTs among all comparisons. Clearly the hormone treatments induced significant transcriptomic changes that led to the vegetative apex into a floral apex. It is worth noting isotig63639, a DE PUT putatively encoding a histone deacetylase. This type of the enzymes allows histone proteins to wrap the DNA more tightly leading to low or no expression of target genes. Isotig63639 was found downregulated in TR1, suggesting that target genes were turned on or expressed at a higher level upon hormone treatments. In addition, while both male and female buds showed decreased expression of ELF3, Flowering time control protein FPA and Zinc finger CCCH domain-containing protein encoding genes, differences in existed between male and female buds during their initiation (Tables 2, 3). Transcript isotig32790 putatively coding for MADS-box transcription factor 50 showed upregulation in male bud initiation and downregulation in female bud initiation. Another upregulated transcript in male bud initiation encoded an Agamous-like MADS-box protein AGL21. In contrast, PUTs for transcriptional corepressor LEUNIG and Flowering time control protein FCA were found upregualted in female bud initiation. There results were consistent with the differences of male and female buds in initiation time, hormone response, and sugar contents.

In summary, the present study identified numerous DE genes among the four stages of bud development, and provided novel information on female and male bud initiation. The results from the in silico DE analyses were reliable because 29 out of the 32 genes chosen for qPCR showed the same upregulation or downregulation trend as the in silico analyses (Table 6, Figure 8). Overall, this study identified a large number of DE genes regulated in the male and female bud transition and provides new insights into the interaction and overlapping genes in the sugar-mediated flowering pathway during the bud-formation stage.

FIGURE 8

Figure 8. Correlation between qPCR and Illumina sequencing data. Correlation between qRT-PCR and Illumina sequencing data of 32 selected genes: 17 upregulated genes and 15 downregulated genes in 16 pairs of amplified RNA samples. Spearman Rank Correlation coefficient = 0.81 (P < 0.05).

Conflict of Interest Statement

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Acknowledgments

This research was supported by The Fundamental Research Funds for the Central Universities of China (No. BLYJ201309), the National Natural Science Foundation of China (31270656), Joint Programs of the Scientific Research and Graduate Training from BMEC, Program for Changjiang Scholars and Innovative Research Team in University (IRT13047) and “111”Project of Beijing Forestry University (B13007).

Supplementary Material

The Supplementary Material for this article can be found online at: http://journal.frontiersin.org/article/10.3389/fpls.2015.00467

References

Alvarez-Buylla, E. R., Liljegren, S. J., Pelaz, S., Gold, S. E., Burgeff, C., Ditta, G. S., et al. (2000). MADS-box gene evolution beyond flowers: expression in pollen, endosperm, guard cells, roots and trichomes. Plant J. 24, 457–466. doi: 10.1046/j.1365-313x.2000.00891.x

PubMed Abstract | CrossRef Full Text

Birol, I., Raymond, A., Jackman, S. D., Pleasance, S., Coope, R., Taylor, G. A., et al. (2013). Assembling the 20 Gb white spruce (Picea glauca) genome from whole-genome shotgun sequencing data. Bioinformatics 29, 1492–1497. doi: 10.1093/bioinformatics/btt178

PubMed Abstract | CrossRef Full Text

Bohlenius, H., Huang, T., Charbonnel-Campaa, L., Brunner, A. M., Jansson, S., Strauss, S. H., et al. (2006). CO/FT regulatory module controls timing of flowering and seasonal growth cessation in trees. Science 312, 1040–1043. doi: 10.1126/science.1126038

PubMed Abstract | CrossRef Full Text | Google Scholar

Brunsfeld, S. J., Soltis, P. S., Soltis, D. E., Gadek, P. A., Quinn, C. J., Strenge, D. D., et al. (1994). Phylogenetic relationships among the genera of Taxodiaceae and Cupressaceae: evidence from rbcL sequences. Syst. Bot. 19, 253–262. doi: 10.2307/2419600

CrossRef Full Text | Google Scholar

Bureau, S., Quilot-Turion, B. N. D., Signoret, V. R., Renaud, C., Maucourt, M., Bancel, D., et al. (2013). Determination of the composition in sugars and organic acids in peach using mid infrared spectroscopy: comparison of prediction results according to data sets and different reference methods. Anal. Chem. 85, 11312–11318. doi: 10.1021/ac402428s

PubMed Abstract | CrossRef Full Text | Google Scholar

Burgeff, C., Liljegren, S. J., Tapia-Lopez, R., Yanofsky, M. F., and Alvarez-Buylla, E. R. (2002). MADS-box gene expression in lateral primordia, meristems and differentiated tissues of Arabidopsis thaliana roots. Planta 214, 365–372. doi: 10.1007/s004250100637

PubMed Abstract | CrossRef Full Text | Google Scholar

Carlsbecker, A., Tandre, K., Johanson, U., Englund, M., and Engstrom, P. (2004). The MADS-box gene DAL1 is a potential mediator of the juvenile-to-adult transition in Norway spruce (Picea abies). Plant J. 40, 546–557. doi: 10.1111/j.1365-313X.2004.02226.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Chao, Y., Zhang, T., Yang, Q., Kang, J., Sun, Y., Gruber, M. Y., et al. (2014). Expression of the alfalfa CCCH-type zinc finger protein gene MsZFN delays flowering time in transgenic Arabidopsis thaliana. Plant Sci. 215, 92–99. doi: 10.1016/j.plantsci.2013.10.012

PubMed Abstract | CrossRef Full Text | Google Scholar

Cho, Y.-H., and Yoo, S.-D. (2011). Signaling role of fructose mediated by FINS1/FBP in Arabidopsis thaliana. PLoS Genet. 7:e1001263. doi: 10.1371/journal.pgen.1001263

PubMed Abstract | CrossRef Full Text | Google Scholar

Colasanti, J., Tremblay, R., Wong, A. Y. M., Coneva, V., Kozaki, A., and Mable, B. K. (2006). The maize INDETERMINATEI flowering time regulator defines a highly conserved zinc finger protein family in higher plants. BMC Genomics 7:158. doi: 10.1186/1471-2164-7-158

PubMed Abstract | CrossRef Full Text | Google Scholar

Conesa, A., Götz, S., García-Gómez, J. M., Terol, J., Talón, M., and Robles, M. (2005). Blast2GO: a universal tool for annotation, visualization and analysis in functional genomics research. Bioinformatics 21, 3674–3676. doi: 10.1093/bioinformatics/bti610

PubMed Abstract | CrossRef Full Text | Google Scholar

Conner, J., and Liu, Z. (2000). LEUNIG, a putative transcriptional corepressor that regulates AGAMOUS expression during flower development. Proc. Natl. Acad. Sci. U.S.A. 97, 12902–12907. doi: 10.1073/pnas.230352397

PubMed Abstract | CrossRef Full Text | Google Scholar

de Folter, S., Immink, R. G. H., Kieffer, M., Parenicova, L., Henz, S. R., Weigel, D., et al. (2005). Comprehensive interaction map of the Arabidopsis MADS box transcription factors. Plant Cell 17, 1424–1433. doi: 10.1105/tpc.105.031831

PubMed Abstract | CrossRef Full Text | Google Scholar

Dorca-Fornell, C., Gregis, V., Grandi, V., Coupland, G., Colombo, L., and Kater, M. M. (2011). The Arabidopsis SOC1-like genes AGL42, AGL71 and AGL72 promote flowering in the shoot apical and axillary meristems. Plant J. 67, 1006–1017. doi: 10.1111/j.1365-313X.2011.04653.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Dörken, V. M. (2012). Proliferating seed cones in Metasequoia glyptostroboidesHu & Cheng (Cupressaceae s.l., Coniferales) elucidate the evolution of seed cones and ovules in Cupressaceae s.l. Feddes Repert. 122, 409–420. doi: 10.1002/fedr.201200004

CrossRef Full Text

Eshghi, S., and Tafazol, E. (2006). Possible role of non-structural carbohydrates in flower induction in strawberry. J. Hortic. Sci. Biotechnol. 81, 854–858.

Google Scholar

Feng, W., Jacob, Y., Veley, K. M., Ding, L., Yu, X., Choe, G., et al. (2011). Hypomorphic alleles reveal FCA-independent roles for FY in the regulation of FLOWERING LOCUS C. Plant Physiol. 155, 1425–1434. doi: 10.1104/pp.110.167817

PubMed Abstract | CrossRef Full Text | Google Scholar

Ferrario, S., Busscher, J., Franken, J., Gerats, T., Vandenbussche, M., Angenent, G. C., et al. (2004). Ectopic expression of the petunia MADS box gene UNSHAVEN accelerates flowering and confers leaf-like characteristics to floral organs in a dominant-negative manner. Plant Cell 16, 1490–1505. doi: 10.1105/tpc.019679

PubMed Abstract | CrossRef Full Text | Google Scholar

Goeman, J. J., van de Geer, S. A., de Kort, F., and van Houwelingen, H. C. (2004). A global test for groups of genes: testing association with a clinical outcome. Bioinformatics 20, 93–99. doi: 10.1093/bioinformatics/btg382

PubMed Abstract | CrossRef Full Text | Google Scholar

Gotz, S., Garcia-Gomez, J. M., Terol, J., Williams, T. D., Nagaraj, S. H., Nueda, M. J., et al. (2008). High-throughput functional annotation and data mining with the Blast2GO suite. Nucleic Acids Res. 36, 3420–3435. doi: 10.1093/nar/gkn176

PubMed Abstract | CrossRef Full Text | Google Scholar

Gramzow, L., and Theissen, G. (2010). A hitchhiker's guide to the MADS world of plants. Genome Biol. 11, 214. doi: 10.1186/gb-2010-11-6-214

PubMed Abstract | CrossRef Full Text | Google Scholar

Griffiths, S., Dunford, R. P., Coupland, G., and Laurie, D. A. (2003). The evolution of CONSTANS-like gene families in barley, rice, and Arabidopsis. Plant Physiol. 131, 1855–1867. doi: 10.1104/pp.102.016188

PubMed Abstract | CrossRef Full Text | Google Scholar

Haas, B. J., Papanicolaou, A., Yassour, M., Grabherr, M., Blood, P. D., Bowden, J., et al. (2013). De novo transcript sequence reconstruction from RNA-seq using the Trinity platform for reference generation and analysis. Nat. Protoc. 8, 1494–1512. doi: 10.1038/nprot.2013.084

PubMed Abstract | CrossRef Full Text | Google Scholar

Hashimoto, K., Goto, S., Kawano, S., Aoki-Kinoshita, K. F., Ueda, N., Hamajima, M., et al. (2006). KEGG as a glycome informatics resource. Glycobiology 16, 63–70. doi: 10.1093/glycob/cwj010

PubMed Abstract | CrossRef Full Text | Google Scholar

Hirayama, T., and Shinozaki, K. (2007). Perception and transduction of abscisic acid signals: keys to the function of the versatile plant hormone ABA. Trends Plant Sci. 12, 343–351. doi: 10.1016/j.tplants.2007.06.013

PubMed Abstract | CrossRef Full Text | Google Scholar

Huang, H.-R., Yan, P.-C., Lascoux, M., and Ge, X.-J. (2012). Flowering time and transcriptome variation in Capsella bursa-pastoris (Brassicaceae). New Phytol. 194, 676–689. doi: 10.1111/j.1469-8137.2012.04101.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Jin, B., Tang, L., Lu, Y., Wang, D., Zhang, M., and Ma, J. (2012). Temporal and spatial characteristics of male cone development in Metasequoia glyptostroboides Hu et Cheng. Plant Signal. Behav. 7, 1687–1694. doi: 10.4161/psb.22898

PubMed Abstract | CrossRef Full Text | Google Scholar

Klintenas, M., Pin, P. A., Benlloch, R., Ingvarsson, P. K., and Nilsson, O. (2012). Analysis of conifer FLOWERING LOCUS T/TERMINAL FLOWER1-like genes provides evidence for dramatic biochemical evolution in the angiosperm FT lineage. New Phytol. 196, 1260–1273. doi: 10.1111/j.1469-8137.2012.04332.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Kumar, R., Tyagi, A. K., and Sharma, A. K. (2011). Genome-wide analysis of auxin response factor (ARF) gene family from tomato and analysis of their role in flower and fruit development. Mol. Genet. Genomics 285, 245–260. doi: 10.1007/s00438-011-0602-7

PubMed Abstract | CrossRef Full Text | Google Scholar

Kunzmann, L., and Mai, D. H. (2011). The first record of fossil Metasequoia (Cupressaceae) from continental Europe. Rev. Palaeobot. Palynol. 164, 247–250. doi: 10.1016/j.revpalbo.2011.01.005

CrossRef Full Text | Google Scholar

Langmead, B., and Salzberg, S. L. (2012). Fast gapped-read alignment with Bowtie 2. Nat. Methods 9, 357–U354. doi: 10.1038/nmeth.1923

PubMed Abstract | CrossRef Full Text | Google Scholar

Lee, J., Oh, M., Park, H., and Lee, I. (2008). SOC1 translocated to the nucleus by interaction with AGL24 directly regulates LEAFY. Plant J. 55, 832–843. doi: 10.1111/j.1365-313X.2008.03552.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Li, P., Wind, J. J., Shi, X., Zhang, H., Hanson, J., Smeekens, S. C., et al. (2011). Fructose sensitivity is suppressed in Arabidopsis by the transcription factor ANAC089 lacking the membrane-bound domain. Proc. Natl. Acad. Sci. U.S.A. 108, 3436–3441. doi: 10.1073/pnas.1018665108

PubMed Abstract | CrossRef Full Text | Google Scholar

Li, Y. Y., Chen, X. Y., Zhang, X., Wu, T. Y., Lu, H. P., and Cai, Y. W. (2005). Genetic differences between wild and artificial populations of Metasequoia glyptostroboides: implications for species recovery. Conserv. Biol. 19, 224–231. doi: 10.1111/j.1523-1739.2005.00025.x

CrossRef Full Text | Google Scholar

Liang, H., and Yin, W.-L. (1994). “Changes of IAA, ABA, and GA1+3 contents in leaf and cone buds of Metasequoia glyptostroboides Hu et Cheng,” in Growth and Development Control and Biotechnology in Woody Plants, eds S. S. Wang and X. N. Jiang (Beijing: China Forestry Publishing House), 215–219.

Lovisetto, A., Guzzo, F., Tadiello, A., Toffali, K., Favretto, A., and Casadoro, G. (2012). Molecular analyses of MADS-Box genes trace back to gymnosperms the invention of fleshy fruits. Mol. Biol. Evol. 29, 409–419. doi: 10.1093/molbev/msr244

PubMed Abstract | CrossRef Full Text | Google Scholar

Lunn, J. E., Feil, R., Hendriks, J. H. M., Gibon, Y., Morcuende, R., Osuna, D., et al. (2006). Sugar-induced increases in trehalose 6-phosphate are correlated with redox activation of ADPglucose pyrophosphorylase and higher rates of starch synthesis in Arabidopsis thaliana. Biochem. J. 397, 139–148. doi: 10.1042/BJ20060083

PubMed Abstract | CrossRef Full Text

Macknight, R., Duroux, M., Laurie, R., Dijkwel, P., Simpson, G., and Dean, C. (2002). Functional significance of the alternative transcript processing of the Arabidopsis floral promoter FCA. Plant Cell 14, 877–888. doi: 10.1105/tpc.010456

PubMed Abstract | CrossRef Full Text | Google Scholar

Marioni, J. C., Mason, C. E., Mane, S. M., Stephens, M., and Gilad, Y. (2008). RNA-seq: an assessment of technical reproducibility and comparison with gene expression arrays. Genome Res. 18, 1509–1517. doi: 10.1101/gr.079558.108

PubMed Abstract | CrossRef Full Text | Google Scholar

Matsoukas, I. G. (2014). Interplay between sugar and hormone signaling pathways modulate floral signal transduction. Front. Genet. 5:218. doi: 10.3389/fgene.2014.00218

PubMed Abstract | CrossRef Full Text | Google Scholar

Matsumura, H., Yoshida, K., Luo, S., Kimura, E., Fujibe, T., Albertyn, Z., et al. (2010). High-throughput superSAGE for digital gene expression analysis of multiple samples using next generation sequencing. PLoS ONE 5:e12010. doi: 10.1371/journal.pone.0012010

PubMed Abstract | CrossRef Full Text | Google Scholar

Michaels, S. D., and Amasino, R. M. (2001). Loss of FLOWERING LOCUS C activity eliminates the late-flowering phenotype of FRIGIDA and autonomous pathway mutations but not responsiveness to vernalization. Plant Cell 13, 935–941. doi: 10.1105/tpc.13.4.935

PubMed Abstract | CrossRef Full Text | Google Scholar

Moghaddam, M. R. B., and Van-den-Ende, W. (2013). Sugars, the clock and transition to flowering. Front. Plant Sci. 4:22. doi: 10.3389/fpls.2013.00022

PubMed Abstract | CrossRef Full Text | Google Scholar

Moon, J., Suh, S. S., Lee, H., Choi, K. R., Hong, C. B., Paek, N. C., et al. (2003). The SOC1 MADS-box gene integrates vernalization and gibberellin signals for flowering in Arabidopsis. Plant J. 35, 613–623. doi: 10.1046/j.1365-313X.2003.01833.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Neale, D. B., Wegrzyn, J. L., Stevens, K. A., Zimin, A. V., Puiu, D., Crepeau, M. W., et al. (2014). Decoding the massive genome of loblolly pine using haploid DNA and novel assembly strategies. Genome Biol. 15, 59. doi: 10.1186/gb-2014-15-3-r59

CrossRef Full Text | Google Scholar

Nelson, E. C. (1998). Metasequoia glyptostroboides, the dawn redwood: some Irish glosses on its discovery and introduction into cultivation. Curtis's Bot. Mag. 15, 77–80. doi: 10.1111/1467-8748.00135

CrossRef Full Text | Google Scholar

Nystedt, B., Street, N. R., Wetterbom, A., Zuccolo, A., Lin, Y.-C., Scofield, D. G., et al. (2013). The Norway spruce genome sequence and conifer genome evolution. Nature 497, 579–584. doi: 10.1038/nature12211

PubMed Abstract | CrossRef Full Text | Google Scholar

Ogata, H., Goto, S., Sato, K., Fujibuchi, W., Bono, H., and Kanehisa, M. (1999). KEGG: kyoto encyclopedia of genes and genomes. Nucleic Acids Res. 27, 29–34. doi: 10.1093/nar/27.1.29

PubMed Abstract | CrossRef Full Text | Google Scholar

Poethig, R. S. (2013). Vegetative phase change and shoot maturation in plants. Curr. Top. Dev. Biol. 105, 125–152. doi: 10.1016/b978-0-12-396968-2.00005-1

PubMed Abstract | CrossRef Full Text | Google Scholar

Robinson, M. D., McCarthy, D. J., and Smyth, G. K. (2010). edgeR: a Bioconductor package for differential expression analysis of digital gene expression data. Bioinformatics 26, 139–140. doi: 10.1093/bioinformatics/btp616

PubMed Abstract | CrossRef Full Text | Google Scholar

Rolland, F., Baena-Gonzalez, E., and Sheen, J. (2006). Sugar sensing and signaling in plants: conserved and novel mechanisms. Annu. Rev. Plant Biol. 57, 675–709. doi: 10.1146/annurev.arplant.57.032905.105441

PubMed Abstract | CrossRef Full Text | Google Scholar

Schmittgen, T. D., and Livak, K. J. (2008). Analyzing real-time PCR data by the comparative CT method. Nat. Protoc. 3, 1101–1108. doi: 10.1038/nprot.2008.73

PubMed Abstract | CrossRef Full Text | Google Scholar

Schomburg, F. M., Patton, D. A., Meinke, D. W., and Amasino, R. M. (2001). FPA, a gene involved in floral induction in Arabidopsis, encodes a protein containing RNA-recognition motifs. Plant cell 13, 1427–1436. doi: 10.1105/tpc.13.6.1427

PubMed Abstract | CrossRef Full Text | Google Scholar

Sharma, R., Rawat, V., and Suresh, C. G. (2014). Genome-wide identification and tissue-specific expression analysis of UDP-glycosyltransferases genes confirm their abundance in cicer arietinum (Chickpea) genome. PLoS ONE 9:e109715. doi: 10.1371/journal.pone.0109715

PubMed Abstract | CrossRef Full Text | Google Scholar

Sitaraman, J., Bui, M., and Liu, Z. (2008). LEUNIG_HOMOLOG and LEUNIG perform partially redundant functions during Arabidopsis embryo and floral development. Plant Physiol. 147, 672–681. doi: 10.1104/pp.108.115923

PubMed Abstract | CrossRef Full Text | Google Scholar

Sterling, C. (1949). Some features in the morphology of Metasequoia. Am. J. Bot. 36, 461–471. doi: 10.2307/2438080

CrossRef Full Text | Google Scholar

Tamura, K., Peterson, D., Peterson, N., Stecher, G., Nei, M., and Kumar, S. (2011). MEGA5: molecular evolutionary genetics analysis using maximum likelihood, evolutionary distance, and maximum parsimony methods. Mol. Biol. Evol. 28, 2731–2739. doi: 10.1093/molbev/msr121

PubMed Abstract | CrossRef Full Text | Google Scholar

Thorsoe, K. S., Bak, S., Olsen, C. E., Imberty, A., Breton, C., and Moller, B. L. (2005). Determination of catalytic key amino acids and UDP sugar donor specificity of the cyanohydrin glycosyltransferase UGT85B1 from Sorghum bicolor. Molecular modeling substantiated by site-specific mutagenesis and biochemical analyses. Plant Physiol. 139, 664–673. doi: 10.1104/pp.105.063842

PubMed Abstract | CrossRef Full Text | Google Scholar

Trevaskis, B., Hemming, M. N., Peacock, W. J., and Dennis, E. S. (2006). HvVRN2 responds to daylength, whereas HvVRN1 is regulated by vernalization and developmental status. Plant Physiol. 140, 1397–1405. doi: 10.1104/pp.105.073486

PubMed Abstract | CrossRef Full Text | Google Scholar

Tsuji, H., Taoka, K.-I., and Shimamoto, K. (2011). Regulation of flowering in rice: two florigen genes, a complex gene network, and natural variation. Curr. Opin. Plant Biol. 14, 45–52. doi: 10.1016/j.pbi.2010.08.016

PubMed Abstract | CrossRef Full Text | Google Scholar

Turnbull, C. (2011). Long-distance regulation of flowering time. J. Exp. Bot. 62, 4399–4413. doi: 10.1093/jxb/err191

PubMed Abstract | CrossRef Full Text | Google Scholar

Uddenberg, D., Reimegard, J., Clapham, D., Almqvist, C., Von Arnold, S., Emanuelsson, O., et al. (2013). Early cone setting in Picea abies acrocona is associated with increased transcriptional activity of a MADS Box transcription factor. Plant Physiol. 161, 813–823. doi: 10.1104/pp.112.207746

PubMed Abstract | CrossRef Full Text | Google Scholar

Veley, K. M., and Michaels, S. D. (2008). Functional redundancy and new roles for genes of the autonomous floral-promotion pathway. Plant Physiol. 147, 682–695. doi: 10.1104/pp.108.118927

PubMed Abstract | CrossRef Full Text | Google Scholar

Wang, J.-W., Park, M. Y., Wang, L.-J., Koo, Y., Chen, X.-Y., Weigel, D., et al. (2011). MiRNA control of vegetative phase change in trees. PLoS Genet. 7:e1002012. doi: 10.1371/journal.pgen.1002012

PubMed Abstract | CrossRef Full Text | Google Scholar

Wang, Z., Gerstein, M., and Snyder, M. (2009). RNA-Seq: a revolutionary tool for transcriptomics. Nat. Rev. Genet. 10, 57–63. doi: 10.1038/nrg2484

PubMed Abstract | CrossRef Full Text | Google Scholar

Wenkel, S., Turck, F., Singer, K., Gissot, L., Le Gourrierec, J., Samach, A., et al. (2006). CONSTANS and the CCAAT box binding complex share a functionally important domain and interact to regulate flowering of Arabidopsis. Plant Cell 18, 2971–2984. doi: 10.1105/tpc.106.043299

PubMed Abstract | CrossRef Full Text | Google Scholar

Winfield, M. O., Lu, C., Wilson, I. D., Coghill, J. A., and Edwards, K. J. (2009). Cold- and light-induced changes in the transcriptome of wheat leading to phase transition from vegetative to reproductive growth. BMC Plant Biol. 9:55. doi: 10.1186/1471-2229-9-55

PubMed Abstract | CrossRef Full Text | Google Scholar

Yu, S., Cao, L., Zhou, C.-M., Zhang, T.-Q., Lian, H., Sun, Y., et al. (2013). Sugar is an endogenous cue for juvenile-to-adult phase transition in plants. Elife 2:e00269. doi: 10.7554/eLife.00269

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhang, P., Li, C., Zhang, P., Jin, C., Pan, D., and Bao, Y. (2014). iTRAQ-based proteomics reveals novel members involved in pathogen challenge in sea cucumber Apostichopus japonicus. PLoS ONE 9:e100492. doi: 10.1371/journal.pone.0100492

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhao, S., and Fernald, R. D. (2005). Comprehensive algorithm for quantitative real-time polymerase chain reaction. J. Comput. Biol. 12, 1047–1064. doi: 10.1089/cmb.2005.12.1047

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhao, Y., Thammannagowda, S., Staton, M., Tang, S., Xia, X., Yin, W., et al. (2013). An EST dataset for Metasequoia glyptostroboides buds: the first EST resource for molecular genomics studies in Metasequoia. Planta 237, 755–770. doi: 10.1007/s00425-012-1783-y

PubMed Abstract | CrossRef Full Text | Google Scholar

Zimin, A., Stevens, K. A., Crepeau, M., Holtz-Morris, A., Koriabine, M., Marcais, G., et al. (2014). Sequencing and assembly of the 22-Gb loblolly pine genome. Genetics 196, 875–890. doi: 10.1534/genetics.113.159715

PubMed Abstract | CrossRef Full Text | Google Scholar